44 static int nCalled = 0;
73 (*TauDummy).AddHiState();
74 (*TauDummy).AddLoState();
75 (*TauDummy).AddLine2Stack();
83 for(
long nelem=ipISO; nelem <
LIMELM; nelem++ )
86 if( nelem < 2 ||
dense.lgElmtOn[nelem] )
94 double EnergyRydGround = 0.;
96 for( ipHi=0; ipHi <
iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
98 double EnergyWN, EnergyRyd;
102 EnergyRyd = HIonPoten/
POW2((
double)
N_(ipHi));
117 iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd = EnergyRyd;
119 EnergyRydGround = EnergyRyd;
120 iso_sp[ipISO][nelem].st[ipHi].energy().set(EnergyRydGround-EnergyRyd);
123 for( ipLo=0; ipLo < ipHi; ipLo++ )
125 EnergyWN =
RYD_INF * (
iso_sp[ipISO][nelem].fb[ipLo].xIsoLevNIonRyd -
126 iso_sp[ipISO][nelem].fb[ipHi].xIsoLevNIonRyd);
134 EnergyWN = -1.0 * EnergyWN;
137 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() = (
realnum)EnergyWN;
139 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() >= 0.);
140 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() >= 0.);
141 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyK() >= 0.);
146 iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() = 0.;
151 iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() =
153 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN()/
155 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng() > 0.);
162 for( ipHi=2; ipHi <
iso_ctrl.nLyman_malloc[ipISO]; ipHi++ )
201 for(
long nelem=ipISO; nelem <
LIMELM; nelem++ )
204 if( nelem < 2 ||
dense.lgElmtOn[nelem] )
206 for( ipLo=
ipH1s; ipLo < (
iso_sp[ipISO][nelem].numLevels_max - 1); ipLo++ )
208 for( ipHi=ipLo + 1; ipHi <
iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
227 iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipEmis() = -1;
229 iso_sp[ipISO][nelem].trans(ipHi,ipLo).AddLine2Stack();
231 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() = Aul;
233 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > 0.);
235 if( ipLo == 0 && ipHi ==
iso_ctrl.nLyaLevel[ipISO] )
237 long redis =
iso_ctrl.ipLyaRedist[ipISO];
241 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() = redis;
247 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() =
iso_ctrl.ipResoRedist[ipISO];
253 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().iRedisFun() =
iso_ctrl.ipSubRedist[ipISO];
256 if(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <=
iso_ctrl.SmallA ||
257 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0.)
259 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() = 0.;
260 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() = 0.;
264 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() =
266 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
267 iso_sp[ipISO][nelem].st[ipHi].
g()));
268 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().gf() > 0.);
271 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() =
273 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN(),
274 iso_sp[ipISO][nelem].st[ipLo].
g()));
275 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() > 0.);
292 if( nelem < 2 ||
dense.lgElmtOn[nelem] )
296 for( ipLo=
ipHe1s1S; ipLo<ipHi; ipLo++ )
321 for(
long nelem=ipISO; nelem <
LIMELM; nelem++ )
324 if( nelem < 2 ||
dense.lgElmtOn[nelem] )
327 iso_sp[ipISO][nelem].st[0].lifetime() = -FLT_MAX;
329 for( ipHi=1; ipHi <
iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
333 for( ipLo=0; ipLo < ipHi; ipLo++ )
335 if(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <=
iso_ctrl.SmallA )
338 iso_sp[ipISO][nelem].st[ipHi].lifetime() +=
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
342 iso_sp[ipISO][nelem].st[ipHi].lifetime() = 1./
iso_sp[ipISO][nelem].st[ipHi].lifetime();
344 for( ipLo=0; ipLo < ipHi; ipLo++ )
346 if(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() <= 0. )
349 if(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <=
iso_ctrl.SmallA )
352 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel() = (
realnum)(
353 (1.f/
iso_sp[ipISO][nelem].st[ipHi].lifetime())/
354 PI4/
iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN());
356 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().dampXvel()> 0.);
369 for(
long nelem = ipISO; nelem <
LIMELM; nelem++ )
372 if( nelem == ipISO ||
dense.lgElmtOn[nelem] )
390 for(
long nelem=ipISO; nelem <
LIMELM; ++nelem )
392 if( nelem < 2 ||
dense.lgElmtOn[nelem] )
394 for(
long ipHi= 1; ipHi <
iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
396 for(
long ipLo=0; ipLo < ipHi; ++ipLo )
398 iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk = 0.;
399 iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up = 0.;
725 for( in = 1L; in <=
iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
727 for( il = 0L; il < in; ++il )
729 iso_sp[ipISO][nelem].st[i].n() = in;
730 iso_sp[ipISO][nelem].st[i].S() = is;
731 iso_sp[ipISO][nelem].st[i].l() = il;
732 iso_sp[ipISO][nelem].st[i].j() = -1;
733 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
738 in =
iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
739 for( level = i; level<
iso_sp[ipISO][nelem].numLevels_max; ++level)
741 iso_sp[ipISO][nelem].st[level].n() = in;
742 iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
743 iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
744 iso_sp[ipISO][nelem].st[level].j() = -1;
746 for( il = 0; il < in; ++il )
748 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
758 ASSERT( (in > 0) && (in < (
iso_sp[ipISO][nelem].n_HighestResolved_max +
iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
761 for( in = 2L; in <=
iso_sp[ipISO][nelem].n_HighestResolved_max +
iso_sp[ipISO][nelem].nCollapsed_max; ++in )
763 for( il = 0L; il < in; ++il )
765 ASSERT(
iso_sp[ipISO][nelem].st[
iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
766 if( in <=
iso_sp[ipISO][nelem].n_HighestResolved_max )
770 ASSERT(
iso_sp[ipISO][nelem].st[
iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
771 ASSERT(
iso_sp[ipISO][nelem].st[
iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].
S() == is );
788 for( in = 1L; in <=
iso_sp[ipISO][nelem].n_HighestResolved_max; ++in )
790 for( il = 0L; il < in; ++il )
792 for( is = 3L; is >= 1L; is -= 2 )
797 if( (il == 1L) && (is == 1L) )
800 if( (in == 1L) && (is == 3L) )
806 iso_sp[ipISO][nelem].st[i].n() = in;
807 iso_sp[ipISO][nelem].st[i].S() = is;
808 iso_sp[ipISO][nelem].st[i].l() = il;
810 iso_sp[ipISO][nelem].st[i].j() = il;
811 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
815 else if( (in == 2) && (il == 1) && (is == 3) )
820 iso_sp[ipISO][nelem].st[i].n() = in;
821 iso_sp[ipISO][nelem].st[i].S() = is;
822 iso_sp[ipISO][nelem].st[i].l() = il;
823 iso_sp[ipISO][nelem].st[i].j() = ij;
824 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
832 iso_sp[ipISO][nelem].st[i].n() = in;
833 iso_sp[ipISO][nelem].st[i].S() = is;
834 iso_sp[ipISO][nelem].st[i].l() = il;
835 iso_sp[ipISO][nelem].st[i].j() = -1L;
836 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = i;
844 iso_sp[ipISO][nelem].st[i].n() = in;
845 iso_sp[ipISO][nelem].st[i].S() = 1L;
846 iso_sp[ipISO][nelem].st[i].l() = 1L;
847 iso_sp[ipISO][nelem].st[i].j() = 1L;
848 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][1][1] = i;
853 in =
iso_sp[ipISO][nelem].n_HighestResolved_max + 1;
854 for( level = i; level<
iso_sp[ipISO][nelem].numLevels_max; ++level)
856 iso_sp[ipISO][nelem].st[level].n() = in;
857 iso_sp[ipISO][nelem].st[level].S() = -LONG_MAX;
858 iso_sp[ipISO][nelem].st[level].l() = -LONG_MAX;
859 iso_sp[ipISO][nelem].st[level].j() = -1;
861 for( il = 0; il < in; ++il )
863 for( is = 1; is <= 3; is += 2 )
865 iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] = level;
876 ASSERT( (in > 0) && (in < (
iso_sp[ipISO][nelem].n_HighestResolved_max +
iso_sp[ipISO][nelem].nCollapsed_max + 1) ) );
879 for( in = 2L; in <=
iso_sp[ipISO][nelem].n_HighestResolved_max +
iso_sp[ipISO][nelem].nCollapsed_max; ++in )
881 for( il = 0L; il < in; ++il )
883 for( is = 3L; is >= 1; is -= 2 )
886 if( (in == 1L) && (is == 3L) )
889 ASSERT(
iso_sp[ipISO][nelem].st[
iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].n() == in );
890 if( in <=
iso_sp[ipISO][nelem].n_HighestResolved_max )
894 ASSERT(
iso_sp[ipISO][nelem].st[
iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].l() == il );
895 ASSERT(
iso_sp[ipISO][nelem].st[
iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is] ].
S() == is );
905 for(
long nelem=ipISO; nelem <
LIMELM; nelem++ )
908 if( nelem < 2 ||
dense.lgElmtOn[nelem] )
910 for( ipLo=
ipH1s; ipLo <
iso_sp[ipISO][nelem].numLevels_max; ipLo++ )
912 iso_sp[ipISO][nelem].st[ipLo].nelem() = (int)(nelem+1);
913 iso_sp[ipISO][nelem].st[ipLo].IonStg() = (int)(nelem+1-ipISO);
915 if(
iso_sp[ipISO][nelem].st[ipLo].j() >= 0 )
917 iso_sp[ipISO][nelem].st[ipLo].g() = 2.f*
iso_sp[ipISO][nelem].st[ipLo].j()+1.f;
919 else if(
iso_sp[ipISO][nelem].st[ipLo].l() >= 0 )
921 iso_sp[ipISO][nelem].st[ipLo].g() = (2.f*
iso_sp[ipISO][nelem].st[ipLo].l()+1.f) *
922 iso_sp[ipISO][nelem].st[ipLo].S();
936 char chConfiguration[11] =
" ";
937 long nCharactersWritten = 0;
942 if(
iso_sp[ipISO][nelem].st[ipLo].n() >
iso_sp[ipISO][nelem].n_HighestResolved_max )
944 nCharactersWritten = sprintf( chConfiguration,
"n=%3li",
945 iso_sp[ipISO][nelem].st[ipLo].n() );
947 else if(
iso_sp[ipISO][nelem].st[ipLo].j() > 0 )
949 nCharactersWritten = sprintf( chConfiguration,
"%3li^%li%c_%li",
950 iso_sp[ipISO][nelem].st[ipLo].n(),
951 iso_sp[ipISO][nelem].st[ipLo].
S(),
953 iso_sp[ipISO][nelem].st[ipLo].j() );
957 nCharactersWritten = sprintf( chConfiguration,
"%3li^%li%c",
958 iso_sp[ipISO][nelem].st[ipLo].n(),
959 iso_sp[ipISO][nelem].st[ipLo].
S(),
963 ASSERT( nCharactersWritten <= 10 );
964 chConfiguration[10] =
'\0';
966 strncpy(
iso_sp[ipISO][nelem].st[ipLo].chConfig(), chConfiguration, 10 );
1107 long int i, j, ipLo, ipHi;
1111 SumAPerN = ((
double*)
MALLOC( (
size_t)(
iso_sp[ipISO][nelem].n_HighestResolved_max +
iso_sp[ipISO][nelem].nCollapsed_max + 1 )*
sizeof(
double )));
1112 memset( SumAPerN, 0, (
iso_sp[ipISO][nelem].n_HighestResolved_max +
iso_sp[ipISO][nelem].nCollapsed_max + 1 )*
sizeof(
double ) );
1115 iso_sp[ipISO][nelem].CascadeProb[0][0] = 1.;
1118 iso_sp[ipISO][nelem].fb[0].SigmaAtot = 0.;
1119 iso_sp[ipISO][nelem].ex[0][0].SigmaCascadeProb = 0.;
1125 for( ipHi=1; ipHi<
iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
1135 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipHi] = 1.;
1136 iso_sp[ipISO][nelem].CascadeProb[ipHi][0] = 0.;
1137 iso_sp[ipISO][nelem].BranchRatio[ipHi][0] = 0.;
1141 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = 0.;
1142 iso_sp[ipISO][nelem].ex[ipHi][ipHi].SigmaCascadeProb = 0.;
1149 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1151 SumAs +=
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1154 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ )
1156 if(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <=
iso_ctrl.SmallA )
1158 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1159 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] = 0.;
1163 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] = 0.;
1164 iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] =
1165 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() / SumAs;
1167 ASSERT(
iso_sp[ipISO][nelem].BranchRatio[ipHi][ipLo] <= 1.0000001 );
1169 SumAPerN[
N_(ipHi)] +=
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul();
1174 ASSERT(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyWN() > 0. ||
1175 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <=
iso_ctrl.SmallA );
1181 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot +=
1183 (
double)
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul(), 2. );
1190 iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot = sqrt(
iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot );
1194 for( i=0; i<ipHi; i++ )
1196 for( ipLo=0; ipLo<=i; ipLo++ )
1198 iso_sp[ipISO][nelem].CascadeProb[ipHi][ipLo] +=
iso_sp[ipISO][nelem].BranchRatio[ipHi][i] *
iso_sp[ipISO][nelem].CascadeProb[i][ipLo];
1204 for( ipLo=0; ipLo<ipHi; ipLo++ )
1206 double SigmaCul = 0.;
1207 for( i=ipLo; i<ipHi; i++ )
1209 if(
iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul() >
iso_ctrl.SmallA )
1212 double SigmaA =
iso_sp[ipISO][nelem].ex[ipHi][i].Error[
IPRAD] *
1213 iso_sp[ipISO][nelem].trans(ipHi,i).Emis().Aul();
1215 pow(SigmaA*
iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*
iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
1216 pow(
iso_sp[ipISO][nelem].fb[ipHi].SigmaAtot*
iso_sp[ipISO][nelem].BranchRatio[ipHi][i]*
1217 iso_sp[ipISO][nelem].CascadeProb[i][ipLo]*
iso_sp[ipISO][nelem].st[ipHi].lifetime(), 2.) +
1218 pow(
iso_sp[ipISO][nelem].
ex[i][ipLo].SigmaCascadeProb*
iso_sp[ipISO][nelem].BranchRatio[ipHi][i], 2.);
1221 SigmaCul = sqrt(SigmaCul);
1222 iso_sp[ipISO][nelem].ex[ipHi][ipLo].SigmaCascadeProb = SigmaCul;
1231 enum {DEBUG_LOC=
false};
1255 fprintf(
ioQQQ,
"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
1256 fprintf(
ioQQQ,
"m\t2\t\t3\t\t4\t\t5\t\t6\n");
1258 for( ipHi=
ipHe2p3P2; ipHi<
iso_sp[ipISO][nelem].numLevels_max-
iso_sp[ipISO][nelem].nCollapsed_max; ipHi++ )
1261 if( (
iso_sp[ipISO][nelem].st[ipHi].l() == 1) && (
iso_sp[ipISO][nelem].st[ipHi].
S() == 3) )
1263 fprintf(
ioQQQ,
"\n%ld\t",
iso_sp[ipISO][nelem].st[ipHi].n());
1266 for( i = ipLo; i<=ipHi; i++)
1268 if( (
iso_sp[ipISO][nelem].st[i].l() == hi_l) && (
iso_sp[ipISO][nelem].st[i].
S() == hi_s) )
1272 Bm +=
iso_sp[ipISO][nelem].CascadeProb[ipHi][i] * (
iso_sp[ipISO][nelem].BranchRatio[i][
ipHe2p3P0] +
1276 Bm +=
iso_sp[ipISO][nelem].CascadeProb[ipHi][i] *
iso_sp[ipISO][nelem].BranchRatio[i][ipLo];
1283 fprintf(
ioQQQ,
"%2.4e\t",Bm);
1289 fprintf(
ioQQQ,
"%2.4e\t",Bm);
1296 fprintf(
ioQQQ,
"\n\n");
1304 for( i=2; i <
iso_sp[ipISO][nelem].n_HighestResolved_max; ++i)
1306 ASSERT( (SumAPerN[i] > SumAPerN[i+1]) ||
opac.lgCaseB );
1310 enum {DEBUG_LOC=
false};
1313 for( i = 2; i<= (
iso_sp[ipISO][nelem].n_HighestResolved_max +
iso_sp[ipISO][nelem].nCollapsed_max); ++i)
1315 fprintf(
ioQQQ,
"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[i]);
1524 double bnl_array[4][3][4][10] = {
1528 {6.13E-01, 2.56E-01, 1.51E-01, 2.74E-01, 3.98E-01, 4.98E-01, 5.71E-01, 6.33E-01, 7.28E-01, 9.59E-01},
1529 {1.31E+00, 5.17E-01, 2.76E-01, 4.47E-01, 5.87E-01, 6.82E-01, 7.44E-01, 8.05E-01, 9.30E-01, 1.27E+00},
1530 {1.94E+00, 7.32E-01, 3.63E-01, 5.48E-01, 6.83E-01, 7.66E-01, 8.19E-01, 8.80E-01, 1.02E+00, 1.43E+00},
1531 {2.53E+00, 9.15E-01, 4.28E-01, 6.16E-01, 7.42E-01, 8.13E-01, 8.60E-01, 9.22E-01, 1.08E+00, 1.56E+00}
1534 {5.63E-01, 2.65E-01, 1.55E-01, 2.76E-01, 3.91E-01, 4.75E-01, 5.24E-01, 5.45E-01, 5.51E-01, 5.53E-01},
1535 {1.21E+00, 5.30E-01, 2.81E-01, 4.48E-01, 5.80E-01, 6.62E-01, 7.05E-01, 7.24E-01, 7.36E-01, 7.46E-01},
1536 {1.81E+00, 7.46E-01, 3.68E-01, 5.49E-01, 6.78E-01, 7.51E-01, 7.88E-01, 8.09E-01, 8.26E-01, 8.43E-01},
1537 {2.38E+00, 9.27E-01, 4.33E-01, 6.17E-01, 7.38E-01, 8.05E-01, 8.40E-01, 8.65E-01, 8.92E-01, 9.22E-01}
1540 {2.97E-01, 2.76E-01, 2.41E-01, 3.04E-01, 3.66E-01, 4.10E-01, 4.35E-01, 4.48E-01, 4.52E-01, 4.53E-01},
1541 {5.63E-01, 5.04E-01, 3.92E-01, 4.67E-01, 5.39E-01, 5.85E-01, 6.10E-01, 6.20E-01, 6.23E-01, 6.23E-01},
1542 {7.93E-01, 6.90E-01, 4.94E-01, 5.65E-01, 6.36E-01, 6.79E-01, 7.00E-01, 7.09E-01, 7.11E-01, 7.11E-01},
1543 {1.04E+00, 8.66E-01, 5.62E-01, 6.31E-01, 7.01E-01, 7.43E-01, 7.63E-01, 7.71E-01, 7.73E-01, 7.73E-01}
1549 {6.70E-02, 2.93E-02, 1.94E-02, 4.20E-02, 7.40E-02, 1.12E-01, 1.51E-01, 1.86E-01, 2.26E-01, 3.84E-01},
1550 {2.39E-01, 1.03E-01, 6.52E-02, 1.31E-01, 2.11E-01, 2.91E-01, 3.61E-01, 4.17E-01, 4.85E-01, 8.00E-01},
1551 {4.26E-01, 1.80E-01, 1.10E-01, 2.09E-01, 3.18E-01, 4.15E-01, 4.93E-01, 5.54E-01, 6.34E-01, 1.04E+00},
1552 {6.11E-01, 2.55E-01, 1.51E-01, 2.74E-01, 3.99E-01, 5.02E-01, 5.80E-01, 6.41E-01, 7.30E-01, 1.21E+00}
1555 {6.79E-02, 3.00E-02, 2.00E-02, 4.30E-02, 7.48E-02, 1.11E-01, 1.44E-01, 1.70E-01, 1.87E-01, 1.96E-01},
1556 {2.40E-01, 1.04E-01, 6.62E-02, 1.32E-01, 2.11E-01, 2.87E-01, 3.51E-01, 3.98E-01, 4.32E-01, 4.58E-01},
1557 {4.26E-01, 1.81E-01, 1.11E-01, 2.10E-01, 3.17E-01, 4.12E-01, 4.89E-01, 5.53E-01, 6.14E-01, 6.84E-01},
1558 {6.12E-01, 2.55E-01, 1.51E-01, 2.73E-01, 3.97E-01, 4.98E-01, 5.77E-01, 6.51E-01, 7.82E-01, 1.18E+00}
1561 {4.98E-02, 3.47E-02, 2.31E-02, 4.54E-02, 7.14E-02, 9.37E-02, 1.08E-01, 1.13E-01, 1.13E-01, 1.11E-01},
1562 {1.75E-01, 1.16E-01, 7.36E-02, 1.36E-01, 2.01E-01, 2.50E-01, 2.76E-01, 2.84E-01, 2.81E-01, 2.77E-01},
1563 {3.38E-01, 1.97E-01, 1.18E-01, 2.13E-01, 3.06E-01, 3.72E-01, 4.06E-01, 4.15E-01, 4.10E-01, 4.04E-01},
1564 {6.01E-01, 2.60E-01, 1.53E-01, 2.76E-01, 3.95E-01, 4.87E-01, 5.45E-01, 5.76E-01, 5.93E-01, 6.05E-01}
1570 {1.77E-01, 3.59E-01, 1.54E-01, 2.75E-01, 3.98E-01, 4.94E-01, 5.51E-01, 5.68E-01, 5.46E-01, 4.97E-01},
1571 {4.09E-01, 7.23E-01, 2.83E-01, 4.48E-01, 5.89E-01, 6.78E-01, 7.22E-01, 7.30E-01, 7.07E-01, 6.65E-01},
1572 {6.40E-01, 1.02E+00, 3.74E-01, 5.49E-01, 6.85E-01, 7.63E-01, 7.98E-01, 8.03E-01, 7.84E-01, 7.53E-01},
1573 {8.70E-01, 1.28E+00, 4.42E-01, 6.17E-01, 7.44E-01, 8.13E-01, 8.42E-01, 8.46E-01, 8.34E-01, 8.13E-01}
1576 {1.78E-01, 3.62E-01, 1.55E-01, 2.73E-01, 3.91E-01, 4.73E-01, 5.10E-01, 5.04E-01, 4.70E-01, 4.32E-01},
1577 {4.08E-01, 7.26E-01, 2.83E-01, 4.45E-01, 5.79E-01, 6.54E-01, 6.78E-01, 6.64E-01, 6.30E-01, 5.98E-01},
1578 {6.37E-01, 1.03E+00, 3.73E-01, 5.46E-01, 6.75E-01, 7.40E-01, 7.57E-01, 7.43E-01, 7.15E-01, 6.92E-01},
1579 {8.65E-01, 1.28E+00, 4.41E-01, 6.14E-01, 7.35E-01, 7.92E-01, 8.05E-01, 7.95E-01, 7.74E-01, 7.59E-01}
1582 {2.07E-01, 3.73E-01, 1.73E-01, 2.85E-01, 4.03E-01, 4.76E-01, 5.06E-01, 5.03E-01, 4.84E-01, 4.63E-01},
1583 {4.32E-01, 7.13E-01, 3.06E-01, 4.54E-01, 5.81E-01, 6.44E-01, 6.59E-01, 6.49E-01, 6.28E-01, 6.11E-01},
1584 {6.40E-01, 9.85E-01, 3.98E-01, 5.53E-01, 6.74E-01, 7.27E-01, 7.36E-01, 7.26E-01, 7.10E-01, 6.98E-01},
1585 {8.38E-01, 1.21E+00, 4.67E-01, 6.20E-01, 7.34E-01, 7.79E-01, 7.87E-01, 7.79E-01, 7.69E-01, 7.63E-01}
1591 {9.31E-02, 3.96E-01, 1.36E-01, 2.74E-01, 3.99E-01, 4.95E-01, 5.52E-01, 5.70E-01, 5.48E-01, 4.96E-01},
1592 {2.25E-01, 8.46E-01, 2.49E-01, 4.46E-01, 5.89E-01, 6.79E-01, 7.23E-01, 7.31E-01, 7.08E-01, 6.64E-01},
1593 {3.59E-01, 1.24E+00, 3.30E-01, 5.47E-01, 6.85E-01, 7.63E-01, 7.98E-01, 8.04E-01, 7.85E-01, 7.53E-01},
1594 {4.93E-01, 1.60E+00, 3.91E-01, 6.15E-01, 7.44E-01, 8.13E-01, 8.42E-01, 8.47E-01, 8.35E-01, 8.12E-01}
1597 {9.32E-02, 3.99E-01, 1.35E-01, 2.72E-01, 3.91E-01, 4.75E-01, 5.14E-01, 5.09E-01, 4.73E-01, 4.31E-01},
1598 {2.25E-01, 8.49E-01, 2.49E-01, 4.44E-01, 5.79E-01, 6.56E-01, 6.81E-01, 6.68E-01, 6.31E-01, 5.96E-01},
1599 {3.58E-01, 1.25E+00, 3.29E-01, 5.44E-01, 6.76E-01, 7.42E-01, 7.60E-01, 7.46E-01, 7.16E-01, 6.91E-01},
1600 {4.92E-01, 1.60E+00, 3.90E-01, 6.12E-01, 7.36E-01, 7.93E-01, 8.07E-01, 7.97E-01, 7.74E-01, 7.58E-01}
1603 {1.13E-01, 4.21E-01, 1.47E-01, 2.83E-01, 4.04E-01, 4.80E-01, 5.13E-01, 5.12E-01, 4.93E-01, 4.71E-01},
1604 {2.52E-01, 8.56E-01, 2.61E-01, 4.50E-01, 5.82E-01, 6.48E-01, 6.66E-01, 6.56E-01, 6.35E-01, 6.16E-01},
1605 {3.85E-01, 1.23E+00, 3.41E-01, 5.49E-01, 6.75E-01, 7.30E-01, 7.41E-01, 7.31E-01, 7.15E-01, 7.02E-01},
1606 {5.14E-01, 1.56E+00, 4.01E-01, 6.15E-01, 7.34E-01, 7.82E-01, 7.90E-01, 7.83E-01, 7.72E-01, 7.65E-01}
1611 double temps[4] = {5000., 10000., 15000., 20000. };
1612 double log_dens[3] = {2., 4., 6.};
1621 ASSERT( (ipTe >=0) && (ipTe < 3) );
1622 ASSERT( (ipDens >=0) && (ipDens < 2) );
1624 for(
long nHi=
iso_sp[ipISO][nelem].n_HighestResolved_max+1; nHi<=
iso_sp[ipISO][nelem].n_HighestResolved_max+
iso_sp[ipISO][nelem].nCollapsed_max; nHi++ )
1626 for(
long lHi=0; lHi<nHi; lHi++ )
1628 for(
long sHi=1; sHi<4; sHi++ )
1630 if( ipISO ==
ipH_LIKE && sHi != 2 )
1632 else if( ipISO ==
ipHE_LIKE && sHi != 1 && sHi != 3 )
1635 double bnl_at_lo_den, bnl_at_hi_den, bnl;
1636 double bnl_max, bnl_min, temp, dens;
1638 long ipL =
MIN2(9,lHi);
1663 temp =
MIN2( temps[3], temp );
1665 dens =
MAX2( log_dens[0], log10(
dense.eden) );
1666 dens =
MIN2( log_dens[2], dens );
1672 if( temp < temps[0] && dens < log_dens[0] )
1673 bnl = bnl_array[ip1][0][0][ipL];
1674 else if( temp < temps[0] && dens >= log_dens[2] )
1675 bnl = bnl_array[ip1][2][0][ipL];
1676 else if( temp >= temps[3] && dens < log_dens[0] )
1677 bnl = bnl_array[ip1][0][3][ipL];
1678 else if( temp >= temps[3] && dens >= log_dens[2] )
1679 bnl = bnl_array[ip1][2][3][ipL];
1682 bnl_at_lo_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
1683 (bnl_array[ip1][ipDens][ipTe+1][ipL] - bnl_array[ip1][ipDens][ipTe][ipL]) + bnl_array[ip1][ipDens][ipTe][ipL];
1685 bnl_at_hi_den = ( temp - temps[ipTe]) / (temps[ipTe+1] - temps[ipTe]) *
1686 (bnl_array[ip1][ipDens+1][ipTe+1][ipL] - bnl_array[ip1][ipDens+1][ipTe][ipL]) + bnl_array[ip1][ipDens+1][ipTe][ipL];
1688 bnl = ( dens - log_dens[ipDens]) / (log_dens[ipDens+1] - log_dens[ipDens]) *
1689 (bnl_at_hi_den - bnl_at_lo_den) + bnl_at_lo_den;
1694 bnl_max =
MAX4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
1695 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
1696 ASSERT( bnl <= bnl_max );
1698 bnl_min =
MIN4( bnl_array[ip1][ipDens][ipTe+1][ipL], bnl_array[ip1][ipDens+1][ipTe+1][ipL],
1699 bnl_array[ip1][ipDens][ipTe][ipL], bnl_array[ip1][ipDens+1][ipTe][ipL] );
1700 ASSERT( bnl >= bnl_min );
1703 iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] = bnl;
1705 ASSERT(
iso_sp[ipISO][nelem].bnl_effective[nHi][lHi][sHi] > 0. );