Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[Fortran] Filon's method Fourier Transform

  1. Jun 29, 2015 #1
    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.

    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: Jun 29, 2015
  2. jcsd
  3. Jun 29, 2015 #2

    jtbell

    User Avatar

    Staff: Mentor

    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.
     
  4. Jun 29, 2015 #3
    I do not really know how I should pass that array as the fourth argument. Could you help me?
     
  5. Jun 29, 2015 #4

    Mark44

    Staff: Mentor

    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.
     
  6. Jun 29, 2015 #5
    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
    Code (Fortran):
    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:
    Code (Fortran):

    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

    Code (Fortran):
      Included at function1.f95:26:

               CHAT(NU) = 2.0 * ( ALPHA * D(NMAX) * SIN ( OMEGA * TMAX ))+ BETA * C
               1
    Error: Unclassifiable statement at (1)
     
    Last edited: Jun 29, 2015
  7. Jun 29, 2015 #6

    Mark44

    Staff: Mentor

    Please show us the code that calls the FILON subroutine.

    Also, the code that is shown in the error message --
    Code (Fortran):
    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.

     
  8. Jun 30, 2015 #7
    I have already fixed that. I just had one additional bracket in this formula. Thank you for help. It works properly now
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: [Fortran] Filon's method Fourier Transform
Loading...