3d-fdtd.txt

三维的有限时域差分源程序

hz(1:ie,1:je)=dahz(1:ie,1:je).*hz(1:ie,1:je)+...
              dbhz(1:ie,1:je).*(ex(1:ie,2:jb)-ex(1:ie,1:je)+...
                                ey(1:ie,1:je)-ey(2:ib,1:je));

hz(is,js)=source(n);


%***********************************************************************
%     Update HZX in PML regions
%***********************************************************************

%     FRONT

hzxbcf(1:iefbc,:)=dahzxbcf(1:iefbc,:).*hzxbcf(1:iefbc,:)-...
  dbhzxbcf(1:iefbc,:).*(eybcf(2:ibfbc,:)-eybcf(1:iefbc,:));

%     BACK

hzxbcb(1:iefbc,:)=dahzxbcb(1:iefbc,:).*hzxbcb(1:iefbc,:)-...
  dbhzxbcb(1:iefbc,:).*(eybcb(2:ibfbc,:)-eybcb(1:iefbc,:));

%     LEFT

hzxbcl(1:iebc-1,:)=dahzxbcl(1:iebc-1,:).*hzxbcl(1:iebc-1,:)-...
  dbhzxbcl(1:iebc-1,:).*(eybcl(2:iebc,:)-eybcl(1:iebc-1,:));
hzxbcl(iebc,:)=dahzxbcl(iebc,:).*hzxbcl(iebc,:)-...
  dbhzxbcl(iebc,:).*(ey(1,:)-eybcl(iebc,:));

%     RIGHT

hzxbcr(2:iebc,:)=dahzxbcr(2:iebc,:).*hzxbcr(2:iebc,:)-...
  dbhzxbcr(2:iebc,:).*(eybcr(3:ibbc,:)-eybcr(2:iebc,:));
hzxbcr(1,:)=dahzxbcr(1,:).*hzxbcr(1,:)-...
  dbhzxbcr(1,:).*(eybcr(2,:)-ey(ib,:));

%***********************************************************************
%     Update HZY in PML regions
%***********************************************************************

%     FRONT

hzybcf(:,1:jebc-1)=dahzybcf(:,1:jebc-1).*hzybcf(:,1:jebc-1)-...
  dbhzybcf(:,1:jebc-1).*(exbcf(:,1:jebc-1)-exbcf(:,2:jebc));
hzybcf(1:iebc,jebc)=dahzybcf(1:iebc,jebc).*hzybcf(1:iebc,jebc)-...
  dbhzybcf(1:iebc,jebc).*(exbcf(1:iebc,jebc)-exbcl(1:iebc,1));
hzybcf(iebc+1:iebc+ie,jebc)=...
  dahzybcf(iebc+1:iebc+ie,jebc).*hzybcf(iebc+1:iebc+ie,jebc)-...
  dbhzybcf(iebc+1:iebc+ie,jebc).*(exbcf(iebc+1:iebc+ie,jebc)-...
                                  ex(1:ie,1));
hzybcf(iebc+ie+1:iefbc,jebc)=...
  dahzybcf(iebc+ie+1:iefbc,jebc).*hzybcf(iebc+ie+1:iefbc,jebc)-...
  dbhzybcf(iebc+ie+1:iefbc,jebc).*(exbcf(iebc+ie+1:iefbc,jebc)-...
                                   exbcr(1:iebc,1));

%     BACK

hzybcb(1:iefbc,2:jebc)=dahzybcb(1:iefbc,2:jebc).*hzybcb(1:iefbc,2:jebc)-...
  dbhzybcb(1:iefbc,2:jebc).*(exbcb(1:iefbc,2:jebc)-exbcb(1:iefbc,3:jbbc));
hzybcb(1:iebc,1)=dahzybcb(1:iebc,1).*hzybcb(1:iebc,1)-...
  dbhzybcb(1:iebc,1).*(exbcl(1:iebc,jb)-exbcb(1:iebc,2));
hzybcb(iebc+1:iebc+ie,1)=...
  dahzybcb(iebc+1:iebc+ie,1).*hzybcb(iebc+1:iebc+ie,1)-...
  dbhzybcb(iebc+1:iebc+ie,1).*(ex(1:ie,jb)-exbcb(iebc+1:iebc+ie,2));
hzybcb(iebc+ie+1:iefbc,1)=...
  dahzybcb(iebc+ie+1:iefbc,1).*hzybcb(iebc+ie+1:iefbc,1)-...
  dbhzybcb(iebc+ie+1:iefbc,1).*(exbcr(1:iebc,jb)-...
                                exbcb(iebc+ie+1:iefbc,2));

%     LEFT

hzybcl(:,1:je)=dahzybcl(:,1:je).*hzybcl(:,1:je)-...
  dbhzybcl(:,1:je).*(exbcl(:,1:je)-exbcl(:,2:jb));

%     RIGHT

hzybcr(:,1:je)=dahzybcr(:,1:je).*hzybcr(:,1:je)-...
  dbhzybcr(:,1:je).*(exbcr(:,1:je)-exbcr(:,2:jb));

%***********************************************************************
%     Visualize fields
%***********************************************************************

if mod(n,4)==0;

timestep=int2str(n);

subplot(3,1,1),pcolor(ex');
shading flat;
caxis([-80.0 80.0]);
axis([1 ie 1 jb]);
colorbar;
axis image;
axis off;
title(['Ex at time step = ',timestep]);

subplot(3,1,2),pcolor(ey');
shading flat;
caxis([-80.0 80.0]);
axis([1 ib 1 je]);
colorbar;
axis image;
axis off;
title(['Ey at time step = ',timestep]);

subplot(3,1,3),pcolor(hz');
shading flat;
caxis([-0.2 0.2]);
axis([1 ie 1 je]);
colorbar;
axis image;
axis off;
title(['Hz at time step = ',timestep]);

nn=n/4;
M(:,nn)=getframe(gcf,rect);

end;

%***********************************************************************
%     END TIME-STEPPING LOOP
%***********************************************************************

end

movie(gcf,M,0,10,rect);






/*      
 * ptfsf-demo-file: a bare-bones 2D code which demonstrates use of the
 *      ptfsf code when the incident field has been written, a priori,
 *      to a data file.
 *
 *      The field at four observations points is recorded to a file.
 *      (If the wrtraw package is installed and used, dumps of the
 *      entire computational domain are made periodically and these
 *      can be used to make 2D color maps of the field.)
 *
 * Copyright (C) 2004 John B. Schneider
 *
 * This code uses the FFTw routines for the Fourier transforms.  See
 * www.fftw.org for that code if you wish to use that code too.
 * Otherwise you will have to replace the calls to the FFTw routines
 * to some other discrete Fourier transform routines.
 *
 * To compile this code, you would use something such as:
 *
 *  gcc -Wall -O2 -c ptfsf.c 
 *  gcc -Wall -O2 ptfsf-demo-file.c -o ptfsf-demo-file ptfsf.o\
 *         -lm -lfftw3 -lfftw3_threads -lpthread 
 *
 * For the GNU C compiler the "-Wall" flag turns on all warnings
 * (always a good idea) and "-O2" gives second-level optimization.
 * You must ensure the included header files are on the search path.
 * If you do not want to use the threaded version of FFTw, you may
 * remove those calls (see the FFTw documentation) and then there is
 * no need to link to the pthread library (which may not be installed
 * on some systems).
 *
 * I have written a suite or routines to generate color-mapped images
 * of fields.  Part of that suite is the wrtraw function.  Since I'm
 * not including that here (contact me if you want it), there is a
 * WRTRAW compiler directive which, if unset, removes all the wrtraw
 * stuff.  So, you can safely ignore that directive.
 *
 *********************************************************************
 * This program is free software; you can redistribute it and/or     *
 * modify it under the terms of the GNU General Public License       *
 * as published by the Free Software Foundation (FSF) version 2      *
 * of the License.                                                   *
 *                                                                   *
 * This program is distributed in the hope that it will be useful,   *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of    *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     *
 * GNU General Public License for more details.                      *
 *                                                                   *
 * You should have received a copy of the GNU General Public License *
 * along with this program; if not, write to the Free Software       *
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA         *
 * 02111-1307, USA.  You may also visit the FSF web site at          *
 * www.fsf.org.  The license under which this software is publish    *
 * is available from www.fsf.org/copyleft/gpl.html or                *
 * www.fsf.org/copyleft/gpl.txt.                                     *
 *********************************************************************
 */ 

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "ptfsf.h"

/* wrtraw stuff */
#ifdef WRTRAW
#include "wrtraw.h"
#endif

/* Size of computational domain. */
#define LIMX 141
#define LIMY 141

/* Macros for allocating and accessing arrays */
#define Ez(I,J) ez[(I)*(LIMY)+(J)]
#define Hx(I,J) hx[(I)*(LIMY-1)+(J)]
#define Hy(I,J) hy[(I)*(LIMY)+(J)]

#define ALLOC_2D(name,nx,ny,type) { \
 name = (type *)calloc((nx)*(ny),sizeof(type)); \
 if (!name) { perror("ALLOC_2D"); \
              fprintf(stderr,"Allocation failed for " #name "\n"); \
              exit(-1); \
             }; \
 }

int main()
{
  const double eta=376.7303662;

  /* Electric and magnetic fields.  The array memory is allocated as a
   * large one-dimensional block, but these arrays are treated as
   * two-dimensional things via the macros given above.
   */
  double cdtds;
  double *ez, *hx, *hy;

  double dteta, doeta; // update equation coefficients
  int idum, jdum, // spatial indices (dummy indices)
    time_end_total,   // total length of simulation
    ntime,        // temporal step
    x_ll, y_ll,   // lower-left corner of TFSF region
    x_ur, y_ur;   // upper-right corner of TFSF region

  /* observation-point stuff  */
  char file_name[80];
  FILE *obs_point;

  /* structure for description of TFSF paramters */
  ptfsf_parameters *params;

  /* wrtraw stuff */
#ifdef WRTRAW
  img_sequence *images=NULL;
  char *basename = "junk";
#endif

  /* Allocate and initialize fields arrays. */
  ALLOC_2D(ez,LIMX,LIMY,double);
  ALLOC_2D(hx,LIMX,LIMY-1,double);
  ALLOC_2D(hy,LIMX-1,LIMY,double);

#ifdef WRTRAW
  images = wrtraw_open2d(0,LIMX-1,1,
             0,LIMY-1,1,
             LIMX,LIMY,
             basename,ez);
  images->verbose = WRTRAW_MAX_INFO;
#endif

  /* Get user-settable parameters. */
  printf("Size of computational domain: %d x %d\n",LIMX,LIMY);

  printf("Enter indices for lower-left corner of TF/SF region: ");
  scanf("%d %d",&x_ll,&y_ll);

  printf("Enter number of time steps overall: ");
  scanf("%d",&time_end_total);

  printf("Enter the output file name: ");
  scanf("%s",file_name);
  obs_point = fopen(file_name,"w");
  if (obs_point == NULL) {
    fprintf(stderr,
        "Observation point output file failed to open.  Terminating...\n");
    exit(-1);
  }
     
  printf("Enter the file name where incident field stored: ");
  scanf("%s",file_name);


  /* 
   * Calculate the needed constants.  Courant number is set here.
   */ 
  cdtds  = 1.0/sqrt(3.0);
  doeta  = cdtds/eta;
  dteta  = cdtds*eta;
      
  /* Initialize perfect total-field/scattered-field code so that it
   * will read incident field from a stored data file. */
  params = ptfsf_init_file(
    x_ll,  y_ll,       // indices of lower-left corner of TF region
    LIMX,  LIMY,       // size of computational domain
    cdtds,             // Courant number
    eta,               // characteristic impedance
    ez,  hx,  hy,      // field arrays
    file_name,      // name of file where incident field stored
    PTFSF_PROGRESS | PTFSF_INFO | PTFSF_ESTIMATE // flags to control behavior
    );

  ptfsf_parameters_display(params);

  x_ur = params->x_ur;
  y_ur = params->y_ur;

  printf("Electric field will be written for these points: "
     "(%d,%d), (%d,%d), (%d,%d)\n",
     x_ll,y_ll,
     (int)(x_ur+x_ll)/2,(int)(y_ur+y_ll)/2,
     x_ur,y_ur);

  printf("Incident field assumed to be non-zero for %d time steps.\n",
     params->time_end);

  /* Do the time stepping. */
  for (ntime=0; ntime<time_end_total; ntime++) {

    if (ntime%10 == 0)
      printf("%d...\n",ntime);

    /* Calculate electric field. */
    for (idum=1; idum<LIMX-1; idum++)
      for (jdum=1; jdum<LIMY-1; jdum++)
    Ez(idum,jdum) +=
      + dteta*(Hy(idum,jdum)-Hy(idum-1,jdum)
                   - (Hx(idum,jdum)-Hx(idum,jdum-1)));

    /* Can uncomment the following to throw in a scatterer to make
     * sure TFSF really transparent.
     */
    // idum = LIMX/2;
    // for (jdum=y_ll+15; jdum<y_ur-15; jdum++)
    //   Ez(idum,jdum) = 0.0;

    /* Electric field is correct after this routines is called, i.e.,
     * total field in the total-field region and scattered field in
     * the scattered-field region.  HOWEVER, the magnetic field is
     * not correct until after the magnetic fields have been
     * updated.
     */
    ptfsf_update_file(ntime);

    /* printf value at observations point to a file */
    fprintf(obs_point,"%i %g %g %g %g\n",ntime,
        Ez(x_ll,y_ll),
        Ez((int)(x_ur+x_ll)/2,(int)(y_ur+y_ll)/2),
        Ez(x_ur,y_ur),
        Ez(x_ur+5,y_ur+5));

    /* Calculate Hx. */
    for(idum=0; idum<LIMX; idum++)
      for(jdum=0; jdum<LIMY-1; jdum++)
    Hx(idum,jdum) -= doeta*(Ez(idum,jdum+1)-Ez(idum,jdum));

    /* Calculate Hy. */
    for(idum=0; idum<LIMX-1; idum++)
      for(jdum=0; jdum<LIMY; jdum++)
    Hy(idum,jdum) += doeta*(Ez(idum+1,jdum)-Ez(idum,jdum));
        
    /* generate the output image */
#ifdef WRTRAW
    if (ntime % 10 == 0)
      wrtraw(images);
#endif
    
  }

  return 0;
}
/* ------------------------- end of main() --------------------------*/