[Fortran] Filon's method Fourier Transform

[Galizius]

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.

   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

 

[jtbell]

Before calling this function, use your function to calculate a series of values at points spaced by some amount DT, and store those values in an array. Pass that array as the fourth argument to this function.

 

[Galizius]

Before calling this function, use your function to calculate a series of values at points spaced by some amount DT, and store those values in an array. Pass that array as the fourth argument to this function.

I do not really know how I should pass that array as the fourth argument. Could you help me?

 

[Mark44]

I do not really know how I should pass that array as the fourth argument. Could you help me?

Then it looks like you're stuck, at least until you can figure out how to pass an array to a subroutine in Fortran.

According to the rules of this forum we are not going to simply do your work for you with no effort shown on your part.

 

[Galizius]

Then it looks like you're stuck, at least until you can figure out how to pass an array to a subroutine in Fortran.

According to the rules of this forum we are not going to simply do your work for you with no effort shown on your part.

Yes, I know and I strongly support such a position.
I have done following code for counting values of my function. I do not really know what to do further. How to force a program to count the transform of that function? I would be really grateful for any further suggestions.

The code is :
Auxiliary function

real(10) function g(z)
implicit none
real(10) :: z, infinity, s,b, I
real(10) :: y
y=huge(1.d0)
infinity = y+y
s=2
b=-1
I=3

if (z>=0 .and. z<s) then
  g=infinity
elseif (z>=s .and. z<I) then
  g=b
elseif (z>=I) then
  g=0
endif
end function g

And the main code:

program function1
implicit none
real(10) :: g, a, x, dt
integer, parameter :: nmax=128
integer :: k
real(kind=10), dimension(nmax) :: f
dt=0.1
a=1
open(1, file='wykresik.dat', action='write', status='replace')

do k=1,nmax
  x=float(k)*dt
  if (x>0) then
  f(k)=exp((-a*g(x))-1)
  else
  f(k)=0
  end if
write(1,*)x,f(k)
enddo
end program function1
include 'g.f95'

Edit:
Okay, I know that I exchange the C as f in the subroutine... But now I am getting following error

  Included at function1.f95:26:

           CHAT(NU) = 2.0 * ( ALPHA * D(NMAX) * SIN ( OMEGA * TMAX ))+ BETA * C
           1
Error: Unclassifiable statement at (1)

 

[Mark44]

Please show us the code that calls the FILON subroutine.

Also, the code that is shown in the error message --

CHAT(NU) = 2.0 * ( ALPHA * D(NMAX) * SIN ( OMEGA * TMAX ))+ BETA * C

-- is different from the code you showed in post #1.
Maybe it was just cut off.

CHAT(NU) = 2.0 * ( ALPHA * C(NMAX) * SIN ( OMEGA * TMAX ) + BETA * CE + GAMMA * CO ) * DT

 

[Galizius]

Please show us the code that calls the FILON subroutine.

Also, the code that is shown in the error message --

CHAT(NU) = 2.0 * ( ALPHA * D(NMAX) * SIN ( OMEGA * TMAX ))+ BETA * C

-- is different from the code you showed in post #1.
Maybe it was just cut off.

I have already fixed that. I just had one additional bracket in this formula. Thank you for help. It works properly now