/*       * ptfsf-demo: a bare-bones 2D code which demonstrates use of the ptfsf code. * *     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.c -o ptfsf-demo 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). * * If you want timing data reported, add the -DTIMING directive to the * second compile command. * * 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"/* Timing stuff. */#ifdef TIMING  #include <sys/time.h>  #include <sys/resource.h>#endif/* 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); \             }; \ }/* The ptfsf_init() function is passed a time-series function * that takes a single argument (the time step) and returns a double * (the incident field at that time step).  I usually use a Ricker * wavelet function with three arguments: the time-step, the Courant * number, and the points per wavelength at the most energetic * frequency.  To use the usual Ricker function, we pass * ptfsf_init() a wrapper which only has one argument and then * calls the usual Ricker routine with the missing arguments supplied. */double ricker_wrapper(double ntime);double ricker(double time, double cdtds, double ppw);/* global variables -- for the sake of getting the wrapper to work. */double cdtds, // Courant number  ppw;        // points per wavelength at most energetic frequencyint 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 *ez, *hx, *hy;  double phi,     // incident angle [degrees]    dteta, doeta; // update equation coefficients  int idum, jdum, // spatial indices (dummy indices)    nend_tfsf,    // time at which incident field assumed to go to zero    nend_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    x_ref, y_ref; // reference point where incident time series exists  /* observation-point stuff  */  char filnam[80];  FILE *obs_point;  /* wrtraw stuff */#ifdef WRTRAW  img_sequence *images=NULL;  char *basename = "junk";#endif  /* timing stuff */#ifdef TIMING  int old_seconds, old_useconds;    struct rusage tp;#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 indices for upper-right corner of TF/SF region: ");  scanf("%d %d",&x_ur,&y_ur);  printf("Enter indices for reference point where incident time series\n"	 "  assumed to be given (should be on or in TFSF boundary): ");  scanf("%d %d",&x_ref,&y_ref);  printf("Enter number of time steps overall: ");  scanf("%d",&nend_total);  printf("Enter time step when the TFSF turns off, i.e., the time step\n"	 "  at which the incident field is essentially zero over the TFSF\n"	 "  boundary (should be no larger than %d): ",nend_total);  scanf("%d",&nend_tfsf);  printf("Enter the points per wavelength at peak of Ricker spectrum: ");  scanf("%lf",&ppw);  printf("Enter the incident angle (should be between 0 and 90) [degrees]: ");  scanf("%lf",&phi);  printf("Enter the output file name: ");  scanf("%s",filnam);  obs_point = fopen(filnam,"w");  if (obs_point == NULL) {    fprintf(stderr,	    "Observation point output file failed to open.  Terminating...\n");    exit(-1);  }  /*    * Calculate the needed constants.   * Courant number is set here.   */   cdtds  = 1.0/sqrt(3.0);  doeta  = cdtds/eta;  dteta  = cdtds*eta;      #ifdef TIMING  getrusage(0,&tp);  old_seconds = tp.ru_utime.tv_sec;  old_useconds = tp.ru_utime.tv_usec;#endif   /* Calculate the incident field using the "perfect" TFSF code. */  ptfsf_init(nend_tfsf,    x_ll,  y_ll,       // indices of lower-left corner of TF region    x_ur, y_ur,        // indices of upper-right corner of TF region    x_ref,  y_ref,     // indices of "reference" point    LIMX,  LIMY,       // size of computational domain    phi,               // incident angle [degrees]		          cdtds,             // Courant number    eta,               // characteristic impedance    ez,  hx,  hy,      // field arrays    ricker_wrapper,    // time-stepping function     PTFSF_PROGRESS | PTFSF_INFO | PTFSF_ESTIMATE // flags to control behavior   );#ifdef TIMING  getrusage(0,&tp);  printf("Calculation of incident field took %.3f seconds.\n",	 tp.ru_utime.tv_sec-old_seconds +	 (tp.ru_utime.tv_usec-old_useconds)/1.e6);#endif  printf("Electric field 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);  /* Do the time stepping. */  for (ntime=0; ntime<nend_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(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      }#ifdef TIMING  getrusage(0,&tp);  printf("Total run time: %.3f seconds.\n",	 tp.ru_utime.tv_sec + tp.ru_utime.tv_usec/1.e6);#endif  return 0;}/* ------------------------- end of main() --------------------------*//* ######################## ricker_wrapper() ####################### *//* A trivial wrapper to be able to call my usual ricker function with * a single argument. */double ricker_wrapper(double ntime) {  return ricker(ntime,cdtds,ppw);}/* ############################ ricker() ########################### *//* Ricker wavelet. */double ricker(double time,  // time step 	      double cdtds, // Courant number	      double ppw    // points/wavelength at most energetic frequency	      ) {  double arg,    arg_max = 70.0, // arguments beyond this value are assumed to yield zero                    //         this allows us to avoid calling exponential                    //         when result will be effectively zero    delay = 2.0;    // delay = multiple of inverse of most energetic frequency,                    //         i.e., multiple of period at that frequency  arg = pow(M_PI*((cdtds*time)/ppw - delay),2);  if (arg > arg_max) {    return 0.0;  } else {    return (1.0 - 2.0*arg) * exp(-arg);  }}/* ------------------------- end of ricker() ----------------------- */