154 long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel;
155 double topoff, TotMinusExplicitResolved,
156 TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.;
157 double RecExplictLevels, TotalRadRecomb, RecCollapsed;
159 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
160 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
161 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
162 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
163 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
164 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
174 RecExplictLevels = 0.;
177 ipFirstCollapsed =
iso_sp[ipISO][nelem].numLevels_local-
iso_sp[ipISO][nelem].nCollapsed_local;
179 ASSERT( ipFirstCollapsed ==
iso_sp[ipISO][nelem].numLevels_local -
iso_sp[ipISO][nelem].nCollapsed_local );
182 TeUsed[ipISO][nelem] =
phycon.te;
184 for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel )
189 if( !
iso_ctrl.lgNoRecombInterp[ipISO] )
203 RecExplictLevels +=
iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[
ipRecRad];
209 Total_DR_Added +=
iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb;
217 for( ipLevel=ipFirstCollapsed; ipLevel<
iso_sp[ipISO][nelem].numLevels_local; ++ipLevel )
225 RecCollapsed +=
iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[
ipRecRad];
233 Total_DR_Added +=
iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb;
236 for( ipLevel=
iso_sp[ipISO][nelem].numLevels_local; ipLevel<
iso_sp[ipISO][nelem].numLevels_max;++ipLevel )
238 iso_sp[ipISO][nelem].fb[ipLevel].DielecRecomb = 0.;
241 for( ipLevel = 0; ipLevel<
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
243 if(
N_(ipLevel) == ThisN )
245 TotRRThisN +=
iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[
ipRecRad];
249 ASSERT(
iso_ctrl.lgRandErrGen[ipISO] || (TotRRThisN<TotRRLastN) || (ThisN<=2L) || (
phycon.te>3E7) || (nelem!=
ipHELIUM) || (ThisN==
iso_sp[ipISO][nelem].n_HighestResolved_local+1) );
252 TotRRLastN = TotRRThisN;
253 TotRRThisN =
iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[
ipRecRad];
258 enum {DEBUG_LOC=
false};
260 static long RUNONCE =
false;
262 if( !RUNONCE && DEBUG_LOC )
264 static long FIRSTTIME =
true;
268 fprintf(
ioQQQ,
"Sum of Radiative Recombination at current iso, nelem, temp = %li %li %.2f\n",
273 fprintf(
ioQQQ,
"%li\t%.2e\n",LastN,TotRRLastN );
281 if(
iso_ctrl.lgNoRecombInterp[ipISO] )
284 TotalRadRecomb = RecCollapsed + RecExplictLevels;
297 else if(
iso_sp[ipISO][nelem].lgLevelsLowered )
300 for( ipLevel = 0; ipLevel<
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
301 TotalRadRecomb +=
iso_sp[ipISO][nelem].fb[ipLevel].RadRecomb[
ipRecRad];
306 TotalRadRecomb =
iso_RRCoef_Te( ipISO, nelem,
iso_sp[ipISO][nelem].numLevels_max -
iso_sp[ipISO][nelem].nCollapsed_max );
308 if(ipISO==0 && nelem==0 )
311 TeUsed[ipISO][nelem] =
phycon.te*0.999;
333 iso_sp[ipISO][nelem].RadRec_caseB = TotalRadRecomb -
iso_sp[ipISO][nelem].fb[0].RadRecomb[
ipRecRad];
337 for(
unsigned i = 0; i <
iso_sp[ipISO][nelem].HighestLevelOpacStack.size(); ++i )
339 long index =
iso_sp[ipISO][nelem].numLevels_max-1;
340 opac.OpacStack[
iso_sp[ipISO][nelem].fb[index].ipOpac-1+i] =
iso_sp[ipISO][nelem].HighestLevelOpacStack[i];
346 if( !
iso_sp[ipISO][nelem].lgLevelsLowered )
351 TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels;
353 topoff = TotMinusExplicitResolved - RecCollapsed;
359 fabs( topoff/TotalRadRecomb ) > 0.01 )
361 fprintf(
ioQQQ,
" PROBLEM negative RR topoff for %li, rel error was %.2e, temp was %.2f\n",
362 nelem, topoff/TotalRadRecomb,
phycon.te );
369 topoff =
MAX2( 0., topoff );
370 double scale_factor = 1. + topoff/
iso_sp[ipISO][nelem].fb[
iso_sp[ipISO][nelem].numLevels_max-1].RadRecomb[
ipRecRad];
371 ASSERT( scale_factor >= 1. );
374 for(
unsigned i = 0; i <
iso_sp[ipISO][nelem].HighestLevelOpacStack.size(); ++i )
376 long index =
iso_sp[ipISO][nelem].numLevels_max-1;
377 opac.OpacStack[
iso_sp[ipISO][nelem].fb[index].ipOpac-1+i] *= scale_factor;
384 if( Total_DR_Added > TotalRadRecomb/100. )
386 if( Total_DR_Added /
ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] > 1.02 )
388 fprintf(
ioQQQ,
" PROBLEM negative DR topoff for %li, tot1, tot2 = %.2e\t%.2e rel error was %.2e, temp was %.2f\n",
391 ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO],
392 Total_DR_Added /
ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - 1.0,
403 iso_sp[ipISO][nelem].fb[
iso_sp[ipISO][nelem].numLevels_max-1].DielecRecomb +=
404 MAX2( 0.,
ionbal.DR_Badnell_rate_coef[nelem][nelem-ipISO] - Total_DR_Added );
415 iso_sp[ipISO][nelem].RadRec_effec = 0.;
417 for( ipLevel=0; ipLevel<
iso_sp[ipISO][nelem].numLevels_local; ++ipLevel )
420 if(
opac.lgCaseB && ipLevel==0 )
440 iso_sp[ipISO][nelem].fb[ipLevel].ConOpacRatio );
453 for( ipLevel=
iso_sp[ipISO][nelem].numLevels_local; ipLevel<
iso_sp[ipISO][nelem].numLevels_max; ++ipLevel )
462 if(
trace.lgTrace &&
trace.lgIsoTraceFull[ipISO] && (nelem ==
trace.ipIsoTrace[ipISO]) )
464 fprintf(
ioQQQ,
" iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem );
467 fprintf(
ioQQQ,
" iso_radiative_recomb recomb effic" );
468 for( ipLevel=0; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
472 fprintf(
ioQQQ,
"\n" );
475 fprintf(
ioQQQ,
" iso_radiative_recomb recomb net effic" );
476 for( ipLevel=0; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
480 fprintf(
ioQQQ,
"\n" );
483 fprintf(
ioQQQ,
" iso_radiative_recomb in optic dep" );
484 for( ipLevel=0; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
488 fprintf(
ioQQQ,
"\n" );
491 fprintf(
ioQQQ,
" iso_radiative_recomb out op depth" );
492 for( ipLevel=0; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
496 fprintf(
ioQQQ,
"\n" );
499 fprintf(
ioQQQ,
" iso_radiative_recomb rad rec coef " );
500 for( ipLevel=0; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
504 fprintf(
ioQQQ,
"\n" );
509 fprintf(
ioQQQ,
" iso_radiative_recomb total rec coef" );
511 fprintf(
ioQQQ,
" case A=" );
514 fprintf(
ioQQQ,
" case B=");
516 fprintf(
ioQQQ,
"\n" );
524 for( ipLevel=0; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
529 " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",
551 fprintf(
save.ioRecom,
"%s %s %2li ",
554 fprintf(
save.ioRecom,
"\n" );
567 0.0009f,0.0037f,0.0003f,0.0003f,0.0003f,0.0018f,
568 0.0009f,0.0050f,0.0007f,0.0003f,0.0001f,0.0007f,
569 0.0045f,0.0071f,0.0005f,0.0005f,0.0004f,0.0005f,0.0004f,0.0009f,
570 0.0045f,0.0071f,0.0005f,0.0005f,0.0004f,0.0005f,0.0004f,0.0005f,0.0004f,0.0009f };
574 0.0100f,0.0060f,0.0080f,0.0080f,0.0080f,0.0200f,
575 0.0200f,0.0200f,0.0200f,0.0600f,0.0600f,0.0080f,
576 0.0200f,0.0200f,0.0070f,0.0100f,0.0100f,0.0020f,0.0030f,0.0070f,
577 0.0300f,0.0300f,0.0100f,0.0200f,0.0200f,0.0200f,0.0200f,0.0010f,0.0004f,0.0090f };
580 for(
long ipHi=0; ipHi<
iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
584 if( ipHi >=
iso_sp[ipISO][nelem].numLevels_max -
iso_sp[ipISO][nelem].nCollapsed_max )
589 level =
iso_sp[ipISO][nelem].QuantumNumbers2Index[5][
MIN2(
L_(ipHi),4) ][
S_(ipHi)];
607 for(
long ipHi=0; ipHi <
iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
609 iso_sp[ipISO][nelem].fb[ipHi].RadEffec = 0.;
612 for(
long ipHigher=ipHi; ipHigher <
iso_sp[ipISO][nelem].numLevels_local; ipHigher++ )
614 ASSERT(
iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] >= 0. );
617 iso_sp[ipISO][nelem].fb[ipHi].RadEffec +=
iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] *
626 enum {DEBUG_LOC=
false};
632 fprintf(
ioQQQ,
"Effective recombination, ipISO=%li, nelem=%li, Te = %e\n", ipISO, nelem,
phycon.te );
633 fprintf(
ioQQQ,
"N\tL\tS\tRadEffec\tLifetime\n" );
635 for(
long i=0; i<maxPrt; i++ )
637 fprintf(
ioQQQ,
"%li\t%li\t%li\t%e\t%e\n",
N_(i),
L_(i),
S_(i),
638 iso_sp[ipISO][nelem].fb[i].RadEffec,
639 MAX2( 0.,
iso_sp[ipISO][nelem].st[i].lifetime() ) );
641 fprintf(
ioQQQ,
"\n" );
648 dprintf(
ioQQQ,
"ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" );
651 for(
long ipHi=0; ipHi <
iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
653 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec = 0.;
656 for(
long ipHigher=ipHi; ipHigher <
iso_sp[ipISO][nelem].numLevels_local; ipHigher++ )
659 ASSERT(
iso_sp[ipISO][nelem].
ex[ipHigher][ipHi].SigmaCascadeProb >= 0. );
662 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec += pow(
iso_sp[ipISO][nelem].
ex[
iso_sp[ipISO][nelem].numLevels_max][ipHigher].Error[
IPRAD] *
663 iso_sp[ipISO][nelem].CascadeProb[ipHigher][ipHi] *
iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[
ipRecRad], 2.) +
664 pow(
iso_sp[ipISO][nelem].
ex[ipHigher][ipHi].SigmaCascadeProb *
iso_sp[ipISO][nelem].fb[ipHigher].RadRecomb[
ipRecRad], 2.);
667 ASSERT(
iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec >= 0. );
668 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec = sqrt(
iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec );
670 for(
long ipLo = 0; ipLo < ipHi; ipLo++ )
672 if( ((
L_(ipLo) ==
L_(ipHi) + 1 ) || (
L_(ipHi) ==
L_(ipLo) + 1 )) )
674 double EnergyInRydbergs =
iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd -
iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd;
676 double emissivity =
iso_sp[ipISO][nelem].fb[ipHi].RadEffec *
iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] *
EN1RYD * EnergyInRydbergs;
677 double sigma_emiss = 0., SigmaBranchRatio = 0.;
679 if( ( emissivity > 2.E-29 ) && (
wavelength < 1.E6 ) && (
N_(ipHi)<=5) )
681 SigmaBranchRatio =
iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] * sqrt(
682 pow( (
double)
iso_sp[ipISO][nelem].
ex[ipHi][ipLo].Error[
IPRAD], 2. ) +
683 pow(
iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*
iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) );
685 sigma_emiss =
EN1RYD * EnergyInRydbergs * sqrt(
686 pow( (
double)
iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec *
iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo], 2.) +
687 pow( SigmaBranchRatio *
iso_sp[ipISO][nelem].fb[ipHi].RadEffec, 2.) );
692 fprintf(
ioQQQ,
"\t%e\t%e\t%e\t%e\t%e\t%e\n",
695 iso_sp[ipISO][nelem].fb[ipHi].RadEffec,
696 iso_sp[ipISO][nelem].fb[ipHi].SigmaRadEffec,
697 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo],
840 double RadRecombReturn;
841 long int i, i1, i2, i3, i4, i5;
842 long int ipLo, nelem;
847 const char* chFilename[
NISO] =
848 {
"h_iso_recomb.dat",
"he_iso_recomb.dat" };
856 if(
iso_ctrl.lgCompileRecomb[ipISO] )
858 iso_ctrl.lgNoRecombInterp[ipISO] =
false;
861 if( !
iso_ctrl.lgNoRecombInterp[ipISO] )
867 if( !
iso_ctrl.lgCompileRecomb[ipISO] )
870 fprintf(
ioQQQ,
" iso_recomb_setup opening %s:", chFilename[ipISO] );
875 ioDATA =
open_data( chFilename[ipISO],
"r" );
879 fprintf(
ioQQQ,
" Defaulting to on-the-fly computation, ");
880 fprintf(
ioQQQ,
" but the code runs much faster if you compile %s!\n", chFilename[ipISO]);
886 for( nelem = ipISO; nelem <
LIMELM; nelem++ )
888 if(
dense.lgElmtOn[nelem] )
899 for( ipLo=0; ipLo <
iso_sp[ipISO][nelem].numLevels_max-
iso_sp[ipISO][nelem].nCollapsed_max; ipLo++ )
907 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
931 fprintf(
ioQQQ,
" iso_recomb_setup could not read first line of %s.\n", chFilename[ipISO]);
935 i1 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
936 i2 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
937 i3 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
938 i4 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
942 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
944 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
950 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
952 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
958 for( nelem = ipISO; nelem <
LIMELM; nelem++ )
960 for( ipLo=0; ipLo <=
NumLevRecomb[ipISO][nelem]; ipLo++ )
966 fprintf(
ioQQQ,
" iso_recomb_setup could not read line %li of %s.\n", i5,
972 i1 = (long)
FFmtRead(chLine,&i3,
sizeof(chLine),&lgEOL);
973 i2 = (long)
FFmtRead(chLine,&i3,
sizeof(chLine),&lgEOL);
975 if( i1!=nelem || i2!=ipLo )
977 fprintf(
ioQQQ,
" iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
979 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
988 if( nelem == ipISO ||
dense.lgElmtOn[nelem] )
996 RRCoef[ipISO][nelem][ipLo][i] = ThisCoef;
1001 fprintf(
ioQQQ,
" iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
1003 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
1012 if( nelem == ipISO ||
dense.lgElmtOn[nelem] )
1014 for( ipLo=
NumLevRecomb[ipISO][nelem]; ipLo <
iso_sp[ipISO][nelem].numLevels_max-
iso_sp[ipISO][nelem].nCollapsed_max; ipLo++ )
1026 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
1028 fprintf(
ioQQQ,
" iso_recomb_setup could not read last line of %s.\n", chFilename[ipISO]);
1032 i1 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
1033 i2 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
1034 i3 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
1035 i4 = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
1040 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
1042 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
1048 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
1050 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
1059 else if(
iso_ctrl.lgCompileRecomb[ipISO] )
1071 fprintf(ioRECOMB,
"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient H-LIke [or HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
1082 for( nelem = ipISO; nelem <
LIMELM; nelem++ )
1095 fprintf(ioRECOMB,
"%li\t%li", nelem, ipLo );
1102 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
1103 fprintf(ioRECOMB,
"\t%f",
RRCoef[ipISO][nelem][ipLo][i] );
1105 fprintf(ioRECOMB,
"\n" );
1111 fprintf(ioRECOMB,
"%li\t%li", nelem,
NumLevRecomb[ipISO][nelem] );
1122 fprintf(ioRECOMB,
"\t%f", log10(
TotalRecomb[ipISO][nelem][i] ) );
1124 fprintf(ioRECOMB,
"\n" );
1127 fprintf(ioRECOMB,
"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient [H-LIke/HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
1140 fprintf(
ioQQQ,
"iso_recomb_setup: compilation complete, %s created.\n", chFilename[ipISO] );
1141 fprintf(
ioQQQ,
"The compilation is completed successfully.\n");