home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / FORIF.FOR < prev    next >
Text File  |  1985-11-29  |  3KB  |  102 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE FORIF
  5. C
  6. C        PURPOSE
  7. C           FOURIER ANALYSIS OF A GIVEN PERIODIC FUNCTION IN THE
  8. C           RANGE 0-2PI
  9. C           COMPUTES THE COEFFICIENTS OF THE DESIRED NUMBER OF TERMS
  10. C           IN THE FOURIER SERIES F(X)=A(0)+SUM(A(K)COS KX+B(K)SIN KX)
  11. C           WHERE K=1,2,...,M TO APPROXIMATE THE COMPUTED VALUES OF A
  12. C           GIVEN FUNCTION SUBPROGRAM
  13. C
  14. C        USAGE
  15. C           CALL FORIF(FUN,N,M,A,B,IER)
  16. C
  17. C        DESCRIPTION OF PARAMETERS
  18. C           FUN-NAME OF FUNCTION SUBPROGRAM TO BE USED FOR COMPUTING
  19. C               DATA POINTS
  20. C           N  -DEFINES THE INTERVAL SUCH THAT 2N+1 POINTS ARE TAKEN
  21. C               OVER THE INTERVAL (0,2PI). THE SPACING IS THUS 2PI/2N+1
  22. C           M  -THE MAXIMUM ORDER OF THE HARMONICS TO BE FITTED
  23. C           A  -RESULTANT VECTOR OF FOURIER COSINE COEFFICIENTS OF
  24. C               LENGTH M+1
  25. C               A SUB 0, A SUB 1,..., A SUB M
  26. C           B  -RESULTANT VECTOR OF FOURIER SINE COEFFICIENTS OF
  27. C               LENGTH M+1
  28. C               B SUB 0, B SUB 1,..., B SUB M
  29. C           IER-RESULTANT ERROR CODE WHERE
  30. C               IER=0  NO ERROR
  31. C               IER=1  N NOT GREATER OR EQUAL TO M
  32. C               IER=2  M LESS THAN 0
  33. C
  34. C        REMARKS
  35. C           M MUST BE GREATER THAN OR EQUAL TO ZERO
  36. C           N MUST BE GREATER THAN OR EQUAL TO M
  37. C           THE FIRST ELEMENT IN VECTOR B IS ZERO IN ALL CASES
  38. C
  39. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  40. C           FUN-NAME OF USER FUNCTION SUBPROGRAM USED FOR COMPUTING
  41. C               DATA POINTS
  42. C           CALLING PROGRAM MUST HAVE FORTRAN EXTERNAL STATEMENT
  43. C           CONTAINING NAMES OF FUNCTION SUBPROGRAMS LISTED IN CALL TO
  44. C           FORIF
  45. C
  46. C        METHOD
  47. C           USES RECURSIVE TECHNIQUE DESCRIBED IN A. RALSTON, H. WILF,
  48. C           'MATHEMATICAL METHODS FOR DIGITAL COMPUTERS', JOHN WILEY
  49. C           AND SONS, NEW YORK, 1960, CHAPTER 24. THE METHOD OF
  50. C           INDEXING THROUGH THE PROCEDURE HAS BEEN MODIFIED TO
  51. C           SIMPLIFY THE COMPUTATION.
  52. C
  53. C     ..................................................................
  54. C
  55.       SUBROUTINE FORIF(FUN,N,M,A,B,IER)
  56.       DIMENSION A(1),B(1)
  57. C
  58. C        CHECK FOR PARAMETER ERRORS
  59. C
  60.       IER=0
  61.    20 IF(M) 30,40,40
  62.    30 IER=2
  63.       RETURN
  64.    40 IF(M-N) 60,60,50
  65.    50 IER=1
  66.       RETURN
  67. C
  68. C        COMPUTE AND PRESET CONSTANTS
  69. C
  70.    60 AN=N
  71.       COEF=2.0/(2.0*AN+1.0)
  72.       CONST=3.141593*COEF
  73.       S1=SIN(CONST)
  74.       C1=COS(CONST)
  75.       C=1.0
  76.       S=0.0
  77.       J=1
  78.       FUNZ=FUN(0.0)
  79.    70 U2=0.0
  80.       U1=0.0
  81.       AI=2*N
  82. C
  83. C        FORM FOURIER COEFFICIENTS RECURSIVELY
  84. C
  85.    75 X=AI*CONST
  86.       U0=FUN(X)+2.0*C*U1-U2
  87.       U2=U1
  88.       U1=U0
  89.       AI=AI-1.0
  90.       IF(AI) 80,80,75
  91.    80 A(J)=COEF*(FUNZ+C*U1-U2)
  92.       B(J)=COEF*S*U1
  93.       IF(J-(M+1)) 90,100,100
  94.    90 Q=C1*C-S1*S
  95.       S=C1*S+S1*C
  96.       C=Q
  97.       J=J+1
  98.       GO TO 70
  99.   100 A(1)=A(1)*0.5
  100.       RETURN
  101.       END
  102.