psvec.txt

matlab编写的关于fdtd编程的步骤和具体方法

!
! in your FDTD code. An example of such a call, as well as a list explaining
! the input variables is given below.
!
! ...........................................................
! The following is an example on how to call the PS Vector Field
!      if (mod(n,20).eq.0 .and. n .lt. 800) then
!        io = n/20 - 1
!        iolast = nstop/20 - 1
!
! The following is for E field components
!        PSFile='Esxup.PS'  ! one file to get the appropriate scale
!        tns = t*1.e+9
!        del1 = dely
!        del2 = delz
!        n1s = 40
!        n1e = 70
!        n2s = 1
!        n2e = nz
!        iskip1 = 1
!        iskip2 = 1
!        Iplane = 2
!        Icut = 13
!        units = 'V/m'
!        grid = .false.
!        totalf =.true.
!        ipost = 13
!        iscratch = 15
!        thick_bound = .true.
!        glob_scale = .true.
!        call PSVF (del1,N1s,N1e,Eys,iskip1,
!     &             del2,N2s,N2e,Ezs,iskip2,
!     &             Nx,Ny,Nz,Iplane,icut,PSfile,ipost,iscratch,totalf,
!     &             tns,delx,dely,delz,io,iolast,units,grid,
!     &             thick_bound,glob_scale)
!
!  The variables are as follows:
!
! parameters:
! del1,del2: FDTD cell size in horizontal and vertical directions (in m)
! Ns1,Ne1:   Starting and ending indices for horizontal vector. The word
!            horizontal reffers to the direction on the postscript plot.
! Ns2,Ne2:   Starting and ending indices for vertical vector. The word
!            vertical reffers to the direction on the postscript plot.
! E1,E2:     FDTD field values. They can be either Electric or Magnetic
!            fields. E1 reffers to the horizontal component, and E2
!            reffers to the vertical component
! Iskip1,Iskip2: Skip values in horizontal and vertical directions
! Nx,Ny,Nz: Actual dimensions of field arrays E1,E2 in FDTD code
! Iplane:   The plane cut (1 for x-y , 2 for y-z , 3 for x-z)
! Icut:     The level of cut in the specified plane cut
! PSFile: File name of the PS output
! Ipost:  Unit number for the PS file (i.e. open(unit=Ipost)). Make sure
!         that this is not an open unit used in the FDTD code.
! Iscrach: Unit number for a scratch file that is internal to the postscript
!          routine. Make sure that this is not an open unit used in the FDTD
!          code.
! totalf: logical variable for the scattered (false) or total (true) field
! tns:    Time in ns
! delx, dely,delz: cell sizes in x,y,z directions (in m)
! io:
! iolast:
! units: This is a string variable that can only be be either 'V/m' or 'A/m'.
!        This option is used to detect if the vector plot is
!        for electric or magnetic field. It is necessary to determine the
!        shift and accurate number of FDTD grid points to be uesd.
! grid:  = .true. if fine grid is desired.
!        = .false. if fine grid is not desired.
!   thick_bound = .true. if a thick boundary is to be plotted
!               = .false. if a thin and dashed boundary is to plotted.
!   glob_scale = true/false for whether or not to create a uniform
!                vector scale for all the plots in the given postscript
!                file. This flag is meaningful
!                only when last_call = .true.
!.....7................................................................2
      subroutine PSVF(del1,N1s,N1e,E1,iskip1,
     &                del2,N2s,N2e,E2,iskip2,
     &                Nx,Ny,Nz,Iplane,icut,PSfile,ipost,iscratch,totalf,
     &                tns,delx,dely,delz,io,iolast,units,grid,
     &                thick_bound,glob_scale)
       
       
! input parameters:
! del1,del2: Steps in FDTD in horizontal and vertical directions
! Ns1,Ne1:   Starting and ending indeces for horizontal vector
! Ns2,Ne2:   Starting and ending indeces for vertical vector
! E1,E2:     Horizontal and vertical vectors
! Iskip1,Iskip2: Skip values in horizontal and vertical directions
! Nx,Ny,Nz: Dimensions of E1,E2 in FDTD code
! Iplane:   Identifier for plane cut (1 for x-y , 2 for y-z , 3 for x-z)
! Icut:     Identifier for the level of cut in specified plane
! PSFile:     File name of the PS output
! Ipost:     Unit number for the PS file

! Output parameters:
! n11,n22: number of grid points in horizontal and vertical directions
! dp1,dp2: steps in horizontal and vertical directions
! vector: e1(1,1,icut),e2(1,1,icut), ...  in sequence
!         up to e1(n11,n22,icut),e2(n11,n22,icut) 
!         for xy plane as an example
! grid:   on/off (true/false)
! totalf: logical variable for the scattered (false) or total (true) field

!      implicit  none   
      parameter  (maxelement=25000)
      character*15 psfile,units,cut_val,time_val
      character*80 title1,title2
      logical   first_call,last_call,new_page,landscape,grid,totalf
      logical   third_comp,glob_scale,thick_bound 
      integer*1 Ipost
      integer   iplane,nx,ny,nz,xpatch,ypatch
      real xwidth,ywidth,psxloc,psyloc,psxsize,psysize
      real vector(maxelement),tns,cut,e1(nx,ny,nz),e2(nx,ny,nz)

      xpatch = (n1e-n1s)/iskip1 + 1
      ypatch = (n2e-n2s)/iskip2 + 1
      xwidth = del1 * iskip1
      ywidth = del2 * iskip2

      maxd = xpatch*ypatch
      if (maxd .gt. maxelement)then
       write(*,*) ' Error in PSVF routine'
       write(*,*) ' Increase the dimension of vector'
       write(*,*) ' from ', maxelement, ' to ', maxd
       stop
      end if

      if (third_comp) then
      k = 1
      jj = n2s - iskip2
      do j = 1,ypatch
        jj = jj + iskip2
        ii = n1s - iskip1
        do i = 1,xpatch
          ii = ii + iskip1
          if (iplane .eq. 1) then           ! x-y plane
            cut=delz*icut
            if(totalf) then
             vector(k)   = e1(ii,jj,icut)+exi(ii,jj,icut)
             vector(k+1) = e2(ii,jj,icut)+eyi(ii,jj,icut)
            else
             vector(k)   = e1(ii,jj,icut)
             vector(k+1) = e2(ii,jj,icut)
            endif
          end if

          if (iplane .eq. 2) then       ! y-z plane
            cut=delx*icut
            vector(k)   = e1(icut,ii,jj)    
            vector(k+1) = e2(icut,ii,jj) 
c            write(7,*) ii,jj,e1(icut,ii,jj),e2(icut,ii,jj)  
          end if
 
          if (iplane .eq. 3) then       ! x-z plane
            cut=dely*icut
            vector(k)   = e1(jj,icut,ii)    ! Ez
            vector(k+1) = e2(jj,icut,ii)    ! Ex
          end if

          k = k + 2
        end do
       end do   

      else

      k = 1
      jj = n2s - iskip2
      do j = 1,ypatch
        jj = jj + iskip2
        ii = n1s - iskip1
        do i = 1,xpatch
          ii = ii + iskip1
          if (iplane .eq. 1) then           ! x-y plane
            cut=delz*icut
            if(totalf) then
             vector(k)   = e1(jj,ii,icut)+exi(jj,ii,icut)
             vector(k+1) = e2(jj,ii,icut)+eyi(jj,ii,icut)
            else
             vector(k)   = e1(jj,ii,icut)
             vector(k+1) = e2(jj,ii,icut)
            endif
          end if

          if (iplane .eq. 2) then       ! y-z plane
            cut=delx*icut
            vector(k)   = e1(icut,jj,ii)    
            vector(k+1) = e2(icut,jj,ii) 
c            write(7,*) ii,jj,e1(icut,jj,ii),e2(icut,jj,ii)  
          end if
 
          if (iplane .eq. 3) then       ! x-z plane
            cut=dely*icut
            vector(k)   = e1(ii,icut,jj)
            vector(k+1) = e2(ii,icut,jj)
          end if

          k = k + 2
        end do
       end do
  
       end if 
! .....
       write(cut_val,100) cut
       cut_val =cut_val//' m'   

      write(time_val,100) tns
      time_val = time_val//' ns' 

 100   format(f7.3)

      if (iplane .eq. 1 ) title1 = 'x-y plane' ! 'BIG TITLE'
      if (iplane .eq. 2 ) title1 = 'y-z plane' ! 'BIG TITLE'
      if (iplane .eq. 3 ) title1 = 'x-z plane' ! 'BIG TITLE'
      title2 = 'Absolute value of field' ! 'small title'

c      psxloc = 4.25    ! for one plot per page
c      psyloc = 5.5
c      psxsize = 7.
c      psysize = 10.

! 4 plots on one page
      psxsize = 3.0
      psysize = 4.0
      if     (mod(io,4).eq.0) then   ! left top 
c          write (*,*) ' mmm ', io,mod(io,4)
        new_page=.true.
        psxloc = 2.125
        psyloc = 8.25
      elseif (mod(io,4).eq.1) then   ! right top
        new_page=.false.
        psxloc = 6.15
        psyloc = 8.25
      elseif (mod(io,4).eq.2) then   ! left bottom
        new_page=.false.
        psxloc = 2.125
        psyloc = 2.75
      elseif (mod(io,4).eq.3) then   !right bottom
        new_page=.false.
        psxloc = 6.15
        psyloc = 2.75
      end if

      if (io .eq. 0) then
        new_page=.true.
        first_call=.true.
        last_call=.false.
      elseif (io .lt. iolast) then
        first_call=.false.
        last_call=.false.
      elseif (io .eq. iolast) then
        first_call=.false.
        last_call=.true.
      end if


      landscape=.false.

      call  psvector(xpatch,ypatch,xwidth,ywidth,vector,iplane,
     x                    third_comp,ipost,iscratch,
     x                    title1,title2,cut_val,time_val,psfile,
     x                    first_call,last_call,new_page,landscape,
     x                    psxloc,psyloc,psxsize,psysize,units,grid,
     x                    thick_bound,glob_scale)

      return
      end
!.....7...............................................................2
!      include 'program2.f'
c------------------------------------------------------------------
c   xpatch = number of vectors in the x direction
c   ypatch = number of vectors in the y direction
c   xwidth = width of a cell in the x direction
c   ywidth = width of a cell in the y direction
c        NOTE: xwidth and ywidth SHOULD be in meters, but they
c              dont have to. The only different is in the legend.
c
c   vector - An array consisting of Ex and Ey as such
c            vector(1) = Ex(1)
c            vector(2) = Ey(1)       -------------------
c            vector(3) = Ex(2)       |     |     |     |       Y