CTITLEMSFFT -- FOURIER TRANSFORM, MULTIVARIATE COMPLEX 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR/DESIGNER R. C. SINGLETON 00000020 CA LANGUAGE S/370 FORTRAN H 00000030 CA WRITTEN OCT. 1968 00000040 CA MODIFIED R. G. CURRIE, 2-28-79. 00000050 CA 00000060 CA 00000070 CA CALL MSFFT (A, B, NTOTS, NSS, NSPANS, ISN) 00000080 CA MULTIVARIATE COMPLEX FOURIER TRANSFORM, COMPUTED IN PLACE 00000090 CA USING MIXED-RADIX FAST FOURIER TRANSFORM ALGORITHM. 00000100 CA BY R.C. SINGLETON, STANFORD RESEARCH INSTITUTE, OCT. 1968. 00000110 CA IEEE TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS, 00000120 CA VOL. AU-17, NO. 2, JUNE 1969, PAGES 93-103. 00000130 CA ARRAYS A AND B ORIGINALLY HOLD THE REAL AND IMAGINARY 00000140 CA COMPONENTS OF THE DATA, AND RETURN THE REAL AND 00000150 CA IMAGINARY COMPONENTS OF THE RESULTING FOURIER COEFFICIENTS. 00000160 CA MULTIVARIATE DATA IS INDEXED ACCORDING TO THE FORTRAN 00000170 CA ARRAY ELEMENT SUCCESSOR FUNCTION, WITHOUT LIMIT 00000180 CA ON THE NUMBER OF IMPLIED MULTIPLE SUBSCRIPTS. 00000190 CA THE SUBROUTINE IS CALLED ONCE FOR EACH VARIATE. 00000200 CA THE CALLS FOR A MULTIVARIATE TRANSFORM MAY BE IN ANY ORDER. 00000210 CA NTOT IS THE TOTAL NUMBER OF COMPLEX DATA VALUES. 00000220 CA N IS THE DIMENSION OF THE CURRENT VARIABLE. 00000230 CA NSPAN/N IS THE SPACING OF CONSECUTIVE DATA VALUES 00000240 CA WHILE INDEXING THE CURRENT VARIABLE. 00000250 CA THE SIGN OF ISN DETERMINES THE SIGN OF THE COMPLEX 00000260 CA EXPONENTIAL, AND THE MAGNITUDE OF ISN IS NORMALLY ONE. 00000270 CA A TRI-VARIATE TRANSFORM WITH A(N1,N2,N3),B(N1,N2,N3) 00000280 CA IS COMPUTED BY 00000290 CA CALL MSFFT (A, B, N1*N2*N3, N1, N1, 1) 00000300 CA CALL MSFFT (A, B, N1*N2*N3, N2, N1*N2, 1) 00000310 CA CALL MSFFT (A, B, N1*N2*N3, N3, N1*N2*N3, 1) 00000320 CA FOR A SINGLE-VARIATE TRANSFORM, 00000330 CA NTOT = N = NSPAN = (NUMBER OF COMPLEX DATA VALUES), E.G. 00000340 CA CALL MSFFT (A, B, N, N, N, 1) 00000350 CA THE DATA MAY ALTERNATIVELY BE STORED IN A SINGLE COMPLEX 00000360 CA ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO 00000370 CA GIVE THE CORRECT INDEXING INCREMENT AND A(2) USED TO 00000380 CA PASS THE INITIAL ADDRESS FOR THE SEQUENCE OF IMAGINARY 00000390 CA VALUES, E.G. 00000400 CA CALL MSFFT (A, A(2), NTOT, N, NSPAN, 2) 00000410 CA ARRAYS AT(MAXF), CK(MAXF), BT(MAXF),SK(MAXF), AND NP(MAXP) 00000420 CA ARE USED FOR TEMPORARY STORAGE. IF THE AVAILABLE STORAGE 00000430 CA IS INSUFFICIENT, THE PROGRAM IS TERMINATED BY A STOP. 00000440 CA MAXF MUST BE .GE. THE MAXIMUM PRIME FACTOR OF N. 00000450 CA MAXP MUST BE .GT. THE NUMBER OF PRIME FACTORS OF N. 00000460 CA IN ADDITION, IF THE SQUARE-FREE PORTION K OF N HAS TWO OR 00000470 CA MORE PRIME FACTORS, THEN MAXP MUST BE .GE. K-1. 00000480 SUBROUTINE MSFFT (A, B, NTOTS, NSS, NSPANS, ISN) 00000490 DIMENSION A(1),B(1) 00000500 C ARRAY STORAGE IN NFAC FOR A MAXIMUM OF 11 FACTORS OF N. 00000510 C IF N HAS MORE THAN ONE SQUARE-FREE FACTOR, THE PRODUCT OF THE 00000520 C SQUARE-FREE FACTORS MUST BE .LE.210 00000530 DIMENSION NFAC(11),NP(209) 00000540 C ARRAY STORAGE FOR MAXIMUM PRIME FACTOR OF 23 00000550 DIMENSION AT(23),CK(23),BT(23),SK(23) 00000560 EQUIVALENCE (I,II) 00000570 C THE FOLLOWING TWO CONSTANTS SHOULD AGREE WITH THE ARRAY DIMENSIONS. 00000580 MAXF=23 00000590 MAXP=209 00000600 C CURRIE CODE TO INCORPORATE FFT DIRECTLY INTO HAGEN'S QUAD 00000610 NTOT=2**NTOTS 00000620 N=2**NSS 00000630 NSPAN=2**NSPANS 00000640 C CURRIE CODE TO INCORPORATE FFT DIRECTLY INTO HAGEN'S QUAD 00000650 IF(N .LT. 2) RETURN 00000660 INC=ISN 00000670 C CURRIE CODE FOR ISN=-2 (INVERSE NORMALIZED) 00000680 IF(ISN.EQ.-2) INC=INC+1 00000690 C CURRIE CODE FOR ISN=-2 (INVERSE NORMALIZED) 00000700 RAD=8.0*ATAN(1.0) 00000710 S72=RAD/5.0 00000720 C72=COS(S72) 00000730 S72=SIN(S72) 00000740 S120=SQRT(0.75) 00000750 IF(ISN .GE. 0) GO TO 10 00000760 S72=-S72 00000770 S120=-S120 00000780 RAD=-RAD 00000790 INC=-INC 00000800 10 NT=INC*NTOT 00000810 KS=INC*NSPAN 00000820 KSPAN=KS 00000830 NN=NT-INC 00000840 JC=KS/N 00000850 RADF=RAD*FLOAT(JC)*0.5 00000860 I=0 00000870 JF=0 00000880 C DETERMINE THE FACTORS OF N 00000890 M=0 00000900 K=N 00000910 GO TO 20 00000920 15 M=M+1 00000930 NFAC(M)=4 00000940 K=K/16 00000950 20 IF(K-(K/16)*16 .EQ. 0) GO TO 15 00000960 J=3 00000970 JJ=9 00000980 GO TO 30 00000990 25 M=M+1 00001000 NFAC(M)=J 00001010 K=K/JJ 00001020 30 IF(MOD(K,JJ).EQ.0) GO TO 25 00001030 J=J+2 00001040 JJ=J**2 00001050 IF(JJ .LE. K) GO TO 30 00001060 IF(K .GT. 4) GO TO 40 00001070 KT=M 00001080 NFAC(M+1)=K 00001090 IF(K .NE. 1) M=M+1 00001100 GO TO 80 00001110 40 IF(K-(K/4)*4.NE.0) GO TO 50 00001120 M=M+1 00001130 NFAC(M)=2 00001140 K=K/4 00001150 50 KT=M 00001160 J=2 00001170 60 IF(MOD(K,J).NE.0) GO TO 70 00001180 M=M+1 00001190 NFAC(M)=J 00001200 K=K/J 00001210 70 J=((J+1)/2)*2+1 00001220 IF(J .LE. K) GO TO 60 00001230 80 IF(KT .EQ. 0) GO TO 100 00001240 J=KT 00001250 90 M=M+1 00001260 NFAC(M)=NFAC(J) 00001270 J=J-1 00001280 IF(J .NE. 0) GO TO 90 00001290 C COMPUTE FOURIER TRANSFORM 00001300 100 SD=RADF/FLOAT(KSPAN) 00001310 CD=2.0*SIN(SD)**2 00001320 SD=SIN(SD+SD) 00001330 KK=1 00001340 I=I+1 00001350 IF(NFAC(I) .NE. 2) GO TO 400 00001360 C TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR) 00001370 KSPAN=KSPAN/2 00001380 K1=KSPAN+2 00001390 210 K2=KK+KSPAN 00001400 AK=A(K2) 00001410 BK=B(K2) 00001420 A(K2)=A(KK)-AK 00001430 B(K2)=B(KK)-BK 00001440 A(KK)=A(KK)+AK 00001450 B(KK)=B(KK)+BK 00001460 KK=K2+KSPAN 00001470 IF(KK .LE. NN) GO TO 210 00001480 KK=KK-NN 00001490 IF(KK .LE. JC) GO TO 210 00001500 IF(KK .GT. KSPAN) GO TO 800 00001510 220 C1=1.0-CD 00001520 S1=SD 00001530 230 K2=KK+KSPAN 00001540 AK=A(KK)-A(K2) 00001550 BK=B(KK)-B(K2) 00001560 A(KK)=A(KK)+A(K2) 00001570 B(KK)=B(KK)+B(K2) 00001580 A(K2)=C1*AK-S1*BK 00001590 B(K2)=S1*AK+C1*BK 00001600 KK=K2+KSPAN 00001610 IF(KK .LT. NT) GO TO 230 00001620 K2=KK-NT 00001630 C1=-C1 00001640 KK=K1-K2 00001650 IF(KK .GT. K2) GO TO 230 00001660 AK=C1-(CD*C1+SD*S1) 00001670 S1=(SD*C1-CD*S1)+S1 00001680 C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION 00001690 C ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE 00001700 C C1=AK 00001710 C1=0.5/(AK**2+S1**2)+0.5 00001720 S1=C1*S1 00001730 C1=C1*AK 00001740 KK=KK+JC 00001750 IF(KK .LT. K2) GO TO 230 00001760 K1=K1+INC+INC 00001770 KK=(K1-KSPAN)/2+JC 00001780 IF (KK .LE. JC+JC) GO TO 220 00001790 GO TO 100 00001800 C TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE) 00001810 320 K1=KK+KSPAN 00001820 K2=K1+KSPAN 00001830 AK=A(KK) 00001840 BK=B(KK) 00001850 AJ=A(K1)+A(K2) 00001860 BJ=B(K1)+B(K2) 00001870 A(KK)=AK+AJ 00001880 B(KK)=BK+BJ 00001890 AK=-0.5*AJ+AK 00001900 BK=-0.5*BJ+BK 00001910 AJ=(A(K1)-A(K2))*S120 00001920 BJ=(B(K1)-B(K2))*S120 00001930 A(K1)=AK-BJ 00001940 B(K1)=BK+AJ 00001950 A(K2)=AK+BJ 00001960 B(K2)=BK-AJ 00001970 KK=K2+KSPAN 00001980 IF(KK .LT. NN) GO TO 320 00001990 KK=KK-NN 00002000 IF(KK .LE. KSPAN) GO TO 320 00002010 GO TO 700 00002020 C TRANSFORM FOR FACTOR OF 4 00002030 400 IF(NFAC(I) .NE. 4) GO TO 600 00002040 KSPNN=KSPAN 00002050 KSPAN=KSPAN/4 00002060 410 C1=1.0 00002070 S1=0.0 00002080 420 K1=KK+KSPAN 00002090 K2=K1+KSPAN 00002100 K3=K2+KSPAN 00002110 AKP=A(KK)+A(K2) 00002120 AKM=A(KK)-A(K2) 00002130 AJP=A(K1)+A(K3) 00002140 AJM=A(K1)-A(K3) 00002150 A(KK)=AKP+AJP 00002160 AJP=AKP-AJP 00002170 BKP=B(KK)+B(K2) 00002180 BKM=B(KK)-B(K2) 00002190 BJP=B(K1)+B(K3) 00002200 BJM=B(K1)-B(K3) 00002210 B(KK)=BKP+BJP 00002220 BJP=BKP-BJP 00002230 IF(ISN .LT. 0) GO TO 450 00002240 AKP=AKM-BJM 00002250 AKM=AKM+BJM 00002260 BKP=BKM+AJM 00002270 BKM=BKM-AJM 00002280 IF(S1) 430,460,430 00002290 430 A(K1)=AKP*C1-BKP*S1 00002300 B(K1)=AKP*S1+BKP*C1 00002310 A(K2)=AJP*C2-BJP*S2 00002320 B(K2)=AJP*S2+BJP*C2 00002330 A(K3)=AKM*C3-BKM*S3 00002340 B(K3)=AKM*S3+BKM*C3 00002350 KK=K3+KSPAN 00002360 IF(KK .LE. NT) GO TO 420 00002370 440 C2=C1-(CD*C1+SD*S1) 00002380 S1=(SD*C1-CD*S1)+S1 00002390 C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION 00002400 C ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE 00002410 C C1=C2 00002420 C1=0.5/(C2**2+S1**2)+0.5 00002430 S1=C1*S1 00002440 C1=C1*C2 00002450 C2=C1**2-S1**2 00002460 S2=2.0*C1*S1 00002470 C3=C2*C1-S2*S1 00002480 S3=C2*S1+S2*C1 00002490 KK=KK-NT+JC 00002500 IF(KK .LE. KSPAN) GO TO 420 00002510 KK=KK-KSPAN+INC 00002520 IF(KK .LE. JC) GO TO 410 00002530 IF(KSPAN .EQ. JC) GO TO 800 00002540 GO TO 100 00002550 450 AKP=AKM+BJM 00002560 AKM=AKM-BJM 00002570 BKP=BKM-AJM 00002580 BKM=BKM+AJM 00002590 IF(S1) 430,460,430 00002600 460 A(K1)=AKP 00002610 B(K1)=BKP 00002620 A(K2)=AJP 00002630 B(K2)=BJP 00002640 A(K3)=AKM 00002650 B(K3)=BKM 00002660 KK=K3+KSPAN 00002670 IF(KK .LE. NT) GO TO 420 00002680 GO TO 440 00002690 C TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE) 00002700 510 C2=C72**2-S72**2 00002710 S2=2.0*C72*S72 00002720 520 K1=KK+KSPAN 00002730 K2=K1+KSPAN 00002740 K3=K2+KSPAN 00002750 K4=K3+KSPAN 00002760 AKP=A(K1)+A(K4) 00002770 AKM=A(K1)-A(K4) 00002780 BKP=B(K1)+B(K4) 00002790 BKM=B(K1)-B(K4) 00002800 AJP=A(K2)+A(K3) 00002810 AJM=A(K2)-A(K3) 00002820 BJP=B(K2)+B(K3) 00002830 BJM=B(K2)-B(K3) 00002840 AA=A(KK) 00002850 BB=B(KK) 00002860 A(KK)=AA+AKP+AJP 00002870 B(KK)=BB+BKP+BJP 00002880 AK=AKP*C72+AJP*C2+AA 00002890 AJ=AKM*S72+AJM*S2 00002900 BK=BKP*C72+BJP*C2+BB 00002910 BJ=BKM*S72+BJM*S2 00002920 A(K1)=AK-BJ 00002930 A(K4)=AK+BJ 00002940 B(K1)=BK+AJ 00002950 B(K4)=BK-AJ 00002960 AK=AKP*C2+AJP*C72+AA 00002970 BK=BKP*C2+BJP*C72+BB 00002980 AJ=AKM*S2-AJM*S72 00002990 BJ=BKM*S2-BJM*S72 00003000 A(K2)=AK-BJ 00003010 A(K3)=AK+BJ 00003020 B(K2)=BK+AJ 00003030 B(K3)=BK-AJ 00003040 KK=K4+KSPAN 00003050 IF(KK .LT. NN) GO TO 520 00003060 KK=KK-NN 00003070 IF(KK .LE. KSPAN) GO TO 520 00003080 GO TO 700 00003090 C TRANSFORM FOR ODD FACTORS 00003100 600 K=NFAC(I) 00003110 KSPNN=KSPAN 00003120 KSPAN=KSPAN/K 00003130 IF(K .EQ. 3) GO TO 320 00003140 IF(K .EQ. 5) GO TO 510 00003150 IF(K .EQ. JF) GO TO 640 00003160 JF=K 00003170 S1=RAD/FLOAT(K) 00003180 C1=COS(S1) 00003190 S1=SIN(S1) 00003200 IF(JF .GT. MAXF) GO TO 998 00003210 CK(JF)=1.0 00003220 SK(JF)=0.0 00003230 J=1 00003240 630 CK(J)=CK(K)*C1+SK(K)*S1 00003250 SK(J)=CK(K)*S1-SK(K)*C1 00003260 K=K-1 00003270 CK(K)=CK(J) 00003280 SK(K)=-SK(J) 00003290 J=J+1 00003300 IF(J .LT. K) GO TO 630 00003310 640 K1=KK 00003320 K2=KK+KSPNN 00003330 AA=A(KK) 00003340 BB=B(KK) 00003350 AK=AA 00003360 BK=BB 00003370 J=1 00003380 K1=K1+KSPAN 00003390 650 K2=K2-KSPAN 00003400 J=J+1 00003410 AT(J)=A(K1)+A(K2) 00003420 AK=AT(J)+AK 00003430 BT(J)=B(K1)+B(K2) 00003440 BK=BT(J)+BK 00003450 J=J+1 00003460 AT(J)=A(K1)-A(K2) 00003470 BT(J)=B(K1)-B(K2) 00003480 K1=K1+KSPAN 00003490 IF(K1 .LT. K2) GO TO 650 00003500 A(KK)=AK 00003510 B(KK)=BK 00003520 K1=KK 00003530 K2=KK+KSPNN 00003540 J=1 00003550 660 K1=K1+KSPAN 00003560 K2=K2-KSPAN 00003570 JJ=J 00003580 AK=AA 00003590 BK=BB 00003600 AJ=0.0 00003610 BJ=0.0 00003620 K=1 00003630 670 K=K+1 00003640 AK=AT(K)*CK(JJ)+AK 00003650 BK=BT(K)*CK(JJ)+BK 00003660 K=K+1 00003670 AJ=AT(K)*SK(JJ)+AJ 00003680 BJ=BT(K)*SK(JJ)+BJ 00003690 JJ=JJ+J 00003700 IF(JJ .GT. JF) JJ=JJ-JF 00003710 IF(K .LT. JF) GO TO 670 00003720 K=JF-J 00003730 A(K1)=AK-BJ 00003740 B(K1)=BK+AJ 00003750 A(K2)=AK+BJ 00003760 B(K2)=BK-AJ 00003770 J=J+1 00003780 IF(J .LT. K) GO TO 660 00003790 KK=KK+KSPNN 00003800 IF(KK .LE. NN) GO TO 640 00003810 KK=KK-NN 00003820 IF(KK .LE. KSPAN) GO TO 640 00003830 C MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4) 00003840 700 IF(I .EQ. M) GO TO 800 00003850 KK=JC+1 00003860 710 C2=1.0-CD 00003870 S1=SD 00003880 720 C1=C2 00003890 S2=S1 00003900 KK=KK+KSPAN 00003910 730 AK=A(KK) 00003920 A(KK)=C2*AK-S2*B(KK) 00003930 B(KK)=S2*AK+C2*B(KK) 00003940 KK=KK+KSPNN 00003950 IF(KK .LE. NT) GO TO 730 00003960 AK=S1*S2 00003970 S2=S1*C2+C1*S2 00003980 C2=C1*C2-AK 00003990 KK=KK-NT+KSPAN 00004000 IF(KK .LE. KSPNN) GO TO 730 00004010 C2=C1-(CD*C1+SD*S1) 00004020 S1=S1+(SD*C1-CD*S1) 00004030 C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION 00004040 C ERROR. IF ROUNDED ARITHMETIC IS USED, THE MAY 00004050 C BE DELETED. 00004060 C1=0.5/(C2**2+S1**2)+0.5 00004070 S1=C1*S1 00004080 C2=C1*C2 00004090 KK=KK-KSPNN+JC 00004100 IF(KK .LE. KSPAN) GO TO 720 00004110 KK=KK-KSPAN+JC+INC 00004120 IF(KK .LE. JC+JC) GO TO 710 00004130 GO TO 100 00004140 C PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES 00004150 C PERMUTATION FOR SQUARE FACTORS OF N 00004160 800 NP(1)=KS 00004170 IF(KT .EQ. 0) GO TO 890 00004180 K=KT+KT+1 00004190 IF(M .LT. K) K=K-1 00004200 J=1 00004210 NP(K+1)=JC 00004220 810 NP(J+1)=NP(J)/NFAC(J) 00004230 NP(K)=NP(K+1)*NFAC(J) 00004240 J=J+1 00004250 K=K-1 00004260 IF(J .LT. K) GO TO 810 00004270 K3=NP(K+1) 00004280 KSPAN=NP(2) 00004290 KK=JC+1 00004300 K2=KSPAN+1 00004310 J=1 00004320 IF(N .NE. NTOT) GO TO 850 00004330 C PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE) 00004340 820 AK=A(KK) 00004350 A(KK)=A(K2) 00004360 A(K2)=AK 00004370 BK=B(KK) 00004380 B(KK)=B(K2) 00004390 B(K2)=BK 00004400 KK=KK+INC 00004410 K2=KSPAN+K2 00004420 IF(K2 .LT. KS) GO TO 820 00004430 830 K2=K2-NP(J) 00004440 J=J+1 00004450 K2=NP(J+1)+K2 00004460 IF(K2 .GT. NP(J)) GO TO 830 00004470 J=1 00004480 840 IF(KK .LT. K2) GO TO 820 00004490 KK=KK+INC 00004500 K2=KSPAN+K2 00004510 IF(K2 .LT. KS) GO TO 840 00004520 IF(KK .LT. KS) GO TO 830 00004530 JC=K3 00004540 GO TO 890 00004550 C PERMUTATION FOR MULTIVARIATE TRANSFORM 00004560 850 K=KK+JC 00004570 860 AK=A(KK) 00004580 A(KK)=A(K2) 00004590 A(K2)=AK 00004600 BK=B(KK) 00004610 B(KK)=B(K2) 00004620 B(K2)=BK 00004630 KK=KK+INC 00004640 K2=K2+INC 00004650 IF(KK .LT. K) GO TO 860 00004660 KK=KK+KS-JC 00004670 K2=K2+KS-JC 00004680 IF(KK .LT. NT) GO TO 850 00004690 K2=K2-NT+KSPAN 00004700 KK=KK-NT+JC 00004710 IF(K2 .LT. KS) GO TO 850 00004720 870 K2=K2-NP(J) 00004730 J=J+1 00004740 K2=NP(J+1)+K2 00004750 IF(K2 .GT. NP(J)) GO TO 870 00004760 J=1 00004770 880 IF(KK .LT. K2) GO TO 850 00004780 KK=KK+JC 00004790 K2=KSPAN+K2 00004800 IF(K2 .LT. KS) GO TO 880 00004810 IF(KK .LT. KS) GO TO 870 00004820 JC=K3 00004830 C CURRIE CODE FOR ISN=-2 (INVERSE NORMALIZED) 00004840 C 890 IF(2*KT+1.GE.M) RETURN 00004850 890 IF(2*KT+1.GE.M) GO TO 952 00004860 C CURRIE CODE FOR ISN=-2 (INVERSE NORMALIZED) 00004870 KSPNN=NP(KT+1) 00004880 C PERMUTATION FOR SQUARE-FREE FACTORS OF N 00004890 J=M-KT 00004900 NFAC(J+1)=1 00004910 900 NFAC(J)=NFAC(J)*NFAC(J+1) 00004920 J=J-1 00004930 IF(J .NE. KT) GO TO 900 00004940 KT=KT+1 00004950 NN=NFAC(KT)-1 00004960 IF(NN .GT. MAXP) GO TO 998 00004970 JJ=0 00004980 J=0 00004990 GO TO 906 00005000 902 JJ=JJ-K2 00005010 K2=KK 00005020 K=K+1 00005030 KK=NFAC(K) 00005040 904 JJ=KK+JJ 00005050 IF(JJ .GE. K2) GO TO 902 00005060 NP(J)=JJ 00005070 906 K2=NFAC(KT) 00005080 K=KT+1 00005090 KK=NFAC(K) 00005100 00005110 J=J+1 00005120 IF(J .LE. NN) GO TO 904 00005130 C DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1 00005140 J=0 00005150 GO TO 914 00005160 910 K=KK 00005170 KK=NP(K) 00005180 NP(K)=-KK 00005190 IF(KK .NE. J) GO TO 910 00005200 K3=KK 00005210 914 J=J+1 00005220 KK=NP(J) 00005230 IF(KK .LT. 0) GO TO 914 00005240 IF(KK .NE. J) GO TO 910 00005250 NP(J)=-J 00005260 IF(J .NE. NN) GO TO 914 00005270 MAXF=INC*MAXF 00005280 C RECORDER A AND B, FOLLOWING THE PERMUTATION CYCLES 00005290 GO TO 950 00005300 924 J=J-1 00005310 IF(NP(J) .LT. 0) GO TO 924 00005320 JJ=JC 00005330 926 KSPAN=JJ 00005340 IF(JJ .GT. MAXF) KSPAN=MAXF 00005350 JJ=JJ-KSPAN 00005360 K=NP(J) 00005370 KK=JC*K+II+JJ 00005380 K1=KK+KSPAN 00005390 K2=0 00005400 928 K2=K2+1 00005410 AT (K2) = A(K1) 00005420 BT(K2)=B(K1) 00005430 K1=K1-INC 00005440 IF(K1 .NE. KK) GO TO 928 00005450 932 K1=KK+KSPAN 00005460 K2=K1-JC*(K+NP(K)) 00005470 K=-NP(K) 00005480 936 A(K1)=A(K2) 00005490 B(K1)=B(K2) 00005500 K1=K1-INC 00005510 K2=K2-INC 00005520 IF(K1 .NE. KK) GO TO 936 00005530 KK=K2 00005540 IF(K .NE. J) GO TO 932 00005550 K1=KK+KSPAN 00005560 K2=0 00005570 940 K2=K2+1 00005580 A(K1)=AT(K2) 00005590 B(K1)=BT(K2) 00005600 K1=K1-INC 00005610 IF(K1 .NE. KK) GO TO 940 00005620 IF(JJ .NE. 0) GO TO 926 00005630 IF(J .NE. 1) GO TO 924 00005640 950 J=K3+1 00005650 NT=NT-KSPNN 00005660 II=NT-INC+1 00005670 IF(NT .GE. 0) GO TO 924 00005680 C CURRIE CODE FOR ISN=-2 (INVERSE NORMALIZED) 00005690 952 IF(ISN.GT.0) RETURN 00005700 IF(ISN.EQ.-1) RETURN 00005710 DO 951 I=1,N 00005720 A(I) = A(I)/FLOAT(N) 00005730 B(I) = -B(I)/FLOAT(N) 00005740 951 CONTINUE 00005750 C CURRIE CODE FOR ISN=-2 (INVERSE NORMALIZED) 00005760 RETURN 00005770 C ERROR FINISH, INSUFFICIENT ARRAY STORAGE 00005780 998 ISN=0 00005790 WRITE(5,999) 00005800 999 FORMAT(44H0ARRAY BOUNDS EXCEEDED WITHIN SUBROUTINE FFT) 00005810 STOP 00005820 END 00005830