c*********************************
c 1-dimensional FDTD code
c Cindy Furse
c University of Utah
c 4-22-94
c*********************************
      parameter(nx_max=300) ! maximum size of FDTD space
      implicit complex(c)
      common/time_domain/ey(nx_max),hz(nx_max)
      real d1(nx_max),c1(nx_max),c2(nx_max) ! constants
      real er_mod(0:200),sig_mod(0:200)
      common/wave1/w ! omega
     &            ,B ! beta
     &            ,F ! frequency (Hz)
     &            ,dt ! time step (seconds)
     &            ,ipulse ! =0 for CW, 1 for RC pulse
     &            ,nstep  ! # of time steps in simulation
      common/time_step/k ! present time step
      common/fdtd_size/nx
      common/incident_location/inc
      common/cell_size/dx

c==============================
c     USER INPUT:
c==============================

c     --pulse(1) or CW(0)--
      ipulse = 0 ! CW
c     ipulse = 1 ! pulse

c     --Size of this FDTD space--
      nx = 200 ! (must be less than nx_max in parameter statement)
  
c     --Frequency--
      F = 60 ! Hz
      pi = acos(-1.)
      w = 2.*pi*F

c     --Cell size--
      dx = 6.e-3 ! meters
      vpo = 2.996e8 ! m/s: velocity of propagation  in air
      dt = dx/(2.*vpo) ! seconds

c     --Location of source--
      inc = 10 ! location of source

c     --Number of time steps in simulation--
      To = 1./F ! period of sine wave
      nTo = nint(To/dt) ! number of time steps in one period
      rcycle = 8.0 ! number of cycles to run
c     nstep = rcycle *nTo
      nstep = 2000

c     --Set up media--
c     (for metal use er1=4.0, sig1=3.27e7 or set BC)
      eo = 8.854e-12
      uo = 1.26e-6
      eta = 377.
      er = 1. ! air
      sig = 0.
      ur = 1.
      vp = vpo/sqrt(er)
      do i=1,nx
         er_mod(i) = 1.
         sig_mod(i) = 0.
      enddo
      nrad = 28
      ncenter = 45
      do i=10+ncenter-nrad,10+ncenter+nrad
         sig_mod(i) = .35
      enddo

c     --Set up model and calculate constants--
      do 5 i=1,nx
         aer = er_mod(i)
         asig = sig_mod(i)
         d1(i) = -dt/(uo*ur*dx)
         c1(i) = (1.-asig*dt/(2.*aer*eo))/(1.+asig*dt/(2.*aer*eo))
         c2(i) = dt/(dx*aer*eo*(1.+asig*dt/(aer*eo*2.)))
5     continue

c     --Constants for adding plane wave source terms--
      hzconst = 1./(1. + sig*dt/(2.*er*eo))
      eyconst = +dt/(uo*ur)

c     --Initialize arrays--
      do 10 i=1,nx
         ey(i) = 0.
         hz(i) = 0.
10    continue

c===============================
c Run FDTD 
c===============================
      do 500 k=1,nstep

c        --SOLVE FOR Ey--
         do 20 i=1,nx
            ey(i) = c1(i)*ey(i) + c2(i)*(hz(i-1)-hz(i))
20       continue

c        --Ey SOURCE--
c        --Force Ey source (like antenna feedpoint)--
c        (reflected fields will reflect off this source)
c        ey(inc) = wave(k*dt,0.))

c        --Add Ey source (plane wave)--
c        (also requires Added Hz source)
c        (reflected fields will not reflect off this source,
c         they pass through unperturbed)
         hzinc = wave((k+0.5)*dt, (0.0)*dx)
         ey(inc) = ey(inc) + hzconst*hzinc

c        --Ey BOUNDARY CONDITIONS--
c         --LEFT side--
c         ey(1) = 0. ! metal wall
          RBC_const = (vp*dt-dx)/(vp*dt+dx) ! Mur 1st order
          ey(1) = eyold_2 + RBC_const*(ey(2) - eyold_1)
          eyold_1 = ey(1)
          eyold_2 = ey(2)

c         --RIGHT side--
c         ey(nx) = 0. ! metal
c         --Mur 1st order --
          ey(nx) = eyold_nx1 + RBC_const*(ey(nx-1) - eyold_nx)
          eyold_nx1 = ey(nx-1)
          eyold_nx = ey(nx)

c        --SOLVE FOR Hz--
         do 30 i=1,nx
           hz(i) = hz(i) + d1(i)*(ey(i+1)-ey(i))
30       continue

c        --Hz SOURCE--
         eyinc = (1./eta)*wave((k)*dt,(0.5)*dx)
         hz(inc) = hz(inc) + eyconst*eyinc

c        --NO Hz BOUNDARY CONDITIONS IN THIS MODEL--

c        --OUTPUT--
         call output_time
        
c        --COMPUTE FREQUENCY-DOMAIN FIELDS--
         call peak2E2U ! 2 equation, 2 unknown
         call peak     ! amax
         call peak_dft ! DFT

500   continue

9999  end
c*******************************
c incident waveform
c*******************************
      function wave(t,x)
      implicit complex(c)
      common/wave1/w ! omega
     &            ,B ! beta
     &            ,F ! frequency (Hz)
     &            ,dt ! time step (seconds)
     &            ,ipulse ! =0 for CW, 1 for RC pulse
     &            ,nstep  ! # of time steps in simulation

      if (ipulse.eq.0) then
c     --CW sine wave--
      wave = sin(w*(t-dt) - B*x)
      else

c     --Raised cosine pulse--
      if (t.le.1./F) then
         wave = 1. - cos(w*t-B*x)
      else
         wave = 0.
      endif

c     --Rectanglular pulse--
c     if (t.le.10) then
c        wave = 1.
c     elseif (t.le.20) then
c        wave = -1.
c     else
c        wave = 0.
c     endif
      endif

      return
      end
c***************************
c Ouptut time domain fields
c***************************
      subroutine output_time
      parameter(nx_max=300) ! maximum size of FDTD space
      implicit complex(c)
      common/time_domain/ey(nx_max),hz(nx_max)
      common/time_step/k ! present time step
      common/fdtd_size/nx

c     --Open file--
      if (k.eq.1) then
         open(70,file='hist1')
         open(71,file='hist2')
         open(72,file='hist3')
         open(73,file='hist4')
      endif

c     --Output
      write(70,*) k,ey(14)
      write(71,*) k,ey(20)
      write(72,*) k,ey(30)
      write(73,*) k,ey(40)

      return
      end
c**********************************
c  Magnitude and phase calculations using two equations and
c  two unknowns.  Algorithm assumes a steady state sinusoid.
c  Numerical noise or DC offsets will cause errors.  If the
c  two time steps chosen are too close to numerically equal,
c  numerical round-off errors will cause errors.  This has
c  been a problem only at ultra-low frequencies, and the
c  solution is to pick n1 and n2 further apart.  Ideal accuracy
c  of phase calculations is when n1 and n2 are 1/4 period apar.t
c  This had not been needed, however, for cellular phone calculations,
c  and is expected that it would be needed only for low-observable
c  radar or other applications requiring very high dynamic range
c  of the simulation.  
c
c  Algorithm developed by Cindy Furse, University of Utah
c  3-8-95
c***********************************
      subroutine peak2E2U
      parameter(nx_max=300)
      implicit complex(c)
      common/time_domain/ey(nx_max),hz(nx_max)
      common/wave1/w ! omega
     &            ,B ! beta
     &            ,F ! frequency (Hz)
     &            ,dt ! time step (seconds)
     &            ,ipulse ! =0 for CW, 1 for RC pulse
     &            ,nstep  ! # of time steps in simulation
      common/time_step/k ! current time step
      common/fdtd_size/nx
      common/incident_location/inc
      common/cell_size/dx
c     --Frequency domain fields
      complex cey(nx_max),chz(nx_max)
      static cey,chz
     
c     --Set up samples--
      n1 = nstep-100 ! first sample
      n2 = nstep   ! second sample
      t1 = n1*dt
      t2 = n2*dt

c     --Save first time step
      if (k.eq.n1) then
         do i=1,nx
            cey(i) = ey(i)
            chz(i) = hz(i)
         enddo
      elseif (k.eq.n2) then
c     --Compute complex value
         do i=1,nx
            cey(i) = camag(cey(i),ey(i),t1,t2,w)
            chz(i) = camag(chz(i),hz(i),t1,t2,w)
         enddo
      endif

c     --Output--
      if (k.eq.n2) then
         open(60,file='peak2E2U.out')
         do i=1,nx
            write(60,*) (i-inc)*dx,cabs(cey(i)),cabs(chz(i))
         enddo
         close(60)
      endif

      if (k.eq.n2) then
         write(6,*) 'E14=',cey(14)
      ceinc = cmplx(-.813145,.5919961)
      ceref = cey(14)-ceinc
      cgam = ceref/ceinc
      write(6,*) 'ceref=',ceref
      write(6,*) 'cgam=',cgam
      write(6,*) 'gam=',cabs(cgam)
      endif

      return
      end
c****************************
c Calculate fields by 2 equations, 2 unknowns
c   t1 = first time (seconds)
c   t2 = second time (seconds)
c   q1 = value of sine wave at time t1
c   q2 = value of sine wave at time t2
c   w  = omega, angular frequency (rads)
c   amag = magnitude of sine wave
c   phase = phase of sine wave (radians)
c Algorithm developed by Cindy Furse, University of Utah
c***************************
      function amag(q1,q2,t1,t2,w)
      double precision top,bot,phase,damag
      
      top = q2*sin(w*t1) - q1*sin(w*t2)
      bot = q1*cos(w*t2) - q2*cos(w*t1)
      phase = atan2(top,bot) ! radians
      damag = abs(q1/sin(w*t1 + phase)) ! magnitude
      amag = damag

      return
      end
c****************************
c Calculate fields by 2 equations, 2 unknowns
c   t1 = first time (seconds)
c   t2 = second time (seconds)
c   q1 = value of sine wave at time t1
c   q2 = value of sine wave at time t2
c   w  = omega, angular frequency (rads)
c   amag = magnitude of sine wave
c   phase = phase of sine wave (radians)
c Algorithm developed by Cindy Furse, University of Utah
c***************************
      complex function camag(q1,q2,t1,t2,w)
      implicit complex(c)
      double precision top,bot,phase,damag
      
      cj = cmplx(0.,1.)
      top = q2*sin(w*t1) - q1*sin(w*t2)
      bot = q1*cos(w*t2) - q2*cos(w*t1)
      phase = atan2(top,bot) ! radians
      damag = abs(q1/sin(w*t1 + phase)) ! magnitude
      camag = damag*(cos(phase) + cj*sin(phase))

      return
      end

c***************************
c Calculate magnitude from max
c***************************
      subroutine peak
      parameter(nx_max=300)
      implicit complex(c)
      common/time_domain/ey(nx_max),hz(nx_max)
      common/wave1/w ! omega
     &            ,B ! beta
     &            ,F ! frequency (Hz)
     &            ,dt ! time step (seconds)
     &            ,ipulse ! =0 for CW, 1 for RC pulse
     &            ,nstep  ! # of time steps in simulation
      common/time_step/k ! current time step
      common/fdtd_size/nx
c     --magnitude fields
      real eyp(nx_max),hzp(nx_max)
      static eyp,hzp

c     --initialize--
      if (k.eq.1) then
         do i=1,nx
              eyp(i) = 0.
            hzp(i) = 0.
         enddo
      endif
     
c     --Run for last cycle only--
      To = 1./F ! period of sine wave
      nTo = nint(To/dt) ! number of time steps in one period
      if (k.ge.nstep-nTo+1) then
         do i=1,nx
            eyp(i) = max(eyp(i),abs(ey(i)))
            hzp(i) = max(hzp(i),abs(hz(i)))
         enddo
      endif

c     --Output--
      if (k.eq.nstep) then
          open(60,file='peak.out')
          do i=1,nx
             write(60,*) i,eyp(i),hzp(i)
          enddo
         close(60)
      endif

      return
      end
c**************************
c Calculate magnitude and phase using DFT
c**************************
      subroutine peak_dft
      parameter(nx_max=300) ! maximum size of FDTD space
      implicit complex(c)
      common/time_domain/ey(nx_max),hz(nx_max)
      common/wave1/w ! omega
     &            ,B ! beta
     &            ,F ! frequency (Hz)
     &            ,dt ! time step (seconds)
     &            ,ipulse ! =0 for CW, 1 for RC pulse
     &            ,nstep  ! # of time steps in simulation
      common/time_step/k ! present time step
      common/fdtd_size/nx
      complex cey(nx_max),chz(nx_max) ! frequency domain fields
      static cey,chz

c     --Initialize--
      cj = cmplx(0.,1.)
      if (k.eq.1) then
         do i=1,nx
            cey(i) = cmplx(0.,0.)
            chz(i) = cmplx(0.,0.)
         enddo
      endif

c     --Run for last cycle only-
      To = 1./F ! period of sine wave
      nTo = nint(To/dt) ! number of time steps in one period
      if (k.ge.nstep-nTo+1) then
         ce = (2./nTo)*cexp(cmplx(0.,-w*dt*k))
         do i=1,nx
            cey(i) = cey(i) + ce*ey(i)
            chz(i) = chz(i) + ce*hz(i)
         enddo
      endif

c     --Output--
      if (k.eq.nstep) then
         open(60,file='peak_dft.out')
         do i=1,nx
            write(60,*) i,cabs(cey(i)),cabs(chz(i))
         enddo
         close(60)
      endif

      return
      end