convconv.c

蒙特卡罗算法。用盟特卡罗算法模拟光在组织钟的传播

double
FlatIntegration(float (*Func) (float))
{
  double      rc = ConvVar.r;
  double      R = ConvVar.in_ptr->beam.R;
  double      b_max = (ConvVar.in_ptr->nr - 0.5) * ConvVar.in_ptr->dr;
  double      a = MAX(0, rc - R);
  double      b = MIN(b_max, rc + R);

  if (a >= b)
    return (0);
  else
    return (qtrap(Func, a, b, ConvVar.in_ptr->eps));
}

/****************************************************************
 ****/
double
GaussIntegration(float (*Func) (float))
{
  double      rc = ConvVar.r;
  double      R = ConvVar.in_ptr->beam.R;
  double      b_max = (ConvVar.in_ptr->nr - 0.5) * ConvVar.in_ptr->dr;
  double      a = MAX(0, rc - GAUSSLIMIT * R);
  double      b = MIN(b_max, rc + GAUSSLIMIT * R);

  if (a >= b)
    return (0);
  else
    return (qtrap(Func, a, b, ConvVar.in_ptr->eps));
}

/****************************************************************
 ****/
void
ConvA_rz(InputStruct * In_Ptr,
	 OutStruct * Out_Ptr)
{
  short       irc, iz;
  double      rc, P = In_Ptr->beam.P, R = In_Ptr->beam.R;

  puts("The convolution may take a little while. Wait...");
  for (irc = 0; irc < In_Ptr->nrc; irc++) {
    rc = (irc + 0.5) * In_Ptr->drc;
    ConvVar.r = rc;
    ConvVar.tree = NULL;	/* init the tree. */
    for (iz = 0; iz < In_Ptr->nz; iz++) {
      ConvVar.iz = iz;
      if (In_Ptr->beam.type == flat)
	Out_Ptr->A_rzc[irc][iz] = 2 * P / (R * R)
	  * FlatIntegration(A_rzFGIntegrand);
      else			/* Gaussian. */
	Out_Ptr->A_rzc[irc][iz] = 4 * P / (R * R)
	  * GaussIntegration(A_rzFGIntegrand);
    }
    FreeTree(ConvVar.tree);
  }

  Out_Ptr->conved.A_rz = 1;
}

/****************************************************************
 ****/
void
ConvRd_ra(InputStruct * In_Ptr,
	  OutStruct * Out_Ptr)
{
  short       irc, ia;
  double      rc, P = In_Ptr->beam.P, R = In_Ptr->beam.R;
  double      b_max = (In_Ptr->nr - 1) * In_Ptr->dr;

  puts("The convolution may take a little while. Wait...");
  for (irc = 0; irc < In_Ptr->nrc; irc++) {
    rc = (irc + 0.5) * In_Ptr->drc;
    ConvVar.r = rc;
    ConvVar.tree = NULL;	/* init the tree. */
    for (ia = 0; ia < In_Ptr->na; ia++) {
      ConvVar.ia = ia;
      if (In_Ptr->beam.type == flat)
	Out_Ptr->Rd_rac[irc][ia] = 2 * P / (R * R)
	  * FlatIntegration(Rd_raFGIntegrand);
      else			/* Gaussian. */
	Out_Ptr->Rd_rac[irc][ia] = 4 * P / (R * R)
	  * GaussIntegration(Rd_raFGIntegrand);
    }
    FreeTree(ConvVar.tree);
  }

  Out_Ptr->conved.Rd_ra = 1;
}

/****************************************************************
 ****/
void
ConvRd_r(InputStruct * In_Ptr, OutStruct * Out_Ptr)
{
  short       irc;
  double      rc, P = In_Ptr->beam.P, R = In_Ptr->beam.R;
  double      b_max = (In_Ptr->nr - 1) * In_Ptr->dr;

  for (irc = 0; irc < In_Ptr->nrc; irc++) {
    rc = (irc + 0.5) * In_Ptr->drc;
    ConvVar.r = rc;
    if (In_Ptr->beam.type == flat)
      Out_Ptr->Rd_rc[irc] = 2 * P / (R * R)
	* FlatIntegration(Rd_rFGIntegrand);
    else			/* Gaussian. */
      Out_Ptr->Rd_rc[irc] = 4 * P / (R * R)
	* GaussIntegration(Rd_rFGIntegrand);
  }

  Out_Ptr->conved.Rd_r = 1;
}

/****************************************************************
 ****/
void
ConvTt_ra(InputStruct * In_Ptr, OutStruct * Out_Ptr)
{
  short       irc, ia;
  double      rc, P = In_Ptr->beam.P, R = In_Ptr->beam.R;
  double      b_max = (In_Ptr->nr - 1) * In_Ptr->dr;

  puts("The convolution may take a little while. Wait...");
  for (irc = 0; irc < In_Ptr->nrc; irc++) {
    rc = (irc + 0.5) * In_Ptr->drc;
    ConvVar.r = rc;
    ConvVar.tree = NULL;	/* init the tree. */
    for (ia = 0; ia < In_Ptr->na; ia++) {
      ConvVar.ia = ia;
      if (In_Ptr->beam.type == flat)
	Out_Ptr->Tt_rac[irc][ia] = 2 * P / (R * R)
	  * FlatIntegration(Tt_raFGIntegrand);
      else			/* Gaussian. */
	Out_Ptr->Tt_rac[irc][ia] = 4 * P / (R * R)
	  * GaussIntegration(Tt_raFGIntegrand);
    }
    FreeTree(ConvVar.tree);
  }

  Out_Ptr->conved.Tt_ra = 1;
}

/****************************************************************
 ****/
void
ConvTt_r(InputStruct * In_Ptr, OutStruct * Out_Ptr)
{
  short       irc;
  double      rc, P = In_Ptr->beam.P, R = In_Ptr->beam.R;
  double      b_max = (In_Ptr->nr - 1) * In_Ptr->dr;

  for (irc = 0; irc < In_Ptr->nrc; irc++) {
    rc = (irc + 0.5) * In_Ptr->drc;
    ConvVar.r = rc;
    if (In_Ptr->beam.type == flat)
      Out_Ptr->Tt_rc[irc] = 2 * P / (R * R)
	* FlatIntegration(Tt_rFGIntegrand);
    else			/* Gaussian. */
      Out_Ptr->Tt_rc[irc] = 4 * P / (R * R)
	* GaussIntegration(Tt_rFGIntegrand);
  }

  Out_Ptr->conved.Tt_r = 1;
}

/****************************************************************
 ****/
void
ShowOutConvMenu(char *in_fname)
{
  printf("Arz = absorption vs r & z [J/cm3]\n");
  printf("Frz = fluence vs r & z [J/cm2]\n");

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

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

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


/****************************************************************
 *	3 numbers each line: r, z, A[r][z].
 ****/
void
WriteA_rzc(InputStruct * In_Ptr,
	   double **A_rzc)
{
  FILE       *file;
  short       ir, iz, nr = In_Ptr->nr, nz = In_Ptr->nz;
  double      r, z, dr = In_Ptr->dr, dz = In_Ptr->dz;
  char        fname[STRLEN];

#if IBMPC
  strcpy(fname, "Arz");
#else
  strcpy(fname, "Arzc");
#endif
  if ((file = GetWriteFile(fname)) == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[J/cm3]\n",
	  "r[cm]", "z[cm]", fname);
  for (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    for (iz = 0; iz < nz; iz++) {
      z = (iz + 0.5) * dz;
      fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n",
	      r, z, A_rzc[ir][iz]);
    }
  }

  fclose(file);
}

/****************************************************************
 *	3 numbers each line: r, z, F[r][z].
 ****/
void
WriteF_rzc(InputStruct * In_Ptr,
	   double **A_rzc)
{
  FILE       *file;
  short       ir, iz, nr = In_Ptr->nr, nz = In_Ptr->nz;
  double      mua, r, z, dr = In_Ptr->dr, dz = In_Ptr->dz;
  char        fname[STRLEN];

#if IBMPC
  strcpy(fname, "Frz");
#else
  strcpy(fname, "Frzc");
#endif
  if ((file = GetWriteFile(fname)) == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[J/cm3]\n",
	  "r[cm]", "z[cm]", fname);
  for (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    for (iz = 0; iz < nz; iz++) {
      z = (iz + 0.5) * dz;
      mua = In_Ptr->layerspecs[IzToLayer(iz, In_Ptr)].mua;
      if (mua > 0.0)
	fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n",
		r, z, A_rzc[ir][iz] / mua);
      else
	fprintf(file, "%-12.4E\t%-12.4E\t%-12.4E\n", r, z, 0.0);
    }
  }

  fclose(file);
}

/****************************************************************
 *	3 numbers each line: r, a, Rd[r][a].
 ****/
void
WriteRd_rac(InputStruct * In_Ptr,
	    double **Rd_rac)
{
  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];

#if IBMPC
  strcpy(fname, "Rra");
#else
  strcpy(fname, "Rrac");
#endif
  if ((file = GetWriteFile(fname)) == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[J/(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, Rd_rac[ir][ia]);
    }
  }

  fclose(file);
}

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

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

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

  fclose(file);
}

/****************************************************************
 *	3 numbers each line:r, a, Tt[r][a]. a = theta.
 ****/
void
WriteTt_rac(InputStruct * In_Ptr,
	    double **Tt_rac)
{
  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];

#if IBMPC
  strcpy(fname, "Tra");
#else
  strcpy(fname, "Trac");
#endif
  if ((file = GetWriteFile(fname)) == NULL)
    return;

  fprintf(file, "%-12s\t%-12s\t%-s[J/(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_rac[ir][ia]);
    }
  }

  fclose(file);
}

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

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

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

  fclose(file);
}

/****************************************************************
 ****/
void
BranchOutConvA(char *Cmd_Str,
	       InputStruct * In_Ptr,
	       OutStruct * Out_Ptr)
{
  switch (toupper(Cmd_Str[1])) {
    case 'R':
    if (toupper(Cmd_Str[2]) == 'Z') {	/* A_rzc. */
      if (!Out_Ptr->conved.A_rz)
	ConvA_rz(In_Ptr, Out_Ptr);
      WriteA_rzc(In_Ptr, Out_Ptr->A_rzc);
    } else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void
BranchOutConvF(char *Cmd_Str,
	       InputStruct * In_Ptr,
	       OutStruct * Out_Ptr)
{
  switch (toupper(Cmd_Str[1])) {
    case 'R':
    if (toupper(Cmd_Str[2]) == 'Z') {	/* F_rzc. */
      if (!Out_Ptr->conved.A_rz)
	ConvA_rz(In_Ptr, Out_Ptr);
      WriteF_rzc(In_Ptr, Out_Ptr->A_rzc);
    } else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

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

  switch (toupper(Cmd_Str[1])) {
  case 'R':
    ch = toupper(Cmd_Str[2]);
    if (ch == '\0') {		/* Rd_rc. */
      if (!Out_Ptr->conved.Rd_r)
	ConvRd_r(In_Ptr, Out_Ptr);
      WriteRd_rc(In_Ptr, Out_Ptr->Rd_rc);
    } else if (ch == 'A') {	/* Rd_rac. */
      if (!Out_Ptr->conved.Rd_ra)
	ConvRd_ra(In_Ptr, Out_Ptr);
      WriteRd_rac(In_Ptr, Out_Ptr->Rd_rac);
    } else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

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

  switch (toupper(Cmd_Str[1])) {
  case 'R':
    ch = toupper(Cmd_Str[2]);
    if (ch == '\0') {		/* Tt_rc. */
      if (!Out_Ptr->conved.Tt_r)
	ConvTt_r(In_Ptr, Out_Ptr);
      WriteTt_rc(In_Ptr, Out_Ptr->Tt_rc);
    } else if (ch == 'A') {	/* Tt_rac. */
      if (!Out_Ptr->conved.Tt_ra)
	ConvTt_ra(In_Ptr, Out_Ptr);
      WriteTt_rac(In_Ptr, Out_Ptr->Tt_rac);
    } else
      puts("...Wrong command");
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 *	The global ConvVar is used.
 ****/
void
BranchOutConvCmd(char *Cmd_Str,
		 InputStruct * In_Ptr,
		 OutStruct * Out_Ptr)
{
  char        ch;

  ConvVar.in_ptr = In_Ptr;
  ConvVar.out_ptr = Out_Ptr;

  switch (toupper(Cmd_Str[0])) {
  case 'A':
    BranchOutConvA(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'F':
    BranchOutConvF(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'R':
    BranchOutConvR(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'T':
    BranchOutConvT(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'H':
    ShowOutConvMenu(In_Ptr->in_fname);
    break;
  case 'Q':
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void