- #1
Galizius
- 14
- 0
I was told to do a Fourier transform of function by using a Filon's method. I have found the code but I don't know how to include any function to the subroutine. I would be grateful for any example of how to use this code.
Fortran:
SUBROUTINE FILONC ( DT, DOM, NMAX, C, CHAT )
C *******************************************************************
C ** CALCULATES THE FOURIER COSINE TRANSFORM BY FILON'S METHOD **
C ** **
C ** A CORRELATION FUNCTION, C(T), IN THE TIME DOMAIN, IS **
C ** TRANSFORMED TO A SPECTRUM CHAT(OMEGA) IN THE FREQUENCY DOMAIN.**
C ** **
C ** REFERENCE: **
C ** **
C ** FILON, PROC ROY SOC EDIN, 49 38, 1928. **
C ** **
C ** PRINCIPAL VARIABLES: **
C ** **
C ** REAL C(NMAX) THE CORRELATION FUNCTION. **
C ** REAL CHAT(NMAX) THE 1-D COSINE TRANSFORM. **
C ** REAL DT TIME INTERVAL BETWEEN POINTS IN C. **
C ** REAL DOM FREQUENCY INTERVAL FOR CHAT. **
C ** INTEGER NMAX NO. OF INTERVALS ON THE TIME AXIS **
C ** REAL OMEGA THE FREQUENCY **
C ** REAL TMAX MAXIMUM TIME IN CORRL. FUNCTION **
C ** REAL ALPHA, BETA, GAMMA FILON PARAMETERS **
C ** INTEGER TAU TIME INDEX **
C ** INTEGER NU FREQUENCY INDEX **
C ** **
C ** USAGE: **
C ** **
C ** THE ROUTINE REQUIRES THAT THE NUMBER OF INTERVALS, NMAX, IS **
C ** EVEN AND CHECKS FOR THIS CONDITION. THE FIRST VALUE OF C(T) **
C ** IS AT T=0. THE MAXIMUM TIME FOR THE CORRELATION FUNCTION IS **
C ** TMAX=DT*NMAX. FOR AN ACCURATE TRANSFORM C(TMAX)=0. **
C *******************************************************************
INTEGER NMAX
REAL DT, DOM, C(0:NMAX), CHAT(0:NMAX)
REAL TMAX, OMEGA, THETA, SINTH, COSTH, CE, CO
REAL SINSQ, COSSQ, THSQ, THCUB, ALPHA, BETA, GAMMA
INTEGER TAU, NU
C *******************************************************************
C ** CHECKS NMAX IS EVEN **
IF ( MOD ( NMAX, 2 ) .NE. 0 ) THEN
STOP ' NMAX SHOULD BE EVEN '
ENDIF
TMAX = REAL ( NMAX ) * DT
C ** LOOP OVER OMEGA **
DO 30 NU = 0, NMAX
OMEGA = REAL ( NU ) * DOM
THETA = OMEGA * DT
C ** CALCULATE THE FILON PARAMETERS **
SINTH = SIN ( THETA )
COSTH = COS ( THETA )
SINSQ = SINTH * SINTH
COSSQ = COSTH * COSTH
THSQ = THETA * THETA
THCUB = THSQ * THETA
IF ( THETA. EQ. 0.0 ) THEN
ALPHA = 0.0
BETA = 2.0 / 3.0
GAMMA = 4.0 / 3.0
ELSE
ALPHA = ( 1.0 / THCUB )
: * ( THSQ + THETA * SINTH * COSTH - 2.0 * SINSQ )
BETA = ( 2.0 / THCUB )
: * ( THETA * ( 1.0 + COSSQ ) -2.0 * SINTH * COSTH )
GAMMA = ( 4.0 / THCUB ) * ( SINTH - THETA * COSTH )
ENDIF
C ** DO THE SUM OVER THE EVEN ORDINATES **
CE = 0.0
DO 10 TAU = 0, NMAX, 2
CE = CE + C(TAU) * COS ( THETA * REAL ( TAU ) )
10 CONTINUE
C ** SUBTRACT HALF THE FIRST AND LAST TERMS **
CE = CE - 0.5 * ( C(0) + C(NMAX) * COS ( OMEGA * TMAX ) )
C ** DO THE SUM OVER THE ODD ORDINATES **
CO = 0.0
DO 20 TAU = 1, NMAX - 1, 2
CO = CO + C(TAU) * COS ( THETA * REAL ( TAU ) )
20 CONTINUE
C ** FACTOR OF TWO FOR THE REAL COSINE TRANSFORM **
CHAT(NU) = 2.0 * ( ALPHA * C(NMAX) * SIN ( OMEGA * TMAX )
: + BETA * CE + GAMMA * CO ) * DT
30 CONTINUE
RETURN
END
Last edited by a moderator: