ptfsf.c

The code assumes a two-dimensional computational domain with TMz polarization (i.e., non-zero field

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

  if (flags & PTFSF_ASCII_DUMP)
    in_file = fopen(file_name,"r");
  else
    in_file = fopen(file_name,"rb");
  if (in_file==NULL) {
    perror("ptfsf_init_file");
    fprintf(stderr,"ptfsf_init_file: "
	    "error opening input file %s.  Terminating...\n",file_name);
    exit(-1);
  }

  /* Read "header" stuff from stored file. */
  if (flags & PTFSF_ASCII_DUMP) {
    fscanf(in_file,"%d",&time_end);
    fscanf(in_file,"%d %d",&x_ur_read,&y_ur_read);
    fscanf(in_file,"%d %d",&x_ref_read,&y_ref_read);
    fscanf(in_file,"%lf %lf %lf",&phi_read,&cdtds_read,&eta_read);
  } else {
    fread(&time_end,sizeof(int),1,in_file);
    fread(&x_ur_read,sizeof(int),1,in_file);
    fread(&y_ur_read,sizeof(int),1,in_file);
    fread(&x_ref_read,sizeof(int),1,in_file);
    fread(&y_ref_read,sizeof(int),1,in_file);
    fread(&phi_read,sizeof(double),1,in_file);
    fread(&cdtds_read,sizeof(double),1,in_file);
    fread(&eta_read,sizeof(double),1,in_file);
  }

  x_ur_read += the_x_ll;
  y_ur_read += the_y_ll;
    
  x_ref_read += the_x_ll;
  y_ref_read += the_y_ll;
  
  if (flags & PTFSF_ASCII_DUMP) {
    sprintf(comp_string,"%g",the_cdtds);
    sscanf(comp_string,"%lf",&cdtds_compare);
    sprintf(comp_string,"%g",the_eta);
    sscanf(comp_string,"%lf",&eta_compare);
  } else {
    cdtds_compare = the_cdtds;
    eta_compare = the_eta;
  }


  if (cdtds_read != cdtds_compare) {
    fprintf(stderr,"ptfsf_init_file: "
	    "Courant number used to generate incident field and current\n"
	    "  Courant number disagree.  Stored = %g, current = %g.\n"
	    "  Make sure these were calculated by *exactly* the same\n"
	    "  method.\n",cdtds_read, cdtds_compare);
    error++;
  }
      
  if (eta_read != eta_compare) {
    fprintf(stderr,"ptfsf_init_file: "
	    "Impedance used to generate incident field and current\n"
	    "  Courant number disagree.  Stored = %g, current = %g.\n"
	    "  Make sure these were calculated by *exactly* the same\n"
	    "  method.\n",eta_read, eta_compare);
    error++;
  }

  /* Check arguments and/or set corresponding global variables. */
  error = argument_check("ptfsf_init_file",
			 the_x_ll, the_y_ll, 
			 x_ur_read, y_ur_read,
                         x_ref_read, y_ref_read,
			 lim_x, lim_y,
			 phi_read,
			 the_cdtds, the_eta,
			 the_ez, the_hx, the_hy);

  if (error != 0) {
    fprintf(stderr,"ptfsf_init: terminating...\n");
    exit(-1);
  }

  /* 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 = NULL;
  params->file_name = file_name;
  params->flags = flags;

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


/* ======================== ptfsf_generate_file() ======================== */
/* ptfsf_generate_file: incident fields are calculated and then sent to 
 *      specified output file
 */
void ptfsf_generate_file(
      int time_end,                // time steps incident field non-zero
      int x_size, int y_size,      // size of TF region
      int x_ref, int y_ref,        // indices of "reference" point
      double phi,                  // incident angle [degrees]		      
      double cdtds,                // Courant number
      double eta,                  // characteristic impedance
      double (*time_func)(double), // time-stepping function 
      char *file_name,             // output file name
      int the_flags                // controlling flags
      )
{
  /* call usual "init" function to do the real work -- just make sure
     arguments are set appropriately */
  ptfsf_init(
    time_end,           // number of time steps desired
    0,  0,              // indices of lower-left corner of TF region
    x_size-1, y_size-1, // indices of upper-right corner of TF region
    x_ref,  y_ref,      // indices of "reference" point
    x_size+1, y_size+1, // size of computational domain -- doesn't matter
                        //   but these values keep checker from complaining
    phi,                // incident angle [degrees]		      
    cdtds,              // Courant number
    eta,                // characteristic impedance
    NULL, NULL, NULL,   // field arrays -- don't have any in this case
    time_func,          // time-stepping function 
    the_flags | PTFSF_NULL_OK // make sure NULL field arrays don't cause problem
   );

  /* dump the incident fields to the specified file */
  ptfsf_dump_file(file_name);

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


/* ========================== ptfsf_update_file() ======================== */
/* ptfsf_update_file: function called between update of E and H if
 *      incident fields stored in a file
 */
void ptfsf_update_file(int n_time) {
  int idum, jdum;
  double e0, e1, h0, h1, field[4];  // holders for the input values

  /* If the current time is beyond the "end time", then there is no
   * correction left to do, so just return.
   */
  if (n_time > time_end)
    return;
  else if (n_time == time_end) {
    fclose(in_file);
    return;
  }

  if (flags & PTFSF_ASCII_DUMP) {
    for (idum=x_ll; idum<=x_ur; idum++) {
      fscanf(in_file,"%lf %lf %lf %lf",&e0,&e1,&h0,&h1);
      // DEBUG
      //if (idum==x_ll)
      //printf("debug: %d %g %g %g %g\n",n_time,e0,e1,h0,h1);
      Ez(idum,y_ll) += cdtds_times_eta*h0;
      Ez(idum,y_ur) -= cdtds_times_eta*h1;
      Hx(idum,y_ll-1) += cdtds_over_eta*e0;
      Hx(idum,y_ur) -= cdtds_over_eta*e1;
    }
    for (jdum=y_ll; jdum<=y_ur; jdum++) {
      fscanf(in_file,"%lf %lf %lf %lf",&e0,&e1,&h0,&h1);
      // DEBUG
      //if (jdum==y_ll)
      //	printf("debug: %d %g %g %g %g\n",n_time,e0,e1,h0,h1);
      Ez(x_ll,jdum) -= cdtds_times_eta*h0;
      Ez(x_ur,jdum) += cdtds_times_eta*h1;
      Hy(x_ll-1,jdum) -= cdtds_over_eta*e0;
      Hy(x_ur,jdum) += cdtds_over_eta*e1;
    }

  } else {
    for (idum=x_ll; idum<=x_ur; idum++) {
      fread(field,sizeof(double),4,in_file);
      Ez(idum,y_ll) += cdtds_times_eta*field[2];
      Ez(idum,y_ur) -= cdtds_times_eta*field[3];
      Hx(idum,y_ll-1) += cdtds_over_eta*field[0];
      Hx(idum,y_ur) -= cdtds_over_eta*field[1];
    }
    for (jdum=y_ll; jdum<=y_ur; jdum++) {
      fread(field,sizeof(double),4,in_file);
      Ez(x_ll,jdum) -= cdtds_times_eta*field[2];
      Ez(x_ur,jdum) += cdtds_times_eta*field[3];
      Hy(x_ll-1,jdum) -= cdtds_over_eta*field[0];
      Hy(x_ur,jdum) += cdtds_over_eta*field[1];
    }

  }

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