CAINDMS2SCG2 -- INITIALIZE FOR FFT 00010000 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CTITLE S2SCG2 -- INITIALIZE FOR FFT 00020000 CA AUTHOR S. NELAN/MIKE MAROLDA 00030000 CA DESIGNER S. NELAN 00040000 CA LANGUAGE FORTRAN 00050000 CA SYSTEM IBM (SEE CRAY) 00060000 CA WRITTEN 11-18-81 00070000 C REVISED 03-15-86 MJM/ESN CREATED FOR IBM VECTOR 00080000 C FACILITY ESSL FROM CRAY CODE. 00090000 C REVISED 09-08-86 ESN CORRECTED LENGTH OF AUX ARRAYS00100000 C REQUIRED FOR SCFT. 00110000 C REVISED 09-11-86 ESN PASS IN LENGTH OF AUX ARRAYS 00120000 C ESSL SUBS IN DOUBLE-WORDS. 00130000 C REVISED 01-07-87 RDK CORRECT NUMBER OF WORDS 00140000 C COPIED TO AUX4 AREA IN COMPLEX 00150000 C FFTS. REVISE NAUX4. 00160000 C REVISED 01-20-87 ESN MERGED IN MFORSP. 00170000 C REVISED 04-28-87 ESN USE GETMN2 TO KEEP ALL TEMP 00180000 C ARRAYS. 00190000 C REVISED 03-05-92 REP. CHANGED TO NOT GIVE ERROR 00200000 C ON SECOND ENTRY TO S2SCG2, JUST 00210000 C RETURN. 00220000 C REVISED 03-06-92 ESN INCREASE THE SIZE OF THE 00221001 C AUX ARRAYS FOR THE 3090J. 00222001 CA 00230000 CA 00240000 CA CALL S2SCG2 (MAXMAG, WORK) 00250000 CA INPUT MAXMAG = MAXIMUM MAGNITUDE FOR FOURIER TRANSFORMS I4 00260000 CA I/0 WORK = BEGINNING ADDRESS FOR WORK ARRAY R4 00270000 CA 00280000 CA 00290000 CA INITIALIZES FOR THE FFT SUBROUTINES. 00300000 CAEND 00310000 C 00320000 CTITLE S2DFT2 -- FAST FOURIER TRANSFORM (REAL) 00330000 CA AUTHOR S. NELAN 00340000 CA DESIGNER S. NELAN 00350000 CA LANGUAGE FORTRAN 00360000 CA SYSTEM IBM (SEE CRAY) 00370000 CA WRITTEN 11-18-81 00380000 C REVISED 00390000 CA 00400000 CA 00410000 CA CALL S2DFT2 (MAG, X, *ERR) 00420000 CA INPUT MAG = MAGNITUDE OF TRANSFORM I4 00430000 CA I/O X = BASE ADDRESS OF INPUT/OUTPUT ARRAY R4 00440000 CA ERR = ERROR RETURN WHEN MAG > MAXMAG FROM S2SCG2 00450000 CA 00460000 CA 00470000 CA THIS ROUTINE PERFORMS A DIRECT FAST FOURIER TRANSFORM BY THE 00480000 CA USE OF THE CFFT2 SUBROUTINE ON THE CRAY. THE INPUT X CONTAINS 00490000 CA NP = 2**MAG REAL VALUES. ON OUTPUT X CONTAINS NY = NP/2 + 1 00500000 CA COMPLEX VALUES CORRESPONDING TO FREQUENCIES FROM ZERO THROUGH 00510000 CA THE NYQUIST. THE J-TH COMPLEX VALUE IS GIVEN BY 00520000 CA 00530000 CA NP 00540000 CA SQRT(2/NP) * SUM( X(K) * EXP(-I*2*PI*(J-1)*(K-1)/NP) ) 00550000 CA K=1 00560000 CA 00570000 CA MAG MUST BE LESS THAN OR EQUAL TO MAXMAG USED IN THE CALL TO 00580000 CA S2SCG2, WHICH MUST BE CALLED PRIOR TO S2DFT2. 00590000 CA 00600000 CAEND 00610000 C 00620000 CTITLE S2DFI2 -- INVERSE FAST FOURIER TRANSFORM (REAL) 00630000 CA AUTHOR S. NELAN 00640000 CA DESIGNER S. NELAN 00650000 CA LANGUAGE FORTRAN 00660000 CA SYSTEM IBM (SEE CRAY) 00670000 CA WRITTEN 11-18-81 00680000 C REVISED 00690000 CA 00700000 CA 00710000 CA CALL S2DFI2 (MAG, X, *ERR) 00720000 CA INPUT MAG = MAGNITUDE OF TRANSFORM I4 00730000 CA I/O X = BASE ADDRESS OF INPUT/OUTPUT ARRAY R4 00740000 CA ERR = ERROR RETURN WHEN MAG > MAXMAG FROM S2SCG2 00750000 CA 00760000 CA 00770000 CA THIS ROUTINE PERFORMS AN INVERSE FAST FOURIER TRANSFORM BY THE 00780000 CA USE OF THE CFTT2 SUBROUTINE ON THE CRAY. THE INPUT X CONTAINS 00790000 CA NY = NP/2 + 1 COMPLEX VALUES CORRESPONDING TO FREQUENCIES FROM 00800000 CA ZERO THROUGH THE NYQUIST, WHERE NP = 2**MAG. ON OUTPUT X 00810000 CA CONTAINS NP REAL VALUES. THE J-TH REAL VALUE IS GIVEN BY 00820000 CA 00830000 CA NY 00840000 CA SQRT(.5/NP) * SUM( REAL( C(K) * E(K) * EXP(I*2*PI*(J-1)*(K-1)/NP) ))00850000 CA K=1 00860000 CA 00870000 CA WHERE C(K) IS THE K-TH COMPLEX VALUE, E(K) = 1 FOR K = 1 OR NY, 00880000 CA E(K) = 2 FOR 1 < K < NY, AND PI = 3.14159265... . MAG MUST BE 00890000 CA LESS THAN OR EQUAL TO MAXMAG USED IN THE CALL TO S2SCG2, WHICH 00900000 CA MUST BE CALLED PRIOR TO S2DFI2. 00910000 CA 00920000 CAEND 00930000 C 00940000 CTITLE S2FFT2 -- FAST FOURIER TRANSFORM (COMPLEX) 00950000 CA AUTHOR S. NELAN 00960000 CA DESIGNER S. NELAN 00970000 CA LANGUAGE FORTRAN 00980000 CA SYSTEM IBM (SEE CRAY) 00990000 CA WRITTEN 11-18-81 01000000 C REVISED 01010000 CA 01020000 CA 01030000 CA CALL S2FFT2 (MAG, X, *ERR) 01040000 CA INPUT MAG = MAGNITUDE OF TRANSFORM I4 01050000 CA I/O X = BASE ADDRESS OF INPUT/OUTPUT ARRAY C4 01060000 CA ERR = ERROR RETURN WHEN MAG > MAXMAG FROM S2SCG2 01070000 CA 01080000 CA 01090000 CA THIS ROUTINE PERFORMS A DIRECT FAST FOURIER TRANSFORM ON THE 01100000 CA COMPLEX INPUT ARRAY X BY THE USE OF THE CFFT2 SUBROUTINE ON THE 01110000 CA CRAY. THE LENGTH OF THE FFT IS 2**MAG. MAG MUST BE LESS THAN 01120000 CA OR EQUAL TO MAXMAG USED IN THE CALL TO S2SCG2. 01130000 CA 01140000 CA 01150000 CAEND 01160000 C 01170000 CTITLE S2FFI2 -- INVERSE FAST FOURIER TRANSFORM (COMPLEX) 01180000 CA AUTHOR S. NELAN 01190000 CA DESIGNER S. NELAN 01200000 CA LANGUAGE FORTRAN 01210000 CA SYSTEM IBM (SEE CRAY) 01220000 CA WRITTEN 11-18-81 01230000 C REVISED 01240000 CA 01250000 CA 01260000 CA CALL S2FFI2 (MAG, X, *ERR) 01270000 CA INPUT MAG = MAGNITUDE OF TRANSFORM I4 01280000 CA I/O X = BASE ADDRESS OF INPUT/OUTPUT ARRAY C4 01290000 CA ERR = ERROR RETURN WHEN MAG > MAXMAG FROM S2SCG2 01300000 CA 01310000 CA 01320000 CA THIS ROUTINE PERFORMS AN INVERSE FAST FOURIER TRANSFORM ON THE 01330000 CA COMPLEX INPUT ARRAY X BY USE OF THE CFFT2 SUBROUTINE ON THE CRAY. 01340000 CA THE LENGTH OF THE FFT IS 2**MAG. MAG MUST BE LESS THAN OR EQUAL 01350000 CA TO MAXMAG USED IN THE CALL TO S2SCG2. 01360000 CA 01370000 CAEND 01380000 C 01390000 CTITLE MFORSP -- RADIX 2 FFT PERFORMED IN PLACE 01400000 CA 01410000 CA AUTHOR D. D. THOMPSON 01420000 CA DESIGNER D. D. THOMPSON 01430000 CA LANGUAGE FORTRAN 01440000 CA SYSTEM CRAY 01450000 CA WRITTEN 1972 01460000 C REVISED 06/81; SINGLE PRECISION FFT BASED 01470000 C ON DOUBLE PRECISION ROUTINE MFORAX. 01480000 C REVISED 04/86; RDK; MOVED TO THE CRAY. 01490000 C REVISED 01/87; ESN; IMPLEMENTED IN S2SCG2. 01500000 CA 01510000 CA CALL MFORSP (N, X, Y, ISIGN) 01520000 CA 01530000 CA IN/OUT ARG TYPE DESCRIPTION 01540000 CA 01550000 CA IN N I4 BASE 2 LOG OF THE FFT LENGTH 01560000 CA IN/OUT X R4 REAL DATA ARRAY 01570000 CA IN/OUT Y R4 IMAGINARY DATA ARRAY 01580000 CA IN ISIGN I4 INPUT FLAG 01590000 CA = -1 FOR DIRECT TRANSFORM 01600000 CA = +1 FOR NORMALIZED INVERSE TRANSFORM 01610000 CA = +2 FOR UNNORMALIZED INVERSE TRANSFORM 01620000 CA 01630000 CA THIS ROUTINE PERFORMS EITHER FORWARD OR REVERSE FFT IN PLACE. 01640000 CA 01650000 CAEND 01660000 C 01670000 C ALL SIGN CONVENTIONS AND SCALING OF RESULTS IS HANDLED WITHIN 01680000 C THE INDIVIDUAL ROUTINES. 01690000 C 01700000 SUBROUTINE S2SCG2 (MAXMAG, WORK) 01710000 C 01720000 C ============================== 01730000 C REAL ARRAYS -- INPUT ARGUMENTS 01740000 C ============================== 01750000 C 01760000 REAL WORK (1) 01770000 REAL X (1) 01780000 REAL Y (1) 01790000 C 01800000 C ======================= 01810000 C INTEGER ARRAYS -- LOCAL 01820000 C ======================= 01830000 C 01840000 INTEGER IPT0M1(15) 01850000 INTEGER IPT0P1(15) 01860000 INTEGER IPT1M1(15) 01870000 INTEGER IPT1P1(15) 01880000 INTEGER IPT2M1(15) 01890000 INTEGER IPT2P1(15) 01900000 INTEGER LAUX1 (15) 01910000 INTEGER LAUX2 (15) 01920000 C 01930000 C ===================================== 01940000 C DATA STATEMENTS (FROM IBM ESSL GUIDE) 01950000 C ===================================== 01960000 C 01970000 C DOUBLE WORD AUXILIARY STORAGE 01980000 C REQUIREMENTS FOR SRCFT 01990000 C DATA LAUX1 / 25000, 25000, 25000, 25000, 25000, 02000000 C * 25000, 25000, 25000, 25000, 25000, 02010000 C * 25000, 25000, 25000, 25000, 47000 / 02020000 C DATA LAUX2 / 20000, 20000, 20000, 20000, 20000, 02030000 C * 20000, 20000, 20000, 20000, 20000, 02040000 C * 20000, 20000, 20000, 20000, 38700 / 02050000 C DOUBLE WORD AUXILIARY STORAGE 02060000 C REQUIREMENTS FOR SCRFT 02070000 C DATA LAUX1 / 25000, 25000, 25000, 25000, 25000, 02080000 C * 25000, 25000, 25000, 25000, 25000, 02090000 C * 25000, 25000, 25000, 25000, 47000 / 02100000 C DATA LAUX2 / 20000, 20000, 20000, 20000, 20000, 02110000 C * 20000, 20000, 20000, 20000, 20000, 02120000 C * 20000, 20000, 20000, 20000, 38700 / 02130000 C DOUBLE WORD AUXILIARY STORAGE 02140000 C REQUIREMENTS FOR SCFT 02150000 C (NOTE THAT THE ARRAY HAS BEEN SHIFTED 02160000 C ONE LOCATION FOR STANDARDIZATION WITH 02170000 C SRCFT AND SCRFT, AND ONE LOCATION FOR 02180000 C SPARC COMPATABILITY) 02190000 C DATA LAUX1 / 20000, 20000, 20000, 20000, 20000, 02200000 C * 20000, 20000, 20000, 20000, 20000, 02210000 C * 20000, 20000, 20000, 38700, 57400 / 02220000 C DATA LAUX2 / 20000, 20000, 20000, 20000, 20000, 02230000 C * 20000, 20000, 20000, 20000, 20000, 02240000 C * 20000, 20000, 20000, 38700, 57400 / 02250000 C MAX DOUBLE WORD AUXILIARY STORAGE 02260000 C REQUIREMENTS FOR SRCFT, SCRFT, SCFT 02270000 DATA LAUX1 / 25000, 25000, 25000, 25000, 25000, 02280000 * 25000, 25000, 25000, 25000, 25000, 02290000 * 25000, 25000, 25000, 38700, 57400 / 02300000 DATA LAUX2 / 20000, 20000, 20000, 20000, 20000, 02310000 * 20000, 20000, 20000, 20000, 20000, 02320000 * 20000, 20000, 20000, 38700, 57400 / 02330000 C 02340000 DATA IPT0M1 / 15*0 / 02350000 DATA IPT0P1 / 15*0 / 02360000 DATA IPT1M1 / 15*0 / 02370000 DATA IPT1P1 / 15*0 / 02380000 DATA IPT2M1 / 15*0 / 02390000 DATA IPT2P1 / 15*0 / 02400000 DATA LASTIX / -999 / 02410000 C 02420000 C ============== 02430000 C INITIALIZATION 02440000 C ============== 02450000 C 02460000 IF (MAXMAG .LT. 1) GO TO 9010 02470000 IF (MAXMAG .GT. 15) GO TO 9010 02480000 C 02490000 IF (LASTIX .NE. -999) GO TO 8000 02500000 LASTIX = 0 02510000 C 02520000 IAUX3 = 1 02530000 NAUX3 = 0 02540000 NAUX4 = 2**MAXMAG 02550000 C 02560000 LENI = 2 * NAUX4 02570000 CALL GETMN2 (WORK, LENI, I, LENO) 02580000 IF (LENI .NE. LENO) CALL XDUMPX 02590000 IAUX4 = I + 1 02600000 C 02610000 GO TO 8000 02620000 C ============================================================== 02630000 C 02640000 ENTRY S2DFT2 (MAG, X, *) 02650000 C 02660000 IF (MAG .GT. MAXMAG) GO TO 9000 02670000 C 02680000 NPTS = 2**MAG 02690000 SCALE = 2.0/SQRT(NPTS*2.0) 02700000 C 02710000 IF (IPT0P1(MAG) .NE. 0) GO TO 200 02720000 LENI = 2 * (LAUX1(MAG)+LAUX2(MAG)) 02730000 CALL GETMN2 (WORK, LENI, I, LENO) 02740000 IF (LENI .NE. LENO) CALL XDUMPX 02750000 IPT0P1(MAG) = I + 1 02760000 IAUX1 = IPT0P1(MAG) 02770000 NAUX1 = LAUX1(MAG) 02780000 IAUX2 = IAUX1 + 2*NAUX1 02790000 NAUX2 = LAUX2(MAG) 02800000 C 02810000 I = LOC(X(1)) 02820000 IF ((I/8)*8 .NE. I) GO TO 110 02830000 CALL SRCFT (1,X,NPTS,X,NPTS,NPTS,1,1,SCALE, 02840000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 02850000 * WORK(IAUX3),NAUX3) 02860000 GO TO 210 02870000 C 02880000 110 DO 120 I = 1, NPTS 02890000 120 WORK(IAUX4+I-1) = X(I) 02900000 CALL SRCFT (1,WORK(IAUX4),NPTS,WORK(IAUX4),NPTS,NPTS, 02910000 * 1,1,SCALE, 02920000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 02930000 * WORK(IAUX3),NAUX3) 02940000 GO TO 415 02950000 C 02960000 200 IAUX1 = IPT0P1(MAG) 02970000 NAUX1 = LAUX1(MAG) 02980000 IAUX2 = IAUX1 + 2*NAUX1 02990000 NAUX2 = LAUX2(MAG) 03000000 I = LOC(X(1)) 03010000 IF ((I/8)*8 .NE. I) GO TO 400 03020000 210 CALL SRCFT (0,X,NPTS,X,NPTS,NPTS,1,1,SCALE, 03030000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 03040000 * WORK(IAUX3),NAUX3) 03050000 GO TO 8000 03060000 C 03070000 400 DO 410 I = 1, NPTS 03080000 410 WORK(IAUX4+I-1) = X(I) 03090000 415 CALL SRCFT (0,WORK(IAUX4),NPTS,WORK(IAUX4),NPTS,NPTS, 03100000 * 1,1,SCALE, 03110000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 03120000 * WORK(IAUX3),NAUX3) 03130000 NPTSP2 = NPTS + 2 03140000 DO 420 I = 1, NPTSP2 03150000 420 X(I) = WORK(IAUX4+I-1) 03160000 C 03170000 GO TO 8000 03180000 C ============================================================== 03190000 C 03200000 ENTRY S2DFI2 (MAG, X, *) 03210000 C 03220000 IF (MAG .GT. MAXMAG) GO TO 9000 03230000 C 03240000 NPTS = 2**MAG 03250000 SCALE = 1.0/SQRT(NPTS*2.0) 03260000 C 03270000 IF (IPT0M1(MAG) .NE. 0) GO TO 1200 03280000 LENI = 2 * (LAUX1(MAG)+LAUX2(MAG)) 03290000 CALL GETMN2 (WORK, LENI, I, LENO) 03300000 IF (LENI .NE. LENO) CALL XDUMPX 03310000 IPT0M1(MAG) = I + 1 03320000 IAUX1 = IPT0M1(MAG) 03330000 NAUX1 = LAUX1(MAG) 03340000 IAUX2 = IAUX1 + 2*NAUX1 03350000 NAUX2 = LAUX2(MAG) 03360000 C 03370000 I = LOC(X(1)) 03380000 IF ((I/8)*8 .NE. I) GO TO 1110 03390000 CALL SCRFT (1, X,NPTS,X,NPTS,NPTS,1,-1,SCALE, 03400000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 03410000 * WORK(IAUX3),NAUX3) 03420000 GO TO 1210 03430000 C 03440000 1110 NPTSP2 = NPTS + 2 03450000 DO 1120 I = 1, NPTSP2 03460000 1120 WORK(IAUX4+I-1) = X(I) 03470000 CALL SCRFT (1,WORK(IAUX4),NPTS,WORK(IAUX4),NPTS,NPTS, 03480000 * 1,-1,SCALE, 03490000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 03500000 * WORK(IAUX3),NAUX3) 03510000 GO TO 1415 03520000 C 03530000 1200 IAUX1 = IPT0M1(MAG) 03540000 NAUX1 = LAUX1(MAG) 03550000 IAUX2 = IAUX1 + 2*NAUX1 03560000 NAUX2 = LAUX2(MAG) 03570000 I = LOC(X(1)) 03580000 IF ((I/8)*8 .NE. I) GO TO 1400 03590000 1210 CALL SCRFT (0, X,NPTS,X,NPTS,NPTS,1,-1,SCALE, 03600000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 03610000 * WORK(IAUX3),NAUX3) 03620000 GO TO 8000 03630000 C 03640000 1400 NPTSP2 = NPTS + 2 03650000 DO 1410 I = 1, NPTSP2 03660000 1410 WORK(IAUX4+I-1) = X(I) 03670000 1415 CALL SCRFT (0,WORK(IAUX4),NPTS,WORK(IAUX4),NPTS,NPTS, 03680000 * 1,-1,SCALE, 03690000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2, 03700000 * WORK(IAUX3),NAUX3) 03710000 DO 1420 I = 1, NPTS 03720000 1420 X(I) = WORK(IAUX4+I-1) 03730000 C 03740000 GO TO 8000 03750000 C ============================================================== 03760000 C 03770000 ENTRY S2FFT2 (MAG, X, *) 03780000 C 03790000 IX = 1 03800000 GO TO 2000 03810000 C ============================================================== 03820000 C 03830000 ENTRY S2FFI2 (MAG, X, *) 03840000 C 03850000 IX = -1 03860000 2000 IF (MAG .GT. MAXMAG) GO TO 9000 03870000 C 03880000 NPTS = 2**MAG 03890000 SCALE = 1.0/SQRT(NPTS*1.0) 03900000 C 03910000 IF (IX .EQ. 1) THEN 03920000 IF (IPT1P1(MAG) .NE. 0) GO TO 2200 03930000 ELSE 03940000 IF (IPT1M1(MAG) .NE. 0) GO TO 2200 03950000 ENDIF 03960000 LENI = 2 * (LAUX1(MAG)+LAUX2(MAG)) 03970000 CALL GETMN2 (WORK, LENI, I, LENO) 03980000 IF (LENI .NE. LENO) CALL XDUMPX 03990000 IF (IX .EQ. 1) THEN 04000000 IPT1P1(MAG) = I + 1 04010000 IAUX1 = IPT1P1(MAG) 04020000 ELSE 04030000 IPT1M1(MAG) = I + 1 04040000 IAUX1 = IPT1M1(MAG) 04050000 ENDIF 04060000 NAUX1 = LAUX1(MAG) 04070000 IAUX2 = IAUX1 + 2*NAUX1 04080000 NAUX2 = LAUX2(MAG) 04090000 C 04100000 I = LOC(X(1)) 04110000 IF ((I/8)*8 .NE. I) GO TO 2110 04120000 CALL SCFT (1,X,1,NPTS,X,1,NPTS,NPTS,1,IX,SCALE, 04130000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2) 04140000 GO TO 2210 04150000 C 04160000 2110 NPTS2 = NPTS * 2 + 0 04170000 DO 2120 I = 1, NPTS2 04180000 2120 WORK(IAUX4+I-1) = X(I) 04190000 CALL SCFT (1,WORK(IAUX4),1,NPTS,WORK(IAUX4),1,NPTS,NPTS, 04200000 * 1,IX,SCALE, 04210000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2) 04220000 GO TO 2415 04230000 C 04240000 2200 IF (IX .EQ. 1) THEN 04250000 IAUX1 = IPT1P1(MAG) 04260000 ELSE 04270000 IAUX1 = IPT1M1(MAG) 04280000 ENDIF 04290000 NAUX1 = LAUX1(MAG) 04300000 IAUX2 = IAUX1 + 2*NAUX1 04310000 NAUX2 = LAUX2(MAG) 04320000 I = LOC(X(1)) 04330000 IF ((I/8)*8 .NE. I) GO TO 2400 04340000 2210 CALL SCFT (0, X,1,NPTS,X,1,NPTS,NPTS,1,IX,SCALE, 04350000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2) 04360000 GO TO 8000 04370000 C 04380000 2400 NPTS2 = NPTS * 2 + 0 04390000 DO 2410 I = 1, NPTS2 04400000 2410 WORK(IAUX4+I-1) = X(I) 04410000 2415 CALL SCFT (0,WORK(IAUX4),1,NPTS,WORK(IAUX4),1,NPTS,NPTS, 04420000 * 1,IX,SCALE, 04430000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2) 04440000 NPTS2 = NPTS * 2 04450000 DO 2420 I = 1, NPTS2 04460000 2420 X(I) = WORK(IAUX4+I-1) 04470000 C 04480000 GO TO 8000 04490000 C ============================================================== 04500000 C 04510000 ENTRY MFORSP (MAG, X, Y, ISIGN) 04520000 C 04530000 IF (MAG .GT. MAXMAG) GO TO 9000 04540000 C 04550000 NPTS = 2**MAG 04560000 SCALE = 1.0 04570000 C 04580000 IX = -ISIGN 04590000 IF (ISIGN .EQ. 2) IX = -1 04600000 C 04610000 DO 3100 I = 1, NPTS 04620000 3100 WORK(IAUX4+2*I-2) = X(I) 04630000 DO 3200 I = 1, NPTS 04640000 3200 WORK(IAUX4+2*I-1) = Y(I) 04650000 C 04660000 IF (IX .EQ. 1) THEN 04670000 IF (IPT2P1(MAG) .NE. 0) GO TO 3400 04680000 ELSE 04690000 IF (IPT2M1(MAG) .NE. 0) GO TO 3400 04700000 ENDIF 04710000 LENI = 2 * (LAUX1(MAG)+LAUX2(MAG)) 04720000 CALL GETMN2 (WORK, LENI, I, LENO) 04730000 IF (LENI .NE. LENO) CALL XDUMPX 04740000 IF (IX .EQ. 1) THEN 04750000 IPT2P1(MAG) = I + 1 04760000 IAUX1 = IPT2P1(MAG) 04770000 ELSE 04780000 IPT2M1(MAG) = I + 1 04790000 IAUX1 = IPT2M1(MAG) 04800000 ENDIF 04810000 NAUX1 = LAUX1(MAG) 04820000 IAUX2 = IAUX1 + 2*NAUX1 04830000 NAUX2 = LAUX2(MAG) 04840000 C 04850000 CALL SCFT (1,WORK(IAUX4),1,NPTS,WORK(IAUX4),1,NPTS,NPTS, 04860000 * 1,IX,SCALE, 04870000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2) 04880000 C 04890000 3400 IF (IX .EQ. 1) THEN 04900000 IAUX1 = IPT2P1(MAG) 04910000 ELSE 04920000 IAUX1 = IPT2M1(MAG) 04930000 ENDIF 04940000 NAUX1 = LAUX1(MAG) 04950000 IAUX2 = IAUX1 + 2*NAUX1 04960000 NAUX2 = LAUX2(MAG) 04970000 CALL SCFT (0,WORK(IAUX4),1,NPTS,WORK(IAUX4),1,NPTS,NPTS, 04980000 * 1,IX,SCALE, 04990000 * WORK(IAUX1),NAUX1,WORK(IAUX2),NAUX2) 05000000 C 05010000 IF (ISIGN .EQ. 1) SCALE = 1.0 / NPTS 05020000 DO 3600 I = 1, NPTS 05030000 3600 X(I) = WORK(IAUX4+2*I-2) * SCALE 05040000 DO 3700 I = 1, NPTS 05050000 3700 Y(I) = WORK(IAUX4+2*I-1) * SCALE 05060000 C 05070000 GO TO 8000 05080000 C 05090000 C ==== 05100000 C EXIT 05110000 C ==== 05120000 C 05130000 8000 RETURN 05140000 C 05150000 C ERROR EXIT 05160000 C 05170000 9000 CONTINUE 05180000 WRITE (6,9500) MAG, MAXMAG 05190000 RETURN 1 05200000 C 05210000 9010 WRITE (6,9510) MAXMAG 05220000 CALL XDUMPX 05230000 GO TO 8000 05240000 C 05250000 C FORMAT STATEMENTS 05260000 C 05270000 9500 FORMAT (1X,'*** INPUT MAG OF',I10,' IS GREATER THAN MAXMAG OF', 05280000 * I10) 05290000 C 05300000 9510 FORMAT (1X,'*** MAXIMUM FFT POWER OF 2 INPUT IS',I10,', WHICH IS',05310000 * ' OUT OF THE ALLOWABLE RANGE ***') 05320000 C 05330000 END 05340000