program ComputeEB
!
!
implicit none
integer,parameter :: iconst = 6 ! integer used to control the resolution
integer,parameter :: Nz = iconst*5 ! half the number of zones in each domain
integer,parameter :: ntot = 3*Nz ! total number of timesteps
integer,parameter :: printnum = iconst ! number of timesteps between data printout
integer,parameter :: c = 3E8 ! speed of prop in a vacuum
REAL, PARAMETER :: pi = 3.1415927
integer :: i,j,k, n, Nzh, numout
double precision :: x(-Nz:Nz), y(-Nz:Nz), z(-Nz:Nz)
double precision :: Ex(-Nz:Nz), Bx(-Nz:Nz), Ey(-Nz:Nz), By(-Nz:Nz), Ez(-Nz:Nz), Bz(-Nz:Nz)
double precision :: Exh(-Nz:Nz), Bxh(-Nz:Nz), Eyh(-Nz:Nz), Byh(-Nz:Nz), Ezh(-Nz:Nz), Bzh(-Nz:Nz)
double precision :: Jx(-Nz:Nz), Jy(-Nz:Nz), Jz(-Nz:Nz)
double precision :: v, L, dt, h, twoh, dth, ed, time, Energy
double precision :: R, rho, omega, Toff, epsilon, mu
character(len=20) :: filenum
!
!
!
! Input parameters and constants
L = 0.1d0 ! length of half the domain
h = L/Nz ! length of each zone
twoh = 2.0d0*h ! 2h
dt = 0.25d0*h/c ! timestep
write(*,*) 'stuff = ', v, h, dt
dth = 0.5d0*dt
Nzh = Nz/2
R = .025d0
rho = 1E10
omega = 1E-3
Toff = 8E-11
epsilon = 8.85E-12
mu = pi*4E-7
!
! Set up the 3D matrix
do i = -Nz, Nz
x(i) = i*h
y(i) = i*h
z(i) = i*h
enddo
!
! Initial conditions
do i = -Nz, Nz
Ex(i) = 0.0d0
Ey(i) = 0.0d0
Ez(i) = 0.0d0
Bx(i) = 0.0d0
By(i) = 0.0d0
Exh(i) = 0.0d0
Eyh(i) = 0.0d0
Ezh(i) = 0.0d0
Bxh(i) = 0.0d0
Byh(i) = 0.0d0
Bzh(i) = 0.0d0
Jx(i) = 0.0d0
Jy(i) = 0.0d0
Jz(i) = 0.0d0
enddo
open(unit=11,file='Ex0')
do i = -Nz, Nz
write(11,111) x(i), Ex(i) ! write data to file
enddo
close(11)
open(unit=11,file='Ey0')
do i = -Nz, Nz
write(11,111) y(i), Ey(i) ! write data to file
enddo
close(11)
open(unit=11,file='Ez0')
do i = -Nz, Nz
write(11,111) z(i), Ez(i) ! write data to file
enddo
close(11)
open(unit=11,file='Bx0')
do i = -Nz, Nz
write(11,111) x(i), Bx(i) ! write data to file
enddo
close(11)
open(unit=11,file='By0')
do i = -Nz, Nz
write(11,111) y(i), By(i) ! write data to file
enddo
close(11)
open(unit=11,file='Bz0')
do i = -Nz, Nz
write(11,111) z(i), Bz(i) ! write data to file
enddo
close(11)
!
! Evolve
time = 0.0d0
numout = 0
do n = 1, ntot ! Compute E and B
do i = -Nz+1, Nz-1
do j = -Nz+1, Nz-1
do k = -Nz+1, Nz-1
! determine if we are within the current
if (time < Toff .AND. sqrt(x(i)**2+y(j)**2+z(k)**2) < R) then
Jx(k) = (-y(j)/2)*rho*omega*(1+cos(pi*sqrt(x(i)**2+y(j)**2+z(k)**2)/R))
Jy(k) = (x(i)/2)*rho*omega*(1+cos(pi*sqrt(x(i)**2+y(j)**2+z(k)**2)/R))
endif
! advance a half timestep
! E
Exh(k) = Ex(k) + (dt/(4*epsilon*mu*h))*(Bz(j+1)-Bz(j-1)-By(k+1)+By(k-1))&
& - (dt/(2*epsilon))*Jx(k)
Eyh(k) = Ey(k) + (dt/(4*epsilon*mu*h))*(Bx(k+1)-Bx(k-1)-Bz(i+1)+Bz(i-1))&
& - (dt/(2*epsilon))*Jy(k)
Ezh(k) = Ex(k) + (dt/(4*epsilon*mu*h))*(By(i+1)-By(i-1)-Bx(j+1)+Bx(j-1))&
& - (dt/(2*epsilon))*Jz(k)
! B
Bxh(k) = Bx(k) - (dt/(4*h))*(Ez(j+1)-Ez(j-1)-Ey(k+1)+Ey(k-1))
Byh(k) = By(k) - (dt/(4*h))*(Ex(k+1)-Ex(k-1)-Ez(i+1)+Ez(i-1))
Bzh(k) = Bz(k) - (dt/(4*h))*(Ey(i+1)-Ey(i-1)-Ex(j+1)+Ex(j-1))
! advance a full timestep
! E
Ex(k) = Ex(k) + (dt/(2*epsilon*mu*h))*(Bzh(j+1)-Bzh(j-1)-Byh(k+1)+Byh(k-1))&
& - (dt/(2*epsilon))*Jx(k)
Ey(k) = Ey(k) + (dt/(2*epsilon*mu*h))*(Bxh(k+1)-Bxh(k-1)-Bzh(i+1)+Bzh(i-1))&
& - (dt/(2*epsilon))*Jy(k)
Ez(k) = Ex(k) + (dt/(2*epsilon*mu*h))*(Byh(i+1)-Byh(i-1)-Bxh(j+1)+Bxh(j-1))&
& - (dt/(2*epsilon))*Jz(k)
! B
Bx(k) = Bx(k) - (dt/(2*h))*(Ezh(j+1)-Ezh(j-1)-Eyh(k+1)+Eyh(k-1))
By(k) = By(k) - (dt/(2*h))*(Exh(k+1)-Exh(k-1)-Ezh(i+1)+Ezh(i-1))
Bz(k) = Bz(k) - (dt/(2*h))*(Eyh(i+1)-Eyh(i-1)-Exh(j+1)+Exh(j-1))
enddo
enddo
enddo
time = time + dt
write(*,*) 'step = ', n, ' time = ', time
! Print out data if n equals a multiple of printnum
if (n .eq. (n/printnum)*printnum) then
numout = numout + 1 ! increment numout (number of print-outs)
write(filenum,*) numout ! write integer numout to character filenum
filenum = adjustl(filenum) ! remove leading spaces
open(unit=33,file='Ex'//trim(filenum)) ! create file with appropriate filename
do i = -Nz, Nz
write(33,111) x(i), Ex(i) ! write data to file
enddo
close(33)
open(unit=33,file='Ey'//trim(filenum)) ! create file with appropriate filename
do i = -Nz, Nz
write(33,111) y(i), Ey(i) ! write data to file
enddo
close(33)
open(unit=33,file='Ez'//trim(filenum)) ! create file with appropriate filename
do i = -Nz, Nz
write(33,111) z(i), Ez(i) ! write data to file
enddo
close(33)
open(unit=33,file='Bx'//trim(filenum)) ! create file with appropriate filename
do i = -Nz, Nz
write(33,111) x(i), Bx(i) ! write data to file
enddo
close(33)
open(unit=33,file='By'//trim(filenum)) ! create file with appropriate filename
do i = -Nz, Nz
write(33,111) y(i), By(i) ! write data to file
enddo
close(33)
open(unit=33,file='Bz'//trim(filenum)) ! create file with appropriate filename
do i = -Nz, Nz
write(33,111) z(i), Bz(i) ! write data to file
enddo
close(33)
endif
enddo
!
111 format(2E20.10)
555 format(2E20.10)
999 format(2E20.10)
end program