/* * ptfsf: Routines for implementation of a "perfect" *   total-field/scattered-field boundary for the Yee FDTD scheme. * * Copyright (C) 2004  John B. Schneider *  ********************************************************************* * 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 <math.h>#include <stdio.h>#include <stdlib.h>#include <fftw3.h>#include "ptfsf.h"#include "fdtdgen.h"/* local functions not visible to outside world */void wave_numbers(int ntot, double cdtds, double phi,		  double *kx_del, double *ky_del);double find_root(int m, int ntot, double cdtds, double cos_p, double sin_p);double dispersion_relation(double x);double dispersion_relation_prime(double x);int argument_check(   char *label,  // name of function using this check (for reporting purposes)   int the_x_ll, int the_y_ll,    // lower-left corner   int the_x_ur, int the_y_ur,    // upper-right corner   int x_ref, int y_ref,          // reference point   int the_lim_x, int the_lim_y,  // size of computational domain   double the_phi,                // incident angle (degrees)   double the_cdtds,              // Courant number   double the_eta,                // impedance    double *the_ez, double *the_hx, double *the_hy // field arrays   );/* Macros for accessing arrays */#define Ez(I,J) ez[(I)*(lim_y)+(J)]#define Hx(I,J) hx[(I)*(lim_y-1)+(J)]#define Hy(I,J) hy[(I)*(lim_y)+(J)]/* First argument indicates position, second argument indicates time-step. */#define Time_series_ez_x0(I,N) time_series_ez_x0[(I)*time_end+(N)]#define Time_series_ez_x1(I,N) time_series_ez_x1[(I)*time_end+(N)]#define Time_series_ez_y0(J,N) time_series_ez_y0[(J)*time_end+(N)]#define Time_series_ez_y1(J,N) time_series_ez_y1[(J)*time_end+(N)]#define Time_series_hx_x0(I,N) time_series_hx_x0[(I)*time_end+(N)]#define Time_series_hx_x1(I,N) time_series_hx_x1[(I)*time_end+(N)]#define Time_series_hy_y0(J,N) time_series_hy_y0[(J)*time_end+(N)]#define Time_series_hy_y1(J,N) time_series_hy_y1[(J)*time_end+(N)]/* Global variables. *//* length = number of samples in FFT space with a default value of     128.  This value will be increased to a suitable multiple of 2 to     ensure it is longer than the length of the time series. */int length=128,  time_end,     // number of time steps incident field non-zero  x_ll, y_ll,   // indices of lower-left TFSF boundary corner  x_ur, y_ur,   // indices of upper-right TFSF boundary corner  x_ref, y_ref, // reference point  lim_x, lim_y, // size of computational domain  flags;        // flags to control some peripheral behavior (see header file)double   eta,                 // impedance  phi,                 // incident angle [radians!!]  cdtds,               // Courant number  cdtds_over_eta,      // Courant number divided by impedance  cdtds_times_eta;     // Courant number times impedancedouble *ez, *hx, *hy;  // electric and magnetic fields/* Arrays which hold the times series at each node adjacent to TFSF boundary */double *time_series_ez_x0=NULL, *time_series_ez_x1=NULL,  *time_series_ez_y0=NULL, *time_series_ez_y1=NULL,  *time_series_hx_x0=NULL, *time_series_hx_x1=NULL,  *time_series_hy_y0=NULL, *time_series_hy_y1=NULL;FILE *in_file; // input file if fields stored in an array/* ===p========================= ptfsf_init() ============================ *//* ptfsf_init: initialization routine -- calculates incident fields */ptfsf_parameters *ptfsf_init(      int the_time_end,           // time steps inc field non-zero      int the_x_ll, int the_y_ll, // indices of lower-left corner of TF region      int the_x_ur, int the_y_ur, // indices of upper-right corner of TF region      int x_ref, int y_ref,       // indices of "reference" point      int the_lim_x, int the_lim_y, // size of computational domain      double the_phi,             // incident angle [degrees]		            double the_cdtds,           // Courant number      double the_eta,             // characteristic impedance      double *the_ez, double *the_hx, double *the_hy, // field arrays      double (*time_func)(double), // time-stepping function       int the_flags                // collection of flags      ){  int fft_flag=0;  int error=0;  int i, j, idum, delay, x_off, y_off, tf_width, tf_height;  double *kx_del, *ky_del, *pi_m_over_nt, *impedance;  ptfsf_parameters *params;  /* FFT and transfer function stuff. */  double *source_temporal, *source_xy_temporal;  fftw_complex *source_spectral, *source_xy_spectral;  fftw_plan planf, planb;  flags = the_flags;  /* Perform error checking and initialize global variables. */  /* Check the arguments that are common to both the case where the     incident field is written to a file or not. */  error = argument_check("ptfsf_init",			 the_x_ll, the_y_ll, 			 the_x_ur, the_y_ur,			 x_ref, y_ref,			 the_lim_x, the_lim_y,			 the_phi,			 the_cdtds, the_eta,			 the_ez, the_hx, the_hy);  /* The following four constants are used various places. */  x_off = x_ref-x_ll;  // offset from bottom to reference point  y_off = y_ref-y_ll;  // offset from left side to reference point  tf_width  = x_ur - x_ll + 1; // width of TF region  tf_height = y_ur - y_ll + 1; // height of TF region  time_end = the_time_end;  if (time_end < 0) {    fprintf(stderr,"ptfsf_init: "	   "illegal number of time steps (must be positive).\n");    error++;  }  /* increase the FFT length by factors of 2 if necessary */  while (time_end > length)    length *= 2;  if (flags & PTFSF_INFO)    printf("ptfsf_init: FFT's will use a length of %d.\n",length);  /* Not much we can check with Courant number or impedance. */  cdtds = the_cdtds;  if (cdtds <= 0.0) {    fprintf(stderr,"ptfsf_init: illegal Courant number\n");    error++;  }    eta = the_eta;  if (eta <= 0.0) {    fprintf(stderr,"ptfsf_init: illegal impedance\n");    error++;  }  cdtds_over_eta = cdtds/eta;  cdtds_times_eta = cdtds*eta;  if (time_func==NULL) {    fprintf(stderr,"ptfsf_init: illegal time-stepping function\n");    error++;  }  if (error != 0) {    fprintf(stderr,"ptfsf_init: terminating...\n");    exit(-1);  }  /* Allocate and initialize miscellaneous arrays. */  ALLOC_1D_STRING(ptfsf_init:,kx_del,length/2+1,double);  ALLOC_1D_STRING(ptfsf_init:,ky_del,length/2+1,double);  ALLOC_1D_STRING(ptfsf_init:,pi_m_over_nt,length/2+1,double);  ALLOC_1D_STRING(ptfsf_init:,impedance,length/2+1,double);  for (idum=0; idum<length/2+1; idum++)    pi_m_over_nt[idum] = M_PI*idum/(double)length;  /* Allocate FFT-related arrays. */  ALLOC_1D_STRING(ptfsf_init:,source_temporal,length,double);  ALLOC_1D_STRING(ptfsf_init:,source_xy_temporal,length,double);  ALLOC_1D_STRING(ptfsf_init:,source_spectral,length/2+1,fftw_complex);  ALLOC_1D_STRING(ptfsf_init:,source_xy_spectral,length/2+1,fftw_complex);  /* Allocate space for source functions */  ALLOC_1D_STRING(ptfsf_init:,time_series_ez_x0,time_end*tf_width,double);  ALLOC_1D_STRING(ptfsf_init:,time_series_ez_x1,time_end*tf_width,double);  ALLOC_1D_STRING(ptfsf_init:,time_series_hx_x0,time_end*tf_width,double);  ALLOC_1D_STRING(ptfsf_init:,time_series_hx_x1,time_end*tf_width,double);  ALLOC_1D_STRING(ptfsf_init:,time_series_ez_y0,time_end*tf_height,double);  ALLOC_1D_STRING(ptfsf_init:,time_series_ez_y1,time_end*tf_height,double);  ALLOC_1D_STRING(ptfsf_init:,time_series_hy_y0,time_end*tf_height,double);  ALLOC_1D_STRING(ptfsf_init:,time_series_hy_y1,time_end*tf_height,double);    /* Find all the wavenumbers */  if (flags & PTFSF_PROGRESS)    printf("ptfsf_init: Generating wavenumbers...\n");  wave_numbers(length, cdtds, phi, kx_del, ky_del);  if (flags & PTFSF_INFO)    printf("ptfsf_init: \n"	   "  Can estimate optimial FFT using flag PTFSF_ESTIMATE.  This\n"	   "    does not require much overhead to perform the initial FFT,\n"	   "    but the FFT will not be optimal.\n"	   "  Can potentially get faster FFTs using the flag PTFSF_MEASURE,\n"	   "    but with more initial overhead incurred.\n"	   "  The FFTw code will search for the best FFT if you use the flag\n"	   "    PTFSF_PATIENT (but it takes a while).  Overhead associated\n"	   "    with this will likely exceed any gains realized by the\n"	   "    faster FFTs (except perhaps in some extreme circumstances).\n"	   );  if (flags & PTFSF_PATIENT) {    fft_flag = FFTW_PATIENT;    if (flags & PTFSF_INFO)      printf("ptfsf_init: "	     "Using flag PTFSF_PATIENT.  This may take a moment...\n");  } else if (flags & PTFSF_MEASURE) {    fft_flag = FFTW_MEASURE;    if (flags & PTFSF_INFO)      printf("ptfsf_init: Using flag PTFSF_MEASURE.\n");  } else {    fft_flag = FFTW_ESTIMATE;    if (flags & PTFSF_INFO)      printf("ptfsf_init: Using flag PTFSF_ESTIMATE.\n");  }  /* Create the FFT plans used by FFTW -- use threaded version. */  if (flags & PTFSF_PROGRESS)    printf("ptfsf_init: Generating the FFTW plans...\n");  fftw_init_threads();  /* Can change the number of threads by changing the argument of following    * function.  "2" is good for dual-processor machines.   */  fftw_plan_with_nthreads(2);  planf = fftw_plan_dft_r2c_1d(length, source_temporal, source_spectral,			       fft_flag);  planb = fftw_plan_dft_c2r_1d(length, source_xy_spectral, source_xy_temporal,			       fft_flag);  if (planf == NULL || planb == NULL ) {    fprintf(stderr,"ptfsf_init: "	    "FFTW plan allocation failed.  Terminating...\n");    exit(-1);  }  /* Create source time-series. */  /* If reference point does not corresponds to lower left corner,     delay the source pulse a sufficient amount to ensure we get     required spectrum at the desired point.  See "Exact TF/SF     boundar" notes, pages 15 for the calculation of this. */   if (x_ref != x_ll || y_ref != y_ll) {    double dist, multiplier, tmp;    multiplier = 9.0/8.0;    tmp = sin(asin(cdtds)/multiplier);    dist = sqrt((x_ref-x_ll)*(x_ref-x_ll)+(y_ref-y_ll)*(y_ref-y_ll))*      cos(atan2((y_ref-y_ll),(x_ref-x_ll))-phi);        delay = dist/cdtds*sqrt(1.0-tmp*tmp)/cos(asin(tmp/cdtds));    if (flags & PTFSF_INFO) {      printf("ptfsf_init: Since reference point not lower-left\n"	     "  corner, pulse is being delayed %d time steps.\n"	     "  This is based on a distance to the phase front of %g.\n",	     delay,dist);      printf("  Delay assumes group velocity %g times speed of light.\n"	     "  This is determined by velocity of fields discretized at\n"	     "  %g times the minumum discretization for propagating waves.\n"	     "  To change this value, change the variable 'multiplier' in\n"	     "  the code and recompile.\n",	     cos(asin(tmp/cdtds))/sqrt(1.0-tmp*tmp),multiplier);      }  } else    delay = 0;  for(idum=0; idum<delay; idum++)    source_temporal[idum] = 0.0;  for(idum=delay; idum<length; idum++)    source_temporal[idum] = time_func(idum-delay);  /* Find transform of source function. */  fftw_execute(planf);  /* Now, for each point, create the necessary time series. */  if (flags & PTFSF_PROGRESS)    printf("ptfsf_init: "	   "Generating Ez source terms along x boundaries...\n");  for (i=0; i<tf_width; i++) {    /* Along bottom boundary. */    /* Multiply source spectrum and transfer function. */    for(idum=0; idum<length/2+1; idum++) {      double ctmp, stmp, arg;      if (kx_del[idum] != 0.0) {	arg = -((i-x_off)*kx_del[idum] - y_off*ky_del[idum]);	ctmp = cos(arg);	stmp = sin(arg);	source_xy_spectral[idum][0] = 	  ctmp*source_spectral[idum][0] - stmp*source_spectral[idum][1];	source_xy_spectral[idum][1] = 	  stmp*source_spectral[idum][0] + ctmp*source_spectral[idum][1];      }	else {	source_xy_spectral[idum][0] = 0.0;	source_xy_spectral[idum][1] = 0.0;      }    }    /* Take inverse FFT of theoretical spectrum to get (unnormalized)       predicted time series.  Store result back in source_in. */    fftw_execute(planb);    /* Store time-series for later use. */    for (idum=0; idum<time_end; idum++)      Time_series_ez_x0(i,idum) = source_xy_temporal[idum]/length;    /* Along top boundary. */    /* Multiply source spectrum and transfer function. */    for(idum=0; idum<length/2+1; idum++) {      double ctmp, stmp, arg;      if (kx_del[idum] != 0.0) {	arg = -((i-x_off)*kx_del[idum] + (y_ur-y_ref)*ky_del[idum]);	ctmp = cos(arg);	stmp = sin(arg);	source_xy_spectral[idum][0] = 	  ctmp*source_spectral[idum][0] - stmp*source_spectral[idum][1];	source_xy_spectral[idum][1] = 	  stmp*source_spectral[idum][0] + ctmp*source_spectral[idum][1];      }	else {	source_xy_spectral[idum][0] = 0.0;	source_xy_spectral[idum][1] = 0.0;      }    }    /* Take inverse FFT of theoretical spectrum to get (unnormalized)       predicted time series.  Store result back in source_in. */    fftw_execute(planb);    /* Store time-series for later use. */    for (idum=0; idum<time_end; idum++)      Time_series_ez_x1(i,idum) = source_xy_temporal[idum]/length;  }  if (flags & PTFSF_PROGRESS)    printf("ptfsf_init: "	   "Generating Hx source terms along x boundaries...\n");  for(idum=0; idum<length/2+1; idum++)    impedance[idum] =       cdtds_over_eta*sin(ky_del[idum]/2.0)/sin(pi_m_over_nt[idum]);  for (i=0; i<tf_width; i++) {    /* Along bottom boundary. */    /* Multiply source spectrum and transfer function. */    for(idum=0; idum<length/2+1; idum++) {      double ctmp, stmp, arg;      if (kx_del[idum] != 0.0) {	arg = -(pi_m_over_nt[idum] +		(i-x_off)*kx_del[idum] - (y_off+0.5)*ky_del[idum]);	ctmp = cos(arg);	stmp = sin(arg);	source_xy_spectral[idum][0] = impedance[idum]*	  (ctmp*source_spectral[idum][0] - stmp*source_spectral[idum][1]);	source_xy_spectral[idum][1] = impedance[idum]*	  (stmp*source_spectral[idum][0] + ctmp*source_spectral[idum][1]);      }	else {	source_xy_spectral[idum][0] = 0.0;	source_xy_spectral[idum][1] = 0.0;      }    }    /* Take inverse FFT of theoretical spectrum to get (unnormalized)       predicted time series.  Store result back in source_in. */    fftw_execute(planb);    /* Store time-series for later use. */    for (idum=0; idum<time_end; idum++)      Time_series_hx_x0(i,idum) = source_xy_temporal[idum]/length;    /* Along top boundary. */    /* Multiply source spectrum and transfer function. */    for(idum=0; idum<length/2+1; idum++) {      double ctmp, stmp, arg;      if (kx_del[idum] != 0.0) {	arg = -(pi_m_over_nt[idum] +		(i-x_off)*kx_del[idum] + (y_ur-y_ref+0.5)*ky_del[idum]);	ctmp = cos(arg);	stmp = sin(arg);	source_xy_spectral[idum][0] = impedance[idum]*	  (ctmp*source_spectral[idum][0] - stmp*source_spectral[idum][1]);	source_xy_spectral[idum][1] = impedance[idum]*	  (stmp*source_spectral[idum][0] + ctmp*source_spectral[idum][1]);      }	else {	source_xy_spectral[idum][0] = 0.0;	source_xy_spectral[idum][1] = 0.0;      }    }    /* Take inverse FFT of theoretical spectrum to get (unnormalized)       predicted time series.  Store result back in source_in. */    fftw_execute(planb);    /* Store time-series for later use. */    for (idum=0; idum<time_end; idum++)      Time_series_hx_x1(i,idum) = source_xy_temporal[idum]/length;  }  if (flags & PTFSF_PROGRESS)    printf("ptfsf_init: "	   "Generating Ez source terms along y boundaries...\n");  for (j=0; j<tf_height; j++) {    /* Along left boundary. */    /* Multiply source spectrum and transfer function. */    for(idum=0; idum<length/2+1; idum++) {      double ctmp, stmp, arg;      if (ky_del[idum] != 0.0) {	arg = -(-x_off*kx_del[idum] + (j-y_off)*ky_del[idum]);	ctmp = cos(arg);	stmp = sin(arg);	source_xy_spectral[idum][0] = 	  ctmp*source_spectral[idum][0] - stmp*source_spectral[idum][1];	source_xy_spectral[idum][1] = 	  stmp*source_spectral[idum][0] + ctmp*source_spectral[idum][1];      } else {	source_xy_spectral[idum][0] = 0.0;	source_xy_spectral[idum][1] = 0.0;      }    }    /* Take inverse FFT of theoretical spectrum to get (unnormalized)       predicted time series.  Store result back in source_in. */    fftw_execute(planb);