/*      
 * ptfsf-file-maker: Generate a file which contains the incident
 *    field suitable for use with the "perfect"
 *    total-field/scattered-field code.
 *
 * 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-file-maker.c -o ptfsf-file-maker 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).
 *
 *********************************************************************
 * 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

/* The tfsf_perfect_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
 * tfsf_perfect_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()
{
  double phi;       // incident angle [degrees]
  int n_end,        // time at which incident field assumed to go to zero
    x_size, y_size, // upper-right corner of TFSF region
    x_ref, y_ref;   // reference point where incident time series exists
  double eta=376.7303662, // impedance
    scale;          // scale factor for impedance

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

  /* timing stuff */
#ifdef TIMING
  int old_seconds, old_useconds;  
  struct rusage tp;
#endif

  /* Get user-settable parameters. */

  printf("Enter horizontal and vertical size of TFSF boundary: ");
  scanf("%d %d",&x_size,&y_size);

  printf("Enter indices for reference point where incident time series\n"
	 "  assumed to be given (should be on or in TFSF boundary).\n"
	 "  Horizontal value should be between 0 and %d, vertical value\n"
	 "  between 0 and %d: ",x_size-1,y_size-1);
  scanf("%d %d",&x_ref,&y_ref);

  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: ");
  scanf("%d",&n_end);

  printf("Using a Ricker wavelet for source function.  Replace the function\n"
	 "  ricker_wrapper() in the source code and recompile if you wish to\n"
	 "  use a different function.\n"
	 "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 multiplier for the Courant number.  Absolute values of\n"
	 "  negative numbers will scale 1/sqrt(3) and positive values will\n"
	 "  scale 1/sqrt(2), e.g., -0.95 yields Courant number of\n"
	 "  0.95/sqrt(3) while 0.5 would yield 0.5/sqrt(2): ");
  scanf("%lf",&cdtds);
  if (cdtds < 0.0)
    cdtds *= -1.0/sqrt(3.0);
  else
    cdtds *= 1.0/sqrt(2.0);
  printf("Using a Courant number of %g.\n",cdtds);

  printf("Impedance currently %.8g.  Enter scale factor for this: ",eta);
  scanf("%lf",&scale);
  eta *= scale;
  printf("Using an impedance of number of %g.\n",eta);

  printf("Enter the output file name: ");
  scanf("%s",file_name);
      
#ifdef TIMING
  getrusage(0,&tp);
  old_seconds = tp.ru_utime.tv_sec;
  old_useconds = tp.ru_utime.tv_usec;
#endif 

  /* have the ptfsf routines do their stuff and generate an output file */
  ptfsf_generate_file(
       n_end,                // time steps incident field non-zero
       x_size, y_size,       // size of TF region
       x_ref, y_ref,         // indices of "reference" point
       phi,                  // incident angle [degrees]		      
       cdtds,                // Courant number
       eta,                  // characteristic impedance
       ricker_wrapper,       // time-stepping function 
       file_name,            // output file name
       PTFSF_PROGRESS | PTFSF_INFO | PTFSF_ESTIMATE // control flags
      );

#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

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


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