I am trying to calculate the harmonic function using the Filon's method. My code looks like that:
program function1
implicit none
real(10) :: g, a, x, b, k, s
! frequency intervals for function and her transform
real(10) :: step, dstep
! number of inervals
integer, parameter :: nmax=10
integer :: j
! function and result of her transformation
real(kind=10), dimension(0:nmax) :: f, wynik
step=0.01
dstep=0.1
!constant value for function
a=1
b=-1
s=2
open(1, file='wykresik.dat', action='write', status='replace')
! loop over all frequencies in the function
do j=0,1000
k=float(j)*step
f=-1*cos(k*s)
enddo
call filonc(step,dstep,nmax,f,wynik)
write(1,*) k, wynik
! filon subroutine call
end program function1
include 'filon.f95'
And the included subroutine is the filons quadrature for high oscilatory functions and it looks like this:
SUBROUTINE FILONC ( DT, DOM, NMAX, D, CHAT )
INTEGER :: NMAX
REAL(kind=10) :: DT, DOM
REAL(kind=10), dimension(0:NMAX) ::D, CHAT
REAL(kind=10) :: TMAX, OMEGA, THETA, SINTH, COSTH, CE, CO
REAL(kind=10) :: SINSQ, COSSQ, THSQ, THCUB, ALPHA, BETA, GAMMA
INTEGER TAU, NU
IF ( MOD ( NMAX, 2 ) .NE. 0 ) then
stop ' NMAX SHOULD BE EVEN '
ENDIF
TMAX = float(NMAX) * DT
DO NU = 0, NMAX
OMEGA = float(NU) * DOM
THETA = OMEGA * DT
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
CE = 0.0
DO TAU = 0, NMAX, 2
CE = CE + D(TAU) * COS ( THETA * float( TAU ) )
ENDDO
CE = CE - 0.5 * ( D(0) + D(NMAX) * COS ( OMEGA * TMAX ) )
CO = 0.0
DO TAU = 1, NMAX - 1, 2
CO = CO + D(TAU) * COS ( THETA * REAL ( TAU ) )
ENDDO
CHAT(NU) = 2.0 * ( ALPHA * D(NMAX) * SIN ( OMEGA * TMAX )+ BETA * CE + GAMMA * CO ) * DT
ENDDO
ENDSUBROUTINE FILONC
I would like to know, why I am getting 11 results for a one k in that code, I have no clue why this is working like that.
This line
write(1,*) k, wynik
writes out the value of the scalar k and the entire rank-1 array wynik. In these two lines
integer, parameter :: nmax=10
...
real(kind=10), dimension(0:nmax) :: f, wynik
wynik is declared to have 11 elements.
So, naturally, the program writes out one value for k and 11 values for wynik
This seems so obvious that I suspect I have misunderstood what you are really asking - in which case please try to be much clearer about what you want to know.