home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / DHARM.FOR < prev    next >
Text File  |  1985-11-29  |  13KB  |  522 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DHARM
  5. C
  6. C        PURPOSE
  7. C           PERFORMS DISCRETE COMPLEX FOURIER TRANSFORMS ON A COMPLEX
  8. C           DOUBLE PRECISION,THREE DIMENSIONAL ARRAY
  9. C
  10. C        USAGE
  11. C           CALL DHARM(A,M,INV,S,IFSET,IFERR)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           A     - A DOUBLE PRECISION VECTOR
  15. C                   AS INPUT, A CONTAINS THE COMPLEX, 3-DIMENSIONAL
  16. C                   ARRAY TO BE TRANSFORMED.  THE REAL PART OF
  17. C                   A(I1,I2,I3) IS STORED IN VECTOR FASHION IN A CELL
  18. C                   WITH INDEX 2*(I3*N1*N2 + I2*N1 + I1) + 1 WHERE
  19. C                   NI = 2**M(I), I=1,2,3 AND I1 = 0,1,...,N1-1 ETC.
  20. C                   THE IMAGINARY PART IS IN THE CELL IMMEDIATELY
  21. C                   FOLLOWING.  NOTE THAT THE SUBSCRIPT I1 INCREASES
  22. C                   MOST RAPIDLY AND I3 INCREASES LEAST RAPIDLY.
  23. C                   AS OUTPUT, A CONTAINS THE COMPLEX FOURIER
  24. C                   TRANSFORM.  THE NUMBER OF CORE LOCATIONS OF
  25. C                   ARRAY A IS 2*(N1*N2*N3)
  26. C           M     - A THREE CELL VECTOR WHICH DETERMINES THE SIZES
  27. C                   OF THE 3 DIMENSIONS OF THE ARRAY A.   THE SIZE,
  28. C                   NI, OF THE I DIMENSION OF A IS 2**M(I), I = 1,2,3
  29. C           INV   - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION
  30. C                   OF DIMENSION ONE FOURTH OF THE QUANTITY
  31. C                   MAX(N1,N2,N3)
  32. C                   LOCATIONS OF A, VIZ., (1/8)*2*N1*N2*N3
  33. C           S     - A DOUBLE PRECISION VECTOR WORK AREA FOR SINE TABLES
  34. C                   WITH DIMENSION THE SAME AS INV
  35. C           IFSET - AN OPTION PARAMETER WITH THE FOLLOWING SETTINGS
  36. C                      0    SET UP SINE AND INV TABLES ONLY
  37. C                      1    SET UP SINE AND INV TABLES ONLY AND
  38. C                           CALCULATE FOURIER TRANSFORM
  39. C                     -1    SET UP SINE AND INV TABLES ONLY AND
  40. C                           CALCULATE INVERSE FOURIER TRANSFORM (FOR
  41. C                           THE MEANING OF INVERSE SEE THE EQUATIONS
  42. C                           UNDER METHOD BELOW)
  43. C                      2    CALCULATE FOURIER TRANSFORM ONLY (ASSUME
  44. C                           SINE AND INV TABLES EXIST)
  45. C                     -2    CALCULATE INVERSE FOURIER TRANSFORM ONLY
  46. C                           (ASSUME SINE AND INV TABLES EXIST)
  47. C           IFERR - ERROR INDICATOR.   WHEN IFSET IS 0,+1,-1,
  48. C                   IFERR = 1 MEANS THE MAXIMUM M(I) IS GREATER THAN
  49. C                   20, I=1,2,3   WHEN IFSET IS 2,-2 , IFERR = 1
  50. C                   MEANS THAT THE SINE AND INV TABLES ARE NOT LARGE
  51. C                   ENOUGH OR HAVE NOT BEEN COMPUTED .
  52. C                   IF ON RETURN IFERR = 0 THEN NONE OF THE ABOVE
  53. C                   CONDITIONS ARE PRESENT
  54. C
  55. C        REMARKS
  56. C           THIS SUBROUTINE IS TO BE USED FOR COMPLEX, DOUBLE PRECISION,
  57. C           3-DIMENSIONAL ARRAYS IN WHICH EACH DIMENSION IS A POWER OF
  58. C           2. THE MAXIMUM M(I) MUST NOT BE LESS THAN 3 OR GREATER THAN
  59. C           20, I = 1,2,3.
  60. C
  61. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  62. C           NONE
  63. C
  64. C        METHOD
  65. C           FOR IFSET = +1, OR +2, THE FOURIER TRANSFORM OF COMPLEX
  66. C           ARRAY A IS OBTAINED.
  67. C
  68. C                  N1-1   N2-1   N3-1                L1   L2   L3
  69. C     X(J1,J2,J3)=SUM    SUM    SUM    A(K1,K2,K3)*W1  *W2  *W3
  70. C                  K1=0   K2=0   K3=0
  71. C
  72. C                  WHERE WI IS THE N(I) ROOT OF UNITY AND L1=K1*J1,
  73. C                        L2=K2*J2, L3=K3*J3
  74. C
  75. C
  76. C           FOR IFSET = -1, OR -2, THE INVERSE FOURIER TRANSFORM A OF
  77. C           COMPLEX ARRAY X IS OBTAINED.
  78. C
  79. C     A(K1,K2,K3)=
  80. C               1      N1-1   N2-1   N3-1                -L1  -L2  -L3
  81. C           -------- *SUM    SUM    SUM    X(J1,J2,J3)*W1  *W2  *W3
  82. C           N1*N2*N3   J1=0   J2=0   J3=0
  83. C
  84. C
  85. C           SEE J.W. COOLEY AND J.W. TUKEY, 'AN ALGORITHM FOR THE
  86. C           MACHINE CALCULATION OF COMPLEX FOURIER SERIES',
  87. C           MATHEMATICS OF COMPUTATIONS, VOL. 19 (APR. 1965), P. 297.
  88. C
  89. C     ..................................................................
  90. C
  91.       SUBROUTINE DHARM(A,M,INV,S,IFSET,IFERR)
  92.       DIMENSION A(1),INV(1),S(1),N(3),M(3),NP(3),W(2),W2(2),W3(2)
  93.       DOUBLE PRECISION A,R,W3,AWI,THETA,ROOT2,S,T,W,W2,FN,AWR
  94.       EQUIVALENCE (N1,N(1)),(N2,N(2)),(N3,N(3))
  95.    10 IF( IABS(IFSET) - 1) 900,900,12
  96.    12 MTT=MAX0(M(1),M(2),M(3)) -2
  97.       ROOT2=DSQRT(2.0D0)
  98.       IF (MTT-MT ) 14,14,13
  99.    13 IFERR=1
  100.       RETURN
  101.    14 IFERR=0
  102.       M1=M(1)
  103.       M2=M(2)
  104.       M3=M(3)
  105.       N1=2**M1
  106.       N2=2**M2
  107.       N3=2**M3
  108.    16 IF(IFSET) 18,18,20
  109.    18 NX= N1*N2*N3
  110.       FN = NX
  111.       DO 19 I = 1,NX
  112.       A(2*I-1) = A(2*I-1)/FN
  113.    19 A(2*I) = -A(2*I)/FN
  114.    20 NP(1)=N1*2
  115.       NP(2)= NP(1)*N2
  116.       NP(3)=NP(2)*N3
  117.       DO 250 ID=1,3
  118.       IL = NP(3)-NP(ID)
  119.       IL1 = IL+1
  120.       MI = M(ID)
  121.       IF (MI)250,250,30
  122.    30 IDIF=NP(ID)
  123.       KBIT=NP(ID)
  124.       MEV = 2*(MI/2)
  125.       IF (MI - MEV )60,60,40
  126. C
  127. C     M IS ODD. DO L=1 CASE
  128.    40 KBIT=KBIT/2
  129.       KL=KBIT-2
  130.       DO 50 I=1,IL1,IDIF
  131.       KLAST=KL+I
  132.       DO 50 K=I,KLAST,2
  133.       KD=K+KBIT
  134. C
  135. C     DO ONE STEP WITH L=1,J=0
  136. C     A(K)=A(K)+A(KD)
  137. C     A(KD)=A(K)-A(KD)
  138. C
  139.       T=A(KD)
  140.       A(KD)=A(K)-T
  141.       A(K)=A(K)+T
  142.       T=A(KD+1)
  143.       A(KD+1)=A(K+1)-T
  144.    50 A(K+1)=A(K+1)+T
  145.       IF (MI - 1)250,250,52
  146.    52 LFIRST =3
  147. C
  148. C     DEF - JLAST = 2**(L-2) -1
  149.       JLAST=1
  150.       GO TO 70
  151. C
  152. C     M IS EVEN
  153.    60 LFIRST = 2
  154.       JLAST=0
  155.    70 DO 240 L=LFIRST,MI,2
  156.       JJDIF=KBIT
  157.       KBIT=KBIT/4
  158.       KL=KBIT-2
  159. C
  160. C     DO FOR J=0
  161.       DO 80 I=1,IL1,IDIF
  162.       KLAST=I+KL
  163.       DO 80 K=I,KLAST,2
  164.       K1=K+KBIT
  165.       K2=K1+KBIT
  166.       K3=K2+KBIT
  167. C
  168. C     DO TWO STEPS WITH J=0
  169. C     A(K)=A(K)+A(K2)
  170. C     A(K2)=A(K)-A(K2)
  171. C     A(K1)=A(K1)+A(K3)
  172. C     A(K3)=A(K1)-A(K3)
  173. C
  174. C     A(K)=A(K)+A(K1)
  175. C     A(K1)=A(K)-A(K1)
  176. C     A(K2)=A(K2)+A(K3)*I
  177. C     A(K3)=A(K2)-A(K3)*I
  178. C
  179.       T=A(K2)
  180.       A(K2)=A(K)-T
  181.       A(K)=A(K)+T
  182.       T=A(K2+1)
  183.       A(K2+1)=A(K+1)-T
  184.       A(K+1)=A(K+1)+T
  185. C
  186.       T=A(K3)
  187.       A(K3)=A(K1)-T
  188.       A(K1)=A(K1)+T
  189.       T=A(K3+1)
  190.       A(K3+1)=A(K1+1)-T
  191.       A(K1+1)=A(K1+1)+T
  192. C
  193.       T=A(K1)
  194.       A(K1)=A(K)-T
  195.       A(K)=A(K)+T
  196.       T=A(K1+1)
  197.       A(K1+1)=A(K+1)-T
  198.       A(K+1)=A(K+1)+T
  199. C
  200.       R=-A(K3+1)
  201.       T = A(K3)
  202.       A(K3)=A(K2)-R
  203.       A(K2)=A(K2)+R
  204.       A(K3+1)=A(K2+1)-T
  205.    80 A(K2+1)=A(K2+1)+T
  206.       IF (JLAST) 235,235,82
  207.    82 JJ=JJDIF   +1
  208. C
  209. C     DO FOR J=1
  210.       ILAST= IL +JJ
  211.       DO 85 I = JJ,ILAST,IDIF
  212.       KLAST = KL+I
  213.       DO 85 K=I,KLAST,2
  214.       K1 = K+KBIT
  215.       K2 = K1+KBIT
  216.       K3 = K2+KBIT
  217. C
  218. C     LETTING W=(1+I)/ROOT2,W3=(-1+I)/ROOT2,W2=I,
  219. C     A(K)=A(K)+A(K2)*I
  220. C     A(K2)=A(K)-A(K2)*I
  221. C     A(K1)=A(K1)*W+A(K3)*W3
  222. C     A(K3)=A(K1)*W-A(K3)*W3
  223. C
  224. C     A(K)=A(K)+A(K1)
  225. C     A(K1)=A(K)-A(K1)
  226. C     A(K2)=A(K2)+A(K3)*I
  227. C     A(K3)=A(K2)-A(K3)*I
  228. C
  229.       R =-A(K2+1)
  230.       T = A(K2)
  231.       A(K2) = A(K)-R
  232.       A(K) = A(K)+R
  233.       A(K2+1)=A(K+1)-T
  234.       A(K+1)=A(K+1)+T
  235. C
  236.       AWR=A(K1)-A(K1+1)
  237.       AWI = A(K1+1)+A(K1)
  238.       R=-A(K3)-A(K3+1)
  239.       T=A(K3)-A(K3+1)
  240.       A(K3)=(AWR-R)/ROOT2
  241.       A(K3+1)=(AWI-T)/ROOT2
  242.       A(K1)=(AWR+R)/ROOT2
  243.       A(K1+1)=(AWI+T)/ROOT2
  244.       T= A(K1)
  245.       A(K1)=A(K)-T
  246.       A(K)=A(K)+T
  247.       T=A(K1+1)
  248.       A(K1+1)=A(K+1)-T
  249.       A(K+1)=A(K+1)+T
  250.       R=-A(K3+1)
  251.       T=A(K3)
  252.       A(K3)=A(K2)-R
  253.       A(K2)=A(K2)+R
  254.       A(K3+1)=A(K2+1)-T
  255.    85 A(K2+1)=A(K2+1)+T
  256.       IF(JLAST-1) 235,235,90
  257.    90 JJ= JJ + JJDIF
  258. C
  259. C     NOW DO THE REMAINING J'S
  260.       DO 230 J=2,JLAST
  261. C
  262. C     FETCH W'S
  263. C     DEF- W=W**INV(J), W2=W**2, W3=W**3
  264.    96 I=INV(J+1)
  265.    98 IC=NT-I
  266.       W(1)=S(IC)
  267.       W(2)=S(I)
  268.       I2=2*I
  269.       I2C=NT-I2
  270.       IF(I2C)120,110,100
  271. C
  272. C     2*I IS IN FIRST QUADRANT
  273.   100 W2(1)=S(I2C)
  274.       W2(2)=S(I2)
  275.       GO TO 130
  276.   110 W2(1)=0.
  277.       W2(2)=1.
  278.       GO TO 130
  279. C
  280. C     2*I IS IN SECOND QUADRANT
  281.   120 I2CC = I2C+NT
  282.       I2C=-I2C
  283.       W2(1)=-S(I2C)
  284.       W2(2)=S(I2CC)
  285.   130 I3=I+I2
  286.       I3C=NT-I3
  287.       IF(I3C)160,150,140
  288. C
  289. C     I3 IN FIRST QUADRANT
  290.   140 W3(1)=S(I3C)
  291.       W3(2)=S(I3)
  292.       GO TO 200
  293.   150 W3(1)=0.
  294.       W3(2)=1.
  295.       GO TO 200
  296. C
  297.   160 I3CC=I3C+NT
  298.       IF(I3CC)190,180,170
  299. C
  300. C     I3 IN SECOND QUADRANT
  301.   170 I3C=-I3C
  302.       W3(1)=-S(I3C)
  303.       W3(2)=S(I3CC)
  304.       GO TO 200
  305.   180 W3(1)=-1.
  306.       W3(2)=0.
  307.       GO TO 200
  308. C
  309. C     3*I IN THIRD QUADRANT
  310.   190 I3CCC=NT+I3CC
  311.       I3CC = -I3CC
  312.       W3(1)=-S(I3CCC)
  313.       W3(2)=-S(I3CC)
  314.   200 ILAST=IL+JJ
  315.       DO 220 I=JJ,ILAST,IDIF
  316.       KLAST=KL+I
  317.       DO 220 K=I,KLAST,2
  318.       K1=K+KBIT
  319.       K2=K1+KBIT
  320.       K3=K2+KBIT
  321. C
  322. C     DO TWO STEPS WITH J NOT 0
  323. C     A(K)=A(K)+A(K2)*W2
  324. C     A(K2)=A(K)-A(K2)*W2
  325. C     A(K1)=A(K1)*W+A(K3)*W3
  326. C     A(K3)=A(K1)*W-A(K3)*W3
  327. C
  328. C     A(K)=A(K)+A(K1)
  329. C     A(K1)=A(K)-A(K1)
  330. C     A(K2)=A(K2)+A(K3)*I
  331. C     A(K3)=A(K2)-A(K3)*I
  332. C
  333.       R=A(K2)*W2(1)-A(K2+1)*W2(2)
  334.       T=A(K2)*W2(2)+A(K2+1)*W2(1)
  335.       A(K2)=A(K)-R
  336.       A(K)=A(K)+R
  337.       A(K2+1)=A(K+1)-T
  338.       A(K+1)=A(K+1)+T
  339. C
  340.       R=A(K3)*W3(1)-A(K3+1)*W3(2)
  341.       T=A(K3)*W3(2)+A(K3+1)*W3(1)
  342.       AWR=A(K1)*W(1)-A(K1+1)*W(2)
  343.       AWI=A(K1)*W(2)+A(K1+1)*W(1)
  344.       A(K3)=AWR-R
  345.       A(K3+1)=AWI-T
  346.       A(K1)=AWR+R
  347.       A(K1+1)=AWI+T
  348.       T=A(K1)
  349.       A(K1)=A(K)-T
  350.       A(K)=A(K)+T
  351.       T=A(K1+1)
  352.       A(K1+1)=A(K+1)-T
  353.       A(K+1)=A(K+1)+T
  354.       R=-A(K3+1)
  355.       T=A(K3)
  356.       A(K3)=A(K2)-R
  357.       A(K2)=A(K2)+R
  358.       A(K3+1)=A(K2+1)-T
  359.   220 A(K2+1)=A(K2+1)+T
  360. C     END OF I AND K LOOPS
  361. C
  362.   230 JJ=JJDIF+JJ
  363. C     END OF J-LOOP
  364. C
  365.   235 JLAST=4*JLAST+3
  366.   240 CONTINUE
  367. C     END OF  L  LOOP
  368. C
  369.   250 CONTINUE
  370. C     END OF  ID  LOOP
  371. C
  372. C     WE NOW HAVE THE COMPLEX FOURIER SUMS BUT THEIR ADDRESSES ARE
  373. C     BIT-REVERSED.  THE FOLLOWING ROUTINE PUTS THEM IN ORDER
  374.       NTSQ=NT*NT
  375.       M3MT=M3-MT
  376.   350 IF(M3MT) 370,360,360
  377. C
  378. C     M3 GR. OR EQ. MT
  379.   360 IGO3=1
  380.       N3VNT=N3/NT
  381.       MINN3=NT
  382.       GO TO 380
  383. C
  384. C     M3 LESS THAN MT
  385.   370 IGO3=2
  386.       N3VNT=1
  387.       NTVN3=NT/N3
  388.       MINN3=N3
  389.   380 JJD3 = NTSQ/N3
  390.       M2MT=M2-MT
  391.   450 IF (M2MT)470,460,460
  392. C
  393. C     M2 GR. OR EQ. MT
  394.   460 IGO2=1
  395.       N2VNT=N2/NT
  396.       MINN2=NT
  397.       GO TO 480
  398. C
  399. C     M2 LESS THAN MT
  400.   470 IGO2 = 2
  401.       N2VNT=1
  402.       NTVN2=NT/N2
  403.       MINN2=N2
  404.   480 JJD2=NTSQ/N2
  405.       M1MT=M1-MT
  406.   550 IF(M1MT)570,560,560
  407. C
  408. C     M1 GR. OR EQ. MT
  409.   560 IGO1=1
  410.       N1VNT=N1/NT
  411.       MINN1=NT
  412.       GO TO 580
  413. C
  414. C     M1 LESS THAN MT
  415.   570 IGO1=2
  416.       N1VNT=1
  417.       NTVN1=NT/N1
  418.       MINN1=N1
  419.   580 JJD1=NTSQ/N1
  420.   600 JJ3=1
  421.       J=1
  422.       DO 880 JPP3=1,N3VNT
  423.       IPP3=INV(JJ3)
  424.       DO 870 JP3=1,MINN3
  425.       GO TO (610,620),IGO3
  426.   610 IP3=INV(JP3)*N3VNT
  427.       GO TO 630
  428.   620 IP3=INV(JP3)/NTVN3
  429.   630 I3=(IPP3+IP3)*N2
  430.   700 JJ2=1
  431.       DO 870 JPP2=1,N2VNT
  432.       IPP2=INV(JJ2)+I3
  433.       DO 860 JP2=1,MINN2
  434.       GO TO (710,720),IGO2
  435.   710 IP2=INV(JP2)*N2VNT
  436.       GO TO 730
  437.   720 IP2=INV(JP2)/NTVN2
  438.   730 I2=(IPP2+IP2)*N1
  439.   800 JJ1=1
  440.       DO 860 JPP1=1,N1VNT
  441.       IPP1=INV(JJ1)+I2
  442.       DO 850 JP1=1,MINN1
  443.       GO TO (810,820),IGO1
  444.   810 IP1=INV(JP1)*N1VNT
  445.       GO TO 830
  446.   820 IP1=INV(JP1)/NTVN1
  447.   830 I=2*(IPP1+IP1)+1
  448.       IF (J-I) 840,850,850
  449.   840 T=A(I)
  450.       A(I)=A(J)
  451.       A(J)=T
  452.       T=A(I+1)
  453.       A(I+1)=A(J+1)
  454.       A(J+1)=T
  455.   850 J=J+2
  456.   860 JJ1=JJ1+JJD1
  457. C
  458.   870 JJ2=JJ2+JJD2
  459. C     END OF JPP2 AND JP3 LOOPS
  460. C
  461.   880 JJ3 = JJ3+JJD3
  462. C     END OF JPP3 LOOP
  463. C
  464.   890 IF(IFSET)891,895,895
  465.   891 DO 892 I = 1,NX
  466.   892 A(2*I) = -A(2*I)
  467.   895 RETURN
  468. C
  469. C     THE FOLLOWING PROGRAM COMPUTES THE SIN AND INV TABLES.
  470. C
  471.   900 MT=MAX0(M(1),M(2),M(3)) -2
  472.       MT = MAX0(2,MT)
  473.   904 IF (MT-18) 906,906,13
  474.   906 IFERR=0
  475.       NT=2**MT
  476.       NTV2=NT/2
  477. C
  478. C     SET UP SIN TABLE
  479. C     THETA=PIE/2**(L+1) FOR L=1
  480.   910 THETA=.7853981633974483
  481. C
  482. C     JSTEP=2**(MT-L+1) FOR L=1
  483.       JSTEP=NT
  484. C
  485. C     JDIF=2**(MT-L) FOR L=1
  486.       JDIF=NTV2
  487.       S(JDIF)=DSIN(THETA)
  488.       DO 950 L=2,MT
  489.       THETA=THETA/2.0D0
  490.       JSTEP2=JSTEP
  491.       JSTEP=JDIF
  492.       JDIF=JSTEP/2
  493.       S(JDIF)=DSIN(THETA)
  494.       JC1=NT-JDIF
  495.       S(JC1)=DCOS(THETA)
  496.       JLAST=NT-JSTEP2
  497.       IF(JLAST - JSTEP) 950,920,920
  498.   920 DO 940 J=JSTEP,JLAST,JSTEP
  499.       JC=NT-J
  500.       JD=J+JDIF
  501.   940 S(JD)=S(J)*S(JC1)+S(JDIF)*S(JC)
  502.   950 CONTINUE
  503. C
  504. C     SET UP INV(J) TABLE
  505. C
  506.   960 MTLEXP=NTV2
  507. C
  508. C     MTLEXP=2**(MT-L). FOR L=1
  509.       LM1EXP=1
  510. C
  511. C     LM1EXP=2**(L-1). FOR L=1
  512.       INV(1)=0
  513.       DO 980 L=1,MT
  514.       INV(LM1EXP+1) = MTLEXP
  515.       DO 970 J=2,LM1EXP
  516.       JJ=J+LM1EXP
  517.   970 INV(JJ)=INV(J)+MTLEXP
  518.       MTLEXP=MTLEXP/2
  519.   980 LM1EXP=LM1EXP*2
  520.   982 IF(IFSET)12,895,12
  521.       END
  522.