C DHAR 10
C ..................................................................DHAR 20
C DHAR 30
C SUBROUTINE DHARM DHAR 40
C DHAR 50
C PURPOSE DHAR 60
C PERFORMS DISCRETE COMPLEX FOURIER TRANSFORMS ON A COMPLEX DHAR 70
C DOUBLE PRECISION,THREE DIMENSIONAL ARRAY DHAR 80
C DHAR 90
C USAGE DHAR 100
C CALL DHARM(A,M,INV,S,IFSET,IFERR) DHAR 110
C DHAR 120
C DESCRIPTION OF PARAMETERS DHAR 130
C A - A DOUBLE PRECISION VECTOR DHAR 140
C AS INPUT, A CONTAINS THE COMPLEX, 3-DIMENSIONAL DHAR 150
C ARRAY TO BE TRANSFORMED. THE REAL PART OF DHAR 160
C A(I1,I2,I3) IS STORED IN VECTOR FASHION IN A CELL DHAR 170
C WITH INDEX 2*(I3*N1*N2 + I2*N1 + I1) + 1 WHERE DHAR 180
C NI = 2**M(I), I=1,2,3 AND I1 = 0,1,...,N1-1 ETC. DHAR 190
C THE IMAGINARY PART IS IN THE CELL IMMEDIATELY DHAR 200
C FOLLOWING. NOTE THAT THE SUBSCRIPT I1 INCREASES DHAR 210
C MOST RAPIDLY AND I3 INCREASES LEAST RAPIDLY. DHAR 220
C AS OUTPUT, A CONTAINS THE COMPLEX FOURIER DHAR 230
C TRANSFORM. THE NUMBER OF CORE LOCATIONS OF DHAR 240
C ARRAY A IS 2*(N1*N2*N3) DHAR 250
C M - A THREE CELL VECTOR WHICH DETERMINES THE SIZES DHAR 260
C OF THE 3 DIMENSIONS OF THE ARRAY A. THE SIZE, DHAR 270
C NI, OF THE I DIMENSION OF A IS 2**M(I), I = 1,2,3 DHAR 280
C INV - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION DHAR 290
C OF DIMENSION ONE FOURTH OF THE QUANTITY DHAR 300
C MAX(N1,N2,N3) DHAR 310
C LOCATIONS OF A, VIZ., (1/8)*2*N1*N2*N3 DHAR 310
C S - A DOUBLE PRECISION VECTOR WORK AREA FOR SINE TABLES DHAR 320
C WITH DIMENSION THE SAME AS INV DHAR 330
C IFSET - AN OPTION PARAMETER WITH THE FOLLOWING SETTINGS DHAR 340
C 0 SET UP SINE AND INV TABLES ONLY DHAR 350
C 1 SET UP SINE AND INV TABLES ONLY AND DHAR 360
C CALCULATE FOURIER TRANSFORM DHAR 370
C -1 SET UP SINE AND INV TABLES ONLY AND DHAR 380
C CALCULATE INVERSE FOURIER TRANSFORM (FOR DHAR 390
C THE MEANING OF INVERSE SEE THE EQUATIONS DHAR 400
C UNDER METHOD BELOW) DHAR 410
C 2 CALCULATE FOURIER TRANSFORM ONLY (ASSUME DHAR 420
C SINE AND INV TABLES EXIST) DHAR 430
C -2 CALCULATE INVERSE FOURIER TRANSFORM ONLY DHAR 440
C (ASSUME SINE AND INV TABLES EXIST) DHAR 450
C IFERR - ERROR INDICATOR. WHEN IFSET IS 0,+1,-1, DHAR 460
C IFERR = 1 MEANS THE MAXIMUM M(I) IS GREATER THAN DHAR 470
C 20, I=1,2,3 WHEN IFSET IS 2,-2 , IFERR = 1 DHAR 480
C MEANS THAT THE SINE AND INV TABLES ARE NOT LARGE DHAR 490
C ENOUGH OR HAVE NOT BEEN COMPUTED . DHAR 500
C IF ON RETURN IFERR = 0 THEN NONE OF THE ABOVE DHAR 510
C CONDITIONS ARE PRESENT DHAR 520
C DHAR 530
C REMARKS DHAR 540
C THIS SUBROUTINE IS TO BE USED FOR COMPLEX, DOUBLE PRECISION,DHAR 550
C 3-DIMENSIONAL ARRAYS IN WHICH EACH DIMENSION IS A POWER OF DHAR 560
C 2. THE MAXIMUM M(I) MUST NOT BE LESS THAN 3 OR GREATER THAN DHAR 570
C 20, I = 1,2,3. DHAR 580
C DHAR 590
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DHAR 600
C NONE DHAR 610
C DHAR 620
C METHOD DHAR 630
C FOR IFSET = +1, OR +2, THE FOURIER TRANSFORM OF COMPLEX DHAR 640
C ARRAY A IS OBTAINED. DHAR 650
C DHAR 660
C N1-1 N2-1 N3-1 L1 L2 L3 DHAR 670
C X(J1,J2,J3)=SUM SUM SUM A(K1,K2,K3)*W1 *W2 *W3 DHAR 680
C K1=0 K2=0 K3=0 DHAR 690
C DHAR 700
C WHERE WI IS THE N(I) ROOT OF UNITY AND L1=K1*J1, DHAR 710
C L2=K2*J2, L3=K3*J3 DHAR 720
C DHAR 730
C DHAR 740
C FOR IFSET = -1, OR -2, THE INVERSE FOURIER TRANSFORM A OF DHAR 750
C COMPLEX ARRAY X IS OBTAINED. DHAR 760
C DHAR 770
C A(K1,K2,K3)= DHAR 780
C 1 N1-1 N2-1 N3-1 -L1 -L2 -L3 DHAR 790
C -------- *SUM SUM SUM X(J1,J2,J3)*W1 *W2 *W3 DHAR 800
C N1*N2*N3 J1=0 J2=0 J3=0 DHAR 810
C DHAR 820
C DHAR 830
C SEE J.W. COOLEY AND J.W. TUKEY, 'AN ALGORITHM FOR THE DHAR 840
C MACHINE CALCULATION OF COMPLEX FOURIER SERIES', DHAR 850
C MATHEMATICS OF COMPUTATIONS, VOL. 19 (APR. 1965), P. 297. DHAR 860
C DHAR 870
C ..................................................................DHAR 880
C DHAR 890
SUBROUTINE DHARM(A,M,INV,S,IFSET,IFERR) DHAR 900
DIMENSION A(1),INV(1),S(1),N(3),M(3),NP(3),W(2),W2(2),W3(2) DHAR 910
DOUBLE PRECISION A,R,W3,AWI,THETA,ROOT2,S,T,W,W2,FN,AWR DHAR 920
EQUIVALENCE (N1,N(1)),(N2,N(2)),(N3,N(3)) DHAR 930
10 IF( IABS(IFSET) - 1) 900,900,12 DHAR 940
12 MTT=MAX0(M(1),M(2),M(3)) -2 DHAR 950
ROOT2=DSQRT(2.0D0) DHAR 960
IF (MTT-MT ) 14,14,13 DHAR 970
13 IFERR=1 DHAR 980
RETURN DHAR 990
14 IFERR=0 DHAR1000
M1=M(1) DHAR1010
M2=M(2) DHAR1020
M3=M(3) DHAR1030
N1=2**M1 DHAR1040
N2=2**M2 DHAR1050
N3=2**M3 DHAR1060
16 IF(IFSET) 18,18,20 DHAR1070
18 NX= N1*N2*N3 DHAR1080
FN = NX DHAR1090
DO 19 I = 1,NX DHAR1100
A(2*I-1) = A(2*I-1)/FN DHAR1110
19 A(2*I) = -A(2*I)/FN DHAR1120
20 NP(1)=N1*2 DHAR1130
NP(2)= NP(1)*N2 DHAR1140
NP(3)=NP(2)*N3 DHAR1150
DO 250 ID=1,3 DHAR1160
IL = NP(3)-NP(ID) DHAR1170
IL1 = IL+1 DHAR1180
MI = M(ID) DHAR1190
IF (MI)250,250,30 DHAR1200
30 IDIF=NP(ID) DHAR1210
KBIT=NP(ID) DHAR1220
MEV = 2*(MI/2) DHAR1230
IF (MI - MEV )60,60,40 DHAR1240
C DHAR1250
C M IS ODD. DO L=1 CASE DHAR1260
40 KBIT=KBIT/2 DHAR1270
KL=KBIT-2 DHAR1280
DO 50 I=1,IL1,IDIF DHAR1290
KLAST=KL+I DHAR1300
DO 50 K=I,KLAST,2 DHAR1310
KD=K+KBIT DHAR1320
C DHAR1330
C DO ONE STEP WITH L=1,J=0 DHAR1340
C A(K)=A(K)+A(KD) DHAR1350
C A(KD)=A(K)-A(KD) DHAR1360
C DHAR1370
T=A(KD) DHAR1380
A(KD)=A(K)-T DHAR1390
A(K)=A(K)+T DHAR1400
T=A(KD+1) DHAR1410
A(KD+1)=A(K+1)-T DHAR1420
50 A(K+1)=A(K+1)+T DHAR1430
IF (MI - 1)250,250,52 DHAR1440
52 LFIRST =3 DHAR1450
C DHAR1460
C DEF - JLAST = 2**(L-2) -1 DHAR1470
JLAST=1 DHAR1480
GO TO 70 DHAR1490
C DHAR1500
C M IS EVEN DHAR1510
60 LFIRST = 2 DHAR1520
JLAST=0 DHAR1530
70 DO 240 L=LFIRST,MI,2 DHAR1540
JJDIF=KBIT DHAR1550
KBIT=KBIT/4 DHAR1560
KL=KBIT-2 DHAR1570
C DHAR1580
C DO FOR J=0 DHAR1590
DO 80 I=1,IL1,IDIF DHAR1600
KLAST=I+KL DHAR1610
DO 80 K=I,KLAST,2 DHAR1620
K1=K+KBIT DHAR1630
K2=K1+KBIT DHAR1640
K3=K2+KBIT DHAR1650
C DHAR1660
C DO TWO STEPS WITH J=0 DHAR1670
C A(K)=A(K)+A(K2) DHAR1680
C A(K2)=A(K)-A(K2) DHAR1690
C A(K1)=A(K1)+A(K3) DHAR1700
C A(K3)=A(K1)-A(K3) DHAR1710
C DHAR1720
C A(K)=A(K)+A(K1) DHAR1730
C A(K1)=A(K)-A(K1) DHAR1740
C A(K2)=A(K2)+A(K3)*I DHAR1750
C A(K3)=A(K2)-A(K3)*I DHAR1760
C DHAR1770
T=A(K2) DHAR1780
A(K2)=A(K)-T DHAR1790
A(K)=A(K)+T DHAR1800
T=A(K2+1) DHAR1810
A(K2+1)=A(K+1)-T DHAR1820
A(K+1)=A(K+1)+T DHAR1830
C DHAR1840
T=A(K3) DHAR1850
A(K3)=A(K1)-T DHAR1860
A(K1)=A(K1)+T DHAR1870
T=A(K3+1) DHAR1880
A(K3+1)=A(K1+1)-T DHAR1890
A(K1+1)=A(K1+1)+T DHAR1900
C DHAR1910
T=A(K1) DHAR1920
A(K1)=A(K)-T DHAR1930
A(K)=A(K)+T DHAR1940
T=A(K1+1) DHAR1950
A(K1+1)=A(K+1)-T DHAR1960
A(K+1)=A(K+1)+T DHAR1970
C DHAR1980
R=-A(K3+1) DHAR1990
T = A(K3) DHAR2000
A(K3)=A(K2)-R DHAR2010
A(K2)=A(K2)+R DHAR2020
A(K3+1)=A(K2+1)-T DHAR2030
80 A(K2+1)=A(K2+1)+T DHAR2040
IF (JLAST) 235,235,82 DHAR2050
82 JJ=JJDIF +1 DHAR2060
C DHAR2070
C DO FOR J=1 DHAR2080
ILAST= IL +JJ DHAR2090
DO 85 I = JJ,ILAST,IDIF DHAR2100
KLAST = KL+I DHAR2110
DO 85 K=I,KLAST,2 DHAR2120
K1 = K+KBIT DHAR2130
K2 = K1+KBIT DHAR2140
K3 = K2+KBIT DHAR2150
C DHAR2160
C LETTING W=(1+I)/ROOT2,W3=(-1+I)/ROOT2,W2=I, DHAR2170
C A(K)=A(K)+A(K2)*I DHAR2180
C A(K2)=A(K)-A(K2)*I DHAR2190
C A(K1)=A(K1)*W+A(K3)*W3 DHAR2200
C A(K3)=A(K1)*W-A(K3)*W3 DHAR2210
C DHAR2220
C A(K)=A(K)+A(K1) DHAR2230
C A(K1)=A(K)-A(K1) DHAR2240
C A(K2)=A(K2)+A(K3)*I DHAR2250
C A(K3)=A(K2)-A(K3)*I DHAR2260
C DHAR2270
R =-A(K2+1) DHAR2280
T = A(K2) DHAR2290
A(K2) = A(K)-R DHAR2300
A(K) = A(K)+R DHAR2310
A(K2+1)=A(K+1)-T DHAR2320
A(K+1)=A(K+1)+T DHAR2330
C DHAR2340
AWR=A(K1)-A(K1+1) DHAR2350
AWI = A(K1+1)+A(K1) DHAR2360
R=-A(K3)-A(K3+1) DHAR2370
T=A(K3)-A(K3+1) DHAR2380
A(K3)=(AWR-R)/ROOT2 DHAR2390
A(K3+1)=(AWI-T)/ROOT2 DHAR2400
A(K1)=(AWR+R)/ROOT2 DHAR2410
A(K1+1)=(AWI+T)/ROOT2 DHAR2420
T= A(K1) DHAR2430
A(K1)=A(K)-T DHAR2440
A(K)=A(K)+T DHAR2450
T=A(K1+1) DHAR2460
A(K1+1)=A(K+1)-T DHAR2470
A(K+1)=A(K+1)+T DHAR2480
R=-A(K3+1) DHAR2490
T=A(K3) DHAR2500
A(K3)=A(K2)-R DHAR2510
A(K2)=A(K2)+R DHAR2520
A(K3+1)=A(K2+1)-T DHAR2530
85 A(K2+1)=A(K2+1)+T DHAR2540
IF(JLAST-1) 235,235,90 DHAR2550
90 JJ= JJ + JJDIF DHAR2560
C DHAR2570
C NOW DO THE REMAINING J'S DHAR2580
DO 230 J=2,JLAST DHAR2590
C DHAR2600
C FETCH W'S DHAR2610
C DEF- W=W**INV(J), W2=W**2, W3=W**3 DHAR2620
96 I=INV(J+1) DHAR2630
98 IC=NT-I DHAR2640
W(1)=S(IC) DHAR2650
W(2)=S(I) DHAR2660
I2=2*I DHAR2670
I2C=NT-I2 DHAR2680
IF(I2C)120,110,100 DHAR2690
C DHAR2700
C 2*I IS IN FIRST QUADRANT DHAR2710
100 W2(1)=S(I2C) DHAR2720
W2(2)=S(I2) DHAR2730
GO TO 130 DHAR2740
110 W2(1)=0. DHAR2750
W2(2)=1. DHAR2760
GO TO 130 DHAR2770
C DHAR2780
C 2*I IS IN SECOND QUADRANT DHAR2790
120 I2CC = I2C+NT DHAR2800
I2C=-I2C DHAR2810
W2(1)=-S(I2C) DHAR2820
W2(2)=S(I2CC) DHAR2830
130 I3=I+I2 DHAR2840
I3C=NT-I3 DHAR2850
IF(I3C)160,150,140 DHAR2860
C DHAR2870
C I3 IN FIRST QUADRANT DHAR2880
140 W3(1)=S(I3C) DHAR2890
W3(2)=S(I3) DHAR2900
GO TO 200 DHAR2910
150 W3(1)=0. DHAR2920
W3(2)=1. DHAR2930
GO TO 200 DHAR2940
C DHAR2950
160 I3CC=I3C+NT DHAR2960
IF(I3CC)190,180,170 DHAR2970
C DHAR2980
C I3 IN SECOND QUADRANT DHAR2990
170 I3C=-I3C DHAR3000
W3(1)=-S(I3C) DHAR3010
W3(2)=S(I3CC) DHAR3020
GO TO 200 DHAR3030
180 W3(1)=-1. DHAR3040
W3(2)=0. DHAR3050
GO TO 200 DHAR3060
C DHAR3070
C 3*I IN THIRD QUADRANT DHAR3080
190 I3CCC=NT+I3CC DHAR3090
I3CC = -I3CC DHAR3100
W3(1)=-S(I3CCC) DHAR3110
W3(2)=-S(I3CC) DHAR3120
200 ILAST=IL+JJ DHAR3130
DO 220 I=JJ,ILAST,IDIF DHAR3140
KLAST=KL+I DHAR3150
DO 220 K=I,KLAST,2 DHAR3160
K1=K+KBIT DHAR3170
K2=K1+KBIT DHAR3180
K3=K2+KBIT DHAR3190
C DHAR3200
C DO TWO STEPS WITH J NOT 0 DHAR3210
C A(K)=A(K)+A(K2)*W2 DHAR3220
C A(K2)=A(K)-A(K2)*W2 DHAR3230
C A(K1)=A(K1)*W+A(K3)*W3 DHAR3240
C A(K3)=A(K1)*W-A(K3)*W3 DHAR3250
C DHAR3260
C A(K)=A(K)+A(K1) DHAR3270
C A(K1)=A(K)-A(K1) DHAR3280
C A(K2)=A(K2)+A(K3)*I DHAR3290
C A(K3)=A(K2)-A(K3)*I DHAR3300
C DHAR3310
R=A(K2)*W2(1)-A(K2+1)*W2(2) DHAR3320
T=A(K2)*W2(2)+A(K2+1)*W2(1) DHAR3330
A(K2)=A(K)-R DHAR3340
A(K)=A(K)+R DHAR3350
A(K2+1)=A(K+1)-T DHAR3360
A(K+1)=A(K+1)+T DHAR3370
C DHAR3380
R=A(K3)*W3(1)-A(K3+1)*W3(2) DHAR3390
T=A(K3)*W3(2)+A(K3+1)*W3(1) DHAR3400
AWR=A(K1)*W(1)-A(K1+1)*W(2) DHAR3410
AWI=A(K1)*W(2)+A(K1+1)*W(1) DHAR3420
A(K3)=AWR-R DHAR3430
A(K3+1)=AWI-T DHAR3440
A(K1)=AWR+R DHAR3450
A(K1+1)=AWI+T DHAR3460
T=A(K1) DHAR3470
A(K1)=A(K)-T DHAR3480
A(K)=A(K)+T DHAR3490
T=A(K1+1) DHAR3500
A(K1+1)=A(K+1)-T DHAR3510
A(K+1)=A(K+1)+T DHAR3520
R=-A(K3+1) DHAR3530
T=A(K3) DHAR3540
A(K3)=A(K2)-R DHAR3550
A(K2)=A(K2)+R DHAR3560
A(K3+1)=A(K2+1)-T DHAR3570
220 A(K2+1)=A(K2+1)+T DHAR3580
C END OF I AND K LOOPS DHAR3590
C DHAR3600
230 JJ=JJDIF+JJ DHAR3610
C END OF J-LOOP DHAR3620
C DHAR3630
235 JLAST=4*JLAST+3 DHAR3640
240 CONTINUE DHAR3650
C END OF L LOOP DHAR3660
C DHAR3670
250 CONTINUE DHAR3680
C END OF ID LOOP DHAR3690
C DHAR3700
C WE NOW HAVE THE COMPLEX FOURIER SUMS BUT THEIR ADDRESSES ARE DHAR3710
C BIT-REVERSED. THE FOLLOWING ROUTINE PUTS THEM IN ORDER DHAR3720
NTSQ=NT*NT DHAR3730
M3MT=M3-MT DHAR3740
350 IF(M3MT) 370,360,360 DHAR3750
C DHAR3760
C M3 GR. OR EQ. MT DHAR3770
360 IGO3=1 DHAR3780
N3VNT=N3/NT DHAR3790
MINN3=NT DHAR3800
GO TO 380 DHAR3810
C DHAR3820
C M3 LESS THAN MT DHAR3830
370 IGO3=2 DHAR3840
N3VNT=1 DHAR3850
NTVN3=NT/N3 DHAR3860
MINN3=N3 DHAR3870
380 JJD3 = NTSQ/N3 DHAR3880
M2MT=M2-MT DHAR3890
450 IF (M2MT)470,460,460 DHAR3900
C DHAR3910
C M2 GR. OR EQ. MT DHAR3920
460 IGO2=1 DHAR3930
N2VNT=N2/NT DHAR3940
MINN2=NT DHAR3950
GO TO 480 DHAR3960
C DHAR3970
C M2 LESS THAN MT DHAR3980
470 IGO2 = 2 DHAR3990
N2VNT=1 DHAR4000
NTVN2=NT/N2 DHAR4010
MINN2=N2 DHAR4020
480 JJD2=NTSQ/N2 DHAR4030
M1MT=M1-MT DHAR4040
550 IF(M1MT)570,560,560 DHAR4050
C DHAR4060
C M1 GR. OR EQ. MT DHAR4070
560 IGO1=1 DHAR4080
N1VNT=N1/NT DHAR4090
MINN1=NT DHAR4100
GO TO 580 DHAR4110
C DHAR4120
C M1 LESS THAN MT DHAR4130
570 IGO1=2 DHAR4140
N1VNT=1 DHAR4150
NTVN1=NT/N1 DHAR4160
MINN1=N1 DHAR4170
580 JJD1=NTSQ/N1 DHAR4180
600 JJ3=1 DHAR4190
J=1 DHAR4200
DO 880 JPP3=1,N3VNT DHAR4210
IPP3=INV(JJ3) DHAR4220
DO 870 JP3=1,MINN3 DHAR4230
GO TO (610,620),IGO3 DHAR4240
610 IP3=INV(JP3)*N3VNT DHAR4250
GO TO 630 DHAR4260
620 IP3=INV(JP3)/NTVN3 DHAR4270
630 I3=(IPP3+IP3)*N2 DHAR4280
700 JJ2=1 DHAR4290
DO 870 JPP2=1,N2VNT DHAR4300
IPP2=INV(JJ2)+I3 DHAR4310
DO 860 JP2=1,MINN2 DHAR4320
GO TO (710,720),IGO2 DHAR4330
710 IP2=INV(JP2)*N2VNT DHAR4340
GO TO 730 DHAR4350
720 IP2=INV(JP2)/NTVN2 DHAR4360
730 I2=(IPP2+IP2)*N1 DHAR4370
800 JJ1=1 DHAR4380
DO 860 JPP1=1,N1VNT DHAR4390
IPP1=INV(JJ1)+I2 DHAR4400
DO 850 JP1=1,MINN1 DHAR4410
GO TO (810,820),IGO1 DHAR4420
810 IP1=INV(JP1)*N1VNT DHAR4430
GO TO 830 DHAR4440
820 IP1=INV(JP1)/NTVN1 DHAR4450
830 I=2*(IPP1+IP1)+1 DHAR4460
IF (J-I) 840,850,850 DHAR4470
840 T=A(I) DHAR4480
A(I)=A(J) DHAR4490
A(J)=T DHAR4500
T=A(I+1) DHAR4510
A(I+1)=A(J+1) DHAR4520
A(J+1)=T DHAR4530
850 J=J+2 DHAR4550
860 JJ1=JJ1+JJD1 DHAR4560
C DHAR4580
870 JJ2=JJ2+JJD2 DHAR4590
C END OF JPP2 AND JP3 LOOPS DHAR4600
C DHAR4610
880 JJ3 = JJ3+JJD3 DHAR4620
C END OF JPP3 LOOP DHAR4630
C DHAR4640
890 IF(IFSET)891,895,895 DHAR4650
891 DO 892 I = 1,NX DHAR4660
892 A(2*I) = -A(2*I) DHAR4670
895 RETURN DHAR4680
C DHAR4690
C THE FOLLOWING PROGRAM COMPUTES THE SIN AND INV TABLES. DHAR4700
C DHAR4710
900 MT=MAX0(M(1),M(2),M(3)) -2 DHAR4720
MT = MAX0(2,MT) DHAR4730
904 IF (MT-18) 906,906,13 DHAR4740
906 IFERR=0 DHAR4770
NT=2**MT DHAR4780
NTV2=NT/2 DHAR4790
C DHAR4800
C SET UP SIN TABLE DHAR4810
C THETA=PIE/2**(L+1) FOR L=1 DHAR4820
910 THETA=.7853981633974483 DHAR4830
C DHAR4840
C JSTEP=2**(MT-L+1) FOR L=1 DHAR4850
JSTEP=NT DHAR4860
C DHAR4870
C JDIF=2**(MT-L) FOR L=1 DHAR4880
JDIF=NTV2 DHAR4890
S(JDIF)=DSIN(THETA) DHAR4900
DO 950 L=2,MT DHAR4910
THETA=THETA/2.0D0 DHAR4920
JSTEP2=JSTEP DHAR4930
JSTEP=JDIF DHAR4940
JDIF=JSTEP/2 DHAR4950
S(JDIF)=DSIN(THETA) DHAR4960
JC1=NT-JDIF DHAR4970
S(JC1)=DCOS(THETA) DHAR4980
JLAST=NT-JSTEP2 DHAR4990
IF(JLAST - JSTEP) 950,920,920 DHAR5000
920 DO 940 J=JSTEP,JLAST,JSTEP DHAR5010
JC=NT-J DHAR5020
JD=J+JDIF DHAR5030
940 S(JD)=S(J)*S(JC1)+S(JDIF)*S(JC) DHAR5040
950 CONTINUE DHAR5050
C DHAR5060
C SET UP INV(J) TABLE DHAR5070
C DHAR5080
960 MTLEXP=NTV2 DHAR5090
C DHAR5100
C MTLEXP=2**(MT-L). FOR L=1 DHAR5110
LM1EXP=1 DHAR5120
C DHAR5130
C LM1EXP=2**(L-1). FOR L=1 DHAR5140
INV(1)=0 DHAR5150
DO 980 L=1,MT DHAR5160
INV(LM1EXP+1) = MTLEXP DHAR5170
DO 970 J=2,LM1EXP DHAR5180
JJ=J+LM1EXP DHAR5190
970 INV(JJ)=INV(J)+MTLEXP DHAR5200
MTLEXP=MTLEXP/2 DHAR5210
980 LM1EXP=LM1EXP*2 DHAR5220
982 IF(IFSET)12,895,12 DHAR5230
END DHAR5240