ptfsf.c

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

    /* 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);
    /* Store time-series for later use. */
    for (idum=0; idum<time_end; idum++)
      Time_series_ez_y0(j,idum) = source_xy_temporal[idum]/length;

    /* Along right 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_ur-x_ref)*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);
    /* Store time-series for later use. */
    for (idum=0; idum<time_end; idum++)
      Time_series_ez_y1(j,idum) = source_xy_temporal[idum]/length;
  }

  if (flags & PTFSF_PROGRESS)
    printf("ptfsf_init: "
	   "Generating Hy source terms along y boundaries...\n");
  for(idum=0; idum<length/2+1; idum++)
    impedance[idum] = 
      -cdtds_over_eta*sin(kx_del[idum]/2.0)/sin(pi_m_over_nt[idum]);
  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 = -(pi_m_over_nt[idum] - 
		(x_off+0.5)*kx_del[idum] + (j-y_off)*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_hy_y0(j,idum) = source_xy_temporal[idum]/length;

    /* Along right 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 = -(pi_m_over_nt[idum]+
		(x_ur-x_ref+0.5)*kx_del[idum]+(j-y_off)*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_hy_y1(j,idum) = source_xy_temporal[idum]/length;
  }

  /* free arrays that were used locally */
  free(kx_del);
  free(ky_del);
  free(pi_m_over_nt);
  free(impedance);
  free(source_temporal);
  free(source_xy_temporal);
  free(source_spectral);
  free(source_xy_spectral);

  /* clean up the resources used by FFTW */
  fftw_destroy_plan(planf);
  fftw_destroy_plan(planb);
  fftw_cleanup_threads();

  /* create and fill the parameters structure */
  params = malloc(sizeof(ptfsf_parameters));
  if (params == NULL) {
    fprintf(stderr,"ptfsf_init: "
	    "could not allocate space for parameter structure.\n"
	    "  Terminating...\n");
    exit(-1);
  }

  params->flags = flags;
  params->time_end = time_end;
  params->x_ll = x_ll;
  params->y_ll = y_ll;
  params->x_ur = x_ur;
  params->y_ur = y_ur;
  params->x_ref = x_ref;
  params->y_ref = y_ref;
  params->lim_x = lim_x;
  params->lim_y = lim_y;
  params->phi = phi*180.0/M_PI;
  params->cdtds = cdtds;
  params->eta = eta;
  params->ez = ez;
  params->hx = hx;
  params->hy = hy;
  params->time_func = time_func;
  params->file_name = NULL;
  params->flags = flags;

  return params;
}
/* -------------------------- end of ptfsf_init() ------------------------ */

/* =========================== ptfsf_update() ============================ */
/* ptfsf_update: function to be called at each time step to correct
 *      the fields on the TFSF boundary -- called between E and H update.
 */
void ptfsf_update(int n_time) {
  int idum, jdum;

  /* If the current time is beyond the "end time", then there is no
   * correction left to do, so just return.  If current time is the
   * end time, free the arrays that were used to store the temporal
   * data.
   */
  if (n_time > time_end)
    return;
  else if (n_time == time_end) {
    free(time_series_ez_x0);
    free(time_series_ez_x1);
    free(time_series_ez_y0);
    free(time_series_ez_y1);
    free(time_series_hx_x0);
    free(time_series_hx_x1);
    free(time_series_hy_y0);
    free(time_series_hy_y1);
    return;
  }

  /* Correction for Ez nodes along TF/SF boundary */
  for (idum=x_ll; idum<=x_ur; idum++) {
    Ez(idum,y_ll) += cdtds_times_eta*Time_series_hx_x0(idum-x_ll,n_time);
    Ez(idum,y_ur) -= cdtds_times_eta*Time_series_hx_x1(idum-x_ll,n_time);
  }
  for (jdum=y_ll; jdum<=y_ur; jdum++) {
    Ez(x_ll,jdum) -= cdtds_times_eta*Time_series_hy_y0(jdum-y_ll,n_time);
    Ez(x_ur,jdum) += cdtds_times_eta*Time_series_hy_y1(jdum-y_ll,n_time);
  }    

  /* Correction for Hy along TF/SF boundary. */
  for (jdum=y_ll; jdum<=y_ur; jdum++) {
    Hy(x_ll-1,jdum) -= cdtds_over_eta*Time_series_ez_y0(jdum-y_ll,n_time);
    Hy(x_ur,jdum) += cdtds_over_eta*Time_series_ez_y1(jdum-y_ll,n_time);
  }

  /* Correction for Hx along TF/SF boundary. */
  for (idum=x_ll; idum<=x_ur; idum++) {
    Hx(idum,y_ll-1) += cdtds_over_eta*Time_series_ez_x0(idum-x_ll,n_time);
    Hx(idum,y_ur) -= cdtds_over_eta*Time_series_ez_x1(idum-x_ll,n_time);
  }


  return;
}
/* ------------------------ end of ptfsf_update() ------------------------ */


/* ========================== ptfsf_init_file() ========================== */
/* ptfsf_init_file: initialization function for case where incident
 *      field stored in data file
 */
ptfsf_parameters *ptfsf_init_file(
    int the_x_ll,  int the_y_ll, // indices of lower-left corner of TF region
    int lim_x, int lim_y,        // size of computational domain
    double the_cdtds,            // Courant number
    double the_eta,              // characteristic impedance
    double *the_ez, double *the_hx, double *the_hy, // field arrays
    char *file_name,             // input file name (contains incident field)
    int the_flags              // flags to control (mostly peripheral) behavior
   )
{
  int error=0, x_ur_read, y_ur_read, x_ref_read, y_ref_read;
  double cdtds_read, eta_read, phi_read;
  /* Want to compare stored and currently defined values for Courant
     number and impedance, but have to be careful is ASCII dump was
     used (since that truncates the precision). */
  double cdtds_compare, eta_compare;
  char comp_string[80];

  ptfsf_parameters *params;

  flags = the_flags;