ptfsf-demo-file.c

The code assumes a two-dimensional computational domain with TMz polarization (i.e., non-zero field

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