/*      
 * 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 frequency


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