ptfsf.c

The tar file contains the following files: ptfsf.c: heart of the perfect TFSF code ptfsf.h: he

    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 ([e-mail])\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. !!!!!!!!!!! */