psvec.txt

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

c            vector(4) = Ey(2)       |  4  |  5  |  6  |     ^
c            vector(5) = Ex(3)       |     |     |     |     |
c            vector(6) = Ey(3)       -------------------     |
c            vector(7) = Ex(4)       |     |     |     |     ----->  X
c            vector(8) = Ey(4)       |  1  |  2  |  3  |
c            vector(9) = Ex(5)       |     |     |     |
c            vector(10)= Ey(5)       -------------------
c   iplane   = 1 if plot is a x-y cut
c            = 2 if plot is a z-y cut
c            = 3 if plot is a x-z cut
c            for any other value, no axis will be drawn
c   third_comp = .true. if the third component of the cut points out
c                of the page
c              = .false. if the third component of the cut points into
c                the page.
c              EXAMPLE if iplane=1 and third_comp=.true. then the
c              z-axis points out of the page.
c              Note: third_comp is meaningful only if 0 < iplane < 4.
c   title1,title2, titles (string variables) to go on top of vector
c                  plot.
c   cut_val = string to indicate the third value.
c        example: if iplane=3 cut_val would be "30 cm" and in the
c                 postscript file it would appear as "y = 30 cm".
c        Special Feature: if the first character is ~, then nothing
c                         will be printed
c   time_val = string to indicate the time.
c        example: "30.5 ns", "200 steps" would yield "t = 30.5 ns"
c                 and "t = 200 steps" respectively.
c        Special Feature: if the first character is ~, then nothing
c                         will be printed
c   postfile = name of postscript file.
c   first_call = .true. if it is the first time that the subroutine
c                has been called (for a given postfile)
c              = .false. if the subroutine has been called beforei
c                 (for a given postfile).
c   last_call  = .true. if it is the last time that the subroutine
c                will be called (for a given postfile).
c              = .false. if the subroutine will be called again
c                 (for a given postfile).
c              last_call is significant for making a uniform vector
c              scale across all plots.
c   new_page   = .true. if a new page is desired (ignored if
c                 first_call is .true.)
c              = .false. if no new page is desired (i.e. user wants
c                 to place multiple plots on same page).
c   landscape  = .true. if landscape orientation is desired
c              = .false. if portrait orientation is desired.
c              NOTE: this flag is effective only when first_page=.true.
c                    or new_page=.true.
c   psxloc = the x location (in inch) of where on the paper to draw
c            the vectors (center alignment).
c   psyloc = the y location (in inch) of where on the paper to draw
c            the vectors (center alignment).
c   psxsize = constraining size (in inch) in x dircetion of the plot.
c   psysize = constraining size (in inch) in y dircetion of the plot.
c   units = string discribing the units of the vector.
c   grid  = .true. if the fine (dashed) grid is desired.
c         = .false. if the find (dashed) grid is not desired.
c   thick_bound = .true. if a thick boundary is to be plotted
c               = .false. if a thin and dashed boundary is to plotted.
c   glob_scale = true/false for whether or not to create a uniform
c                vector scale for all the plots in the given postscript
c                file. This flag is meaningful
c                only when last_call = .true.
c------------------------------------------------------------------
      subroutine psvector(xpatch,ypatch,xwidth,ywidth,vector,iplane,
     x                    third_comp,ipost,iscratch,
     x                    title1,title2,cut_val,time_val,
     x                    postfile,first_call,last_call,new_page,
     x                    landscape,psxloc,psyloc,psxsize,psysize,
     x                    units,grid,thick_bound,glob_scale)
      implicit  none
      integer   maxelement,iplane,ipost,iscratch
      parameter  (maxelement=25000)
      integer   xpatch,ypatch
      real      xwidth,ywidth,psxloc,psyloc,psxsize,psysize
      real      vector(maxelement)
      character*15 postfile,units,cut_val,time_val
      character*80 title1,title2
      logical   first_call,last_call,new_page,landscape,grid
      logical   glob_scale,third_comp,thick_bound

      integer loop,pagenumber
      real    bigx,bigy,dummy1,dummy2,locscale,sizescale,pi180
      real    dummy3
      integer pulse,const
      character*80 trash
      real    min_e_scale,min_h_scale
      character*15 units2

      pi180 = 4*atan(1.0)/180.0
      pagenumber = 0
      


c check for under-dimensioning.
      if (xpatch*ypatch*2.gt.maxelement) then
         print *, 'DIMENSIONING ERROR IN psvector !!!!!!!!!!!'
         print *, 'PLEASE INCREASE (maxelement) to  at least',
     &            2*xpatch*ypatch
         stop
      end if

c make sure units = 'A/m' or 'V/m' this determines whether data is
c magnetic or electric field
      if ((units.ne.'V/m').and.(units.ne.'A/m')) then
         print *, 'Currently, psvector, is designed to work only',
     &            'with (A/m) or (V/m)'
         print *, 'You have units = ',units
         stop
      end if

c determine the scale for the drawing. const strips
c out one cell from the x and y direction. By setting const to 0 for both the
c electric and magnetic fields, it is effectively defunct.

      const = 0
      if (units.eq.'V/m') then
         if (thick_bound) then
            const = 1
         end if
      end if

      if ((psxsize/(real(xpatch-const)*xwidth)).lt.
     -              (psysize/(real(ypatch-const)*ywidth))) then
         locscale = psxsize/(real(xpatch-const)*xwidth)
      else
         locscale = psysize/(real(ypatch-const)*ywidth)
      end if

c determine the largest field in x and in y

      bigx = 0.0
      bigy = 0.0
      do 10 loop = 1, 2*xpatch*ypatch,2
         if (abs(vector(loop)).gt.bigx) bigx=abs(vector(loop))
         if (abs(vector(loop+1)).gt.bigy) bigy=abs(vector(loop+1))
 10   continue

c determine the scale for the vectors on the drawing
c be careful to avoid a division by 0 error.

      if ((bigx.eq.0.0).and.(bigy.eq.0.0)) then
         sizescale = 1.E+30
      else if (bigx.eq.0.0) then
         sizescale = locscale/bigy *.85 * ywidth
      else if (bigy.eq.0.0) then
         sizescale = locscale/bigx *.85 * xwidth
      else
         if (xwidth/bigx.lt.ywidth/bigy) then
            sizescale = locscale/bigx *.85 * xwidth
         else
            sizescale = locscale/bigy *.85 * ywidth
         end if
      end if

 99      format(a)

c open postscript file, and write vital subroutines.

      if (first_call) then
         open (unit=ipost,file=postfile,status='UNKNOWN')
         write(ipost,99)'%!PS-Adobe-2.0'
         write(ipost,99)'%%BoundingBox: 0 0 612 792'
         write(ipost,99)
         write(ipost,99)'/inch {72 mul} def'
         write(ipost,99)
         write(ipost,99)'/landscape {90 rotate 0 inch -8.5 inch '//
     x              'translate} def'
         write(ipost,99)
         write(ipost,99)'/grid %xctr,yctr,#x_pulses,#y_pulses,'//
     x            'width_x_pulse,width_y_pulse,'
         write(ipost,99)'   %border_line_width,cell_line_width,'//
     x             'grid_flag'
         write(ipost,99)'{ 9 dict begin'
         write(ipost,99)'    newpath'
         write(ipost,99)'%----- declare variables ------'
         write(ipost,99)'    /grid_flag exch def'
         write(ipost,99)'    /cell_line_width exch def'
         write(ipost,99)'    /border_line_width exch def'
         write(ipost,99)'    /width_y_pulse exch inch def'
         write(ipost,99)'    /width_x_pulse exch inch def'
         write(ipost,99)'    /num_y_pulses exch def'
         write(ipost,99)'    /num_x_pulses exch def'
         write(ipost,99)'    inch exch inch exch'
         write(ipost,99)'    newpath'
         write(ipost,99)'    moveto'
         write(ipost,99)'%----- draw the border ------'
         write(ipost,99)'    width_x_pulse num_x_pulses mul -2 div'
         write(ipost,99)'    width_y_pulse num_y_pulses mul -2 div'
         write(ipost,99)'    rmoveto'
         write(ipost,99)'    0 width_y_pulse num_y_pulses mul rlineto'
         write(ipost,99)'    width_x_pulse num_x_pulses mul 0 rlineto'
         write(ipost,99)'    0 width_y_pulse num_y_pulses '//
     x               ' mul -1 mul rlineto'
         write(ipost,99)'    closepath'
         write(ipost,99)'    border_line_width setlinewidth'
         write(ipost,99)'    currentpoint  stroke  moveto'
         write(ipost,99)'%------ draw horizontal lines -------'
         write(ipost,99)' grid_flag 1 eq {'
         write(ipost,99)'  gsave  [.5 2] 0 setdash'
         write(ipost,99)'  cell_line_width setlinewidth'
         write(ipost,99)'  gsave'
         write(ipost,99)'    1 1 num_y_pulses 1 sub'
         write(ipost,99)'    { pop'
         write(ipost,99)'      0 width_y_pulse rmoveto '
         write(ipost,99)'      gsave'
         write(ipost,99)'         width_x_pulse num_x_pulses '//
     x               'mul 0 rlineto stroke'
         write(ipost,99)'      grestore'
         write(ipost,99)'    } for'
         write(ipost,99)'  grestore'
         write(ipost,99)'%------ draw vertical lines --------'
         write(ipost,99)'    1 1 num_x_pulses 1 sub'
         write(ipost,99)'    { pop'
         write(ipost,99)'      width_x_pulse 0 rmoveto'
         write(ipost,99)'      gsave'
         write(ipost,99)'         0 width_y_pulse num_y_pulses '//
     x               'mul rlineto stroke'
         write(ipost,99)'      grestore'
         write(ipost,99)'    } for   grestore'
         write(ipost,99)'} if'
         write(ipost,99)'  end'
         write(ipost,99)'} def'
         write(ipost,99)
         write(ipost,99)'/vector   % x,y,length,angle. (distances'//
     x                  ' are in inches, angles in degrees)'
         write(ipost,99)'{ 4 dict begin'
         write(ipost,99)'    newpath'
         write(ipost,99)'    /angle exch def'
         write(ipost,99)'    /length exch inch def'
         write(ipost,99)'    gsave'
         write(ipost,99)'    inch exch inch exch'
         write(ipost,99)'    moveto '
         write(ipost,99)'    angle rotate'
c removing the following commented line will make the vector center aligned.
c         write(ipost,99)'    length -2 div 0 rmoveto'
         write(ipost,99)'    length .6 mul 0 rlineto'
         write(ipost,99)'    .015 length mul setlinewidth'
         write(ipost,99)'    currentpoint stroke moveto'
         write(ipost,99)'    0 -.15 length mul rlineto'
         write(ipost,99)'    length .4 mul length .15 mul rlineto'
         write(ipost,99)'    length -.4 mul length .15 mul rlineto'
         write(ipost,99)'    closepath'
         write(ipost,99)'    0 setgray fill'
         write(ipost,99)'    grestore'
         write(ipost,99)'  end'
         write(ipost,99)'} def'
         write(ipost,99)
         write(ipost,99)'/draw_axis    %xpos,ypos,iplane,3rd_axis'
         write(ipost,99)'              %3rd_axis = 0/1 then the third'//
     x                               ' axis points in/out of page'
         write(ipost,99)'{7 dict begin'
         write(ipost,99)'/third_axis exch def'
         write(ipost,99)'/iplane exch def'
         write(ipost,99)'/ypos exch def'
         write(ipost,99)'/xpos exch def'
         write(ipost,99)'gsave'
         write(ipost,99)'/len 0.3 def'
         write(ipost,99)
         write(ipost,99)'  xpos ypos translate'
         write(ipost,99)'  newpath  % make the corner'
         write(ipost,99)'  len .015 mul  setlinewidth'
         write(ipost,99)'  2 0 moveto'
         write(ipost,99)'  0 0 lineto'
         write(ipost,99)'  0 2 lineto stroke'
         write(ipost,99)'  0  0 len 90 vector'
         write(ipost,99)'  0  0 len 0 vector'
         write(ipost,99)'  iplane 1 eq{'
         write(ipost,99)'    third_axis 1 eq'
         write(ipost,99)'      {/str1 (x) def'
         write(ipost,99)'       /str2 (y) def}'
         write(ipost,99)'      {/str1 (y) def'
         write(ipost,99)'       /str2 (x) def} ifelse} if'
         write(ipost,99)'  iplane 2 eq{'
         write(ipost,99)'    third_axis 1 eq'
         write(ipost,99)'      {/str1 (y) def'
         write(ipost,99)'       /str2 (z) def}'
         write(ipost,99)'      {/str1 (z) def'
         write(ipost,99)'       /str2 (y) def} ifelse} if'