! compile: ifort cheby.f90 -o cheby -fpp -lfftw3 -fpp

! Programm zum Berechnen von Ableitungen mit Cheby-Polynomen.

PROGRAM cheby

IMPLICIT NONE

#INCLUDE 'fftw3.f'

INTEGER, PARAMETER :: dm=128
INTEGER :: i
INTEGER*8 :: plan_forw, plan_back
DOUBLE PRECISION :: pi,x
DOUBLE PRECISION, DIMENSION(dm) :: feld, coeff, trigo_feld, gridp


pi=DACOS(-1.0d0)

! FFTW plaene >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#ifdef zeros
! falls nullstellen als kollokationspunkte
CALL dfftw_plan_r2r_1d(plan_forw,dm,feld,coeff,FFTW_REDFT10,FFTW_ESTIMATE)
CALL dfftw_plan_r2r_1d(plan_back,dm,coeff,feld,FFTW_REDFT01,FFTW_ESTIMATE)
#else
! falls extrema und randpunkte als kollokationspunkte
CALL dfftw_plan_r2r_1d(plan_forw,dm,feld,coeff,FFTW_REDFT00,FFTW_ESTIMATE)
CALL dfftw_plan_r2r_1d(plan_back,dm,coeff,feld,FFTW_REDFT00,FFTW_ESTIMATE)
#endif
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

! Feld initialisieren >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
OPEN(9,FILE="vorher.dat")

DO i=1,dm

#ifdef zeros
  ! nullstellen
  gridp(i)=DCOS(pi/DBLE(2*dm)*(2*i-1))
#else
  ! extrema und randpunkte
  gridp(i)=DCOS(pi/DBLE(dm-1)*(i-1))
#endif

  !feld(i)=DEXP(-gridp(i)**2/0.025d0)
  !feld(i)=DCOS(pi*gridp(i))+1
  feld(i)=DSIN(pi*gridp(i))

  WRITE(9,*)  gridp(i), feld(i)

END DO

CLOSE(9)
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

! FFT in den Fourierraum
CALL dfftw_execute_r2r(plan_forw,feld,coeff)
#ifdef zeros
 coeff=coeff*DBLE(2*dm)**-1.0d0
#else
 coeff=coeff*DBLE(2*(dm-1))**-1.0d0
#endif


! Ableitung berechnen
CALL derivative(coeff,dm)
CALL derivative(coeff,dm)


! FFT in den Ortsraum
CALL dfftw_execute_r2r(plan_back,coeff,feld)


! Feld initialisieren >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
OPEN(9,FILE="nacher.dat")
DO i=1,dm
  WRITE(9,*)  gridp(i), feld(i)
END DO
CLOSE(9)
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

END PROGRAM


!******************************************************************************
!******************************************************************************


SUBROUTINE derivative(coeff,dm)

IMPLICIT NONE

INTEGER :: dm
INTEGER :: i
DOUBLE PRECISION :: pi
DOUBLE PRECISION, DIMENSION(dm) :: coeff,coeff_backup

pi=DACOS(-1.0d0)

 coeff_backup=coeff

 coeff(dm)=0.0d0
 coeff(dm-1)=0.0d0

DO i=1,dm-2
  coeff(dm-i-1)=2.0d0*(dm-i-1)*coeff_backup(dm-i)+coeff(dm-i+1)
END DO



END SUBROUTINE


!******************************************************************************
!******************************************************************************

! diese Subroutine wird nicht benutzt
! nur wichtig, falls die Ableitung nicht ueber Ruecktrafo der 
! Koeffizienten bestimmt wird

SUBROUTINE trigo_derivative(coeff,feld_out,dm)

IMPLICIT NONE

INTEGER :: dm
INTEGER :: i,k
DOUBLE PRECISION :: pi, theta
DOUBLE PRECISION, DIMENSION(dm) :: coeff, feld_out

pi=DACOS(-1.0d0)
feld_out=0

DO i=2,dm
  DO k=1,dm

    theta=pi*(2*i-1)/DBLE(2*dm)

    feld_out(i)=feld_out(i)+2.0*coeff(k)*(k-1)*DSIN((k-1)*theta )/DSIN(theta)

  END DO
END DO

END SUBROUTINE