*  Dave Rossi 05/02/99  
*  fdtd2d is a 2d fdtd program based upon xyfdtd and fdtda. It has
*  the ability to model lossless dielectrics, perfect conductors, excited
*  by defineable electric field sources.  Regions for materials may be
*  defined as rectangles,discs, or triangular patches.  The triangular  
*  patch capability is a practical addition to modeling real geometry.  
*  E-field sources can have step profiles across a rectangular, disc, or 
*  triangular area.  For rectangular and disc shaped sources the E-field
*  profile can be shaped as a cosine, or 1+cosine function in either the
*  X direction, the Y direction, or both X and Y simultaneously (discs). 
*  The E-field xcitation can be either sinusoudal continuous wave, or a
*  sinusoidal half wave impulse.  The program does not include fourier 
*  transform or far field transform capabilities, but those type of 
*  analysis capability programs should be readily available.
*  Results are able to be written to output files for viewing in programs
*  such as GNUPLOT, or plotmtv.  If plotmtv is used for viewing X,Y areas
*  with a Z value of Ez or power flow, the the string $data=curve3d
*  should be placed as the first line of the file.  The output files have 
*  been constructed as generic ascii files to allow more choices of  
*  visualization programs.  The output files for area plots and material 
*  viewing plots can be quite large (10 MB).  The output can be split up
*  by designating sub-regions to be viewed.  Up to 99 files are available
*  for output by the program.  This limit may be severly truncated by 
*  settings of your operating system, or by the compiler program you use.
*  I have found some insidious problems under Linux with unprintable
*  characters getting into the problem definition files.  In addition,
*  the program converts and translates with F2C and the GCC compiler
*  slightly differently under Windows95 and Linux. The problem comes with 
*  character representations.  For use under Linux the <sine> definition
*  in the input file td2d.prb may need to be in a single quote <'sine'>  
*  to be correctly recognized by the program. I have compiled and 
*  run fdtd2d under Red Hat Linux 5.2 with F2C and the GCC compiler,
*  under Windows95 with the CYGWIN runtime DLL and F2C and under 
*  DOS/Windows 95 with F2C and the DJ Delorie version of the GCC 
*  compiler using CWSDPMI as the resident memory manager.  The program is 
*  all lower case and the variables are often longer than 6 characters
*  in violation of rigid FORTRAN standards.  The variable names were meant 
*  to be self explanatory, but can be truncated to 6 characters and will 
*  remain unique.  If you compiler is case sensitive, the program can be 
*  UPCASED with a utility, or in a word processor.
*  Use this program at your own risk!  There are no warranties expressed
*  or implied of its accuracy or suitability for any application.       

       program fdtd2d
       include 'fdtd2com.f'


      call gettheinput
      call checktheinput
      call initialize
      call assignproblemvalues
      call mirrorconstants
      t=0
      print*,'Program Calculations Starting'
      call openplotfiles

* Field feeding

*      pi=3.14159265359
*	speedc=3e8
      rk0=2*pi*freq/c
      betas=7.27e6
      betaa=7.95e6
      gamma5s=sqrt(betas**2-epsrmat(1)*rk0**2)
      gamma5a=sqrt(betaa**2-epsrmat(1)*rk0**2)
      rk2s=sqrt(epsrmat(2)*rk0**2-betas**2)
      rk2a=sqrt(epsrmat(2)*rk0**2-betaa**2)
      write(*,*) gamma5s, gamma5a, rk2s, rk2a
      do 50 im=1,90
      if (im .lt. 30) then
         ez(im,100)=0.25*exp(gamma5s*(im-30)*sizex)
     *    -0.25*exp(gamma5a*(im-30)*sizex)
      elseif ((im .ge. 30) .and. (im .lt. 40)) then
         ez(im,100)=4*cos(rk2s*(im-35)*sizex)
     *                 +4*sin(rk2s*(im-35)*sizex)
     *                 -4*cos(rk2a*(im-35)*sizex)
     *                 -4*sin(rk2a*(im-35)*sizex)
      elseif ((im. ge. 40) .and. (im .lt. 50)) then
         ez(im,100)=2.54*exp(-gamma5s*(im-45)*sizex)
     *            + 2.54*exp(gamma5s*(im-45)*sizex)
     *            + -3*exp(-gamma5a*(im-45)*sizex)
     *            + 3*exp(gamma5a*(im-45)*sizex)
      elseif ((im. ge. 50) .and. (im .lt. 60)) then
         ez(im,100)=4*cos(rk2s*(im-55)*sizex)
     *            -4*sin(rk2s*(im-55)*sizex)
     *            +4*cos(rk2a*(im-55)*sizex)
     *            -4*sin(rk2a*(im-55)*sizex)
      elseif ((im. ge. 60) .and. (im .lt. 90)) then
         ez(im,100)=0.25*exp(-gamma5s*(im-60)*sizex)
     *            + 0.25*exp(-gamma5a*(im-60)*sizex)
      endif
50      continue
*  Variables to add are: betas and betaa


      do 100 n=1,numofdt
        call calcez
        t=t+(dt/2)
        call murboundary
        call calchxy
        t=t+(dt/2)
        counter = n
        call writeplotfiles
        print*,(counter*100/numofdt),'% completed.'

 100   continue
       call closeplotfiles 
      end


************************************************************************
          subroutine gettheinput
          include 'fdtd2com.f'
        open(9,file='echotd2d.prb',status='unknown')
        open(8,file='td2d.prb',status='old')

      read(8,*)freq,amplitude
      write(9,*)freq,amplitude

      read(8,*)wave
      write(9,*)wave

      read(8,*)nxused,sizex,nyused,sizey
      write(9,*)nxused,sizex,nyused,sizey

      read(8,*)numofmats
      write(9,*)numofmats

      do n=1,numofmats
         read(8,*)mater(n),epsrmat(n),murmat(n)
         write(9,*)mater(n),epsrmat(n),murmat(n)
      end do
      read(8,*)numofregions
      write(9,*)numofregions
      do n=1,numofregions
         read(8,*)matdesig(n),shape(n),ri(n,1),rj(n,1),ri(n,2),rj(n,2),
     >            patchtype(n),rr(n)
      write(9,*)matdesig(n),'  ',shape(n),ri(n,1),rj(n,1),ri(n,2),
     >           rj(n,2),' ',patchtype(n),rr(n)

      end do

      read(8,*)dt,numofdt
      write(9,*)dt,numofdt

      read(8,*)numofplotfiles
      write(9,*)numofplotfiles
      do n=1,numofplotfiles
      read(8,*)plottype(n),fieldtype(n),pli(n,1),plj(n,1),pli(n,2),
     >        plj(n,2),tstep(n)
      write(9,*)plottype(n),' ',fieldtype(n),pli(n,1),plj(n,1),
     >          pli(n,2),plj(n,2),tstep(n)
      end do
      close(8)
      close(9)
      print*,'Input File Read and Echoed in file named <Echotd2d.prb>'

      return
      end
************************************************************************

          subroutine checktheinput
          include 'fdtd2com.f'

      ipass = 1
      do 10 i=1,numofmats
       if((sizex.gt.(c/(sqrt(epsrmat(i)*murmat(i))*10*freq))).or.
     >      (sizey.gt.(c/(sqrt(epsrmat(i)*murmat(i))*10*freq)))) then
        ipass = 0
        print*,'The grid size of material ',i,': is too coarse'
       end if

 10    continue

       if ((dt.gt.(sizex/(2*c))).or.(dt.gt.(sizey/(2*c)))) then
        ipass = 0
        print*,'The time step size is too large'
       end if

       if (ipass.eq.1) then
        print*,'Simple stability check passed'
       else
       print*,'Failed Input Check!!'
       end if

      return
      end

************************************************************************
          subroutine initialize
          include 'fdtd2com.f'

* take care of special material numbers 100 & 200

      epsrmat(100)=1.0
      murmat(100)=1.0
      epsrmat(200)=1.0
      murmat(200)=1.0

* now do regular specified materials

      do 20 i=1,nxused
      do 10 j=1,nyused
        epsr(i,j)=1.0
        mur(i,j)=1.0
        flgsrc(i,j)=0
        metlflg(i,j)=0

 10   continue
 20   continue

* calculate constants for Ez,Hx, and Hy

      dtoedx=dt/(sizex*eps0)
      dtoedy=dt/(sizey*eps0)
      dtomdx=dt/(sizex*amu0)
      dtomdy=dt/(sizey*amu0)


* calculate MUR coeficients in 2D

      cxd=(c*dt-sizex)/(c*dt+sizex)
      cyd=(c*dt-sizey)/(c*dt+sizey)
      cxu=cxd
      cyu=cyd
      do 40 i=1,nxused
       do 30 j=1,nyused
        ez(i,j)=0
        hx(i,j)=0
        hy(i,j)=0
  30    continue
  40   continue

      do 60 i=1,nxused
       do 50 j=1,4
       ezy1(i,j)=0
       ezy2(i,j)=0
  50   continue
  60  continue

      do 80 i=1,4
       do 70 j=1,nyused
       ezx1(i,j)=0
       ezx2(i,j)=0
  70   continue
  80  continue

      print*,'Initialization complete'
      return
      end

************************************************************************
          subroutine assignproblemvalues
          include 'fdtd2com.f'


      do 150 n=1,numofregions
       if ((shape(n).eq.'R').or.(shape(n).eq.'r')) then
        do 20 i=ri(n,1),ri(n,2)
        do 10 j=rj(n,1),rj(n,2)
         epsr(i,j)=epsrmat(matdesig(n))
         mur(i,j)= murmat(matdesig(n))
         if (matdesig(n).eq.200) then
          flgsrc(i,j)=matdesig(n)+ n
         end if
         if(matdesig(n).eq.100) then
          metlflg(i,j)=1
         else 
          metlflg(i,j)=0
         end if
 10     continue
 20     continue
       end if
       if ((shape(n).eq.'D').or.(shape(n).eq.'d')) then
        do 40 i=nint(ri(n,1)-rr(n)),nint(ri(n,1)+rr(n))
        do 30 j=nint(rj(n,1)-rr(n)),nint(rj(n,1)+rr(n))
         if (((i.gt.0).and.(i.lt.nxused)).and.((j.gt.0).and.
     >      (j.lt.nyused)).and.( ((ri(n,1)-i)**2+(rj(n,1)-j)**2).le.
     >       (rr(n)**2))) then
           epsr(i,j)=epsrmat(matdesig(n))
           mur(i,j)= murmat(matdesig(n))
           if(matdesig(n).eq.200) then
            flgsrc(i,j)=matdesig(n)+ n
           end if
           if(matdesig(n).eq.100) then
            metlflg(i,j)=1
           else 
            metlflg(i,j)=0
           end if
         end if 
 30     continue
 40     continue
       end if 

***********  Right triangular patches

      if ((shape(n).eq.'P').or.(shape(n).eq.'p')) then

*****  Start with Patch type LR (Lower Right)

         if (patchtype(n).eq.'LR') then
           j=rj(n,1)
           i=ri(n,1)
           epsr(i,j)=epsrmat(matdesig(n))
           mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then
              flgsrc(i,j)=matdesig(n)+ n 
             end if
             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if

          do 60 i=(ri(n,1)+1),ri(n,2)
          do 50 j=rj(n,1),rj(n,2)
           if( (((j-rj(n,1))*1.0)/(i-ri(n,1))).le. 
     >       (((rj(n,2)-rj(n,1))*1.0)/(ri(n,2)-ri(n,1))) ) then
            epsr(i,j)=epsrmat(matdesig(n))
            mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then
              flgsrc(i,j)=matdesig(n)+ n 
             end if
             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if
           end if
 50       continue
 60       continue
         end if

********    Patch type UL (Upper Left)


       if (patchtype(n).eq.'UL') then
           i=ri(n,1)
           do 65 j=rj(n,1),rj(n,2)
           epsr(i,j)=epsrmat(matdesig(n))
           mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then
              flgsrc(i,j)=matdesig(n) + n
             end if
             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if
 65       continue

          do 80 i=(ri(n,1)+1),ri(n,2)
          do 70 j=rj(n,1),rj(n,2)
           if( (((j-rj(n,1))*1.0)/(i-ri(n,1))).ge. 
     >       (((rj(n,2)-rj(n,1))*1.0)/(ri(n,2)-ri(n,1))) ) then
            epsr(i,j)=epsrmat(matdesig(n))
            mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then

              flgsrc(i,j)=matdesig(n)+ n
             end if

             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if
           end if
 70       continue
 80       continue
       end if

********    Patch type LL (Lower Left)

       if (patchtype(n).eq.'LL') then
           i=ri(n,1)
           do 85 j=rj(n,1),rj(n,2)
           epsr(i,j)=epsrmat(matdesig(n))
           mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then
              flgsrc(i,j)=matdesig(n) + n 
             end if
             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if
 85       continue

          do 100 i=(ri(n,1)+1),ri(n,2)
          do 90 j=rj(n,1),rj(n,2)
           if( (((j-rj(n,2))*1.0)/(i-ri(n,1))).le. 
     >       (((rj(n,1)-rj(n,2))*1.0)/(ri(n,2)-ri(n,1))) ) then
            epsr(i,j)=epsrmat(matdesig(n))
            mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then
              flgsrc(i,j)=matdesig(n) + n
             end if
             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if
           end if
 90       continue
 100       continue
       end if

********    Patch type UR (Upper Right)

       if (patchtype(n).eq.'UR') then
           i=ri(n,1)
           j=rj(n,2)
           epsr(i,j)=epsrmat(matdesig(n))
           mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then
              flgsrc(i,j)=matdesig(n) + n
             end if
             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if

          do 120 i=(ri(n,1)+1),ri(n,2)
          do 110 j=rj(n,1),rj(n,2)
           if( (((j-rj(n,2))*1.0)/(i-ri(n,1))).ge. 
     >       (((rj(n,1)-rj(n,2))*1.0)/(ri(n,2)-ri(n,1))) ) then
            epsr(i,j)=epsrmat(matdesig(n))
            mur(i,j)= murmat(matdesig(n))
             if(matdesig(n).eq.200) then
              flgsrc(i,j)=matdesig(n)+ n
             end if
             if(matdesig(n).eq.100) then
              metlflg(i,j)=1
             else 
              metlflg(i,j)=0
             end if
           end if
 110       continue
 120       continue
       end if

      end if


 150  continue
      print*,'Got through Assignproblemvalues'
      return
      end



************************************************************************

          subroutine openplotfiles
          include 'fdtd2com.f'
          character*2 num2str
          external num2str

      do 10 i=1,numofplotfiles
       if ((plottype(i).eq.'t').or.(plottype(i).eq.'T')) then  
        filename(i)='P'//num2str(i)//fieldtype(i)//'.plt'
        print*,'Opening time sequence file',filename(i)
        open((i+10),file=filename(i),status='unknown')  
       end if

       if ((plottype(i).eq.'r').or.(plottype(i).eq.'R')) then
        filename(i)='P'//num2str(i)//'t'//fieldtype(i)//'.plt'
        print*,'Opening area map file ',filename(i)
        print*,' for timestep',tstep(i)
        open((i+10),file=filename(i),status='unknown')     
       end if

       if ((plottype(i).eq.'m').or.(plottype(i).eq.'M')) then
        filename(i)='P'//num2str(i)//'mat.plt'
        print*,'Opening area material map file ',filename(i)
        open((i+10),file=filename(i),status='unknown')  
       end if

       if ((plottype(i).eq.'f').or.(plottype(i).eq.'F')) then
        filename(i)='P'//num2str(i)//fieldtype(i)//'flo.plt'
        print*,'Opening area power flow file ',filename(i)
        print*,' for timestep',tstep(i)
        open((i+10),file=filename(i),status='unknown')
       end if





 10   continue
      return
      end 



************************************************************************


          subroutine calcez
          include 'fdtd2com.f'
      

      do 20 i=2,nxused-1
      do 10 j=2,nyused-1
      if (metlflg(i,j).eq.0) then
       if (flgsrc(i,j).lt.201) then
       ez(i,j)= (ez(i,j)+dtoedx*(((hy(i,j))/epsr(i,j))
     >  -((hy(i-1,j))/epsr(i-1,j)))- dtoedy*(((hx(i,j))/epsr(i,j))
     >  -((hx(i,j-1))/epsr(i,j-1))))
       else
        ez(i,j)=srcez(i,j)
       end if
      end if

       if (metlflg(i,j).eq.1) then
        ez(i,j)=0
       end if  

 10   continue

 20   continue

      return
      end
************************************************************************
          subroutine murboundary
          include 'fdtd2com.f'

* murzx part

      j=2
      ez(1,j)=ezx1(2,j)+cxd*(ez(2,j)-ezx1(1,j))
      ez(nxused,j)=ezx1(3,j)+cxu*(ez(nxused-1,j)-ezx1(4,j))

      j=nyused-1
      ez(1,j)=ezx1(2,j)+cxd*(ez(2,j)-ezx1(1,j))
      ez(nxused,j)=ezx1(3,j)+cxu*(ez(nxused-1,j)-ezx1(4,j))

      do 10  j=3,nyused-2
       ez(1,j)=ezx1(2,j)+cxd*(ez(2,j)-ezx1(1,j))
       ez(nxused,j)=ezx1(3,j)+cxu*(ez(nxused-1,j)-ezx1(4,j))

 10   continue

      do 20 j=2,nyused-1
       ezx2(1,j)=ezx1(1,j)
       ezx2(2,j)=ezx1(2,j)
       ezx2(3,j)=ezx1(3,j)
       ezx2(4,j)=ezx1(4,j)
       ezx1(1,j)=ez(1,j)
       ezx1(2,j)=ez(2,j)
       ezx1(3,j)=ez((nxused-1),j)
       ezx1(4,j)=ez(nxused,j)

 20    continue


*  murzy part
 
      i=2
      ez(i,1)=ezy1(i,2)+cyd*(ez(i,2)-ezy1(i,1))
      ez(1,nyused)=ezy1(i,3)+cyd*(ez(i,nyused-1)-ezy1(i,4))

      i=nxused-1
      ez(i,1)=ezy1(i,2)+cyd*(ez(i,2)-ezy1(i,1))
      ez(i,nyused)=ezy1(i,3)+cyd*(ez(i,nyused-1)-ezy1(i,4))

      do 30  i=3,nxused-2
       ez(i,1)=ezy1(i,2)+cyd*(ez(i,2)-ezy1(i,1))
       ez(i,nyused)=ezy1(i,3)+cyd*(ez(i,nyused-1)-ezy1(i,4))


 30   continue

      do 40 i=2,nxused-1
       ezy2(i,1)=ezy1(i,1)
       ezy2(i,2)=ezy1(i,2)
       ezy2(i,3)=ezy1(i,3)
       ezy2(i,4)=ezy1(i,4)
       ezy1(i,1)=ez(i,1)
       ezy1(i,2)=ez(i,2)
       ezy1(i,3)=ez(i,nyused-1)
       ezy1(i,4)=ez(i,nyused)

 40    continue
 

      return
      end
************************************************************************
          subroutine calchxy
          include 'fdtd2com.f'

      do 20 i=1,nxused
      do 10 j=1,nyused-1
      hx(i,j)=hx(i,j)-dtomdy*((ez(i,j+1)/mur(i,j+1))-(ez(i,j)/mur(i,j)))
 10   continue

 20   continue

      do 40 i=1,nxused-1
      do 30 j=1,nyused
      hy(i,j)=hy(i,j)+dtomdx*((ez(i+1,j)/mur(i+1,j))-(ez(i,j)/mur(i,j)))
 30   continue
 40   continue

 
      return
      end
************************************************************************
          subroutine writeplotfiles
          include 'fdtd2com.f'
          character*2 num2str
          external num2str

      do 100 i=1,numofplotfiles


**Time series Ez,Hx,Hy fields, or XX,YY, Poynting flow or XY total power flow

       if ((plottype(i).eq.'t').or.(plottype(i).eq.'T')) then
        filename(i)='P'//num2str(i)//fieldtype(i)//'.plt'
         if (fieldtype(i).eq.'Ez') then
          write((i+10),*)counter,ez(pli(i,1),plj(i,1))
         end if
         if (fieldtype(i).eq.'Hx') then
          write((i+10),*)counter,hx(pli(i,1),plj(i,1))
         end if
         if (fieldtype(i).eq.'Hy') then
          write((i+10),*)counter,hy(pli(i,1),plj(i,1))
         end if
         if (fieldtype(i).eq.'XX') then
          write((i+10),*)counter,ez(pli(i,1),plj(i,1))*
     >                   hy(pli(i,1),plj(i,1))
         end if
         if (fieldtype(i).eq.'YY') then
          write((i+10),*)counter,ez(pli(i,1),plj(i,1))*
     >                   hx(pli(i,1),plj(i,1))
         end if
         if (fieldtype(i).eq.'XY') then
          write((i+10),*)counter,sqrt((ez(pli(i,1),plj(i,1))*
     >     hy(pli(i,1),plj(i,1)))**2 + (ez(pli(i,1),plj(i,1))*
     >      hx(pli(i,1),plj(i,1)))**2)
         end if
       end if

* Rectangular area outputs

       if (((plottype(i).eq.'R').or.(plottype(i).eq.'r')).and.
     >    (counter.eq.tstep(i))) then
        filename(i)='P'//num2str(i)//'t'//fieldtype(i)//'.plt'
        print*,'Writing slice data to ',filename(i)
        print*,'at timestep ',tstep(i)
        print*,t,' seconds from start'
        do 20 ii=pli(i,1),pli(i,2)
        do 10 jj=plj(i,1),plj(i,2)
         if (fieldtype(i).eq.'Ez') then
          write((i+10),*)ii,jj,ez(ii,jj)
         end if
         if (fieldtype(i).eq.'Hx') then
          write((i+10),*)ii,jj,hx(ii,jj)
         end if
         if (fieldtype(i).eq.'Hy') then
          write((i+10),*)ii,jj,hy(ii,jj)
         end if
 10     continue
 20     continue
       end if

***  Materials inspection area plots


       if (((plottype(i).eq.'M').or.(plottype(i).eq.'m')).and.
     >     (counter.eq.1)) then
        filename(i)='P'//num2str(i)//'mat.plt'
        print*,'Writing material plot to ',filename(i)
        do 40 ii=pli(i,1),pli(i,2)
        do 30 jj=plj(i,1),plj(i,2)
         k=0
          if (flgsrc(ii,jj).ge.200) then
           k=flgsrc(ii,jj)-199 
          end if
         write((i+10),*)ii,jj,(  (epsr(ii,jj)*mur(ii,jj)*
     >          (1-metlflg(ii,jj))) + k )
 
 30     continue
 40     continue

       end if

***  Power flow area plots XX, YY, Poynting or XY total power



       if (((plottype(i).eq.'F').or.(plottype(i).eq.'f')).and.
     >    (counter.eq.tstep(i))) then
        filename(i)='P'//num2str(i)//fieldtype(i)//'flo.plt'
        print*,'Writing power flow data to ',filename(i)
        print*,'at timestep ',tstep(i)
        print*,t,' seconds from start'
        do 60 ii=pli(i,1),pli(i,2)
        do 50 jj=plj(i,1),plj(i,2)
         if (fieldtype(i).eq.'XX') then
          write((i+10),*)ii,jj,(ez(ii,jj)*hy(ii,jj))
         end if
         if (fieldtype(i).eq.'YY') then
          write((i+10),*)ii,jj,(ez(ii,jj)*hx(ii,jj))
         end if
         if (fieldtype(i).eq.'XY') then
          write((i+10),*)ii,jj,sqrt((ez(ii,jj)*hy(ii,jj))**2 
     >                   +(ez(ii,jj)*hx(ii,jj))**2)
         end if

 50     continue
 60     continue

       end if

 100   continue
      return
      end 
************************************************************************
 
          subroutine closeplotfiles
          include 'fdtd2com.f'
          character*2 num2str
          external num2str

      do 10 i=1,numofplotfiles
        filename(i)=' '//num2str(i)
        print*,'Closing file',filename(i)
        close(i+10)  


 10   continue
      return
      end 


************************************************************************



         real function srcez(i,j)
          include 'fdtd2com.f'

      if (wave.eq.'sine') then
       srcez=swave(i,j)

        if (flgsrc(i,j).gt.200) then
         k=flgsrc(i,j)-200
         if (patchtype(k).eq.'SX') then
          srcez=(srcez/2)*(1+cos(2*pi*(i-((ri(k,2)+ri(k,1))/2.0))/
     >               (ri(k,2)-ri(k,1))))
         end if
         if (patchtype(k).eq.'SY')then
         srcez=(srcez/2)*( 1+cos(2*pi*(j-((rj(k,2)+rj(k,1))/2.0))/
     >               (rj(k,2)-rj(k,1))))
         end if
         if (patchtype(k).eq.'SB') then
         srcez=(srcez/4)*(1+cos(2*pi*(i-((ri(k,2)+ri(k,1))/2.0))/
     >             (ri(k,2)-ri(k,1))))*(1+cos(2*pi*(j-((rj(k,2)-
     >               rj(k,1))/2))/(rj(k,2)-rj(k,1))))
         end if
 
         if (patchtype(k).eq.'CX') then
          srcez=(srcez)*(cos(pi*(i-((ri(k,2)+ri(k,1)))/2.0)/
     >               (ri(k,2)-ri(k,1))))
         end if
         if (patchtype(k).eq.'CY')then
         srcez=(srcez)*(cos(pi*(j-((rj(k,2)+rj(k,1))/2.0))/
     >               (rj(k,2)-rj(k,1))))
         end if
         if (patchtype(k).eq.'CB') then
          srcez=(srcez)*(cos(pi*(i-((ri(k,2)+ri(k,1))/2.0))/
     >              (ri(k,2)-ri(k,1))))*(cos(pi*(j-((rj(k,2)+
     >               rj(k,1))/2))/(rj(k,2)-rj(k,1))))
         end if

        end if
       
      else 
       srcez=pulse(i,j)
        if (flgsrc(i,j).gt.200) then
         k=flgsrc(i,j)-200
         if (patchtype(k).eq.'SX') then
          srcez=(srcez/2)*(1+cos(2*pi*(i-((ri(k,2)+ri(k,1))/2.0))/
     >               (ri(k,2)-ri(k,1))))
         end if
         if (patchtype(k).eq.'SY')then
         srcez=(srcez/2)*( 1+cos(2*pi*(j-((rj(k,2)+rj(k,1))/2.0))/
     >               (rj(k,2)-rj(k,1))))
         end if
         if (patchtype(k).eq.'SB') then
         srcez=(srcez/4)*(1+cos(2*pi*(i-((ri(k,2)+ri(k,1))/2.0))/
     >             (ri(k,2)-ri(k,1))))*(1+cos(2*pi*(j-((rj(k,2)-
     >               rj(k,1))/2.0))/(rj(k,2)-rj(k,1))))
         end if
 
         if (patchtype(k).eq.'CX') then
          srcez=(srcez)*(cos(pi*(i-((ri(k,2)+ri(k,1)))/2)/
     >               (ri(k,2)-ri(k,1))))
         end if
         if (patchtype(k).eq.'CY')then
         srcez=(srcez)*(cos(pi*(j-((rj(k,2)+rj(k,1))/2))/
     >               (rj(k,2)-rj(k,1))))
         end if
         if (patchtype(k).eq.'CB') then
          srcez=(srcez)*(cos(pi*(i-((ri(k,2)+ri(k,1))/2))/
     >              (ri(k,2)-ri(k,1))))*(cos(pi*(j-((rj(k,2)+
     >               rj(k,1))/2))/(rj(k,2)-rj(k,1))))
         end if

        end if




      end if
      return
      end

************************************************************************
          real function swave(i,j)
          include 'fdtd2com.f'

      if(flgsrc(i,j).ge.200) then
        swave=amplitude*sin(2*pi*freq*t)
       else swave=0
      end if

      return

      end

************************************************************************
          real function pulse(i,j)
          include 'fdtd2com.f'

      if ((flgsrc(i,j).ge.200).and.(t.le.(1/(2*freq)))) then
       pulse=amplitude*sin(2*pi*freq*t)
       else pulse=0

      end if
      return
      end

************************************************************************
          character*2 function num2str(i)
      num2str='XX'
      do 10 k=0,9
       if(  (i.ge.10*k).and.(i.lt.10*(k+1)) ) then
        num2str= char(k+48)//char(mod(i,10)+48)
       end if
 10   continue
      return
      end

************************************************************************
          subroutine mirrorconstants
          include 'fdtd2com.f'

      print*,'sizex,sizey,nxused,nyused,'
      print*,sizex,sizey,nxused,nyused
      print*,'numofmats,numofregions,numofplotfiles'
      print*,numofmats,numofregions,numofplotfiles
      print*,'freq,amplitude,dt,numofdt'
      print*,freq,amplitude,dt,numofdt 

      return
      end

************************************************************************
















