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