/*       * 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() --------------------------*/