SUBROUTINE GSFT(T,F,NPTS,NSPEC,AMP,PHI) C C This subroutine does a Fourier series transform on C the set of points (T,F). C C Input: C C T Independent variable array. C F Array of points for which the transform is to be taken. C NPTS Number of points. C C Output: C C NSPEC Number of harmonics to calculate. C AMP Array of amplitudes. C PHI Array of phase angles (radians). C C Note: If An are the coefficients of the cosine terms and C Bn are the coefficients of the sine terms then, AMP(n) = C SQRT( An*An + Bn*Bn ) & PHI(n) = ATAN2( An, Bn ). C DIMENSION T(NPTS), F(NPTS), AMP(NSPEC), PHI(NSPEC) REAL n2PiP DATA PI/3.14159265358979/, ROOTWO/1.41421356237310/ C PERIOD = T(NPTS) - T(1) DO 2000 n=1,NSPEC n2PiP = n*2.*Pi/PERIOD An = F(1)*( T(2)-T(1) ) + F(NPTS)*( T(NPTS)-T(NPTS-1) ) Bn = 0. NPTSM1 = NPTS-1 C DO 1000 j=2,NPTSM1 THETA = n2PiP*( T(j)-T(1) ) DELTAT = T(j+1) - T(j-1) An = An + F(j)*COS( THETA )*DELTAT 1000 Bn = Bn + F(j)*SIN( THETA )*DELTAT An = An/PERIOD Bn = Bn/PERIOD C AMP(n) = SQRT( An*An + Bn*Bn )/ROOTWO PHI(n) = 0. IF( AMP(n) .EQ. 0. ) GO TO 2000 PHI(n) = ATAN2( An, Bn ) 2000 CONTINUE C RETURN END