/* ======================== ptfsf_binary_to_ascii() ====================== *//* ptfsf_binary_to_ascii: convert a binary incident-field file to an *      ASCII one */void ptfsf_binary_to_ascii(char *file_name_in, char *file_name_out) {  int idum, jdum, n;  /* holders for the input values */  int time_end_read, x_ur_read, y_ur_read, x_ref_read, y_ref_read;  double phi_read, cdtds_read, eta_read, field[4];  FILE *in, *out;  in=fopen(file_name_in,"rb");  out=fopen(file_name_out,"w");  if (in==NULL || out==NULL) {    perror("ptfsf_binary_to_ascii");    fprintf(stderr,"ptfsf_binary_to_ascii: "	    "error opening input or output file.  Terminating...\n");    exit(-1);  }  fread(&time_end_read,sizeof(int),1,in);  fread(&x_ur_read,sizeof(int),1,in);  fread(&y_ur_read,sizeof(int),1,in);  fread(&x_ref_read,sizeof(int),1,in);  fread(&y_ref_read,sizeof(int),1,in);  fread(&phi_read,sizeof(double),1,in);  fread(&cdtds_read,sizeof(double),1,in);  fread(&eta_read,sizeof(double),1,in);  fprintf(out,"%d\n",time_end_read);  fprintf(out,"%d %d\n",x_ur_read,y_ur_read);  fprintf(out,"%g\n",phi_read);  fprintf(out,"%g %g\n",cdtds_read,eta_read);  for (n=0; n<time_end_read; n++) {    if ((n+1)%100 == 0)       printf("ptfsf_binary_to_ascii: converted through time step %d...\n",n+1);    for (idum=0; idum<=x_ur; idum++) {      fread(field,sizeof(double),4,in);      fprintf(out,"%g %g %g %g\n",field[0],field[1],field[2],field[3]);    }    fprintf(out,"\n");    for (jdum=0; jdum<=y_ur; jdum++) {      fread(field,sizeof(double),4,in);      fprintf(out,"%g %g %g %g\n",field[0],field[1],field[2],field[3]);    }    fprintf(out,"\n");  }  printf("ptfsf_binary_to_ascii: done converting\n");  return;}/* -------------------- end of ptfsf_binary_to_ascii() ------------------- *//* ======================== ptfsf_dump_file() ============================ *//* ptfsf_dump_file: dump the incident fields to a file.  Dump is done *   such that all fields are written for a given time step, then the time *   step is advanced, their dumped, advance, dump... */void ptfsf_dump_file(char *file_name){  int n_time, idum, jdum;  FILE *out;  /* Make sure incident fields have been calculated already or haven't   *  been calculated and freed.   */  if (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) {    fprintf(stderr,"ptfsf_dump_file: "	    "time-series arrays have not been calculated or have\n"	    "  been calculated and freed already.  Output file not\n"	    "  generated.\n");    return;  }  if (flags & PTFSF_ASCII_DUMP)    out = fopen(file_name,"w");  else    out = fopen(file_name,"wb");  if (out == NULL) {    perror("ptfsf_dump_file");    fprintf(stderr,"ptfsf_dump_file: "	    " failed to open file %s.\n",file_name);    return;  }  /* The time-series arrays, because of the nature of the FFT's, store   * the entire temporal history for a point continuously.   * Unfortunately that works against us when we're trying to write   * all the fields out for a given times step (which is to say we   * can't just write entire blocks of contiguous memory).   */  if (flags & PTFSF_ASCII_DUMP) {     printf("ptfsf_dump_file: Writing an ASCII output file.\n"	   "  USING BINARY FORMAT MUCH BETTER!  Consider unsetting\n"	   "  PTFST_ASCII_DUMP flag.\n");    /* Write the number of time steps. */    fprintf(out,"%d\n",time_end);    /* Write size (sort of) of TFSF boundary. */    fprintf(out,"%d %d\n",x_ur,y_ur);    /* Write incident angle.  Radians at this point. */    fprintf(out,"%g\n",phi);    /* Write Courant number and impedance. */    fprintf(out,"%g %g\n",cdtds,eta);    /* Write all the fields for each time step. */    for (n_time=0; n_time<time_end; n_time++) {      /* Fields along horizontal boundaries. */      for (idum=0; idum<=x_ur; idum++)	fprintf(out,"%g %g %g %g\n",		Time_series_ez_x0(idum,n_time),		Time_series_ez_x1(idum,n_time),		Time_series_hx_x0(idum,n_time),		Time_series_hx_x1(idum,n_time));      fprintf(out,"\n");      /* Fields along vertical boundaries. */      for (jdum=0; jdum<=y_ur; jdum++)	fprintf(out,"%g %g %g %g\n",		Time_series_ez_y0(jdum,n_time),		Time_series_ez_y1(jdum,n_time),		Time_series_hy_y0(jdum,n_time),		Time_series_hy_y1(jdum,n_time));      fprintf(out,"\n");    }  } else {    printf("ptfsf_dump_file: Writing a binary output file.\n");    /* Write the number of time steps. */    fwrite(&time_end,sizeof(int),1,out);    /* Write size (well, sort of) of TFSF boundary. */    fwrite(&x_ur,sizeof(int),1,out);    fwrite(&y_ur,sizeof(int),1,out);    /* Write reference point. */    fwrite(&x_ref,sizeof(int),1,out);    fwrite(&y_ref,sizeof(int),1,out);    /* Write incident angle.  Radians at this point. */    fwrite(&phi,sizeof(double),1,out);    /* Write Courant number and impedance. */    fwrite(&cdtds,sizeof(double),1,out);    fwrite(&eta,sizeof(double),1,out);    /* Write all the fields for each time step. */    for (n_time=0; n_time<time_end; n_time++) {      /* Fields along horizontal boundaries. */      for (idum=0; idum<=x_ur; idum++) {	fwrite(&Time_series_ez_x0(idum,n_time),sizeof(double),1,out);	fwrite(&Time_series_ez_x1(idum,n_time),sizeof(double),1,out);	fwrite(&Time_series_hx_x0(idum,n_time),sizeof(double),1,out);	fwrite(&Time_series_hx_x1(idum,n_time),sizeof(double),1,out);      }      /* Fields along vertical boundaries. */      for (jdum=0; jdum<=y_ur; jdum++) {	fwrite(&Time_series_ez_y0(jdum,n_time),sizeof(double),1,out);	fwrite(&Time_series_ez_y1(jdum,n_time),sizeof(double),1,out);	fwrite(&Time_series_hy_y0(jdum,n_time),sizeof(double),1,out);	fwrite(&Time_series_hy_y1(jdum,n_time),sizeof(double),1,out);      }    }  }    fclose(out);  return;}/* ----------------------- end of ptfsf_dump_file() ---------------------- *//* ====================== ptfsf_parameter_display() ====================== *//* ptsfs_parameters_display: display some of the parameters values *      (should add more to this!) */void ptfsf_parameters_display(ptfsf_parameters *params){  printf("TFSF simulation parameters:\n");  printf("  Incident field assumed zero after %d time steps.\n",	 params->time_end);  printf("  Indices of lower-left corner: (%d,%d)\n",	 params->x_ll,params->y_ll);  printf("  Indices of upper-right corner: (%d,%d)\n",	 params->x_ur,params->y_ur);  printf("  Indices of reference point: (%d,%d)\n",	 params->x_ref,params->y_ref);  printf("  Size of computational domain: %d x %d\n",	 params->lim_x,params->lim_y);  printf("  Incident angle: %g\n",params->phi);  printf("  Courant number: %g\n",params->cdtds);  printf("  Impedance: %g\n",params->eta);  return ;}/* ------------------- end of ptfsf_parameter_display() ------------------ *//*   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   !!!!!!  Remaining functions are not visible to the "outside world." !!!!!!   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*//* ========================== argument_check() =========================== *//* argument_check: check the supplied arguments for reasonable values and/or *      set the corresponding global variable */int argument_check(   char *label, // label to use for error messages   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   ){  int error=0;  /* Printer some header nonsense. */  if (!(flags & PTFSF_SILENT)) {    printf("%s: \n"	   "  \"Perfect\" total-field/scattered-field FDTD code.\n\n"	   "  Copyright (C) 2004  John Schneider (schneidj@eecs.wsu.edu)\n\n"	   "  This code is released under the GNU General Public License.\n"	   "  See <http://www.fsf.org/copyleft/gpl.html> for details.\n\n",	   label);    printf("  You must call ptfsf_update() or ptfsf_update_file() between\n"	   "  the electric and the magnetic field updates.  After calling\n"	   "  it, the electric field will be correct throughout the\n"	   "  computational domain.  After the subsequent update of the\n"	   "  magnetic fields, they will be correct throughout the\n"	   "  computational domain.\n\n");    printf("  Note: This version of the code uses the FFTw threaded library\n"	   "  with two threads.  This will be less than optimal on a single\n"	   "  processor machine (but should run fine).  Currently two\n"	   "  threads are used which is ideal for dual-processor machines.\n"	   "  Change the argument of fftw_plan_with_nthreads() in the\n"	   "  ptfsf.c source code if you wish to change this number.\n"	   "  See <http://www.fftw.org> for details concerning the FFTw\n"	   "  routines\n\n");    printf("  To turn off this message, make sure the PTFSF_SILENT flag is\n"	   "  set when you call the ptfsf routines.\n");  }  /* Perform error checking and initialize global variables. */  lim_x = the_lim_x;  lim_y = the_lim_y;  if (lim_x < 0 || lim_y <  0) {    fprintf(stderr,"%s: "	    "field array sizes must have non-negative size.\n",label);    error++;  }  x_ll = the_x_ll;  y_ll = the_y_ll;  if (x_ll < 0 || x_ll >= lim_x ||      y_ll < 0 || y_ll >= lim_y) {    fprintf(stderr,"%s: lower-left corner location illegal.\n",label);    error++;  }  x_ur = the_x_ur;  y_ur = the_y_ur;  if (x_ur < x_ll || x_ur >= lim_x ||      y_ur < y_ll || y_ur >= lim_y) {    fprintf(stderr,"%s: upper-right corner location illegal.\n",label);    error++;  }  if (x_ref < x_ll || x_ref > x_ur ||      y_ref < y_ll || y_ref > x_ur) {    fprintf(stderr,"ptfsf_init: reference point not in TF region.\n");    error++;  }  ez = the_ez;  hx = the_hx;  hy = the_hy;  if ((ez==NULL || hx==NULL || hy==NULL) && !(flags & PTFSF_NULL_OK)) {    fprintf(stderr,"%s: one or more field pointers are NULL\n",label);    error++;  }  phi = the_phi;  if (phi==0.0 || phi==90.0) {    fprintf(stderr,"ptfsf_init: "	    "there are more efficient methods for grid-aligned propagation\n");    fprintf(stderr,"ptfsf_init: "	    "this code not suitable for that that case\n");    error++;  }  if (phi>90.0 || phi<0.0) {    fprintf(stderr,"ptfsf_init: "	    "propagation angle should be between 0 and 90 degrees\n");    error++;  }  /* convert angle to radians */  phi *= M_PI/180.0;  /* Not much we can check with Courant number or impedance. */  cdtds = the_cdtds;  if (cdtds <= 0.0) {    fprintf(stderr,"%s: illegal Courant number\n",label);    error++;  }    eta = the_eta;  if (eta <= 0.0) {    fprintf(stderr,"%s: illegal impedance\n",label);    error++;  }  cdtds_over_eta = cdtds/eta;  cdtds_times_eta = cdtds*eta;  return error;}/* ----------------------- end of argument_check() ----------------------- *//* =========================== wave_numbers() ============================ *//* wave_numbers: generate the wave numbers at all the necessary frequencies */void wave_numbers(int ntot, double cdtds, double phi,		  double *kx_del, double *ky_del){  int i, m_max;  double cos_phi, sin_phi, k_ratio;  cos_phi = cos(phi);  sin_phi = sin(phi);  /* Calculate maximum allowed spectral index that yields a real     solution and chop things off there -- some of the higher     frequencies which may yield real wave number may not be able to     be coupled into the grid.  I'm not sure about that.  Perhaps     something to look into later. */  m_max = ntot/M_PI*asin(cdtds);  for (i=0; i<m_max; i++) {    k_ratio = find_root(i,ntot,cdtds,cos_phi,sin_phi);    kx_del[i] = k_ratio*i*2.0*M_PI*cos_phi/cdtds/ntot;    ky_del[i] = k_ratio*i*2.0*M_PI*sin_phi/cdtds/ntot;  }  for (i=m_max; i<ntot/2+1; i++) {    kx_del[i] = 0.0;    ky_del[i] = 0.0;  }  return;}/* ------------------------ end of wave_numbers() ------------------------ *//* Global variables for the root-finding functions. */double arg_x, arg_y, right_side;/* ============================== find_root() ============================ *//* find_root: Newton's method used to solve dispersion relation. *      Only returns non-zero results for wavenumbers that are purely real. *//* Error tolerance. */#define TOLERANCE (1.0e-11)double find_root(int m, int ntot, double cdtds, double cos_phi, double sin_phi){  double x=0.99;  arg_x = M_PI*m*cos_phi/cdtds/ntot;  arg_y = M_PI*m*sin_phi/cdtds/ntot;  right_side = sin(M_PI*m/(double)ntot)/cdtds;  right_side = right_side*right_side;  /* Check if a real root can be found.  Worst case is grid-aligned     propagation, so do it for that.  */  if (fabs(sin(M_PI*m/(double)ntot)/cdtds) >= 1.0)    return 0.0;  while (fabs(dispersion_relation(x)) > TOLERANCE)    x = x - dispersion_relation(x)/dispersion_relation_prime(x);  return x;}/* ------------------------- end of find_root() -------------------------- *//* ======================== dispersion_relation() ======================== */double dispersion_relation(double x)/* dispersion_relation: the 2D dispersion relation */{  double tmp_x, tmp_y;  tmp_x = sin(x*arg_x);  tmp_y = sin(x*arg_y);  return tmp_x*tmp_x + tmp_y*tmp_y - right_side;}/* --------------------- end of dispersion_relation() -------------------- *//* ===================== dispersion_relation_prime() ===================== *//* dispersion_relation_prime: Derivative of the dispersion relation. */double dispersion_relation_prime(double x){  return sin(2.0*x*arg_x)*arg_x + sin(2.0*x*arg_y)*arg_y;}/* ---------------- end of dispersion_relation_prime() ------------------- *//* !!!!!!!!!!!!! Do not add externally visible functions here. !!!!!!!!!!! */