24#define LIM_H2_POP_LOOP 10
27#define H2_DISS_ALLISON_DALGARNO 6e-19f
70 double source_so_far = 0.;
71 double source_so_far_s = 0.;
72 double sink_so_far = 0.;
75 double pop_tot_s = 0.;
111 pop_tot_s +=
states[ipHi].Pop();
117 pop_tot +=
states[ipHi].Pop();
122 double sink_tot =
mole.sink_rate_tot(
sp) * pop_tot;
124 double sink_left = sink_tot - sink_so_far;
127 sink_left /= pop_tot;
128 if( sink_left >= 0. )
141 double sink_tot_s =
mole.sink_rate_tot(
sp_star) * pop_tot_s;
143 double sink_left_s = sink_tot_s - sink_so_far_s;
146 sink_left_s /= pop_tot_s;
149 if( sink_left_s >= 0. )
159 double source_tot =
mole.source_rate_tot(
sp);
160 double source_left = source_tot - source_so_far;
161 double source_tot_s =
mole.source_rate_tot(
sp_star);
162 double source_left_s = source_tot_s - source_so_far_s;
163 if( source_left+source_left_s >= 0. )
165 double rpop_lte = 1.;
166 double rpop_lte_s = 0.;
167 if (
hmi.lgLeiden_Keep_ipMH2s)
170 double pop_lte_s = 0.;
173 long iElec =
states[ipHi].n();
174 long iVib =
states[ipHi].v();
175 long iRot =
states[ipHi].J();
181 rpop_lte = 1./
SDIV(pop_lte);
182 rpop_lte_s = 1./
SDIV(pop_lte_s);
186 long iElec =
states[ipHi].n();
187 long iVib =
states[ipHi].v();
188 long iRot =
states[ipHi].J();
203 DEBUG_ENTRY(
"diatomics::H2_X_coll_rate_evaluate()" );
230 fprintf(
ioQQQ,
" Collider densities are:");
245 for(
long ipLo=0; ipLo<ipHi; ++ipLo )
248 double colldown = 0.;
288 ASSERT( (*tr).Emis().Aul() > 0. );
290 (*tr).Emis().ipFine() =
ipFineCont( (*tr).EnergyRyd());
311 ASSERT( (*tr).ipCont() > 0 );
312 drive += (*tr).Emis().pump() * (*tr).Emis().PopOpc() * (*tr).EnergyErg();
337 ASSERT( (*tr).ipCont() > 0 );
338 if( (*(*tr).Hi()).Pop() > smallfloat && (*tr).Emis().PopOpc() > smallfloat )
346 " H2_RadPress returns, radiation pressure is %.2e\n",
364 energy += st->Pop() * st->energy();
383 (*tr).outline_resonance();
433 ASSERT( (*tr).ipCont() > 0 );
483 double rot_cooling , dCoolDT;
500 " The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n",
548 fprintf(
ioQQQ,
" H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n");
620 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
625 ASSERT( (*tr).ipCont() > 0 );
632 AulEscp[ihi][ilo] = (*tr).Emis().Aul()*(
633 (*tr).Emis().Pesc() +
634 (*tr).Emis().Pelec_esc());
635 AulDest[ihi][ilo] = (*tr).Emis().Aul()*(*tr).Emis().Pdest();
636 AulPump[ilo][ihi] = (*tr).Emis().pump();
650 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
655 ASSERT( (*tr).ipCont() > 0 );
661 (*(*tr).Hi()).Pop() * (
662 (*tr).Emis().Aul()*( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() ) +
663 (*tr).Emis().pump() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g() );
696 if(lgDeBug)fprintf(
ioQQQ,
"DEBUG H2_Level_low_matrix, ilo=%li",ilo);
712 if(lgDeBug)fprintf(
ioQQQ,
"\n");
724 if(lgDeBug)fprintf(
ioQQQ,
"\t%.1e",ratein);
727 double rateout = ratein *
735 if(lgDeBug)fprintf(
ioQQQ,
"\n");
771 enum {DEBUG_LOC=
false};
772 if( DEBUG_LOC || lgDeBug)
774 fprintf(
ioQQQ,
"DEBUG H2 matexcit");
777 fprintf(
ioQQQ,
"\t%li",ilo );
791 fprintf(
ioQQQ,
"AulEscp[n][]\\[][n] = Aul*Pesc\n");
794 fprintf(
ioQQQ,
"\t%li",ilo );
799 fprintf(
ioQQQ,
"%li", ihi);
807 fprintf(
ioQQQ,
"AulPump [n][]\\[][n]\n");
810 fprintf(
ioQQQ,
"\t%li",ilo );
815 fprintf(
ioQQQ,
"%li", ihi);
823 fprintf(
ioQQQ,
"CollRate_levn [n][]\\[][n]\n");
826 fprintf(
ioQQQ,
"\t%li",ilo );
831 fprintf(
ioQQQ,
"%li", ihi);
838 fprintf(
ioQQQ,
"SOURCE");
843 fprintf(
ioQQQ,
"\nSINK");
888 fprintf(
ioQQQ,
"\n DEBUG H2_Level_lowJ dense_total: %.3e matrix rel pops\n", *
dense_total);
889 fprintf(
ioQQQ,
"v\tJ\tpop\n");
894 fprintf(
ioQQQ,
"%3li\t%3li\t%.3e\t%.3e\t%.3e\n",
902 fprintf(
ioQQQ,
" H2_Level_low_matrix called atom_levelN which returned negative populations.\n");
922 double old_solomon_rate=-1.;
923 long int n_pop_oscil = 0;
930 lgOrthoParaRatioConv;
931 double quant_old=-1.,
934 bool lgH2_pops_oscil=
false,
935 lgH2_pops_ever_oscil=
false;
938 double PopChgMax_relative=0. , PopChgMaxOld_relative=0., PopChgMax_total=0., PopChgMaxOld_total=0.;
939 long int iRotMaxChng_relative , iVibMaxChng_relative,
940 iRotMaxChng_total , iVibMaxChng_total,
942 double popold_relative , popnew_relative , popold_total , popnew_total;
947 double converge_pops_relative=1e-2,
948 converge_pops_total=1e-3,
949 converge_ortho_para=1e-2;
951 double dens_rel_to_lim_react =
mole.species[
sp->index].xFracLim;
956 "\n***************H2_LevelPops %s call %li this iteration, zone is %.2f, H2/H:%.e Te:%e ne:%e\n",
960 dens_rel_to_lim_react,
967 static long int nzone_prt=-1;
968 if(
nzone!=nzone_prt )
971 fprintf(
ioQQQ,
"DEBUG zone %li species %s rel_to_lim:%.3e Te:%.3e *ne:%.3e n(%s):%.3e\n",
974 dens_rel_to_lim_react,
997 " H2_LevelPops %s pops too small, not computing, set to LTE and return, H2/H is %.2e and H2_to_H_limit is %.2e.",
999 dens_rel_to_lim_react,
1046 fprintf(
ioQQQ,
"%s 1st call - using LTE level pops\n",
label.c_str() );
1052 long iElec = (*st).n();
1053 long iVib = (*st).v();
1054 long iRot = (*st).J();
1089 double pop_total = 0.;
1092 long iElec = (*st).n();
1093 long iVib = (*st).v();
1095 pop_total += (*st).Pop();
1106 "%s H2_renorm_chemistry is %.4e, *dense_total is %.4e pops_per_elec[0] is %.4e\n",
1115 long iElec = (*st).n();
1116 long iVib = (*st).v();
1117 long iRot = (*st).J();
1125 " H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n",
1133 converge_pops_relative *= 2.;
1134 converge_pops_total *= 3.;
1135 converge_ortho_para *= 3.;
1138 if( !
conv.nTotalIoniz )
1149 lgConv_h2_soln =
false;
1151 long loop_h2_pops = 0;
1162 lgConv_h2_soln =
true;
1207 fprintf(
ioQQQ,
" Rel pop(e=%li)" ,iElecHi);
1219 long iElec = (*st).n();
1220 if( iElec > 0 )
continue;
1221 long iVib = (*st).v();
1236 fprintf(
ioQQQ,
"\n");
1238 fprintf(
ioQQQ,
" Rel pop(0,J)");
1241 long iElec = (*st).n();
1242 if( iElec > 0 )
continue;
1243 long iVib = (*st).v();
1244 if( iVib > 0 )
continue;
1247 fprintf(
ioQQQ,
"\n");
1251 double pop_total = 0.;
1253 pop_total += (*st).Pop();
1259 fprintf(
ioQQQ,
"DEBUG H2 %ld %ld %g %g %g %g %g\n",
1269 (*st).Pop() *= H2_renorm_conserve;
1270 long iElec = (*st).n();
1271 long iVib = (*st).v();
1280 double sum_pops_matrix = 0.;
1283 sum_pops_matrix +=
states[i].Pop();
1292 PopChgMaxOld_relative = PopChgMax_relative;
1293 PopChgMaxOld_total = PopChgMax_total;
1294 PopChgMax_relative = 0.;
1295 PopChgMax_total = 0.;
1296 iRotMaxChng_relative =-1;
1297 iVibMaxChng_relative = -1;
1298 iRotMaxChng_total =-1;
1299 iVibMaxChng_total = -1;
1300 popold_relative = 0.;
1301 popnew_relative = 0.;
1323 long iElec = (*st).n();
1324 long iVib = (*st).v();
1325 long iRot = (*st).J();
1334 SDIV(pop) > fabs(PopChgMax_relative) &&
1342 PopChgMax_relative =
1344 iRotMaxChng_relative = iRot;
1345 iVibMaxChng_relative = iVib;
1347 popnew_relative = pop;
1355 if( fabs(rel_change) > fabs(PopChgMax_total) )
1357 PopChgMax_total = rel_change;
1358 iRotMaxChng_total = iRot;
1359 iVibMaxChng_total = iVib;
1384 else if( rel_change> 0.1 )
1403 double H2_renorm_conserve_init = *
dense_total/sumold;
1410 long iElec = (*st).n();
1411 long iVib = (*st).v();
1412 long iRot = (*st).J();
1426 long iElec = (*st).n();
1427 long iVib = (*st).v();
1428 long iRot = (*st).J();
1429 const double& pop = (*st).Pop();
1472 if( loop_h2_pops>2 && (
1474 (PopChgMax_relative*PopChgMaxOld_relative<0. ) ) )
1476 lgH2_pops_oscil =
true;
1477 if( loop_h2_pops > 6 )
1480 lgH2_pops_ever_oscil =
true;
1486 lgH2_pops_oscil =
false;
1490 lgH2_pops_ever_oscil =
false;
1500 lgConv_h2_soln =
true;
1501 lgPopsConv_total =
true;
1502 lgPopsConv_relative =
true;
1504 lgSolomonConv =
true;
1505 lgOrthoParaRatioConv =
true;
1509 if( fabs(PopChgMax_relative)>converge_pops_relative )
1512 lgConv_h2_soln =
false;
1513 lgPopsConv_relative =
false;
1516 quant_old = PopChgMaxOld_relative;
1518 quant_new = PopChgMax_relative;
1520 strcpy( chReason ,
"rel pops changed" );
1525 else if( fabs(PopChgMax_total)>converge_pops_total)
1527 lgConv_h2_soln =
false;
1528 lgPopsConv_total =
false;
1531 quant_old = PopChgMaxOld_total;
1533 quant_new = PopChgMax_total;
1535 strcpy( chReason ,
"tot pops changed" );
1546 lgConv_h2_soln =
false;
1547 lgOrthoParaRatioConv =
false;
1550 strcpy( chReason ,
"ortho/para ratio changed" );
1554 else if( !
thermal.lgTemperatureConstant &&
1556 conv.HeatCoolRelErrorAllowed/10.
1563 lgConv_h2_soln =
false;
1567 strcpy( chReason ,
"heating changed" );
1577 else if(
rfield.lgInducProcess &&
1580 hmi.H2_frac_abund_set==0 &&
1584 conv.EdenErrorAllowed/5.)
1586 lgConv_h2_soln =
false;
1587 lgSolomonConv =
false;
1588 quant_old = old_solomon_rate;
1590 strcpy( chReason ,
"Solomon rate changed" );
1594 if( !lgConv_h2_soln )
1605 fprintf(
ioQQQ,
" %s loop %3li no conv oscl?%c why:%s ",
1608 TorF(lgH2_pops_ever_oscil),
1610 if( !lgPopsConv_relative )
1611 fprintf(
ioQQQ,
" PopChgMax_relative:%.4e v:%li J:%li old:%.4e new:%.4e",
1613 iVibMaxChng_relative,
1614 iRotMaxChng_relative ,
1617 else if( !lgPopsConv_total )
1618 fprintf(
ioQQQ,
" PopChgMax_total:%.4e v:%li J:%li old:%.4e new:%.4e",
1624 else if( !lgHeatConv )
1625 fprintf(
ioQQQ,
" heat:%.4e old:%.4e new:%.4e",
1630 else if( !lgSolomonConv )
1632 else if( !lgOrthoParaRatioConv )
1633 fprintf(
ioQQQ,
" current, old, older ratios are %.4e %.4e %.4e",
1637 fprintf(
ioQQQ,
"\n");
1642 if(
trace.nTrConvg >= 5 )
1645 " H2 5lev %li Conv?%c",
1647 TorF(lgConv_h2_soln) );
1649 if( fabs(PopChgMax_relative)>0.1 )
1650 fprintf(
ioQQQ,
" pops, rel chng %.3e",PopChgMax_relative);
1652 fprintf(
ioQQQ,
" rel heat %.3e rel chng %.3e H2 heat/cool %.2e",
1658 " Oscil?%c Ever Oscil?%c",
1659 TorF(lgH2_pops_oscil) ,
1660 TorF(lgH2_pops_ever_oscil) );
1661 fprintf(
ioQQQ,
"\n");
1667 "H2 loop\t%li\tkase pop chng\t%i\tchem renorm fac\t%.4e\tortho/para ratio:\t%.3e\tfrac of pop in matrix: %.3f\n",
1675 if( iVibMaxChng_relative>=0 && iRotMaxChng_relative>=0 && PopChgMax_relative>1e-10 )
1677 "end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n",
1679 PopChgMax_relative ,
1682 iVibMaxChng_relative , iRotMaxChng_relative
1689 if( !lgConv_h2_soln && !
conv.lgSearch )
1691 conv.lgConvPops =
false;
1692 lgPopsConverged =
false;
1693 old_val = quant_old;
1694 new_val = quant_new;
1699 ASSERT( (*st).Pop() >= 0. );
1705 (*tr).Coll().cool() = 0.;
1706 (*tr).Coll().heat() = 0.;
1708 (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
1711 (*tr).Emis().phots() = (*tr).Emis().Aul() * ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc()) * (*(*tr).Hi()).Pop();
1714 (*tr).Emis().xIntensity() = (*tr).Emis().phots() * (*tr).EnergyErg();
1723 double popTimesE = (*st).Pop() * (*st).energy().WN();
1757 static long ip_cut_off = -1;
1758 if( ip_cut_off < 0 )
1761 ip_cut_off =
ipoint( 1.14 );
1766 double flux_accum_photodissoc_BigH2_H2s = 0;
1768 long ip_H2_level =
ipoint( 1.07896 - 2.5 /
EVRYD);
1769 for(
long i= ip_H2_level; i < ip_cut_off; ++i )
1771 flux_accum_photodissoc_BigH2_H2s += (
rfield.flux[0][i-1] +
rfield.ConInterOut[i-1]+
1778 long iElec = (*st).n();
1779 if( iElec > 0 )
continue;
1780 long iVib = (*st).v();
1781 long iRot = (*st).J();
1782 const double &pop = (*st).Pop();
1784 const double mass_stat_factor = 3.634e-5/(2*2);
1801 if( arg_ratio > 0. )
1805 H2_stat[0][iVib][iRot] * mass_stat_factor;
1812 double flux_accum_photodissoc_BigH2_H2g = 0;
1814 ip_H2_level =
ipoint( 1.07896 - (*st).energy().Ryd() );
1816 for(
long i= ip_H2_level; i < ip_cut_off; ++i )
1818 flux_accum_photodissoc_BigH2_H2g += (
rfield.flux[0][i-1] +
rfield.ConInterOut[i-1]+
1833 if( arg_ratio > 0. )
1836 H2_stat[0][iVib][iRot] * mass_stat_factor;
1880 fprintf(
ioQQQ,
" H2_LevelPops exit1 %8.2f loops:%3li H2/H:%.3e Sol dis old %.3e new %.3e Sol dis star %.3e g-to-s %.3e photodiss star %.3e\n",
1883 dens_rel_to_lim_react,
1913 const double FRAC = 0.99999;
1915 double sum_pop = 0.;
1918 const bool PRT =
false;
1919 if( PRT ) fprintf(
ioQQQ,
"DEBUG pops ");
1930 if( PRT ) fprintf(
ioQQQ,
"\n");
1943 DEBUG_ENTRY(
"diatomics::SolveExcitedElectronicLevels()" );
1949 double CosmicRayHILyaExcitationRate = (
hmi.lgLeidenCRHack ) ?
secondaries.x12tot : 0.;
1954 long iElecLo = (*Lo).n();
1955 long iVibLo = (*Lo).v();
1956 long iRotLo = (*Lo).J();
1958 long iElecHi = (*Hi).n();
1959 if( iElecHi < 1 )
continue;
1960 long iVibHi = (*Hi).v();
1961 long iRotHi = (*Hi).J();
1968 double rate_up = (*tr).Emis().pump() + CosmicRayHILyaExcitationRate * (*tr).Coll().col_str();
1970 (*tr).Emis().Aul() * ( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() ) +
1971 rate_up * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
1977 rate_in[iElecHi][iVibHi][iRotHi] +=
H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_up;
1988 ASSERT( rate_up >= 0. && rate_down >= 0. );
1996 long iElec = (*st).n();
1997 long iVib = (*st).v();
1998 long iRot = (*st).J();
2027 fprintf(
ioQQQ,
" Pop(e=%li):",iElec);
2030 fprintf(
ioQQQ,
"\n");
2037 if( (*Lo).n() != 0 )
continue;
2039 if( (*Hi).n() < 1 )
continue;
2042 double rate = (*Hi).Pop() *
2043 ((*tr).Emis().Aul() * ( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() ) +
2044 (*tr).Emis().pump() * (*Lo).g() / (*Hi).g());
2054 DEBUG_ENTRY(
"diatomics::SolveSomeGroundElectronicLevels()");
2070 double H2boltz =
H2_Boltzmann[iElecHi][iVibHi][iRotHi];
2072 for(
long ipLo=0; ipLo<ipHi; ++ipLo )
2081 H2stat /
H2_stat[0][iVibLo][iRotLo] *
2130 else if( nEner == 1 )
2147 double pump_from_below = 0.;
2148 for(
long ipLo = 0; ipLo<nEner; ++ipLo )
2157 if( ( abs(iRotLo-iRot) == 2 || iRotLo == iRot ) && (iVibLo <= iVib) && (*tr).ipCont() > 0 )
2159 double rateone = (*tr).Emis().Aul() * ( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() );
2162 double CosmicRayHILyaExcitationRate = (
hmi.lgLeidenCRHack ) ?
secondaries.x12tot : 0.;
2163 double pump_up = (*tr).Emis().pump() + CosmicRayHILyaExcitationRate * (*tr).Coll().col_str();
2165 rateone += pump_up * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
2185#if defined(__ICC) && defined(__i386)
2186#pragma optimization_level 1
2208 long iElec = (*st).n();
2209 long iVib = (*st).v();
2210 long iRot = (*st).J();
2233 double H2boltzHi =
H2_Boltzmann[iElecHi][iVibHi][iRotHi];
2234 double H2popHi =
states[ipHi].Pop();
2235 double ewnHi =
states[ipHi].energy().WN();
2237 for(
long ipLo=0; ipLo<ipHi; ++ipLo )
2241 double rate_dn_heat = 0.;
2250 double rate_up_cool = rate_dn_heat *
states[ipLo].Pop() *
2252 H2statHi /
H2_stat[iElecLo][iVibLo][iRotLo] *
2255 rate_dn_heat *= H2popHi;
2261 double conversion = (ewnHi -
states[ipLo].energy().WN() ) *
ERG1CM;
2262 double heatone = rate_dn_heat * conversion;
2263 double coolone = rate_up_cool * conversion;
2265 double oneline = heatone - coolone;
2275 (rate_up_cool==0 && rate_dn_heat==0) ||
2276 (
states[ipHi].energy().WN() >
states[ipLo].energy().WN()) );
2283 " DEBUG H2 heat fnzone\t%.2f\trenorm\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n",
2296 " H2_Cooling Ctot\t%.4e\t HeatDiss \t%.4e\t HeatDexc \t%.4e\n" ,
2346 fprintf(
ioQQQ,
" iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n");
2366 if( iRot <0 || iVib >
nVib_hi[iElec] || iRot >
nRot_hi[iElec][iVib])
2368 fprintf(
ioQQQ,
" iVib and iRot must lie within X, returning -2.\n");
2369 fprintf(
ioQQQ,
" iVib must be <= %li and iRot must be <= %li.\n",
2386 if( strcmp(chLabel,
"ZERO") == 0 )
2393 else if( strcmp(chLabel,
"ADD ") == 0 )
2398 long iElec = (*st).n();
2399 if( iElec > 0 )
continue;
2400 long iVib = (*st).v();
2401 long iRot = (*st).J();
2412 else if( strcmp(chLabel,
"PRIN") != 0 )
2414 fprintf(
ioQQQ,
" H2_Colden does not understand the label %s\n",
2444 (*tr).Emis().ots() = (*(*tr).Hi()).Pop() * (*tr).Emis().Aul() * (*tr).Emis().Pdest();
2457 double sumpop1 = 0.;
2458 double sumpopA1 = 0.;
2459 double sumpopcollH2O_deexcit = 0.;
2460 double sumpopcollH2p_deexcit = 0.;
2461 double sumpopcollH_deexcit = 0.;
2463 double sumpopcollH2O_excit = 0.;
2464 double sumpopcollH2p_excit = 0.;
2465 double sumpopcollH_excit = 0.;
2470 long iElecHi = (*stHi).n();
2471 if( iElecHi > 0 )
continue;
2472 long iVibHi = (*stHi).v();
2473 long iRotHi = (*stHi).J();
2474 double ewnHi = (*stHi).energy().WN();
2477 long iVibLo = (*stLo).v();
2478 long iRotLo = (*stLo).J();
2479 double ewnLo2 = (*stLo).energy().WN();
2490 double popHi = (*(*tr).Hi()).Pop();
2491 double popLo = (*(*tr).Lo()).Pop();
2499 sumpopcollH2O_deexcit += popHi *
CollRateCoeff[ihi][ilo][2];
2500 sumpopcollH2p_deexcit += popHi *
CollRateCoeff[ihi][ilo][3];
2502 double temp = popLo *
2514 sumpopA1 += popHi * (*tr).Emis().Aul();
2540 double H2_sum_excit_elec_den = 0.;
2547 return H2_sum_excit_elec_den;
void atom_levelN(long int nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double ***AulEscp, double ***col_str, double ***AulDest, double ***AulPump, double ***CollRate, const double source[], const double sink[], bool lgCollRateDone, double *cooltl, double *coolder, const char *chLabel, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE, multi_arr< double, 2 > *Cool, multi_arr< double, 2 > *dCooldT)
sys_float sexp(sys_float x)
bool fp_equal(sys_float x, sys_float y, int n=3)
NORETURN void TotalInsanity(void)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
TransitionProxy::iterator iterator
double Average_collH_deexcit
multi_arr< realnum, 2 > H2_X_coll_rate
multi_arr< double, 2 > H2_col_rate_out
double Average_collH_dissoc_g
multi_arr< bool, 3 > H2_lgOrtho
double H2_renorm_chemistry
double photodissoc_BigH2_H2s
double Average_collH2_deexcit
void H2_Calc_Average_Rates(void)
double pops_per_elec[N_ELEC]
double rate_grain_op_conserve
vector< double > stat_levn
double Average_collH_excit
void H2_CollidRateEvalAll(void)
void H2_X_coll_rate_evaluate(void)
void SolveExcitedElectronicLevels(void)
double Average_collH_dissoc_s
multi_arr< realnum, 3 > H2_dissprob
multi_arr< realnum, 2 > H2_X_formation
valarray< long > ipElec_H2_energy_sort
void H2_Solomon_rate(void)
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef
realnum GetXColden(long iVib, long iRot)
multi_arr< long int, 2 > ipTransitionSort
double ortho_para_current
valarray< long > ipVib_H2_energy_sort
const double *const dense_total
void SolveSomeGroundElectronicLevels(void)
void CalcPhotoionizationRate(void)
multi_arr< realnum, 3 > H2_stat
double Average_collH2_dissoc_s
double Solomon_dissoc_rate_g
double H2_DissocEnergies[N_ELEC]
multi_arr< double, 2 > H2_col_rate_in
multi_arr< double, 3 > H2_populations_LTE
double photodissoc_BigH2_H2g
valarray< long > ipRot_H2_energy_sort
multi_arr< realnum, 3 > CollRateCoeff
void H2_X_sink_and_source(void)
void H2_Colden(const char *chLabel)
double GetExcitedElecDensity(void)
multi_arr< double, 3 > H2_Boltzmann
const double ENERGY_H2_STAR
multi_arr< double, 2 > H2_X_rate_from_elec_excited
double Average_collH2_dissoc_g
void H2_zero_pops_too_low(void)
multi_arr< int, 2 > H2_ipPhoto
long int nCall_this_iteration
void H2_Level_low_matrix(realnum abundance)
multi_arr< bool, 2 > lgH2_radiative
multi_arr< double, 2 > pops_per_vib
void H2_RT_tau_reset(void)
valarray< realnum > H2_X_source
void H2_LevelPops(bool &lgPopsConverged, double &old_value, double &new_value)
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef_H2
multi_arr< realnum, 2 > H2_X_colden_LTE
multi_arr< realnum, 2 > H2_X_colden
multi_arr< realnum, 2 > H2_X_Hmin_back
double Solomon_dissoc_rate_s
valarray< long > nRot_hi[N_ELEC]
multi_arr< double, 3 > H2_old_populations
double Average_collH2_excit
multi_arr< double, 3 > H2_rad_rate_out
TransitionList::iterator rad_end
multi_arr< realnum, 3 > H2_disske
multi_arr< double, 2 > H2_rad_rate_in
double H2_InterEnergy(void)
multi_arr< double, 2 > H2_X_rate_to_elec_excited
long int nLevels_per_elec[N_ELEC]
long int nzone_nlevel_set
double rate_grain_J1_to_J0
valarray< realnum > H2_X_sink
multi_arr< long int, 3 > ipEnergySort
multi_arr< double, 3 > Cont_Dissoc_Rate
ProxyIterator< qStateConstProxy, qStateConstProxy > const_iterator
ProxyIterator< qStateProxy, qStateConstProxy > iterator
long ipFineCont(double energy_ryd)
long ipoint(double energy_ryd)
long ipLineEnergy(double energy, const char *chLabel, long ipIonEnergy)
multi_arr< double, 2 >::const_iterator md2ci
multi_arr< double, 2 >::iterator md2i
multi_arr< realnum, 3 >::const_iterator mr3ci
void ConvFail(const char chMode[], const char chDetail[])
UNUSED const realnum BIGFLOAT
realnum GetDopplerWidth(realnum massAMU)
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
t_mole_global mole_global
molezone * findspecieslocal(const char buf[])
void mole_update_species_cache(void)
static realnum collider_density_total_not_H2
double cdH2_colden(long iVib, long iRot)
#define H2_DISS_ALLISON_DALGARNO
static realnum collider_density[N_X_COLLIDER]
UNUSED const double EVRYD
UNUSED const double EN1EV
UNUSED const double ERG1CM
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
void RT_line_one(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
void RT_line_one_tau_reset(const TransitionProxy &t)
void RT_line_one_tauinc(const TransitionProxy &t, long int mas_species, long int mas_ion, long int mas_hi, long int mas_lo, realnum DopplerWidth)
void RT_OTS_AddLine(double ots, long int ip)
t_secondaries secondaries
static double ** CollRate