convconv.c

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

OutputConvData(InputStruct * In_Ptr,
	       OutStruct * Out_Ptr)
{
  char        cmd_str[STRLEN];

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

/****************************Contours***************************/
/****************************************************************
 *	Absorption density to fluence. A = F/mua;
 ****/
void
A2Fconv(InputStruct * In_Ptr, double **A_rz)
{
  short       nz = In_Ptr->nz, nrc = In_Ptr->nrc;
  short       ir, iz;
  double      mua;

  for (ir = 0; ir < nrc; ir++)
    for (iz = 0; iz < nz; iz++) {
      mua = In_Ptr->layerspecs[IzToLayer(iz, In_Ptr)].mua;
      if (mua > 0.0)
	A_rz[ir][iz] /= mua;
    }
}

/****************************************************************
 *	Fluence to absorption density. F = A*mua;
 ****/
void
F2Aconv(InputStruct * In_Ptr, double **A_rz)
{
  short       nz = In_Ptr->nz, nrc = In_Ptr->nrc;
  short       ir, iz;
  double      mua;

  for (ir = 0; ir < nrc; ir++)
    for (iz = 0; iz < nz; iz++) {
      mua = In_Ptr->layerspecs[IzToLayer(iz, In_Ptr)].mua;
      if (mua > 0.0)
	A_rz[ir][iz] *= mua;
    }
}

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

/****************************************************************
 ****/
void
BranchContConvCmd(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':
    if (!Out_Ptr->conved.A_rz)
      ConvA_rz(In_Ptr, Out_Ptr);
    IsoPlot(Out_Ptr->A_rzc, In_Ptr->nrc - 1, In_Ptr->nz - 1,
	    In_Ptr->drc, In_Ptr->dz);
    break;
  case 'F':
    if (!Out_Ptr->conved.A_rz)
      ConvA_rz(In_Ptr, Out_Ptr);
    A2Fconv(In_Ptr, Out_Ptr->A_rzc);
    IsoPlot(Out_Ptr->A_rzc, In_Ptr->nrc - 1, In_Ptr->nz - 1,
	    In_Ptr->drc, In_Ptr->dz);
    F2Aconv(In_Ptr, Out_Ptr->A_rzc);
    break;
  case 'R':
    if (!Out_Ptr->conved.Rd_ra)
      ConvRd_ra(In_Ptr, Out_Ptr);
    IsoPlot(Out_Ptr->Rd_rac, In_Ptr->nrc - 1, In_Ptr->na - 1,
	    In_Ptr->drc, In_Ptr->da);
    break;
  case 'T':
    if (!Out_Ptr->conved.Tt_ra)
      ConvTt_ra(In_Ptr, Out_Ptr);
    IsoPlot(Out_Ptr->Tt_rac, In_Ptr->nrc - 1, In_Ptr->na - 1,
	    In_Ptr->drc, In_Ptr->da);
    break;
  case 'H':
    ShowContConvMenu(In_Ptr->in_fname);
    break;
  case 'Q':
    break;
  default:
    puts("...Wrong command");
  }
}
/****************************************************************
 ****/
void
ContourConvData(InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  char        cmd_str[STRLEN];

  if (!Out_Ptr->allocated)
    puts("...No data to output");
  else if (In_Ptr->beam.type == pencil)
    puts("...No incident beam specified");
  else
    do {
      printf("\n> Contour output of convolved data (h for help) => ");
      do
	gets(cmd_str);
      while (!strlen(cmd_str));	/* avoid null string. */
      BranchContConvCmd(cmd_str, In_Ptr, Out_Ptr);
    } while (toupper(cmd_str[0]) != 'Q');
}


/****************************Scanning***************************/
/****************************************************************
 ****/
void
ShowScanConvMenu(char *in_fname)
{
  printf("Ar = absorption vs r @ fixed z [J/cm3]\n");
  printf("Az = absorption vs z @ fixed r [J/cm3]\n");
  printf("Fr = fluence vs r @ fixed z [J/cm2]\n");
  printf("Fz = fluence vs z @ fixed r [J/cm2]\n");
  printf("Rr = diffuse reflectance vs r @ fixed angle [J/(cm2 sr)]\n");
  printf("Ra = diffuse reflectance vs angle @ fixed r [J/(cm2 sr)]\n");
  printf("Tr = transmittance vs r @ fixed angle [J/(cm2 sr)]\n");
  printf("Ta = transmittance vs angle @ fixed r [J/(cm2 sr)]\n");
  printf("Q  = quit\n");
  printf("* input filename: %s \n", in_fname);
}

/****************************************************************
 *	Ext is either "Ars" or "Frs".
 ****/
void
ScanConvA_r(char *Ext, InputStruct * In_Ptr, double **A_rzc)
{
  short       irc, iz, nrc = In_Ptr->nrc, nz = In_Ptr->nz;
  double      r, z, drc = In_Ptr->drc, dz = In_Ptr->dz;
  FILE       *file;

  file = GetWriteFile(Ext);
  if (file == NULL)
    return;

  printf("z grid separation is %-10.4lg cm.\n", dz);
  printf("Input fixed z index (0 - %2hd): ", nz - 1);
  iz = GetShort(0, nz - 1);
  fprintf(file, "%-12s\t%-s@z=%-9.3lg\n", "r[cm]", Ext, dz * (iz + 0.5));
  for (irc = 0; irc < nrc; irc++) {
    r = (irc + 0.5) * drc;
    fprintf(file, "%-12.4E\t%-12.4E\n", r, A_rzc[irc][iz]);
  }

  fclose(file);
}

/****************************************************************
 *	Ext is either "Azs" or "Fzs".
 ****/
void
ScanConvA_z(char *Ext, InputStruct * In_Ptr, double **A_rzc)
{
  short       irc, iz, nrc = In_Ptr->nrc, nz = In_Ptr->nz;
  double      r, z, drc = In_Ptr->drc, dz = In_Ptr->dz;
  FILE       *file;

  file = GetWriteFile(Ext);
  if (file == NULL)
    return;

  printf("r grid separation is %-10.4lg cm.\n", drc);
  printf("Input fixed r index (0 - %2hd): ", nrc - 1);
  irc = GetShort(0, nrc - 1);
  fprintf(file, "%-12s\t%-s@r=%-9.3lg\n", "z[cm]", Ext, drc * (irc + 0.5));
  for (iz = 0; iz < nz; iz++) {
    z = (iz + 0.5) * dz;
    fprintf(file, "%-12.4E\t%-12.4E\n", z, A_rzc[irc][iz]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void
ScanConvRd_r(InputStruct * In_Ptr, double **Rd_rac)
{
  short       irc, ia, nrc = In_Ptr->nrc, na = In_Ptr->na;
  double      r, a, drc = In_Ptr->drc, da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

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

  printf("Angle grid separation is %-10.4lg rad.\n", da);
  printf("Input fixed angle index (0 - %2hd): ", na - 1);
  ia = GetShort(0, na - 1);

  fprintf(file, "%-12s\t%-s@a=%-9.3lg\n", "r[cm]", fname, da * (ia + 0.5));
  for (irc = 0; irc < nrc; irc++) {
    r = (irc + 0.5) * drc;
    fprintf(file, "%-12.4E\t%-12.4E\n", r, Rd_rac[irc][ia]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void
ScanConvRd_a(InputStruct * In_Ptr, double **Rd_rac)
{
  short       irc, ia, nrc = In_Ptr->nrc, na = In_Ptr->na;
  double      r, a, drc = In_Ptr->drc, da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

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

  printf("r grid separation is %-10.4lg cm.\n", drc);
  printf("Input fixed r index (0 - %2hd): ", nrc - 1);
  irc = GetShort(0, nrc - 1);

  fprintf(file, "%-12s\t%-s@r=%-9.3lg\n", "a[rad]", fname, drc * (irc + 0.5));
  for (ia = 0; ia < na; ia++) {
    a = (ia + 0.5) * da;
    fprintf(file, "%-12.4E\t%-12.4E\n", a, Rd_rac[irc][ia]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void
ScanConvTt_r(InputStruct * In_Ptr, double **Tt_rac)
{
  short       irc, ia, nrc = In_Ptr->nrc, na = In_Ptr->na;
  double      r, a, drc = In_Ptr->drc, da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

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

  printf("Angle grid separation is %-10.4lg rad.\n", da);
  printf("Input fixed angle index (0 - %2hd): ", na - 1);
  ia = GetShort(0, na - 1);

  fprintf(file, "%-12s\t%-s@a=%-9.3lg\n", "r[cm]", fname, da * (ia + 0.5));
  for (irc = 0; irc < nrc; irc++) {
    r = (irc + 0.5) * drc;
    fprintf(file, "%-12.4E\t%-12.4E\n", r, Tt_rac[irc][ia]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void
ScanConvTt_a(InputStruct * In_Ptr, double **Tt_rac)
{
  short       irc, ia, nrc = In_Ptr->nrc, na = In_Ptr->na;
  double      r, a, drc = In_Ptr->drc, da = In_Ptr->da;
  FILE       *file;
  char        fname[STRLEN];

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

  printf("r grid separation is %-10.4lg cm.\n", drc);
  printf("Input fixed r index (0 - %2hd): ", nrc - 1);
  irc = GetShort(0, nrc - 1);

  fprintf(file, "%-12s\t%-s@r=%-9.3lg\n", "a[rad]", fname, drc * (irc + 0.5));
  for (ia = 0; ia < na; ia++) {
    a = (ia + 0.5) * da;
    fprintf(file, "%-12.4E\t%-12.4E\n", a, Tt_rac[irc][ia]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void
BranchScanConvA(char *Cmd_Str,
		InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  char        fname[STRLEN];

  if (!Out_Ptr->conved.A_rz)
    ConvA_rz(In_Ptr, Out_Ptr);

  switch (toupper(Cmd_Str[1])) {
  case 'R':
#if IBMPC
    strcpy(fname, "Ars");
    ScanConvA_r(fname, In_Ptr, Out_Ptr->A_rzc);
#else
    strcpy(fname, "Arsc");
    ScanConvA_r(fname, In_Ptr, Out_Ptr->A_rzc);
#endif
    break;
  case 'Z':
#if IBMPC
    strcpy(fname, "Azs");
    ScanConvA_z(fname, In_Ptr, Out_Ptr->A_rzc);
#else
    strcpy(fname, "Azsc");
    ScanConvA_z(fname, In_Ptr, Out_Ptr->A_rzc);
#endif
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void
BranchScanConvF(char *Cmd_Str,
		InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  char        fname[STRLEN];

  if (!Out_Ptr->conved.A_rz)
    ConvA_rz(In_Ptr, Out_Ptr);
  A2Fconv(In_Ptr, Out_Ptr->A_rzc);

  switch (toupper(Cmd_Str[1])) {
  case 'R':
#if IBMPC
    strcpy(fname, "Frs");
    ScanConvA_r(fname, In_Ptr, Out_Ptr->A_rzc);
#else
    strcpy(fname, "Frsc");
    ScanConvA_r(fname, In_Ptr, Out_Ptr->A_rzc);
#endif
    break;
  case 'Z':
#if IBMPC
    strcpy(fname, "Fzs");
    ScanConvA_z(fname, In_Ptr, Out_Ptr->A_rzc);
#else
    strcpy(fname, "Fzsc");
    ScanConvA_z(fname, In_Ptr, Out_Ptr->A_rzc);
#endif
    break;
  default:
    puts("...Wrong command");
  }

  F2Aconv(In_Ptr, Out_Ptr->A_rzc);
}

/****************************************************************
 ****/
void
BranchScanConvR(char *Cmd_Str,
		InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  if (!Out_Ptr->conved.Rd_ra)
    ConvRd_ra(In_Ptr, Out_Ptr);
  switch (toupper(Cmd_Str[1])) {
  case 'R':
    ScanConvRd_r(In_Ptr, Out_Ptr->Rd_rac);
    break;
  case 'A':
    ScanConvRd_a(In_Ptr, Out_Ptr->Rd_rac);
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void
BranchScanConvT(char *Cmd_Str,
		InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  if (!Out_Ptr->conved.Tt_ra)
    ConvTt_ra(In_Ptr, Out_Ptr);
  switch (toupper(Cmd_Str[1])) {
  case 'R':
    ScanConvTt_r(In_Ptr, Out_Ptr->Tt_rac);
    break;
  case 'A':
    ScanConvTt_a(In_Ptr, Out_Ptr->Tt_rac);
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void
BranchScanConvCmd(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':
    BranchScanConvA(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'F':
    BranchScanConvF(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'R':
    BranchScanConvR(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'T':
    BranchScanConvT(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'H':
    ShowScanConvMenu(In_Ptr->in_fname);
    break;
  case 'Q':
    break;
  default:
    puts("...Wrong command");
  }
}

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

  if (!Out_Ptr->allocated)
    puts("...No data to output");
  else if (In_Ptr->beam.type == pencil)
    puts("...No incident beam specified");
  else
    do {
      printf("\n> Scans of convolved data (h for help) => ");
      do
	gets(cmd_str);
      while (!strlen(cmd_str));	/* avoid null string. */
      BranchScanConvCmd(cmd_str, In_Ptr, Out_Ptr);
    } while (toupper(cmd_str[0]) != 'Q');
}