ptfsf.c

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

  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() ---------------------- */


/* ======================== 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. */