convo.c

卷积程序

  fclose(file);
}

/****************************************************************
 *	2 numbers each line: a, Rd[a].
 ****/
void 
WriteRd_a(InputStruct * In_Ptr,
	  double *Rd_a)
{
  short       ia, na = In_Ptr->na;
  double      da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Ra");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-s[1/sr]\n", "a[rad]", fname);
  for (ia = 0; ia < na; ia++)
    fprintf(file, "%-12.4E\t%-12.4E\n", (ia + 0.5) * da,
	    Rd_a[ia]);

  fclose(file);
}

/****************************************************************
 *	3 numbers each line:r, a, Tt[r][a]. a = theta.
 ****/
void 
WriteTt_ra(InputStruct * In_Ptr,
	   double **Tt_ra)
{
  short       ir, ia, nr = In_Ptr->nr, na = In_Ptr->na;
  double      r, a, dr = In_Ptr->dr, da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Tra");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[1/(cm2sr)]\n",
	  "r[cm]", "a[rad]", fname);
  for (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    for (ia = 0; ia < na; ia++) {
      a = (ia + 0.5) * da;
      fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n",
	      r, a, Tt_ra[ir][ia]);
    }
  }

  fclose(file);
}

/****
 *	2 numbers each line: r, Tt[r].
 ****/
void 
WriteTt_r(InputStruct * In_Ptr,
	  double *Tt_r)
{
  short       ir, nr = In_Ptr->nr;
  double      dr = In_Ptr->dr;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Tr");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-s[1/cm2]\n", "r[cm]", fname);
  for (ir = 0; ir < nr; ir++)
    fprintf(file, "%-12.4E\t%-12.4E\n", (ir + 0.5) * dr,
	    Tt_r[ir]);

  fclose(file);
}

/****************************************************************
 *	2 numbers each line: theta, Tt[theta].
 ****/
void 
WriteTt_a(InputStruct * In_Ptr,
	  double *Tt_a)
{
  short       ia, na = In_Ptr->na;
  double      da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "Ta");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fprintf(file, "%-12s\t%-s[1/sr]\n", "a[rad]", fname);
  for (ia = 0; ia < na; ia++)
    fprintf(file, "%-12.4E\t%-12.4E\n", (ia + 0.5) * da,
	    Tt_a[ia]);

  fclose(file);
}

/****************************************************************
 *	Write output in M. Keijzer's format so that the file can be
 *	read by the convolution program written by Keijzer in Pascal.
 ****/
void 
WriteKFormat(InputStruct * In_Ptr,
	     OutStruct * Out_Ptr)
{
  short       i, j;
  double      dz, dr;
  FILE       *file;
  char        fname[STRLEN];

  strcpy(fname, "K");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  fputs("output.filename\n", file);
  fprintf(file, "%12hd layers\n", In_Ptr->num_layers);
  fprintf(file, "%12s %12s %12s %12s %12s %12s\n",
	  "layer", "mua", "mus", "g", "nt", "thickness");

  for (i = 1; i <= In_Ptr->num_layers; i++)
    fprintf(file,
	    "%12hd %12.6lf %12.6lf %12.6lf %12.6lf %12.6lf\n",
	    i, In_Ptr->layerspecs[i].mua,
	    In_Ptr->layerspecs[i].mus,
	    In_Ptr->layerspecs[i].g, In_Ptr->layerspecs[i].n,
	    In_Ptr->layerspecs[i].z1 -
	    In_Ptr->layerspecs[i].z0);
  fprintf(file, "%12.6lf index of refraction above\n",
	  In_Ptr->layerspecs[0].n);
  fprintf(file, "%12.6lf index of refraction below\n",
	  In_Ptr->layerspecs[i].n);
  fprintf(file, "\n");

  fprintf(file, "%12ld photons\n", In_Ptr->num_photons);
  fprintf(file, "%12.6lf critical weight\n", In_Ptr->Wth);
  fprintf(file, "%12.6lf depth of boxes micron\n",
	  In_Ptr->dz * 1e4);
  fprintf(file, "%12.6lf width of boxes micron\n",
	  In_Ptr->dr * 1e4);
  fprintf(file, "%12hd number of boxes in z \n", In_Ptr->nz);
  fprintf(file, "%12hd number of boxes in r \n", In_Ptr->nr);
  fprintf(file, "\n");

  fprintf(file,
	  "%12.6lf Total reflection (including direct R)\n",
	  Out_Ptr->Rsp + Out_Ptr->Rd);
  for (i = 1; i <= In_Ptr->num_layers; i++)
    fprintf(file, "%12.6lf Absorbed in layer %12hd\n",
	    Out_Ptr->A_l[i], i);
  fprintf(file, "%12.6lf Total transmission\n", Out_Ptr->Tt);
  fprintf(file, "\n");

  fprintf(file, "Reflectance and Transmission in [cm-2]\n");
  fprintf(file, "Absorption in z-layers in [cm-1]\n");
  fprintf(file, "Absorption in z/r-boxes in [cm-3]\n");

  fprintf(file, "z/r [cm] Layer");
  dr = In_Ptr->dr;
  for (i = 0; i < In_Ptr->nr; i++)
    fprintf(file, "%12.6lf ", i * dr);
  fprintf(file, "\n");

  fprintf(file, "Refl. ");
  for (i = 0; i < In_Ptr->nr; i++)
    fprintf(file, "%12.6lf ", Out_Ptr->Rd_r[i]);
  fprintf(file, "\n");

  dz = In_Ptr->dz;
  for (i = 0; i < In_Ptr->nz; i++) {
    fprintf(file, "%12.6lf %12.6lf", i * dz, Out_Ptr->A_z[i]);
    for (j = 0; j < In_Ptr->nr; j++)
      fprintf(file, "%12.6lf ", Out_Ptr->A_rz[j][i]);
    fprintf(file, "\n");
  }

  fprintf(file, "Transm.");
  for (i = 0; i < In_Ptr->nr; i++)
    fprintf(file, "%12.6lf ", Out_Ptr->Tt_r[i]);
  fprintf(file, "\n");

  fclose(file);
}

/****************************************************************
 ****/
void 
ShowOutMenu(char *in_fname)
{
  printf("I   = Input parameters of mcml\n");
  printf("3   = reflectance, absorption, and transmittance\n");

  printf("AL  = absorption vs layer [-]\n");
  printf("Az  = absorption vs z [1/cm]\n");
  printf("Arz = absorption vs r & z [1/cm3]\n");
  printf("Fz  = fluence vs z [-]\n");
  printf("Frz = fluence vs r & z [1/cm2]\n");

  printf("Rr  = diffuse reflectance vs radius r [1/cm2]\n");
  printf("Ra  = diffuse reflectance vs angle alpha [1/sr]\n");
  printf("Rra = diffuse reflectance vs radius and angle [1/(cm2 sr)]\n");

  printf("Tr  = transmittance vs radius r [1/cm2]\n");
  printf("Ta  = transmittance vs angle alpha [1/sr]\n");
  printf("Tra = transmittance vs radius and angle [1/(cm2 sr)]\n");

  printf("K   = Keijzer's format\n");
  printf("Q   = Quit to main menu\n");
  printf("* input filename: %s \n", in_fname);
}

/****************************************************************
 ****/
void 
BranchOutA(char *Cmd_Str,
	   InputStruct * In_Ptr,
	   OutStruct * Out_Ptr)
{
  switch (toupper(Cmd_Str[1])) {
    case 'L':			/* A_l. */
    WriteA_layer(In_Ptr->num_layers, Out_Ptr->A_l);
    break;
  case 'Z':
    if (toupper(Cmd_Str[2]) == '\0')	/* A_z. */
      WriteA_z(In_Ptr, Out_Ptr->A_z);
    else
      puts("...Wrong command");
    break;
  case 'R':
    if (toupper(Cmd_Str[2]) == 'Z')	/* A_rz. */
      WriteA_rz(In_Ptr, Out_Ptr->A_rz);
    else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void 
BranchOutF(char *Cmd_Str,
	   InputStruct * In_Ptr,
	   OutStruct * Out_Ptr)
{
  switch (toupper(Cmd_Str[1])) {
    case 'Z':
    if (toupper(Cmd_Str[2]) == '\0')	/* F_z. */
      WriteF_z(In_Ptr, Out_Ptr->A_z);
    else
      puts("...Wrong command");
    break;
  case 'R':
    if (toupper(Cmd_Str[2]) == 'Z')	/* F_rz. */
      WriteF_rz(In_Ptr, Out_Ptr->A_rz);
    else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void 
BranchOutR(char *Cmd_Str,
	   InputStruct * In_Ptr,
	   OutStruct * Out_Ptr)
{
  char        ch;

  switch (toupper(Cmd_Str[1])) {
  case 'A':			/* Rd_a. */
    WriteRd_a(In_Ptr, Out_Ptr->Rd_a);
    break;
  case 'R':
    ch = toupper(Cmd_Str[2]);
    if (ch == '\0')		/* Rd_r. */
      WriteRd_r(In_Ptr, Out_Ptr->Rd_r);
    else if (ch == 'A')		/* Rd_ra. */
      WriteRd_ra(In_Ptr, Out_Ptr->Rd_ra);
    else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void 
BranchOutT(char *Cmd_Str,
	   InputStruct * In_Ptr,
	   OutStruct * Out_Ptr)
{
  char        ch;

  switch (toupper(Cmd_Str[1])) {
  case 'A':			/* Tt_a. */
    WriteTt_a(In_Ptr, Out_Ptr->Tt_a);
    break;
  case 'R':
    ch = toupper(Cmd_Str[2]);
    if (ch == '\0')		/* Tt_r. */
      WriteTt_r(In_Ptr, Out_Ptr->Tt_r);
    else if (ch == 'A')		/* Tt_ra. */
      WriteTt_ra(In_Ptr, Out_Ptr->Tt_ra);
    else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void 
BranchOutCmd(char *Cmd_Str,
	     InputStruct * In_Ptr,
	     OutStruct * Out_Ptr)
{
  char        ch;

  switch (toupper(Cmd_Str[0])) {
  case 'I':
    WriteInParm(In_Ptr);
    break;
  case '3':
    WriteRAT(Out_Ptr);
    break;
  case 'K':
    WriteKFormat(In_Ptr, Out_Ptr);
    break;
  case 'A':
    BranchOutA(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'F':
    BranchOutF(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'R':
    BranchOutR(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'T':
    BranchOutT(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'H':
    ShowOutMenu(In_Ptr->in_fname);
    break;
  case 'Q':
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void 
OutputOrigData(InputStruct * In_Ptr,
	       OutStruct * Out_Ptr)
{
  char        cmd_str[STRLEN];

  if (!Out_Ptr->allocated)
    puts("...No data to output");
  else
    do {
      printf("\n> Output mcml data (h for help) => ");
      do
	gets(cmd_str);
      while (!strlen(cmd_str));	/* avoid null string. */
      BranchOutCmd(cmd_str, In_Ptr, Out_Ptr);
    } while (toupper(cmd_str[0]) != 'Q');
}

/****************************Contours***************************/
/****************************************************************
 ****/
void 
ShowContOrigMenu(char *in_fname)
{
  printf("A = absorption vs r & z [1/cm3]\n");
  printf("F = fluence vs r & z [1/cm2]\n");
  printf("R = diffuse reflectance vs radius and angle [1/(cm2 sr)]\n");
  printf("T = transmittance vs radius and angle [1/(cm2 sr)]\n");
  printf("Q   = Quit to main menu\n");
  printf("* input filename: %s \n", in_fname);
}

/****************************************************************
 *	Absorption density to fluence. A = F/mua;
 ****/
void