/*
 * 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 impedance
double *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;