program circle C..... C Computes the fields about a circular PEC that is illuminated C by a TM line source. The program uses "wrtraw" to create C individual raw frames. C C This code is for illustrative purposes only. It was thrown C together in about 15 minutes for the sake of illustrating the C creation of an animation using FDTD code, but don't trust the C FDTD part of this code (even though I hope it's correct, I'm not C willing to spend any time checking it. C C Other than lowecase letters, this should be conforming Fortran C 77 code and should compile with 'f77 circle.f -o circle' C C John B. Schneider C School of Elec. Eng. and Comp. Sci. C Washington State University C Pullman, WA 99164-2752 C schneidj@eecs.wsu.edu C..... integer limx, limy, limxy parameter (limx=500, limy=500) parameter (limxy=500*500) real ez(limx,limy), hx(limx,limy), hy(limx,limy) real doeta, dspace, dtds, dteta, deltim, dtmds, eta, & pconst, pi, rad, xmin, xmax integer ixsrc, iysrc integer idum, iecho, ierr, nhop, iin, iout, ippw, & ihop, nstart, nend, iend, istart, ilead, & jdum, jhop, jend, jstart, & nframe, ntime character base*80 logical inside(limx,limy) common /fields/ ez, hx, hy common /consts/ eta, pi common /iounit/ ierr, iin, iout common /locate/ inside common /minmax/ xmin, xmax common /misccn/ doeta, dteta common /params/ dtds, dtmds, pconst, ippw common /source/ ixsrc, iysrc C.....Initialize field values do 520 idum=1,limx do 510 jdum=1,limy ez(idum,jdum) = 0.e0 hy(idum,jdum) = 0.e0 hx(idum,jdum) = 0.e0 510 continue 520 continue ilead = limx C..... C-----ARRAYS----- C.....Real. C ez = z component of electric field C hx = x component of magnetic field C hx = y component of magnetic field C-----VARIABLES----- C.....Real. C dspace = grid point separation C dtds = ratio of the space increment to the time increment C deltim = distance that light travels in time delta t C dtmds = (deltim-dspace)/(deltim+dspace) C dum1 = dummy variable used as a holder C eta = characteristic impedance of free space C pconst = phase constant C pi = 3.1415... C ippw = wavelength of the incident wave in number of spatial C increments. May specify the points per meter (or provide C some other unitless measure of the discretization). C xmax = maximum of total field C.....Integer. C idum = dummy variable (often used to indicate x location) C iecho = time steps between echoing the current time step C nhop = offset between output time steps C nend = ending time step for which grid values are written C nstart = starting time step for which grid values are written C iend = end value of the rows that are written as output C istart = start value of the rows that are written as output C jend = end values of the columns that are written as output C jstart = start value of the columns that are written as output C limx = number of x grid points (parameter) C limy = number of y grid points (parameter) C ntime = do loop counter (related to current time step) C.....Character. C base = the base name of output files C.....Logical. C inside = true if Ez point is within circle C-----COMMON BLOCKS----- C /consts/ contains constants C /dimens/ contains values that pertain to the grid dimensions C /minmax/ min and max values of the output C /params/ contains the ratio of the spatial dimension to the C time dimension, the wavelength, etc. C-----SUBROUTINES----- C ezstep Updates Ez. C hxstep Updates Hx. C hystep Updates Hy. C wrtfrm Writes the requested information to an output file. C..... C.....Number of time steps between echoing current time step to user. 30 write(ierr,*)'Enter the number of time steps between echoing', & ' the current time step.' write(ierr,*)'(0 = no echo)' read(iin,*)iecho if (iecho .lt. 0) goto 30 C.....Scatterer information. write(6,*)'Enter the radius of scatterer.' read(5,*)rad call ezloc(rad) C.....Incident field info. 40 write(ierr,*)'Enter the number of grid spaces in one wavelength', & ' of the incident' write(ierr,*)'field.' read(iin,*)ippw write(6,*)'Circle is centered at (',limx/2,',',limx/2,')' write(6,*)'Circle has a radius of ',rad write(6,*)'Enter x and y for the source.' read(5,*)ixsrc,iysrc C.....Get time steps information. 50 write(ierr,*)'Enter the time at which output will begin, the', & ' total number of time' write(ierr,*)'steps and the number of time steps between', & ' succesive writing of' write(ierr,*)'output.' read(iin,*)nstart,nend,nhop if (nstart .lt. 0 .or. & nstart .gt. nend .or. & nend .le. 0 .or. & nhop .gt. (nend - nstart)) then write(ierr,1020) 1020 format(1x,'Try again.') goto 50 endif C.....Determine the range over which output is desired. write(iout,1010)'Number of rows in grid: ',limx 1010 format(1x,a,i4) write(iout,1010)'Number of columns in grid:',limy 80 write(ierr,*)'For the output, enter the starting and ending', & ' rows, and modulus.' write(ierr,*)'(To obtain a single row, enter same value', & ' twice).' read(iin,*)istart,iend,ihop if (istart .le. 0 .or. & istart .gt. iend .or. & iend .gt. limx) then write(ierr,1020) goto 80 endif 100 write(ierr,*)'For the output, enter the starting and ending', & ' columns and modulus.' write(ierr,*)'(To obtain a single column, enter same value', & ' twice.)' read(iin,*)jstart,jend,jhop if (jstart .le. 0 .or. & jstart .gt. jend .or. & jend .gt. limy) then write(ierr,1020) goto 100 endif C.....Get output file information. 700 write(ierr,*)'Enter the base name for the output files.' read(iin,'(a)')base C.....Calculate the needed constants. deltim = 1.e0 dspace = 1.4143e0 dtds = deltim/dspace doeta = dtds/eta dteta = dtds*eta dtmds = (deltim-dspace)/(deltim+dspace) pconst = 2.e0*pi/dble(ippw) C.....Initialize the max value to zero. xmax = 0.e0 C.....Do the time stepping. do 530 ntime = 1,nend call ezstep(ntime,ez,hy,hx) call hxstep(ntime,ez,hy,hx) call hystep(ntime,ez,hy,hx) C........Write "permanent" output if appropriate. if (ntime .ge. nstart .and. & mod(ntime-nstart,nhop) .eq. 0) then nframe = (ntime-nstart)/nhop + 1 call wrtraw(nframe,istart,iend,ihop,jstart,jend,jhop, & base,ez,ilead,xmax) endif C........Echo time step if appropriate. if (iecho .ne. 0) then if (mod(ntime,iecho) .eq. 0) & write(iout,*)'Time = ',ntime endif 530 continue write(6,*)'Maximum value of field = ',xmax 2000 stop end C-------------------------- end of main program -------------------------- C########################### subroutine ezstep ########################### subroutine ezstep(ntime,ez,hy,hx) C..... C This subroutine calculates the next values of the pressure C given the current values. C..... integer limx, limy parameter (limx=500, limy=500) real ez(limx,*), hx(limx,*), hy(limx,*) real doeta, dtds, dteta, dtmds, ezhold, & eta, pconst, etime, pi, ricker logical inside(limx,limy) integer ixsrc, iysrc integer idum, ippw, jdum, ntime common /consts/ eta, pi common /locate/ inside common /misccn/ doeta, dteta common /params/ dtds, dtmds, pconst, ippw common /source/ ixsrc, iysrc etime = real(ntime) C.....First, store the corrent field at the source location. ezhold = ez(ixsrc,iysrc) C..... C Calculate Ez -- ignore boundaries. C..... do 520 idum=2,limx do 510 jdum=2,limy if (inside(idum,jdum)) then ez(idum,jdum) = 0.e0 else ez(idum,jdum) = ez(idum,jdum) & + dteta & * ((hy(idum,jdum)-hy(idum-1,jdum)) & - (hx(idum,jdum)-hx(idum,jdum-1))) endif 510 continue 520 continue C.....Calculate the field at the source location. ez(ixsrc,iysrc) = ezhold + ricker(etime) & + dteta*((hy(ixsrc,iysrc)-hy(ixsrc-1,iysrc)) & - (hx(ixsrc,iysrc)-hx(ixsrc,iysrc-1))) return end C----------------------- end of subroutine ezstep ------------------------ C########################### subroutine hxstep ########################### subroutine hxstep(ntime,ez,hy,hx) C..... C This subroutine calculates the next values of the Hx component C of the magnetic field. C..... integer limx, limy parameter (limx=500, limy=500) real ez(limx,*), hx(limx,*), hy(limx,*) real doeta, dtds, dteta, dtmds, eta, pconst, pi integer idum, ippw, jdum, ntime common /consts/ eta, pi common /misccn/ doeta, dteta common /params/ dtds, dtmds, pconst, ippw C..... C Hx in the interior -- ignoring boundary conditions. C..... do 520 jdum=1,limy-1 do 510 idum=1,limx hx(idum,jdum) = hx(idum,jdum) & - doeta*(ez(idum,jdum+1)-ez(idum,jdum)) 510 continue 520 continue return end C----------------------- end of subroutine hxstep ------------------------ C########################### subroutine hystep ########################### subroutine hystep(ntime,ez,hy,hx) C..... C This subroutine calculates the next value of the Hy component of C magnetic field. C..... integer limx, limy parameter (limx=500, limy=500) real ez(limx,*), hx(limx,*), hy(limx,*) real doeta, dtds, dteta, dtmds, eta, pconst, pi integer idum, ippw, jdum, ntime common /consts/ eta, pi common /misccn/ doeta, dteta common /params/ dtds, dtmds, pconst, ippw C..... C Calculate Hy -- ignore boundary conditions. C..... do 520 idum=1,limx-1 do 510 jdum=1,limy hy(idum,jdum) = hy(idum,jdum) & + doeta*(ez(idum+1,jdum)-ez(idum,jdum)) 510 continue 520 continue return end C----------------------- end of subroutine hystep ------------------------ C######################### block data statements ######################### block data blkdat integer limx, limy, limxy parameter (limx=500, limy=500) parameter (limxy=500*500) real eta, pi integer ierr, iin, iout common /consts/ eta, pi common /iounit/ ierr, iin, iout data eta/376.7343096e0/, pi/3.1415926e0/ data ierr/6/, iin/5/, iout/6/ C-----VARIABLES----- C.....real. C eta = characteristic impedance free space C pi = 3.1415... end C---------------------- end of block data statements ---------------------- C############################ function ricker ############################# real function ricker(time) real arg, doeta, dteta, dtds, dtmds, eta, pconst, pi, time integer ippw common /consts/ eta, pi common /misccn/ doeta, dteta common /params/ dtds, dtmds, pconst, ippw arg = pi**2*(dtds*time/dble(ippw)-1)**2 ricker = 100.e0*(1.e0-2.e0*arg)*exp(-arg) return end C--------------------------end of function ricker ------------------------- C########################### subroutine prloc ########################### subroutine ezloc(rad) C..... C This subroutine determines if a pressure point is inside or C outside of the wedge. C..... integer limx, limy parameter (limx=500, limy=500) real rad, dist logical inside(limx,limy) integer ierr, idum, iin, iout, jdum, ixcntr, iycntr common /locate/ inside common /iounit/ ierr, iin, iout ixcntr = limx/2 iycntr = limy/2 C.....To the left of wedge. do 520 jdum=1,limy do 510 idum=1,limx dist = sqrt(real(idum-ixcntr)**2+real(jdum-iycntr)**2) if (dist .gt. rad) then inside(idum,jdum) = .false. else inside(idum,jdum) = .true. endif 510 continue 520 continue return end C----------------------- end of subroutine ezloc ------------------------ C########################## subroutine wrtraw ########################### C.....Subroutine to write raw output files. subroutine wrtraw(nframe,istart,iend,ihop,jstart,jend,jhop, & base,field,ilead,xmax) implicit real (a-h,o-z) real field(ilead,*) character filnam*1024, base*(*), frmt*30 C..... C Local variables: C frmt = internal file used to construct output file name C ilen = number of non-blank characters in the base name C idum,jdum = dummy loop counters C ipnt = record pointer for output file C ndigit = number of digits in current frame number C..... C.....Determine the location of last character in the base name. ilen = index(base,' ') - 1 C.....Calculate the number of digits in the current frame number. ndigit = log10(real(nframe)) + 1 C.....Construct the proper format for the desired output file name. write(frmt,'(a2,i4,a5,i9,a1)')'(a',ilen,',a1,i',ndigit,')' C.....Put this together to obtain the proper output file name. write(filnam,frmt)base,'.',nframe C.....Open the output file and write header info. open(68,file=filnam,status='new',form='unformatted', & access='direct',recl=4,err=2000) write(68,rec=1)real((iend-istart)/ihop + 1) write(68,rec=2)real((jend-jstart)/jhop + 1) C..... C Loop through writing of data. The first entry is for the upper C left corner, the next is for the next column over, etc. C..... ipnt=2 do 20 jdum=jend,jstart,-jhop do 10 idum=istart,iend,ihop ipnt=ipnt + 1 if(abs(field(idum,jdum)).gt.xmax)xmax=abs(field(idum,jdum)) write(68,rec=ipnt)field(idum,jdum) 10 continue 20 continue close(68) return 2000 write(6,*)'Error opening output file in wrtraw.' stop end C---------------------- end of subroutine wrtraw ------------------------