convo.c

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

  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 
A2F(InputStruct * In_Ptr, double **A_rz)
{
  short       nz = In_Ptr->nz, nr = In_Ptr->nr;
  double      dr = In_Ptr->dr, dz = In_Ptr->dz;
  short       ir, iz;
  double      mua;

  for (ir = 0; ir < nr; 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 
F2A(InputStruct * In_Ptr, double **A_rz)
{
  short       nz = In_Ptr->nz, nr = In_Ptr->nr;
  double      dr = In_Ptr->dr, dz = In_Ptr->dz;
  short       ir, iz;
  double      mua;

  for (ir = 0; ir < nr; 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 
BranchContOrigCmd(char *Cmd_Str,
		  InputStruct * In_Ptr,
		  OutStruct * Out_Ptr)
{
  char        ch;

  switch (toupper(Cmd_Str[0])) {
  case 'A':
    IsoPlot(Out_Ptr->A_rz, In_Ptr->nr - 1, In_Ptr->nz - 1,
	    In_Ptr->dr, In_Ptr->dz);
    break;
  case 'F':
    A2F(In_Ptr, Out_Ptr->A_rz);
    IsoPlot(Out_Ptr->A_rz, In_Ptr->nr - 1, In_Ptr->nz - 1,
	    In_Ptr->dr, In_Ptr->dz);
    F2A(In_Ptr, Out_Ptr->A_rz);
    break;
  case 'R':
    IsoPlot(Out_Ptr->Rd_ra, In_Ptr->nr - 1, In_Ptr->na - 1,
	    In_Ptr->dr, In_Ptr->da);
    break;
  case 'T':
    IsoPlot(Out_Ptr->Tt_ra, In_Ptr->nr - 1, In_Ptr->na - 1,
	    In_Ptr->dr, In_Ptr->da);
    break;
  case 'H':
    ShowContOrigMenu(In_Ptr->in_fname);
    break;
  case 'Q':
    break;
  default:
    puts("...Wrong command");
  }
}
/****************************************************************
 ****/
void 
ContourOrigData(InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  char        cmd_str[STRLEN];

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

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

/****************************************************************
 *	Ext is either "Ars" or "Frs".
 ****/
void 
ScanOrigA_r(char *Ext, InputStruct * In_Ptr, double **A_rz)
{
  short       ir, iz, nr = In_Ptr->nr, nz = In_Ptr->nz;
  double      r, z, dr = In_Ptr->dr, 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 (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    fprintf(file, "%-12.4E\t%-12.4E\n", r, A_rz[ir][iz]);
  }

  fclose(file);
}

/****************************************************************
 *	Ext is either "Azs" or "Fzs".
 ****/
void 
ScanOrigA_z(char *Ext, InputStruct * In_Ptr, double **A_rz)
{
  short       ir, iz, nr = In_Ptr->nr, nz = In_Ptr->nz;
  double      r, z, dr = In_Ptr->dr, dz = In_Ptr->dz;
  FILE       *file;

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

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

  fclose(file);
}

/****************************************************************
 ****/
void 
ScanOrigRd_r(InputStruct * In_Ptr, double **Rd_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, "Rrs");
  file = GetWriteFile(fname);
  if (file == 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 (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    fprintf(file, "%-12.4E\t%-12.4E\n", r, Rd_ra[ir][ia]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void 
ScanOrigRd_a(InputStruct * In_Ptr, double **Rd_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, "Ras");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  printf("r grid separation is %-10.4lg cm.\n", dr);
  printf("Input fixed r index (0 - %2hd): ", nr - 1);
  ir = GetShort(0, nr - 1);
  fprintf(file, "%-12s\t%-s@r=%-9.3lg\n", "a[rad]", fname, dr * (ir + 0.5));
  for (ia = 0; ia < na; ia++) {
    a = (ia + 0.5) * da;
    fprintf(file, "%-12.4E\t%-12.4E\n", a, Rd_ra[ir][ia]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void 
ScanOrigTt_r(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, "Trs");
  file = GetWriteFile(fname);
  if (file == 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 (ir = 0; ir < nr; ir++) {
    r = (ir + 0.5) * dr;
    fprintf(file, "%-12.4E\t%-12.4E\n", r, Tt_ra[ir][ia]);
  }

  fclose(file);
}

/****************************************************************
 ****/
void 
ScanOrigTt_a(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, "Tas");
  file = GetWriteFile(fname);
  if (file == NULL)
    return;

  printf("r grid separation is %-10.4lg cm.\n", dr);
  printf("Input fixed r index (0 - %2hd): ", nr - 1);
  ir = GetShort(0, nr - 1);
  fprintf(file, "%-12s\t%-s@r=%-9.3lg\n", "a[rad]", fname, dr * (ir + 0.5));
  for (ia = 0; ia < na; ia++) {
    a = (ia + 0.5) * da;
    fprintf(file, "%-12.4E\t%-12.4E\n", a, Tt_ra[ir][ia]);
  }

  fclose(file);
}

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

  switch (toupper(Cmd_Str[1])) {
  case 'R':
    strcpy(fname, "Ars");
    ScanOrigA_r(fname, In_Ptr, Out_Ptr->A_rz);
    break;
  case 'Z':
    strcpy(fname, "Azs");
    ScanOrigA_z(fname, In_Ptr, Out_Ptr->A_rz);
    break;
  default:
    puts("...Wrong command");
  }
}

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

  A2F(In_Ptr, Out_Ptr->A_rz);

  switch (toupper(Cmd_Str[1])) {
  case 'R':
    strcpy(fname, "Frs");
    ScanOrigA_r(fname, In_Ptr, Out_Ptr->A_rz);
    break;
  case 'Z':
    strcpy(fname, "Fzs");
    ScanOrigA_z(fname, In_Ptr, Out_Ptr->A_rz);
    break;
  default:
    puts("...Wrong command");
  }

  F2A(In_Ptr, Out_Ptr->A_rz);
}

/****************************************************************
 ****/
void 
BranchScanOrigR(char *Cmd_Str,
		InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  switch (toupper(Cmd_Str[1])) {
    case 'R':
    ScanOrigRd_r(In_Ptr, Out_Ptr->Rd_ra);
    break;
  case 'A':
    ScanOrigRd_a(In_Ptr, Out_Ptr->Rd_ra);
    break;
  default:
    puts("...Wrong command");
  }
}

/****************************************************************
 ****/
void 
BranchScanOrigT(char *Cmd_Str,
		InputStruct * In_Ptr,
		OutStruct * Out_Ptr)
{
  switch (toupper(Cmd_Str[1])) {
    case 'R':
    ScanOrigTt_r(In_Ptr, Out_Ptr->Tt_ra);
    break;
  case 'A':
    ScanOrigTt_a(In_Ptr, Out_Ptr->Tt_ra);
    break;
  default:
    puts("...Wrong command");
  }
}

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

  switch (toupper(Cmd_Str[0])) {
  case 'A':
    BranchScanOrigA(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'F':
    BranchScanOrigF(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'R':
    BranchScanOrigR(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'T':
    BranchScanOrigT(Cmd_Str, In_Ptr, Out_Ptr);
    break;
  case 'H':
    ShowScanOrigMenu(In_Ptr->in_fname);
    break;
  case 'Q':
    break;
  default:
    puts("...Wrong command");
  }
}

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

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