22 #define PRT_POPS false
24 #define LIM_H2_POP_LOOP 100
27 #define H2_DISS_ALLISON_DALGARNO 6e-19f
73 { 36118.11, 118375.6, 118375.6, 118375.6, 118375.6,133608.6,133608.6 };
129 fprintf(
ioQQQ,
" Collider densities are:");
206 ipHi = nLevels_per_elec[0];
207 while( (--ipHi) > 0 )
267 for( ipLo=0; ipLo<ipHi; ++ipLo )
281 ASSERT( H2CollRate[nColl]*collider_density[nColl] >= 0. );
307 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
317 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
325 long int lim_elec_lo = 0;
326 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
331 if( iElecLo==iElecHi )
333 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
336 if( iElecLo==iElecHi && iVibHi==iVibLo )
339 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
344 ASSERT(
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul > 0. );
345 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont =
347 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN *
WAVNRYD ,
349 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ipFine =
351 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN * WAVNRYD );
366 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
382 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
390 long int lim_elec_lo = 0;
391 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
396 if( iElecLo==iElecHi )
398 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
401 if( iElecLo==iElecHi && iVibHi==iVibLo )
406 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
412 ASSERT( H2L->ipCont > 0 );
413 h2_drive += H2L->Emis->pump*H2L->EnergyErg*H2L->Emis->PopOpc;
429 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
446 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
454 long int lim_elec_lo = 0;
455 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
460 if( iElecLo==iElecHi )
462 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
465 if( iElecLo==iElecHi && iVibHi==iVibLo )
470 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
474 ASSERT( H2L->ipCont > 0 );
476 if( H2L->Hi->Pop > smallfloat && H2L->Emis->PopOpc > smallfloat )
492 " H2_RadPress returns, radiation pressure is %.2e\n",
502 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
520 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
524 double H2popHi =
H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop;
530 long int lim_elec_lo = 0;
531 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
536 if( iElecLo==iElecHi )
538 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
541 if( iElecLo==iElecHi && iVibHi==iVibLo )
546 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
552 energy += H2popHi*H2L->EnergyErg;
568 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
577 for( iElecHi=0; iElecHi<1; ++iElecHi )
579 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
587 long int lim_elec_lo = 0;
588 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
593 if( iElecLo==iElecHi )
595 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
598 if( iElecLo==iElecHi && iVibHi==iVibLo )
601 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
607 outline( &
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
621 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
632 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
640 long int lim_elec_lo = 0;
641 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
646 if( iElecLo==iElecHi )
648 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
651 if( iElecLo==iElecHi && iVibHi==iVibLo )
655 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
665 RT_line_one( &
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo], lgDoEsc, lgUpdateFineOpac,
false );
676 enum {DEBUG_LOC=
false};
683 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
694 for( iVibLo=0; iVibLo<=
h2.
nVib_hi[iElecLo]; ++iVibLo )
697 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
704 ASSERT(
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul > 0. );
706 sumpop +=
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop;
707 sumpopA +=
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->Pop*
708 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul;
715 fprintf(
ioQQQ,
"DEBUG sumpop = %.3e sumpopA= %.3e A=%.3e\n", sumpop, sumpopA,
716 sumpopA/
SDIV(sumpop) );
726 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
747 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
755 long int lim_elec_lo = 0;
756 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
761 if( iElecLo==iElecHi )
763 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
766 if( iElecLo==iElecHi && iVibHi==iVibLo )
770 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
776 ASSERT(
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont > 0 );
778 -9 , iRotHi , iVibLo , iRotLo);
794 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
804 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
812 long int lim_elec_lo = 0;
813 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
818 if( iElecLo==iElecHi )
820 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
823 if( iElecLo==iElecHi && iVibHi==iVibLo )
826 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
849 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
859 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
867 long int lim_elec_lo = 0;
868 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
873 if( iElecLo==iElecHi )
875 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
878 if( iElecLo==iElecHi && iVibHi==iVibLo )
881 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
910 static double **AulEscp ,
925 static long int levelAsEval=-1;
927 static bool lgFirst=
true;
941 double rot_cooling , dCoolDT;
942 static long int ndimMalloced = 0;
943 double rateout , ratein;
960 " The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n",
971 excit = (
double *)
MALLOC(
sizeof(
double)*(size_t)(ndimMalloced) );
972 stat_levn = (
double *)
MALLOC(
sizeof(
double)*(size_t)(ndimMalloced) );
973 pops = (
double *)
MALLOC(
sizeof(
double)*(size_t)(ndimMalloced) );
974 create = (
double *)
MALLOC(
sizeof(
double)*(size_t)(ndimMalloced) );
975 destroy = (
double *)
MALLOC(
sizeof(
double)*(size_t)(ndimMalloced) );
976 depart = (
double *)
MALLOC(
sizeof(
double)*(size_t)(ndimMalloced) );
978 AulPump = ((
double **)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double *)));
979 CollRate_levn = ((
double **)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double *)));
980 AulDest = ((
double **)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double *)));
981 AulEscp = ((
double **)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double *)));
982 col_str = ((
double **)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double *)));
983 for( i=0; i<(ndimMalloced); ++i )
985 AulPump[i] = ((
double *)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double )));
986 CollRate_levn[i] = ((
double *)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double )));
987 AulDest[i] = ((
double *)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double )));
988 AulEscp[i] = ((
double *)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double )));
989 col_str[i] = ((
double *)
MALLOC((
size_t)(ndimMalloced)*
sizeof(
double )));
992 for( j=0; j < ndimMalloced; j++ )
999 for( j=0; j < ndimMalloced; j++ )
1007 stat_levn[j] =
H2_stat[0][iVib][iRot];
1012 for( j=0; j < ndimMalloced-1; j++ )
1015 ASSERT( excit[j+1] > excit[j] );
1024 fprintf(
ioQQQ,
" H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n");
1046 ASSERT( levelAsEval <= ndimMalloced );
1065 CollRate_levn[j][i] = 0.;
1099 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
1104 ASSERT(
H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].ipCont > 0 );
1111 AulEscp[ihi][ilo] =
H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Aul*(
1112 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pesc +
1113 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pdest +
1114 H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->Pelec_esc);
1115 AulDest[ilo][ihi] = 0.;
1116 AulPump[ilo][ihi] =
H2Lines[0][iVibHi][iRotHi][0][iVib][iRot].Emis->pump;
1133 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
1137 ASSERT(
H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].ipCont > 0 );
1144 (
H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Aul*
1145 (
H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pesc +
1146 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pelec_esc +
1147 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pdest)+
H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->pump *
1148 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Lo->g/
1149 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Hi->g);
1152 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->pump;
1158 create[ilo] += ratein;
1161 destroy[ilo] += rateout;
1185 if(lgDeBug)fprintf(
ioQQQ,
"DEBUG H2_Level_low_matrix, ilo=%li",ilo);
1195 if(lgDeBug)fprintf(
ioQQQ,
"\t%.1e",CollRate_levn[ihi][ilo]);
1198 CollRate_levn[ilo][ihi] = CollRate_levn[ihi][ilo]*
1205 if(lgDeBug)fprintf(
ioQQQ,
"\n");
1221 if(lgDeBug)fprintf(
ioQQQ,
"\t%.1e",ratein);
1230 destroy[ilo] += rateout;
1255 CollRate_levn[1][0] +=
1273 enum {DEBUG_LOC=
false};
1274 if( DEBUG_LOC || lgDeBug)
1276 fprintf(
ioQQQ,
"DEBUG H2 matexcit");
1279 fprintf(
ioQQQ,
"\t%li",ilo );
1281 fprintf(
ioQQQ,
"\n");
1284 fprintf(
ioQQQ,
"\t%.2e",excit[ihi] );
1286 fprintf(
ioQQQ,
"\n");
1289 fprintf(
ioQQQ,
"\t%.2e",stat_levn[ihi] );
1291 fprintf(
ioQQQ,
"\n");
1293 fprintf(
ioQQQ,
"AulEscp[n][]\\[][n] = Aul*Pesc\n");
1296 fprintf(
ioQQQ,
"\t%li",ilo );
1298 fprintf(
ioQQQ,
"\n");
1301 fprintf(
ioQQQ,
"%li", ihi);
1304 fprintf(
ioQQQ,
"\t%.2e",AulEscp[ilo][ihi] );
1306 fprintf(
ioQQQ,
"\n");
1309 fprintf(
ioQQQ,
"AulPump [n][]\\[][n]\n");
1312 fprintf(
ioQQQ,
"\t%li",ilo );
1314 fprintf(
ioQQQ,
"\n");
1317 fprintf(
ioQQQ,
"%li", ihi);
1320 fprintf(
ioQQQ,
"\t%.2e",AulPump[ihi][ilo] );
1322 fprintf(
ioQQQ,
"\n");
1325 fprintf(
ioQQQ,
"CollRate_levn [n][]\\[][n]\n");
1328 fprintf(
ioQQQ,
"\t%li",ilo );
1330 fprintf(
ioQQQ,
"\n");
1333 fprintf(
ioQQQ,
"%li", ihi);
1336 fprintf(
ioQQQ,
"\t%.2e",CollRate_levn[ihi][ilo] );
1338 fprintf(
ioQQQ,
"\n");
1340 fprintf(
ioQQQ,
"SOURCE");
1343 fprintf(
ioQQQ,
"\t%.2e",create[ihi]);
1345 fprintf(
ioQQQ,
"\nSINK");
1348 fprintf(
ioQQQ,
"\t%.2e",destroy[ihi]);
1350 fprintf(
ioQQQ,
"\n");
1398 fprintf(
ioQQQ,
"\n DEBUG H2_Level_lowJ hmi.H2_total: %.3e matrix rel pops\n",
hmi.
H2_total);
1399 fprintf(
ioQQQ,
"v\tJ\tpop\n");
1405 fprintf(
ioQQQ,
"%3li\t%3li\t%.3e\t%.3e\t%.3e\n",
1413 fprintf(
ioQQQ,
" H2_Level_low_matrix called atom_levelN which returned negative H2_populations.\n");
1422 static double TeUsedColl=-1.f;
1423 double H2_renorm_conserve=0.,
1424 H2_renorm_conserve_init=0. ,
1428 double old_solomon_rate=-1.;
1429 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
1431 long int n_pop_oscil = 0;
1433 bool lgConv_h2_soln,
1435 lgPopsConv_relative,
1438 lgOrthoParaRatioConv;
1439 double quant_old=-1.,
1442 bool lgH2_pops_oscil=
false,
1443 lgH2_pops_ever_oscil=
false;
1446 long int iElec , iVib , iRot,ip;
1447 double sum_pops_matrix;
1450 static double ortho_para_old=0. , ortho_para_older=0. , ortho_para_current=0.;
1454 double PopChgMax_relative=0. , PopChgMaxOld_relative=0., PopChgMax_total=0., PopChgMaxOld_total=0.;
1455 long int iRotMaxChng_relative , iVibMaxChng_relative,
1456 iRotMaxChng_total , iVibMaxChng_total,
1457 nXLevelsMatrix_save;
1458 double popold_relative , popnew_relative , popold_total , popnew_total;
1462 double flux_accum_photodissoc_BigH2_H2g, flux_accum_photodissoc_BigH2_H2s;
1463 long int ip_H2_level;
1466 double converge_pops_relative=0.1 ,
1467 converge_pops_total=1e-3,
1468 converge_ortho_para=1e-2;
1480 "\n***************H2_LevelPops call %li this iteration, zone is %.2f, H2/H:%.e Te:%e ne:%e\n",
1490 static long int nzone_prt=-1;
1491 if(
nzone!=nzone_prt )
1494 fprintf(
ioQQQ,
"DEBUG zone %li H2/H:%.3e Te:%.3e *ne:%.3e n(H2):%.3e\n",
1516 " H2_LevelPops pops too small, not computing, set to LTE and return, H2/H is %.2e and mole.H2_to_H_limit is %.2e.",
1549 if( fabs(1. - TeUsedColl /
phycon.
te ) > 0.005 )
1562 fprintf(
ioQQQ,
"H2 1st call - using LTE level pops\n");
1567 for( iVib=0; iVib<=
h2.
nVib_hi[iElec]; ++iVib )
1584 ortho_para_older = ortho_para_current;
1585 ortho_para_old = ortho_para_current;
1593 for( iVib=0; iVib<=
h2.
nVib_hi[iElec]; ++iVib )
1614 for( iVib=0; iVib<=
h2.
nVib_hi[iElec]; ++iVib )
1637 "H2 H2_renorm_chemistry is %.4e, hmi.H2_total is %.4e pops_per_elec[0] is %.4e\n",
1644 for( iVib=0; iVib<=
h2.
nVib_hi[iElec]; ++iVib )
1656 " H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n",
1663 converge_pops_relative *= 2.;
1664 converge_pops_total *= 3.;
1665 converge_ortho_para *= 3.;
1676 lgConv_h2_soln =
false;
1680 static long int nzoneEval=-1;
1681 if(
nzone != nzoneEval )
1697 double rate_in , rate_out;
1698 static double old_HeatH2Dexc_BigH2=0., HeatChangeOld=0. , HeatChange=0.;
1714 for( iVib=0; iVib<=
h2.
nVib_hi[iElec]; ++iVib )
1742 fprintf(
ioQQQ,
" Pop(e=%li):",iElecHi);
1747 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
1767 for( iVibLo=0; iVibLo<=
h2.
nVib_hi[iElecLo]; ++iVibLo )
1772 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1777 ASSERT( H2L->ipCont > 0 );
1783 double cr_cross_section = H2L->Coll.cs;
1787 pow(
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng*1e-8,3)*
1788 (
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g/
1789 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g)*
1790 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Aul*
1791 log(3.)*
HPLANCK/(160.f*pow(
PI,3)*0.5*1e-8*1.6e-12);
1792 ASSERT( fabs(cr_cross_section-
H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cs)/
1793 SDIV(cr_cross_section) < 1e-7 );
1799 double rate_one = H2L->Emis->pump
1807 +factor*cr_cross_section;
1827 rate_one = H2L->Emis->Aul*
1830 H2L->Emis->Pelec_esc +
1838 rate_out += rate_one;
1840 ASSERT( rate_in >= 0. && rate_out >= 0. );
1851 double H2popHi = rate_in /
SDIV( rate_out );
1859 for( iVibLo=0; iVibLo<=
h2.
nVib_hi[iElecLo]; ++iVibLo )
1865 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
1869 ASSERT( H2L->ipCont > 0 );
1877 (H2L->Emis->Pesc + H2L->Emis->Pelec_esc + H2L->Emis->Pdest) +
1879 H2L->Emis->pump * H2L->Lo->g/H2gHi);
1905 fprintf(
ioQQQ,
"\n");
1913 PopChgMaxOld_relative = PopChgMax_relative;
1914 PopChgMaxOld_total = PopChgMax_total;
1915 PopChgMax_relative = 0.;
1916 PopChgMax_total = 0.;
1919 iRotMaxChng_relative =-1;
1920 iVibMaxChng_relative = -1;
1921 iRotMaxChng_total =-1;
1922 iVibMaxChng_total = -1;
1923 popold_relative = 0.;
1924 popnew_relative = 0.;
1940 double col_rate_in = 0.;
1941 double col_rate_out = 0.;
1942 for( ipLo=0; ipLo<nEner; ++ipLo )
1954 H2stat /
H2_stat[0][iVibLo][iRotLo] *
1957 col_rate_in += collup;
1960 for( ipHi=nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi )
1969 H2_stat[0][iVibHi][iRotHi] / H2stat *
1979 col_rate_in += colldn;
1993 nEner = nLevels_per_elec[0];
2003 if( nEner+1 < nLevels_per_elec[0] )
2024 else if( nEner == 1 )
2056 for( ipHi = nEner+1; ipHi<nLevels_per_elec[0]; ++ipHi )
2066 if( ( abs(iRotHi-iRot) ==2 || iRotHi==iRot ) && (iVib <= iVibHi) &&
2067 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].ipCont > 0 )
2072 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Aul*
2073 (
H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pesc +
2074 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pelec_esc +
2075 H2Lines[iElecHi][iVibHi][iRotHi][iElec][iVib][iRot].Emis->Pdest);
2087 for( ipLo = 0; ipLo<nEner; ++ipLo )
2095 if( ( abs(iRotLo-iRot) == 2 || iRotLo == iRot ) && (iVibLo <= iVib) &&
2096 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].ipCont > 0 )
2099 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Aul*
2100 (
H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pesc +
2101 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pelec_esc +
2102 H2Lines[iElec][iVib][iRot][iElecLo][iVibLo][iRotLo].Emis->Pdest);
2129 fprintf(
ioQQQ,
" Rel pop(e=%li)" ,iElecHi);
2135 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
2155 fprintf(
ioQQQ,
"\n");
2157 fprintf(
ioQQQ,
" Rel pop(0,J)");
2164 fprintf(
ioQQQ,
"\n");
2170 sum_pops_matrix = 0.;
2202 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
2219 for( iVib=0; iVib<=
h2.
nVib_hi[iElec]; ++iVib )
2240 PopChgMax_relative =
2244 iRotMaxChng_relative = iRot;
2245 iVibMaxChng_relative = iVib;
2256 if( fabs(rel_change) > fabs(PopChgMax_total) )
2258 PopChgMax_total = rel_change;
2259 iRotMaxChng_total = iRot;
2260 iVibMaxChng_total = iVib;
2277 if( rel_change > 3. &&
2288 else if( rel_change> 0.1 )
2316 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
2332 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
2357 ortho_para_older = ortho_para_old;
2358 ortho_para_old = ortho_para_current;
2372 HeatChangeOld = HeatChange;
2375 static long int loop_h2_oscil=-1;
2379 (HeatChangeOld*HeatChange<0. ) ||
2380 (PopChgMax_relative*PopChgMaxOld_relative<0. ) ) )
2382 lgH2_pops_oscil =
true;
2386 lgH2_pops_ever_oscil =
true;
2389 frac_new_oscil *= 0.8f;
2390 frac_new_oscil =
MAX2( frac_new_oscil , 0.1f);
2396 lgH2_pops_oscil =
false;
2400 frac_new_oscil = 1.f;
2401 lgH2_pops_ever_oscil =
false;
2414 lgConv_h2_soln =
true;
2415 lgPopsConv_total =
true;
2416 lgPopsConv_relative =
true;
2418 lgSolomonConv =
true;
2419 lgOrthoParaRatioConv =
true;
2423 if( fabs(PopChgMax_relative)>converge_pops_relative )
2426 lgConv_h2_soln =
false;
2427 lgPopsConv_relative =
false;
2430 quant_old = PopChgMaxOld_relative;
2432 quant_new = PopChgMax_relative;
2434 strcpy( chReason ,
"rel pops changed" );
2439 else if( fabs(PopChgMax_total)>converge_pops_total)
2441 lgConv_h2_soln =
false;
2442 lgPopsConv_total =
false;
2445 quant_old = PopChgMaxOld_total;
2447 quant_new = PopChgMax_total;
2449 strcpy( chReason ,
"tot pops changed" );
2456 else if( fabs(ortho_para_current-ortho_para_old) /
SDIV(ortho_para_current)> converge_ortho_para )
2460 lgConv_h2_soln =
false;
2461 lgOrthoParaRatioConv =
false;
2462 quant_old = ortho_para_old;
2463 quant_new = ortho_para_current;
2464 strcpy( chReason ,
"ortho/para ratio changed" );
2477 lgConv_h2_soln =
false;
2481 strcpy( chReason ,
"heating changed" );
2500 lgConv_h2_soln =
false;
2501 lgSolomonConv =
false;
2502 quant_old = old_solomon_rate;
2504 strcpy( chReason ,
"Solomon rate changed" );
2508 if( !lgConv_h2_soln )
2519 fprintf(
ioQQQ,
" loop %3li no conv oscl?%c why:%s ",
2521 TorF(lgH2_pops_ever_oscil),
2523 if( !lgPopsConv_relative )
2524 fprintf(
ioQQQ,
" PopChgMax_relative:%.4e v:%li J:%li old:%.4e new:%.4e",
2526 iVibMaxChng_relative,
2527 iRotMaxChng_relative ,
2530 else if( !lgPopsConv_total )
2531 fprintf(
ioQQQ,
" PopChgMax_total:%.4e v:%li J:%li old:%.4e new:%.4e",
2537 else if( !lgHeatConv )
2538 fprintf(
ioQQQ,
" heat:%.4e old:%.4e new:%.4e",
2543 else if( !lgSolomonConv )
2545 else if( !lgOrthoParaRatioConv )
2546 fprintf(
ioQQQ,
" current, old, older ratios are %.4e %.4e %.4e",
2547 ortho_para_current , ortho_para_old, ortho_para_older );
2550 fprintf(
ioQQQ,
"\n");
2569 " H2 5lev %li Conv?%c",
2571 TorF(lgConv_h2_soln) );
2573 if( fabs(PopChgMax_relative)>0.1 )
2574 fprintf(
ioQQQ,
" pops, rel chng %.3e",PopChgMax_relative);
2576 fprintf(
ioQQQ,
" rel heat %.3e rel chng %.3e H2 heat/cool %.2e",
2582 " Oscil?%c Ever Oscil?%c",
2583 TorF(lgH2_pops_oscil) ,
2584 TorF(lgH2_pops_ever_oscil) );
2585 if( lgH2_pops_ever_oscil )
2586 fprintf(
ioQQQ,
" frac_new_oscil %.4f",frac_new_oscil);
2587 fprintf(
ioQQQ,
"\n");
2593 "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",
2601 if( iVibMaxChng_relative>=0 && iRotMaxChng_relative>=0 )
2603 "end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n",
2605 PopChgMax_relative ,
2608 iVibMaxChng_relative , iRotMaxChng_relative
2623 " H2_LevelPops: H2_populations not converged in %li tries; due to %s, old, new are %.4e %.4e, iteration %li zone %.2f.\n",
2637 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
2641 long int lim_elec_lo = 0;
2645 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop = H2popHi;
2652 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
2657 if( iElecLo==iElecHi )
2659 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
2662 if( iElecLo==iElecHi && iVibHi==iVibLo )
2667 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
2672 H2L->Coll.cool = 0.;
2673 H2L->Coll.heat = 0.;
2676 H2L->Lo->Pop - H2popHi * H2L->Lo->g / H2gHi;
2679 H2L->Emis->phots = H2L->Emis->Aul *
2680 (H2L->Emis->Pesc + H2L->Emis->Pelec_esc) *
2684 H2L->Emis->xIntensity =
2693 H2L->Emis->ColOvTot = 1.;
2699 H2L->Emis->ColOvTot = 0.;
2739 static long ip_cut_off = -1;
2740 if( ip_cut_off < 0 )
2743 ip_cut_off =
ipoint( 1.14 );
2748 flux_accum_photodissoc_BigH2_H2s = 0;
2750 for( i= ip_H2_level; i < ip_cut_off; ++i )
2757 for( iVib=0; iVib<=
h2.
nVib_hi[iElec]; ++iVib )
2771 H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2s;
2792 if( arg_ratio > 0. )
2796 (
H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5;
2803 flux_accum_photodissoc_BigH2_H2g = 0;
2807 for( i= ip_H2_level; i < ip_cut_off; ++i )
2814 H2_populations[iElec][iVib][iRot] * flux_accum_photodissoc_BigH2_H2g;
2836 if( arg_ratio > 0. )
2839 (
H2_stat[0][iVib][iRot]/(2.*2.))*3.634e-5);
2871 double sumpop1 = 0.;
2872 double sumpopA1 = 0.;
2873 double sumpopcollH2O_deexcit = 0.;
2874 double sumpopcollH2p_deexcit = 0.;
2875 double sumpopcollH_deexcit = 0.;
2877 double sumpopcollH2O_excit = 0.;
2878 double sumpopcollH2p_excit = 0.;
2879 double sumpopcollH_excit = 0.;
2884 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[0]; ++iVibHi )
2887 for( iRotHi=
h2.
Jlowest[0]; iRotHi<=nr1; ++iRotHi )
2889 double ewnHi =
energy_wn[0][iVibHi][iRotHi];
2890 for( iVibLo=0; iVibLo<=
h2.
nVib_hi[0]; ++iVibLo )
2894 for( iRotLo=
h2.
Jlowest[0]; iRotLo<=nr2; ++iRotLo )
2896 double ewnLo2 =
energy_wn[0][iVibLo][iRotLo];
2910 sumpopcollH2O_deexcit +=
H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][2];
2911 sumpopcollH2p_deexcit +=
H2_populations[0][iVibHi][iRotHi]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][3];
2914 sumpopcollH_excit +=
H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0]*
H2_stat[0][iVibHi][iRotHi] /
H2_stat[0][iVibLo][iRotLo] *
2916 sumpopcollH2O_excit +=
H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][2]*
H2_stat[0][iVibHi][iRotHi] /
H2_stat[0][iVibLo][iRotLo] *
2918 sumpopcollH2p_excit +=
H2_populations[0][iVibLo][iRotLo]*H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][3]*
H2_stat[0][iVibHi][iRotHi] /
H2_stat[0][iVibLo][iRotLo] *
2926 sumpop1 +=
H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Hi->Pop;
2927 sumpopA1 +=
H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Hi->Pop*
2928 H2Lines[0][iVibHi][iRotHi][0][iVibLo][iRotLo].Emis->Aul;
2959 fprintf(
ioQQQ,
" H2_LevelPops exit2 Sol dissoc %.2e (TH85 %.2e)",
2966 fprintf(
ioQQQ,
" H2g Sol %.2e H2s Sol %.2e",
2971 fprintf(
ioQQQ,
" H2g->H2s %.2e (TH85 %.2e)",
2977 fprintf(
ioQQQ,
" H2 con diss %.2e (TH85 %.2e)\n",
2983 fprintf(
ioQQQ,
" H2_LevelPops exit1 %8.2f loops:%3li H2/H:%.3e Sol dis old %.3e new %.3e",
2989 fprintf(
ioQQQ,
"\n");
3014 # define FRAC 0.99999
3016 double sum_pop = 0.;
3020 if(
PRT ) fprintf(
ioQQQ,
"DEBUG pops ");
3048 #if defined(__ICC) && defined(__i386)
3049 #pragma optimization_level 1
3055 const char *chRoutine)
3057 long int iElecHi , iElecLo , iVibHi , iVibLo , iRotHi , iRotLo;
3063 double Big1_heat , Big1_cool,
3064 H2_X_col_cool , H2_X_col_heat;
3065 long int ipVib_big_heat_hi,ipVib_big_heat_lo ,ipRot_big_heat_hi ,
3066 ipRot_big_heat_lo ,ipVib_big_cool_hi,ipVib_big_cool_lo ,
3067 ipRot_big_cool_hi , ipRot_big_cool_lo;
3069 # ifdef DEBUG_DIS_DEAT
3071 long int iElecBig , iVibBig , iRotBig;
3076 enum {DEBUG_H2_COLL_X_HEAT=
false };
3082 static long int nCount=-1;
3083 long int nCountDebug = 930,
3101 # ifdef DEBUG_DIS_DEAT
3110 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
3119 # ifdef DEBUG_DIS_DEAT
3120 if( heatone > heatbig )
3131 # ifdef DEBUG_DIS_DEAT
3132 fprintf(
ioQQQ,
"DEBUG H2 dis heat\t%.2f\t%.3f\t%li\t\t%li\t%li\n",
3153 long int nBug1 = 200 , nBug2 = 201;
3158 nBug1 = (long)
fudge(0);
3159 nBug2 = (long)
fudge(1);
3162 if( DEBUG_H2_COLL_X_HEAT && (nCount == nBug1 || nCount==nBug2) )
3167 ioBAD = fopen(
"firstpop.txt" ,
"w" );
3169 else if( nCount==nBug2 )
3171 ioBAD = fopen(
"secondpop.txt" ,
"w" );
3178 fprintf(ioBAD ,
"%li\t%li\t%.2e\n",
3190 ipVib_big_heat_hi = -1;
3191 ipVib_big_heat_lo = -1;
3192 ipRot_big_heat_hi = -1;
3193 ipRot_big_heat_lo = -1;
3194 ipVib_big_cool_hi = -1;
3195 ipVib_big_cool_lo = -1;
3196 ipRot_big_cool_hi = -1;
3197 ipRot_big_cool_lo = -1;
3209 double H2boltzHi =
H2_Boltzmann[iElecHi][iVibHi][iRotHi];
3211 double ewnHi =
energy_wn[iElecHi][iVibHi][iRotHi];
3213 for( ipLo=0; ipLo<ipHi; ++ipLo )
3215 double coolone , oneline;
3231 rate_up_cool = rate_dn_heat *
H2_populations[iElecLo][iVibLo][iRotLo] *
3233 H2statHi /
H2_stat[iElecLo][iVibLo][iRotLo] *
3236 rate_dn_heat *= H2popHi;
3242 double conversion = (ewnHi -
energy_wn[iElecLo][iVibLo][iRotLo]) *
ERG1CM;
3243 heatone = rate_dn_heat * conversion;
3244 coolone = rate_up_cool * conversion;
3246 oneline = heatone - coolone;
3250 H2_X_col_cool += coolone;
3251 H2_X_col_heat += heatone;
3253 if( 0 && DEBUG_H2_COLL_X_HEAT && (nCount == 692 || nCount==693) )
3255 static FILE *ioBAD=NULL;
3256 if(ipHi == 1 && ipLo == 0)
3260 ioBAD = fopen(
"firstheat.txt" ,
"w" );
3262 else if( nCount==693 )
3264 ioBAD = fopen(
"secondheat.txt" ,
"w" );
3266 fprintf(ioBAD,
"DEBUG start \n");
3268 fprintf(ioBAD,
"DEBUG BAD DAY %li %li %li %li %.3e %.3e\n",
3269 iVibHi , iRotHi, iVibLo , iRotLo , heatone , coolone );
3270 if( ipHi==nLevels_per_elec[iElecHi]-1 &&
3282 if( DEBUG_H2_COLL_X_HEAT )
3284 if( oneline < Big1_cool )
3286 Big1_cool = oneline;
3287 ipVib_big_cool_hi = iVibHi;
3288 ipVib_big_cool_lo = iVibLo;
3289 ipRot_big_cool_hi = iRotHi;
3290 ipRot_big_cool_lo = iRotLo;
3292 else if( oneline > Big1_heat )
3294 Big1_heat = oneline;
3295 ipVib_big_heat_hi = iVibHi;
3296 ipVib_big_heat_lo = iVibLo;
3297 ipRot_big_heat_hi = iRotHi;
3298 ipRot_big_heat_lo = iRotLo;
3304 (rate_up_cool==0 && rate_dn_heat==0) ||
3311 if( DEBUG_H2_COLL_X_HEAT )
3316 if(nCount>nCountDebug )
3319 "DEBUG H2_Cooling A %li %15s, Te %.3e net Heat(Xcol) %.2e "
3320 "heat %.2e cool %.2e H+/0 "
3321 "%.2e n(H2)%.3e Sol rat %.3e grn J1->0%.2e frac heat 1 line "
3322 "%.2e Hi(v,j)%li %li Lo(v,J)%li %li frac cool 1 line %.2e "
3323 "Hi(v,j)%li %li Lo(v,J)%li %li POP(J=1,13)",
3334 ipVib_big_heat_hi , ipRot_big_heat_hi , ipVib_big_heat_lo , ipRot_big_heat_lo ,
3336 ipVib_big_cool_hi , ipRot_big_cool_hi , ipVib_big_cool_lo , ipRot_big_cool_lo );
3338 for( iRotLo=0; iRotLo<14; ++iRotLo )
3340 fprintf(
ioQQQ,
"\t%.2e" ,
3343 fprintf(
ioQQQ,
"\t%li\n",nCount);
3345 fprintf(
ioQQQ,
"DEBUG H2_Cooling B heat Coll Rate (lrg col) dn,up" );
3346 double HeatNet = 0.;
3349 fprintf(
ioQQQ,
"\t%.2e" ,
3350 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3352 fprintf(
ioQQQ,
"\t%.2e" ,
3354 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3357 H2_stat[0][ipVib_big_heat_hi][ipRot_big_heat_hi] /
H2_stat[0][ipVib_big_heat_lo][ipRot_big_heat_lo] *
3358 H2_Boltzmann[0][ipVib_big_heat_hi][ipRot_big_heat_hi] /
3361 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3364 H2_CollRate[ipVib_big_heat_hi][ipRot_big_heat_hi][ipVib_big_heat_lo][ipRot_big_heat_lo][nColl]*
3367 H2_stat[0][ipVib_big_heat_hi][ipRot_big_heat_hi] /
H2_stat[0][ipVib_big_heat_lo][ipRot_big_heat_lo] *
3368 H2_Boltzmann[0][ipVib_big_heat_hi][ipRot_big_heat_hi] /
3373 fprintf(
ioQQQ ,
" HeatNet %.2e",HeatNet);
3374 fprintf(
ioQQQ,
"\n");
3375 fprintf(
ioQQQ,
"DEBUG H2_Cooling C cool Coll Rate (lrg col) dn,up" );
3376 double CoolNet = 0.;
3379 fprintf(
ioQQQ,
"\t%.2e" ,
3380 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3382 fprintf(
ioQQQ,
"\t%.2e" ,
3384 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3387 H2_stat[0][ipVib_big_cool_hi][ipRot_big_cool_hi] /
H2_stat[0][ipVib_big_cool_lo][ipRot_big_cool_lo] *
3388 H2_Boltzmann[0][ipVib_big_cool_hi][ipRot_big_cool_hi] /
3391 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3394 H2_CollRate[ipVib_big_cool_hi][ipRot_big_cool_hi][ipVib_big_cool_lo][ipRot_big_cool_lo][nColl]*
3397 H2_stat[0][ipVib_big_cool_hi][ipRot_big_cool_hi] /
H2_stat[0][ipVib_big_cool_lo][ipRot_big_cool_lo] *
3398 H2_Boltzmann[0][ipVib_big_cool_hi][ipRot_big_cool_hi] /
3403 fprintf(
ioQQQ ,
" CoolNet %.2e",CoolNet);
3404 fprintf(
ioQQQ,
"\n");
3411 " DEBUG H2 heat fnzone\t%.2f\trenrom\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n",
3423 enum {DEBUG_LOC=
false };
3424 if( DEBUG_H2_COLL_X_HEAT && DEBUG_LOC &&
3437 iElecHi = iElecLo = 0;
3438 iVibHi = iVibLo = 0;
3441 rate_dn_heat = rate_up_cool = 0.;
3447 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl]*
3454 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl]*
3455 collider_density[nColl]*
3457 H2_stat[iElecHi][iVibHi][iRotHi] /
H2_stat[iElecLo][iVibLo][iRotLo] *
3462 fprintf(
ioQQQ,
"DEBUG H2_cooling D pop %li ov %li\t%.3e\tdn up 31\t%.3e\t%.3e\n",
3467 if( nCount>= nCountStop )
3472 enum {DEBUG_LOC=
false};
3475 static long nzdone=-1 , nzincre;
3482 fprintf(
ioQQQ,
" H2 nz\t%.2f\tnzinc\t%li\tTe\t%.4e\tH2\t%.3e\tcXH\t%.2e\tdcXH/dt%.2e\tDish\t%.2e \n",
3509 enum {DEBUG_LOC=
false};
3512 fprintf(
ioQQQ,
"DEBUG H2_cooling E %15s %c vib deex %li Te %.3e net heat %.3e cool %.3e heat %.3e\n",
3521 if( 0 && nCount > nCountStop )
3529 " H2_Cooling Ctot\t%.4e\t HeatH2Dish_BigH2 \t%.4e\t HeatH2Dexc_BigH2 \t%.4e\n" ,
3576 fprintf(
ioQQQ,
" iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n");
3587 fprintf(
ioQQQ,
" iVib and iRot must lie within X, returning -2.\n");
3588 fprintf(
ioQQQ,
" iVib must be <= %li and iRot must be <= %li.\n",
3605 long int iVib , iRot;
3613 if( strcmp(chLabel,
"ZERO") == 0 )
3616 for( iVib = 0; iVib <=
h2.
nVib_hi[0]; ++iVib )
3627 else if( strcmp(chLabel,
"ADD ") == 0 )
3630 for( iVib = 0; iVib <=
h2.
nVib_hi[0]; ++iVib )
3645 else if( strcmp(chLabel,
"PRIN") != 0 )
3647 fprintf(
ioQQQ,
" H2_Colden does not understand the label %s\n",
3665 long int iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo;
3674 long int lim_elec_hi = 0;
3675 for( iElecHi=0; iElecHi<=lim_elec_hi; ++iElecHi )
3677 for( iVibHi=0; iVibHi<=
h2.
nVib_hi[iElecHi]; ++iVibHi )
3681 double H2popHi =
H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi->Pop;
3682 long int lim_elec_lo = 0;
3683 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
3688 if( iElecLo==iElecHi )
3690 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
3693 if( iElecLo==iElecHi && iVibHi==iVibLo )
3696 for( iRotLo=
h2.
Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
3699 if( iElecHi==0 &&
lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
3702 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots =
3703 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul * H2popHi *
3704 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Pdest;
3708 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->ots,
3709 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].ipCont );