!  Idealized, D-dimensional harmonic oscillator integrated using
!  the modified Euler method.
!  File:  harmnd.f90

!  Generate random initial condtions with abs(X) <= 1.0 and conduct
!  trials for various time steps.  D is number of dimensions, N = 2*D.
!  K is the spring constant, the mass is taken to be 1.0 .
!  Experiments use time interval [0:t1]; each calculation uses M*S steps,
!  reporting every S steps.  The state vector is organized as
!
!	x(1) x(2) . . . x(D) v(1) v(2) . . . v(D)
!	 1    2          D   D+1  D+2        2*D = N

      subroutine harmonic(D, N, K, t1)
      integer D, N
      real    K, t1

      real, parameter :: eps = 1E-10,  pi2 = 3.14159265/2.0

      INCLUDE 'meuler.h'

      INTERFACE                         ! Needs this to know how to
      function Fharm(X, t, N)           ! pass Fharm to Meuler
      integer N
      real, dimension(N) :: X, Fharm
      real t						! Not actually used
      END INTERFACE

      INTERFACE
      subroutine prstat(X, N, t)
      integer N
      real, dimension(N) :: X
      real t
      END INTERFACE

      EXTERNAL Fharm

      real,    dimension(N) ::  X0, X, A, X1
      logical, dimension(N) ::  WA, WB
      real t, hs, W
      integer M, S, i, seed(3)

      integer Dc			! Common data for convenient
      real    Kc			! interaction
      common /harmfe/ Kc, Dc

      if ( K < eps ) then
	K = 1.0
	print *, "Unsuitable K, using ", K
      end if  
      W  = sqrt(K)
      Kc = K						! Common
      Dc = D
      print '(" 3 integer seeds for RNG > "$)'
      read *, seed
      call random_seed(PUT=seed)
      call Random_Number(X0)
      X0 = 2.0*X0 - 1.0					! -1.0 <= X0 <= 1.0
!  Calculate amplitude and phase for analytic solution
!  Eliminating slightly complex WHERE's since f90 compiler
!  can't seem to deal with them.
      do i = 1, D
        if ( abs(X0(D+i)) < eps ) then          ! velocity = 0
          A(i) = X0(i);  A(D+i) = pi2           ! amplitude/phase
        else if ( abs(X0(i)) < eps ) then       ! position = 0
          A(i) = X0(D+i)/W;  A(D+i) = 0.0
        else                                    ! neither 0
          A(D+i) = atan2( X0(i), X0(D+i)/W )    ! phase
          A(i)   = X0(i)/sin(A(D+i))            ! amplitude
        end if
        X1( i ) =   A(i)*sin(W*t1 + A(D+i))       ! Analytic solution
        X1(D+i) = W*A(i)*cos(W*t1 + A(D+i))
      end do

!      WA(1:D) = ( abs(X0(D+1:N)) < eps )		! velocity = 0
!      WHERE ( WA(1:D) )
!	A( 1 :D) = X0(1:D)				! amplitude `a'
!	A(D+1:N) = pi2					! phase `b'
!     END WHERE
!     WB(1:D) = (.not. WA(1:D))
!     WB(1:D) = (WB(1:D) .and. (abs(X0(1:D)) < eps))	! position = 0
!     WHERE ( WB(1:D) )
!	A( 1 :D) = X0(D+1:N)/W
!	A(D+1:N) = 0.0
!       WA(1:D) = .TRUE.
!    END WHERE   					! neither 0
!     WHERE ( .not. WA(1:D) )
!	A(D+1:N) = atan2( X0(1:D), X0(D+1:N)/W )	! phase
!	A( 1 :D) = X0(1:D)/sin(A(D+1:N))		! amplitude
!     END WHERE
!     X1( 1 :D) =   A(1:D)*sin(W*t1 + A(D+1:N))		! Analytic solution
!     X1(D+1:N) = W*A(1:D)*cos(W*t1 + A(D+1:N))

      call prstat(X0, N, 0.0)
!  Run experiments
      print *, "Specify M <= 0 to exit."
      do while ( .true. )
      print '(" M S > "$)'
      read *, M, S
      if ( M <= 0 ) go to 99				! Exit
      X  = X0						! Initial conditions
      h  = t1/(M*S)					! Time step
      hs = t1/M						! Time for s steps
      t  = 0.0
      do i = 1, M					! S steps and report
	call mEuler(X, Fharm, t, h, S, N)		! Integrate
	t = hs*i					! New time
        call prstat(X, N, t)				! Report
      end do ! i = 1, M
      print *, "Analytic solution:"
      call prstat(X1, N, t1)
      print 21, maxval(abs(X - X1))			! Error at t1
 21   FORMAT (" Max abs error: ", E14.5)
      end do ! while ( .true. )
 99   continue						! For exit
      end subroutine harmonic

!  Function for N/2-dimensional harmonic oscillator

      function Fharm(X, t, N)
      integer N
      real, dimension(N) :: X, Fharm
      real t						! Not actually used

      real    Kc					! Spring constant
      integer Dc					! Dimension = N/2
      common /harmfe/ Kc, Dc				! For Fharm

      Fharm(  1 :Dc) = X(Dc+1:N)			! Velocities
      Fharm(Dc+1:N ) = -Kc*X(1:Dc)			! Accelerations
      end function Fharm

!  Print current status

      subroutine prstat(X, N, t)
      integer N
      real, dimension(N) :: X
      real t

      integer Dc			! Common data for convenient
      real    Kc			! interaction
      common /harmfe/ Kc, Dc

      print *, "Simulated time:  ", t
      print 21, X( 1   : Dc)				! Positions
      print 21, X(Dc+1 :  N)				! Velocities
 21   FORMAT (1X, 5(F12.7:))
      end subroutine prstat

!  Obtain parameters

      program harm

      integer Dc			! Common data for convenient
      real    Kc			! interaction
      common /harmfe/ Kc, Dc

      integer D, N
      real K, t1

      print *, "D-dimensional harmonic oscillator.  Specify D <= 0"
      print *, "to quit."
      do while ( .true. )
	print '(" D K t1 > "$)'
        read *, D, K, t1
        if ( D <= 0 ) go to 99				! Exit
	N = 2*D						! Array size
	call harmonic(D, N, K, t1)
      end do
 99   continue						! For exit
      end program harm
