C ******************************************************************* ** FICHE F.37. ROUTINES TO CALCULATE FOURIER TRANSFORMS. ** C ** THREE SEPARATE ROUTINES FOR DIFFERENT APPLICATIONS. ** C ******************************************************************* 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 SUBROUTINE LADO ( DT, NMAX, C, CHAT ) C ******************************************************************* C ** CALCULATES THE FOURIER COSINE TRANSFORM BY LADO'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 ** LADO, J COMPUT PHYS, 8 417, 1971. ** 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 BETWEEN POINTS IN CHAT.** C ** INTEGER NMAX NO. OF INTERVALS ON THE TIME AXIS ** C ** ** C ** USAGE: ** C ** ** C ** THE CORRELATION FUNCTION IS REQUIRED AT HALF INTEGER ** C ** INTERVALS, I.E. C(T), T=(TAU-0.5)*DT FOR TAU=1 .. NMAX. ** C ** THE COSINE TRANSFORM IS RETURNED AT HALF INTERVALS, I.E. ** C ** CHAT(OMEGA), OMEGA=(NU-0.5)*DOM FOR NU = 1 .. NMAX. ** C ******************************************************************* INTEGER NMAX REAL DT, C(NMAX), CHAT(NMAX) INTEGER TAU, NU REAL TAUH, NUH, NMAXH, PI, SUM PARAMETER ( PI = 3.1415927 ) C ******************************************************************* NMAXH = REAL ( NMAX ) - 0.5 C ** LOOP OVER OMEGA ** DO 20 NU = 1, NMAX NUH = REAL ( NU ) - 0.5 SUM = 0.0 C ** LOOP OVER T ** DO 10 TAU = 1, NMAX TAUH = REAL ( TAU ) - 0.5 SUM = SUM + C(TAU) * COS ( TAUH * NUH * PI / NMAXH ) 10 CONTINUE C ** FACTOR OF TWO FOR THE REAL COSINE TRANSFORM ** CHAT(NU) = 2.0 * DT * SUM 20 CONTINUE RETURN END SUBROUTINE FILONS ( DR, DK, NMAX, H, HHAT ) C ******************************************************************* C ** FOURIER SINE TRANSFORM BY FILON'S METHOD ** C ** ** C ** A SPATIAL CORRELATION FUNCTION, H(R), IS TRANSFORMED TO ** C ** HHAT(K) IN RECIPROCAL SPACE. ** C ** ** C ** REFERENCE: ** C ** ** C ** FILON, PROC ROY SOC EDIN, 49 38, 1928. ** C ** ** C ** PRINCIPAL VARIABLES: ** C ** ** C ** REAL KVEC THE WAVENUMBER ** C ** REAL RMAX MAXIMUM DIST IN CORREL. FUNCTION ** C ** REAL ALPHA, BETA, GAMMA FILON PARAMETERS ** C ** REAL H(NMAX) THE CORRELATION FUNCTION ** C ** REAL HHAT(NMAX) THE 3-D TRANSFORM ** C ** REAL DR INTERVAL BETWEEN POINTS IN H ** C ** REAL DK INTERVAL BETWEEN POINTS IN HHAT ** C ** INTEGER NMAX NO. OF INTERVALS ** 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 H(R) ** C ** IS AT R=0. THE MAXIMUM R FOR THE CORRELATION FUNCTION IS ** C ** RMAX=DR*NMAX. FOR AN ACCURATE TRANSFORM H(RMAX)=0. ** C ******************************************************************* INTEGER NMAX REAL DR, DK, H(0:NMAX), HHAT(0:NMAX) REAL RMAX, K, THETA, SINTH, COSTH REAL SINSQ, COSSQ, THSQ, THCUB, ALPHA, BETA, GAMMA REAL SE, SO, FOURPI, R INTEGER IR, IK C ******************************************************************* C ** CHECKS NMAX IS EVEN ** IF ( MOD ( NMAX, 2 ) .NE. 0 ) THEN STOP ' NMAX SHOULD BE EVEN ' ENDIF FOURPI = 16.0 * ATAN ( 1.0 ) RMAX = REAL ( NMAX ) * DR C ** LOOP OVER K ** DO 30 IK = 0, NMAX K = REAL ( IK ) * DK THETA = K * DR 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 ** THE INTEGRAND IS H(R) * R FOR THE 3-D TRANSFORM ** C ** DO THE SUM OVER THE EVEN ORDINATES ** SE = 0.0 DO 10 IR = 0, NMAX, 2 R = REAL ( IR ) * DR SE = SE + H(IR) * R * SIN ( K * R ) 10 CONTINUE C ** SUBTRACT HALF THE FIRST AND LAST TERMS ** C ** HERE THE FIRST TERM IS ZERO ** SE = SE - 0.5 * ( H(NMAX) * RMAX * SIN ( K * RMAX ) ) C ** DO THE SUM OVER THE ODD ORDINATES ** SO = 0.0 DO 20 IR = 1, NMAX - 1, 2 R = REAL ( IR ) * DR SO = SO + H(IR) * R * SIN ( K * R ) 20 CONTINUE HHAT(IK) = ( - ALPHA * H(NMAX) * RMAX * COS ( K * RMAX) : + BETA * SE + GAMMA * SO ) * DR C ** INCLUDE NORMALISING FACTOR ** HHAT(IK) = FOURPI * HHAT(IK) / K 30 CONTINUE RETURN END