!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
!
MODULE EB
      USE ODE, ONLY: wp
      IMPLICIT NONE
      PUBLIC
      INTEGER, PARAMETER :: n_eq = 6 ! First 3 are velocities, last 3 are coordinates
      ! Here are the physical parameters of the problem with default values:
      REAL (KIND=wp) :: B0 = 1.0, E0 = 1.0, theta = 0.0, omega = 0.0, phi = 0.0, dB = 0.0, dE = 0.0
CONTAINS
      SUBROUTINE EB_force (x, y, dy)
         INTEGER, PARAMETER :: wp = KIND (0.0D0)
         REAL (KIND=wp), INTENT (IN) :: x
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (IN) :: y
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (OUT) :: dy
         REAL (KIND=wp) :: B, E
         B = B0 * (1.0_wp+dB*Cos(omega*x))
         E = E0 * (1.0_wp+dE*Cos(omega*x+phi))
     	 ! Use Lentz's law to derive this:  	
         dy = (/ B * y (2) + E * Sin (theta), - B * y (1), E * Cos (theta), y (1), y (2), y (3) /)
      END SUBROUTINE EB_force
END MODULE EB
!
PROGRAM ODE_EB
      USE EB
      USE ODE
      USE SimpleGraphics
!
      INTEGER, PARAMETER :: sp = KIND (0.0E0), n_pts = 1000
      REAL (KIND=wp), DIMENSION (n_eq, n_pts) :: vr ! The velocities and positions together as y
      REAL (KIND=wp), DIMENSION (3) :: v0 ! Inital velocity
      REAL (KIND=wp), DIMENSION (n_pts) :: t
      REAL (KIND=wp) :: t_min = 0.0_wp, t_max = 50.0_wp
      INTEGER :: i
      CHARACTER (50), DIMENSION (2) :: title ! The title of the plot
!
      ! I only set the non-default physical parameters
      ! B0=1.0_wp
      ! E0=0.05_wp
      ! theta=2*ATAN(1.0_wp) ! theta=pi/2
      ! dB=0.1_wp
      ! omega=1.0_wp
      WRITE (*,*) "Enter the values of E0,B0,dE,dB,theta, omega and phi"
      READ (*,*) E0, B0, dE, dB, theta, omega, phi
      WRITE (*,*) "Enter the initial velocity v0 and the end of the time interval t_max"
      READ (*,*) v0, t_max
!
      CALL ODESolve (f=EB_force, n_eq=n_eq, x_range= (/ t_min, t_max /), n_points=n_pts, &
                     y0= (/ v0, 0.0_wp, 0.0_wp, 0.0_wp /), x=t, y=vr, relerr=1D-5)
!
      ! This is an example of using Plot3D plot:
      title (1) = "Orbit of a charged particle in crossed constant EB fields" ! First line of the title
      title (2) = "A. Donev, 11/5/00" ! Second line of title
      CALL InitGraphics (file="Orbit.png", file_type="PNG", plot_title=title, x_label="x", y_label="y", z_label="z")
      CALL Plot3D(x=REAL(vr(4,:),sp),y=REAL(vr(5,:),sp),z=REAL(vr(6,:),sp), &
                  plot_spec="L-R",new_plot=.TRUE.) ! Just a line plot
      CALL EndGraphics ()
!
END PROGRAM ODE_EB
!