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