MODULE define_DDEs

  IMPLICIT NONE

! Define the number of ODEs.
  INTEGER, PARAMETER :: NEQN=2,NLAGSD=0

! Define the problem dependent parameters.
  DOUBLE PRECISION, PARAMETER :: &
    EPS=0.05D0,LAMBDA=2.0D0,KAPPA=4.0D0,OMEGA=0.375D0

CONTAINS

  SUBROUTINE ODES(T,Y,DY)

! Define the derivatives for ODEAVG.

    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY
    INTENT(IN) :: T,Y
    INTENT(OUT) :: DY

    DY(1) = Y(2)
    DY(2) = - EPS*LAMBDA*Y(2) - Y(1) - EPS*KAPPA*Y(1)**3 &
            + EPS*COS((1.0D0+EPS*OMEGA)*T)

    RETURN
  END SUBROUTINE ODES
  
  SUBROUTINE DDES(T,Y,Z,DY)
! Define the derivatives for pure ODE check solution.
    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(NEQN) :: Y,DY
    DOUBLE PRECISION, DIMENSION(NEQN,:) :: Z
    INTENT(IN) :: T,Y,Z
    INTENT(OUT) :: DY
    CALL ODES(T,Y,DY)
    RETURN
  END SUBROUTINE DDES

END MODULE define_DDEs

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

PROGRAM holmes

! M.H. Holmes, Introduction to Perturbation Methods,
! Springer, New York, 1995, solves the ODE
!
!   y'' + epsilon*lambda*y' + y + ...
!         epsilon*kappa*y^3 = epsilon*cos((1+epsilon*omega)t)
! 
! with parameter values that results in a highly oscillatory
! solution with slowly increasing magnitude.  Here the solution
! is compared to the RMS average computed with ODEAVG.
! The DDEs are defined in the module define_DDEs. The problem is
! solved here with ODEAVG/DDE_SOLVER and its output is written
! to a file holmes.dat. The auxilary function holmes.m imports
! the data into Matlab and plots it.

  USE define_DDEs
  USE DDE_SOLVER_M

  IMPLICIT NONE

  ! The quantities passed to the solver in NVAR are
  !   NEQN  = number of ODEs
  !   AVG_OPT = averaging option
  ! NEQN is defined in the module define_DDEs as a
  ! PARAMETER so it can be used for dimensioning
  ! arrays here.

  INTEGER, DIMENSION(2) :: NVAR
  INTEGER, DIMENSION(2) :: NVARD = (/NEQN,NLAGSD/)
  DOUBLE PRECISION, DIMENSION(NEQN) :: Y0
  DOUBLE PRECISION, DIMENSION(1) :: DELTA,BETA
  DOUBLE PRECISION, DIMENSION(NEQN) :: ABSERR,RELERR

  TYPE(DDE_SOL) :: SOL,SOLODE
  !   SOL%NPTS         -- NPTS,number of output points.
  !   SOL%T(NPTS)      -- values of independent variable, T.  
  !   SOL%Y(NPTS,NEQN) -- values of dependent variable, Y,
  !                       corresponding to values of SOL%T. 
  TYPE(DDE_OPTS) :: OPTS
  DOUBLE PRECISION :: T0,TF
  INTEGER :: AVG_OPT(2),I,J,MKIND,MAVG
  
! Open the output files used by the holmes.m plotting script.
  OPEN(UNIT=6, FILE='holmes.dat')
  OPEN(UNIT=7, FILE='holmesd.dat')
  OPEN(UNIT=8, FILE='holmesode.dat')

! Define the averaging option to be used.
  PRINT *,' Please enter 1 or 2, the number of averages to be'
  PRINT *,' computed, followed by <return>.'
  READ  *, MAVG
! PRINT *,' Please enter 1 for solution averaging or 2 for'
! PRINT *,' RMS amplitude averaging, followed by <return>.'
! READ  *, MKIND
  PRINT *,' '
  MKIND = 2
  AVG_OPT(1) = MAVG
  AVG_OPT(2) = MKIND

! Define the moving average DELTA.
  DELTA(1) = 35.0D0

! Save DELTA and AVG_OPT for use in the Matlab plotting
! script holmes.m.
  WRITE(7,*) DELTA(1), AVG_OPT(1)

! Define the integration parameters.
  NVAR = (/NEQN,AVG_OPT(1)/)
  T0 = 0.0D0
  TF = 250.0D0
  Y0 = (/0.0D0,0.0D0/)
  ABSERR(1:NEQN) = 1.0D-6
  RELERR(1:NEQN) = 1.0D-3

! Supply the integration options for the pure ODE check solution.
  BETA(1) = 0.0D0
  OPTS = DDE_SET(RE_VECTOR=RELERR,AE_VECTOR=ABSERR)
  SOLODE = DDE_SOLVER(NVARD,DDES,BETA,Y0,TSPAN=(/T0,TF/),OPTIONS=OPTS) 
  IF (SOLODE%FLAG == 0) THEN
!    Print integration statistics.
     CALL PRINT_STATS(SOLODE)  
     ! If the solver was successful, write the solution
     ! to a file for subsequent plotting in Matlab.
     DO I = 1, SOLODE%NPTS
        WRITE(UNIT=8,FMT='(3D12.4)') SOLODE%T(I),(SOLODE%Y(I,J),J=1,NEQN)
     END DO
  ELSE
     PRINT *,' Abnormal return from DDE_SOLVER with FLAG = ',&
     SOLODE%FLAG
     STOP
  END IF

! Supply the integration options for the ODEAVG solution.
  OPTS = DDE_SET(RE_VECTOR=RELERR,AE_VECTOR=ABSERR)
! Perform the integration.
  SOL = ODEAVG(AVG_OPT,DELTA,NVAR,ODES,Y0,TSPAN=(/T0,TF/), &
               OPTIONS=OPTS) 
  IF (SOL%FLAG == 0) THEN
!    Print integration statistics.
     CALL PRINT_STATS(SOL)    
     ! If the solver was successful, write the solution
     ! to a file for subsequent plotting in Matlab.
  
     DO I = 1, SOL%NPTS
        WRITE(UNIT=6,FMT='(3D12.4)') SOL%T(I),(SOL%Y(I,J),J=1,NEQN)
     END DO
     PRINT *,' Normal return from DDE_SOLVER with results'
     PRINT *,' written to the file holmes.dat.'
     PRINT *,' '
     PRINT *,' These results can be accessed in Matlab'
     PRINT *,' and plotted by'
     PRINT *,' '
     PRINT *,' >> [t,y] = holmes.m'
     PRINT *,' '
  ELSE
     PRINT *,' Abnormal return from DDE_SOLVER with FLAG = ',&
     SOL%FLAG
     STOP
  END IF

  STOP
END PROGRAM holmes