home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / itrpapsm / drharm.for < prev    next >
Encoding:
Text File  |  1985-11-29  |  4.0 KB  |  128 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DRHARM
  5. C
  6. C        PURPOSE
  7. C           FINDS THE FOURIER COEFFICIENTS OF ONE DIMENSIONAL DOUBLE
  8. C           PRECISION REAL DATA
  9. C
  10. C        USAGE
  11. C           CALL DRHARM(A,M,INV,S,IFERR)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           A     - A DOUBLE PRECISION VECTOR
  15. C                   AS INPUT, CONTAINS ONE DIMENSIONAL REAL DATA. A IS
  16. C                   2*N+4 CORE LOCATIONS, WHERE N = 2**M. 2*N REAL
  17. C                   NUMBERS ARE PUT INTO THE FIRST 2*N CORE LOCATIONS
  18. C                   OF A
  19. C                   AS OUTPUT, A CONTAINS THE FOURIER COEFFICIENTS
  20. C                   A0/2,B0=0,A1,B1,A2,B2,...,AN/2,BN=0 RESPECTIVELY IN
  21. C                   THE FIRST 2N+2 CORE LOCATIONS OF A
  22. C           M     - AN INTEGER WHICH DETERMINES THE SIZE OF THE VECTOR
  23. C                   A. THE SIZE OF A IS 2*(2**M) + 4
  24. C           INV   - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION OF
  25. C                   DIMENSION ONE EIGHTH THE NUMBER OF REAL INPUT, VIZ.,
  26. C                   (1/8)*2*(2**M)
  27. C           S     - A DOUBLE PRECISION VECTOR WORK AREA FOR SINE TABLES
  28. C                   WITH DIMENSION THE SAME AS INV
  29. C           IFERR - A RETURNED VALUE OF 1 MEANS THAT M IS LESS THAN 3 OR
  30. C                   GREATER THAN 20. OTHERWISE IFERR IS SET = 0
  31. C
  32. C        REMARKS
  33. C           THIS SUBROUTINE GIVES THE FOURIER COEFFICIENTS OF 2*(2**M)
  34. C           REAL POINTS. SEE SUBROUTINE DHARM FOR THREE DIMENSIONAL,
  35. C           DOUBLE PRECISION, COMPLEX FOURIER TRANSFORMS.
  36. C
  37. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  38. C           DHARM
  39. C
  40. C        METHOD
  41. C           THE FOURIER COEFFICIENTS A0,B0=0,A1,B1,...,AN,BN=0 ARE
  42. C           OBTAINED FOR INPUT XJ, J=0,1,2,...,2N-1 FOR THE FOLLOWING
  43. C           EQUATION (PI = 3.14159...)
  44. C
  45. C                 N-1                                               J
  46. C     XJ=(1/2)A0+SUM (AK*COS(PI*J*K/N)+BK*SIN(PI*J*K/N))+(1/2)AN(-1)
  47. C                 K=1
  48. C
  49. C           SEE REFERENCE UNDER SUBROUTINE DHARM
  50. C
  51. C     ..................................................................
  52. C
  53.       SUBROUTINE DRHARM(A,M,INV,S,IFERR)
  54.       DIMENSION A(1),L(3),INV(1),S(1)
  55.       DOUBLE PRECISION A,SI,AP1IM,FN,CO,CIRE,AP2IM,S,SS,DEL,CIIM,AP1RE,
  56.      1 CNIRE,SC,SIS,AP2RE,CNIIM
  57.       IFSET=1
  58.       L(1)=M
  59.       L(2)=0
  60.       L(3)=0
  61.       NTOT=2**M
  62.       NTOT2 = 2*NTOT
  63.       FN = NTOT
  64.       DO   3 I = 2,NTOT2,2
  65.    3  A(I) = -A(I)
  66.       DO   6 I = 1,NTOT2
  67.    6  A(I) = A(I)/FN
  68.       CALL DHARM(A,L,INV,S,IFSET,IFERR)
  69. C
  70. C     MOVE LAST HALF OF A(J)S DOWN ONE SLOT AND ADD A(N) AT BOTTOM TO
  71. C     GIVE ARRAY FOR A1PRIME AND A2PRIME CALCULATION
  72. C
  73.    21 DO  52 I=1,NTOT,2
  74.       J0=NTOT2+2-I
  75.       A(J0)=A(J0-2)
  76.    52 A(J0+1)=A(J0-1)
  77.       A(NTOT2+3)=A(1)
  78.       A(NTOT2+4)=A(2)
  79. C
  80. C     CALCULATE A1PRIMES AND STORE IN FIRST N SLOTS
  81. C     CALCULATE A2PRIMES AND STORE IN SECOND N SLOTS IN REVERSE ORDER
  82.       K0=NTOT+1
  83.       DO 104 I=1,K0,2
  84.       K1=NTOT2-I+4
  85.       AP1RE=.5*(A(I)+A(K1))
  86.       AP2RE=-.5*(A(I+1)+A(K1+1))
  87.       AP1IM=.5*(-A(I+1)+A(K1+1))
  88.       AP2IM=-.5*(A(I)-A(K1))
  89.       A(I)=AP1RE
  90.       A(I+1)=AP1IM
  91.       A(K1)=AP2RE
  92.   104 A(K1+1)=AP2IM
  93.       NTO = NTOT/2
  94.   110 NT=NTO+1
  95.       DEL=3.141592653589793/DFLOAT(NTOT)
  96.       SS=DSIN(DEL)
  97.       SC=DCOS(DEL)
  98.       SI=0.0
  99.       CO=1.0
  100. C
  101. C     COMPUTE C(J)S FOR J=0 THRU J=N
  102.   114 DO 116 I=1,NT
  103.       K6=NTOT2-2*I+5
  104.       AP2RE=A(K6)*CO+A(K6+1)*SI
  105.       AP2IM=-A(K6)*SI+A(K6+1)*CO
  106.       CIRE=.5*(A(2*I-1)+AP2RE)
  107.       CIIM=.5*(A(2*I)+AP2IM)
  108.       CNIRE=.5*(A(2*I-1)-AP2RE)
  109.       CNIIM=.5*(A(2*I)-AP2IM)
  110.       A(2*I-1)=CIRE
  111.       A(2*I)=CIIM
  112.       A(K6)=CNIRE
  113.       A(K6+1)=-CNIIM
  114.       SIS=SI
  115.       SI=SI*SC+CO*SS
  116.   116 CO=CO*SC-SIS*SS
  117. C
  118. C     SHIFT C(J)S FOR J=N/2+1 TO J=N UP ONE SLOT
  119.       DO 117 I=1,NTOT,2
  120.       K8=NTOT+4+I
  121.       A(K8-2)=A(K8)
  122.   117 A(K8-1)=A(K8+1)
  123.       DO 500 I=3,NTOT2,2
  124.       A(I) = 2. * A(I)
  125.   500 A(I + 1) = -2. * A(I + 1)
  126.       RETURN
  127.       END
  128.