! compile: ifort cheby_wlg.f90 -o cheby_wlg -fpp -lfftw3 -ldislin -fpp

! Programm zu Simulation der Waermeleitungsgleichung mit Hilfe von Cheby-Polynomen.
! Je nach Anfangsbedingungen koennen Oszillationen auftreten, viel Spass beim Debuggen!

MODULE para

IMPLICIT NONE

INTEGER, PARAMETER :: dm=128, screen_every=100
INTEGER*8 :: plan_forw, plan_back
LOGICAL :: pitfall=.TRUE.
DOUBLE PRECISION :: dt=1e-7, kappa=1.0, pi
INTEGER :: n_iter=1e7

REAL :: mini, maxi


END MODULE


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


PROGRAM cheby_wlg

USE para

IMPLICIT NONE

#INCLUDE 'fftw3.f'


INTEGER :: i_iter,i
REAL, PARAMETER :: pong=ACOS(-1.0)
DOUBLE PRECISION, DIMENSION(dm) :: gridp
DOUBLE PRECISION, DIMENSION(dm) :: feld, k1, k2, k3, k4


pi=DACOS(-1.0d0)

! dislin
CALL METAFL ('XWIN')    ! Dislin Initialisierung
CALL PAGE(2500,2500)
CALL SCRMOD('REVERS')
CALL SCLMOD('FULL')
CALL DISINI()


! initialisieren >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
OPEN(9,FILE="ic.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)-0.75)**2/0.025d0)
  feld(i)=DCOS(pi*gridp(i))+1
  !feld(i)=DSIN(pi*gridp(i))

  IF(gridp(i).GT.-0.7.AND.gridp(i).LT.0.9) feld(i)=1.0

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

END DO
CLOSE(9)

! RB beachten
feld(1)=1.0d0
feld(dm)=0.0d0

mini=MINVAL(feld)
maxi=MAXVAL(feld)
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


DO i_iter=0,n_iter


  IF(MOD(i_iter,screen_every).EQ.0) THEN 

    OPEN(9,FILE="feld.dat")
    DO i=1,dm
      WRITE(9,*)  gridp(i),feld(i)
    END DO
    CLOSE(9)

    CALL display(feld,gridp)


  END IF

  CALL rhs(feld,k1)
  CALL rhs(feld+dt/2.0d0*k1,k2)
  CALL rhs(feld+dt/2.0d0*k2,k3)
  CALL rhs(feld+dt*k3,k4)

  feld=feld+dt/6.0d0*(k1+2.0d0*k2+2.0d0*k3+K4)


END DO



CALL DISFIN()

END PROGRAM 


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


SUBROUTINE rhs(feld_in,feld_out)

USE para

IMPLICIT NONE

#INCLUDE 'fftw3.f'

DOUBLE PRECISION, DIMENSION(dm) :: feld_in,feld_out, coeff

IF(pitfall) THEN

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

  pitfall=.FALSE.
END IF

!RB beachten
feld_in(1)=1.0d0
feld_in(dm)=0.0d0


CALL dfftw_execute_r2r(plan_forw,feld_in,coeff)
#ifdef zeros
 coeff=coeff*DBLE(2*dm)**-1.0d0
#else
 coeff=coeff*DBLE(2*(dm-1))**-1.0d0
#endif

CALL derivative(coeff)
CALL derivative(coeff)


CALL dfftw_execute_r2r(plan_back,coeff,feld_out)

feld_out=kappa*feld_out

END SUBROUTINE


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


SUBROUTINE derivative(coeff)

USE para

IMPLICIT NONE

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


 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


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


SUBROUTINE display(feld_in,gridp)

USE para
USE dislin

IMPLICIT NONE

INTEGER :: i
DOUBLE PRECISION, DIMENSION(dm) :: feld_in,gridp

mini=MINVAL(feld_in)
maxi=MAXVAL(feld_in)

!mini=0.0
!maxi=1.0

CALL ERASE()
CALL AXSSCL('LIN','XY')
CALL NAME('','X')
CALL NAME('','Y')
CALL GRAF(-1.0, 1.0, -1.0, 0.5 , mini, maxi ,mini, (maxi-mini)/10.0)
CALL COLOR('RED')
CALL CURVE(REAL(gridp),REAL(feld_in),dm)
CALL TITLE()
CALL COLOR('WHITE')
CALL ENDGRF()




END SUBROUTINE