cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
mole_h2.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*H2_ContPoint set the ipCont struc element for the H2 molecule, called by ContCreatePointers */
4 /*H2_Accel radiative acceleration due to H2 */
5 /*H2_RadPress rad pressure due to h2 lines called in PresTotCurrent */
6 /*H2_InterEnergy internal energy of H2 called in PresTotCurrent */
7 /*H2_RT_diffuse do emission from H2 - called from RT_diffuse */
8 /*H2_itrzn - average number of H2 pop evaluations per zone */
9 /*H2_RTMake do RT for H2 - called from RT_line_all */
10 /*H2_RT_tau_inc increment optical depth for the H2 molecule, called from RT_tau_inc */
11 /*H2_LineZero initialize optical depths in H2, called from RT_tau_init */
12 /*H2_RT_tau_reset the large H2 molecule, called from RT_tau_reset */
13 /*H2_Colden maintain H2 column densities within X */
14 /*H2_LevelPops do level populations for H2, called by Hydrogenic */
15 /*H2_Level_low_matrix evaluate CO rotation cooling */
16 /*H2_cooling evaluate cooling and heating due to H2 molecule */
17 /*H2_X_coll_rate_evaluate find collisional rates within X */
18 /*cdH2_colden return column density in H2, negative -1 if cannot find state,
19  * header is cddrive */
20 /*H2_DR choose next zone thickness based on H2 big molecule */
21 /* turn this flag on to do minimal debug print of pops */
22 #define PRT_POPS false
23 /* this is limit to number of loops over H2 pops before giving up */
24 #define LIM_H2_POP_LOOP 10
25 /* this is a typical dissociation cross section (cm2) for H2 + Hnu -> 2H + ke */
26 /* >>chng 05 may 11, had been 2.5e-19 */
27 #define H2_DISS_ALLISON_DALGARNO 6e-19f
28 #include "cddefines.h"
29 #include "cddrive.h"
30 #include "atoms.h"
31 #include "conv.h"
32 #include "secondaries.h"
33 #include "pressure.h"
34 #include "trace.h"
35 #include "hmi.h"
36 #include "rt.h"
37 #include "radius.h"
38 #include "ipoint.h"
39 #include "phycon.h"
40 #include "thermal.h"
41 #include "dense.h"
42 #include "h2.h"
43 #include "mole.h"
44 #include "rfield.h"
45 #include "doppvel.h"
46 #include "lines_service.h"
47 
48 static realnum collider_density[N_X_COLLIDER];
49 static realnum collider_density_total_not_H2;
50 
51 void diatomics::H2_X_sink_and_source( void )
52 {
53  DEBUG_ENTRY( "diatomics::H2_X_sink_and_source()" );
54 
55  /* this is total density of all colliders, is only used for collisional dissociation
56  * rates for H2 are not included here, will be added separately*/
57  collider_density_total_not_H2 = collider_density[0] +
58  collider_density[1] + collider_density[4] +
59  dense.eden;
60 
61  for( long ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
62  {
63  H2_X_source[ipHi] = 0.;
64  H2_X_sink[ipHi] = 0.;
65  }
66 
67  double source_so_far = 0.;
68  double source_so_far_s = 0.;
69  double sink_so_far = 0.;
70  double pop_tot = 0.;
71  double sink_so_far_s = spon_diss_tot * H2_den_s;
72  double pop_tot_s = 0.;
73 
74  for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
75  {
76  /* array of energy sorted indices within X */
77  long iVibHi = ipVib_H2_energy_sort[ipHi];
78  long iRotHi = ipRot_H2_energy_sort[ipHi];
79 
80  /* count formation from grains and H- as a collisional formation process */
81  /* cm-3 s-1, evaluated in mole_H2_form */
82  H2_X_source[ipHi] += H2_X_formation[iVibHi][iRotHi];
83 
84  /*>>chng 05 sep 18, GS, H2 + e = H- + H*, H2_X_Hmin_back has units s-1 */
85  H2_X_sink[ipHi] += H2_X_Hmin_back[iVibHi][iRotHi];
86 
87  /* this represents collisional dissociation into continuum of X,
88  * rates are just guesses */
89  H2_X_sink[ipHi] += collider_density_total_not_H2 *
90  H2_coll_dissoc_rate_coef[iVibHi][iRotHi] * lgColl_deexec_Calc;
91 
92  /*>>chng 05 jul 20, GS, collisional dissociation with H2g and H2s are added here*/
93  H2_X_sink[ipHi] += hmi.H2_total*
94  H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] * lgColl_deexec_Calc;
95 
96  /* rate (s-1) out of this level */
97  if( mole_global.lgStancil )
98  {
99  H2_X_sink[ipHi] += Cont_Dissoc_Rate[0][iVibHi][iRotHi];
100  }
101  else
102  H2_X_sink[ipHi] += rfield.flux_accum[H2_ipPhoto[iVibHi][iRotHi]-1]*H2_DISS_ALLISON_DALGARNO;
103 
104  if ( states[ipHi].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
105  {
106  source_so_far_s += H2_X_source[ipHi];
107  sink_so_far_s += H2_X_sink[ipHi] * states[ipHi].Pop();
108  pop_tot_s += states[ipHi].Pop();
109  }
110  else
111  {
112  source_so_far += H2_X_source[ipHi];
113  sink_so_far += H2_X_sink[ipHi] * states[ipHi].Pop();
114  pop_tot += states[ipHi].Pop();
115  }
116  }
117 
118  // cm-3 s-1
119  double sink_tot = mole.sink_rate_tot(sp) * pop_tot;
120  // cm-3 s-1
121  double sink_left = sink_tot - sink_so_far;
122  // divide by population in X to get units s-1
123  ASSERT( pop_tot > 1e-10 * (*dense_total) );
124  sink_left /= pop_tot;
125  if( sink_left >= 0. )
126  {
127  for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
128  {
129  if( states[ipHi].energy().WN() <= ENERGY_H2_STAR || !hmi.lgLeiden_Keep_ipMH2s )
130  {
131  H2_X_sink[ipHi] += sink_left;
132  }
133  }
134  }
135 
136  // cm-3 s-1
137  fixit("kill the second term (sp_star) when H2* is killed in chemistry");
138  double sink_tot_s = mole.sink_rate_tot(sp_star) * pop_tot_s;
139  // cm-3 s-1
140  double sink_left_s = sink_tot_s - sink_so_far_s;
141  // divide by population in X to get units s-1
142  if( pop_tot_s > 1e-30 * (*dense_total) )
143  sink_left_s /= pop_tot_s;
144  else
145  sink_left_s = 0.;
146  if( sink_left_s >= 0. )
147  {
148  for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
149  {
150  if( states[ipHi].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
151  H2_X_sink[ipHi] += sink_left_s;
152  }
153  }
154 
155  fixit("kill the second term (sp_star) when H2* is killed in chemistry");
156  double source_tot = mole.source_rate_tot(sp);
157  double source_left = source_tot - source_so_far;
158  double source_tot_s = mole.source_rate_tot(sp_star);
159  double source_left_s = source_tot_s - source_so_far_s;
160  if( source_left+source_left_s >= 0. )
161  {
162  double rpop_lte = 1.;
163  double rpop_lte_s = 0.;
164  if (hmi.lgLeiden_Keep_ipMH2s)
165  {
166  double pop_lte = 0.;
167  double pop_lte_s = 0.;
168  for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
169  {
170  long iElec = states[ipHi].n();
171  long iVib = states[ipHi].v();
172  long iRot = states[ipHi].J();
173  if( states[ipHi].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
174  pop_lte_s += H2_populations_LTE[iElec][iVib][iRot];
175  else
176  pop_lte += H2_populations_LTE[iElec][iVib][iRot];
177  }
178  rpop_lte = 1./SDIV(pop_lte);
179  rpop_lte_s = 1./SDIV(pop_lte_s);
180  }
181  for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
182  {
183  long iElec = states[ipHi].n();
184  long iVib = states[ipHi].v();
185  long iRot = states[ipHi].J();
186  if( states[ipHi].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
187  H2_X_source[ipHi] += source_left_s * H2_populations_LTE[iElec][iVib][iRot]*rpop_lte_s;
188  else
189  H2_X_source[ipHi] += source_left * H2_populations_LTE[iElec][iVib][iRot]*rpop_lte;
190  }
191  }
192 
193  return;
194 }
195 
196 /*H2_X_coll_rate_evaluate find collisional rates within X -
197  * this is one time upon entry into H2_LevelPops */
198 void diatomics::H2_X_coll_rate_evaluate( void )
199 {
200  DEBUG_ENTRY( "diatomics::H2_X_coll_rate_evaluate()" );
201 
202  /* set collider density
203  * the colliders are:
204  * [0] = H
205  * [1], [5] = He (old and new cs data)
206  * [2] = H2 ortho
207  * [3] = H2 para
208  * [4] = H+ + H3+ */
209  /* atomic hydrogen */
210  collider_density[0] = dense.xIonDense[ipHYDROGEN][0];
211  /* all ortho h2 */
212  /* He - H2 */
213  collider_density[1] = dense.xIonDense[ipHELIUM][0];
214  /* H2 - H2(ortho) */
215  collider_density[2] = h2.ortho_density_f;
216  /* all para H2 */
217  collider_density[3] = h2.para_density_f;
218  /* protons - ionized hydrogen */
219  collider_density[4] = dense.xIonDense[ipHYDROGEN][1];
220  /* H3+ - assume that H3+ has same rates as proton */
221  collider_density[4] += (realnum)findspecieslocal("H3+")->den;
222 
223  ASSERT( fp_equal( hmi.H2_total_f ,collider_density[2]+collider_density[3]) );
224 
225  if( nTRACE >= n_trace_full )
226  {
227  fprintf(ioQQQ," Collider densities are:");
228  for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
229  {
230  fprintf(ioQQQ,"\t%.3e", collider_density[nColl]);
231  }
232  fprintf(ioQQQ,"\n");
233  }
234 
235  H2_X_coll_rate.zero();
236 
237  for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
238  {
239  if( lgColl_deexec_Calc )
240  {
241  /* excitation within X due to thermal particles */
242  for( long ipLo=0; ipLo<ipHi; ++ipLo )
243  {
244  /* collisional interactions with upper levels within X */
245  double colldown = 0.;
246  mr3ci CollRate = CollRateCoeff.begin(ipHi, ipLo);
247  for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
248  {
249  /* downward collision rate, units s-1 */
250  colldown += CollRate[nColl]*collider_density[nColl];
251  ASSERT( CollRate[nColl]*collider_density[nColl] >= 0. );
252  }
253  /* rate in from upper level, units cm-3 s-1 */
254  H2_X_coll_rate[ipHi][ipLo] += colldown;
255  }/* end loop over ipLo */
256  }
257  }
258 
259  return;
260 }
261 
262 /*H2_itrzn - average number of H2 pop evaluations per zone */
263 double diatomics::H2_itrzn( void )
264 {
265  if( lgEnabled && nH2_zone>0 )
266  {
267  return( (double)nH2_pops / (double)nH2_zone );
268  }
269  else
270  {
271  return 0.;
272  }
273 }
274 
275 /* set the ipCont struc element for the H2 molecule, called by ContCreatePointers */
276 void diatomics::H2_ContPoint( void )
277 {
278  if( !lgEnabled )
279  return;
280 
281  DEBUG_ENTRY( "H2_ContPoint()" );
282 
283  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
284  {
285  ASSERT( (*tr).Emis().Aul() > 0. );
286  (*tr).ipCont() = ipLineEnergy( (*tr).EnergyRyd(), label.c_str(), 0 );
287  (*tr).Emis().ipFine() = ipFineCont( (*tr).EnergyRyd());
288  }
289  return;
290 }
291 
292 /* ===================================================================== */
293 /* radiative acceleration due to H2 called in rt_line_driving */
294 double diatomics::H2_Accel(void)
295 {
296  /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */
297  if( !lgEnabled /*|| !nCall_this_zone*/ )
298  return 0.;
299 
300  DEBUG_ENTRY( "H2_Accel()" );
301 
302  /* this routine computes the line driven radiative acceleration */
303 
304  double drive = 0.;
305 
306  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
307  {
308  ASSERT( (*tr).ipCont() > 0 );
309  drive += (*tr).Emis().pump() * (*tr).Emis().PopOpc() * (*tr).EnergyErg();
310  }
311 
312  return drive;
313 }
314 
315 /* ===================================================================== */
316 /* rad pressure due to H2 lines called in PresTotCurrent */
317 double diatomics::H2_RadPress(void)
318 {
319  /* will be used to check on size of opacity, was capped at this value */
320  realnum smallfloat=SMALLFLOAT*10.f;
321 
322  /* radiation pressure sum is expensive - do not evaluate if we did not
323  * bother evaluating large molecule */
324  if( !lgEnabled || !nCall_this_zone )
325  return 0.;
326 
327  DEBUG_ENTRY( "H2_RadPress()" );
328 
329  realnum doppler_width = GetDopplerWidth( mass_amu );
330  double press = 0.;
331 
332  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
333  {
334  ASSERT( (*tr).ipCont() > 0 );
335  if( (*(*tr).Hi()).Pop() > smallfloat && (*tr).Emis().PopOpc() > smallfloat )
336  {
337  press += PressureRadiationLine( *tr, doppler_width );
338  }
339  }
340 
341  if(nTRACE >= n_trace_full)
342  fprintf(ioQQQ,
343  " H2_RadPress returns, radiation pressure is %.2e\n",
344  press );
345  return press;
346 }
347 
348 #if 0
349 /* ===================== */
350 /* internal energy of H2 */
351 double diatomics::H2_InterEnergy(void)
352 {
353  /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */
354  if( !lgEnabled /*|| !nCall_this_zone*/ )
355  return 0.;
356 
357  DEBUG_ENTRY( "H2_InterEnergy()" );
358 
359  double energy = 0.;
360  for( qList::iterator st = trans.states.begin(); st != trans.states.end(); ++st )
361  energy += st->Pop() * st->energy();
362 
363  return energy;
364 }
365 #endif
366 
367 /*H2_RT_diffuse do emission from H2 - called from RT_diffuse */
368 void diatomics::H2_RT_diffuse(void)
369 {
370  if( !lgEnabled || !nCall_this_zone )
371  return;
372 
373  DEBUG_ENTRY( "H2_RT_diffuse()" );
374 
375  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
376  {
377  qList::iterator Hi = (*tr).Hi();
378  if( (*Hi).n() > 0 )
379  continue;
380  (*tr).outline_resonance();
381  }
382 
383  return;
384 }
385 
386 /* RT for H2 lines */
387 void diatomics::H2_RTMake( linefunc line_one )
388 {
389  if( !lgEnabled )
390  return;
391 
392  DEBUG_ENTRY( "H2_RTMake()" );
393 
394  realnum doppler_width = GetDopplerWidth( mass_amu );
395  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
396  {
397  /* >>chng 03 jun 18, added 4th parameter in call to this routine - says to not
398  * include self-shielding of line across this zone. This introduces a dr dependent
399  * variation in the line pumping rate, which made H2 abundance fluctuate due to
400  * Solomon process having slight dr-caused mole. */
401  line_one( *tr, false, 0.f, doppler_width );
402  }
403 
404  return;
405 }
406 
407 /* increment optical depth for the H2 molecule, called from RT_tau_inc which is called by cloudy,
408  * one time per zone */
409 void diatomics::H2_RT_tau_inc(void)
410 {
411  /* >>chng 05 jan 26, now use LTE populations for small H2 abundance case, since electronic
412  * lines become self-shielding surprisingly quickly */
413  if( !lgEnabled /*|| !nCall_this_zone*/ )
414  return;
415 
416  DEBUG_ENTRY( "H2_RT_tau_inc()" );
417 
418  /* remember largest and smallest chemistry renormalization factor -
419  * if both networks are parallel will be unity,
420  * but only do this after we have stable solution */
421  if( nzone > 0 && nCall_this_iteration>2 )
422  {
423  renorm_max = MAX2( H2_renorm_chemistry , renorm_max );
424  renorm_min = MIN2( H2_renorm_chemistry , renorm_min );
425  }
426 
427  realnum doppler_width = GetDopplerWidth( mass_amu );
428  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
429  {
430  ASSERT( (*tr).ipCont() > 0 );
431  RT_line_one_tauinc( *tr,-9, -9, -9, -9, doppler_width );
432  }
433 
434  return;
435 }
436 
437 
438 /* initialize optical depths in H2, called from RT_tau_init */
439 void diatomics::H2_LineZero( void )
440 {
441  if( !lgEnabled )
442  return;
443 
444  DEBUG_ENTRY( "H2_LineZero()" );
445 
446  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
447  {
448  (*tr).Zero();
449  }
450 
451  return;
452 }
453 
454 /* the large H2 molecule, called from RT_tau_reset */
455 void diatomics::H2_RT_tau_reset( void )
456 {
457  if( !lgEnabled )
458  return;
459 
460  DEBUG_ENTRY( "H2_RT_tau_reset()" );
461 
462  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
463  {
464  RT_line_one_tau_reset( *tr );
465  }
466 
467  return;
468 }
469 
470 /*H2_Level_low_matrix evaluate lower populations within X */
471 void diatomics::H2_Level_low_matrix(
472  /* total abundance within matrix */
473  realnum abundance )
474 {
475  /* will need to MALLOC space for these but only on first call */
476  bool lgDoAs;
477  int nNegPop;
478  bool lgDeBug,
479  lgZeroPop;
480  double rot_cooling , dCoolDT;
481 
482  DEBUG_ENTRY( "H2_Level_low_matrix()" );
483 
484  /* option to not use the matrix */
485  if( nXLevelsMatrix <= 1 )
486  {
487  return;
488  }
489 
490  if( lgFirst )
491  {
492  /* check that not more levels than there are in X */
493  if( nXLevelsMatrix > nLevels_per_elec[0] )
494  {
495  /* number is greater than number of levels within X */
496  fprintf( ioQQQ,
497  " The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n",
498  nLevels_per_elec[0]);
499  cdEXIT(EXIT_FAILURE);
500  }
501  /* will never do this again */
502  lgFirst = false;
503  /* remember how much space we malloced in case ever called with more needed */
504  /* >>chng 05 jan 19, allocate max number of levels
505  ndimMalloced = nXLevelsMatrix;*/
506  ndimMalloced = nLevels_per_elec[0];
507  /* allocate the 1D arrays*/
508  excit.resize( ndimMalloced );
509  stat_levn.resize( ndimMalloced );
510  pops.resize( ndimMalloced );
511  create.resize( ndimMalloced );
512  destroy.resize( ndimMalloced );
513  depart.resize( ndimMalloced );
514  /* create space for the 2D arrays */
515  AulPump = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
516  CollRate_levn = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
517  AulDest = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
518  AulEscp = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
519  AulPump[0] = ((double *)MALLOC((size_t)(ndimMalloced*ndimMalloced)*sizeof(double )));
520  CollRate_levn[0] = ((double *)MALLOC((size_t)(ndimMalloced*ndimMalloced)*sizeof(double )));
521  AulDest[0] = ((double *)MALLOC((size_t)(ndimMalloced*ndimMalloced)*sizeof(double )));
522  AulEscp[0] = ((double *)MALLOC((size_t)(ndimMalloced*ndimMalloced)*sizeof(double )));
523  for( long i=1; i< ndimMalloced; ++i )
524  {
525  AulPump[i] = AulPump[i-1]+ndimMalloced;
526  CollRate_levn[i] = CollRate_levn[i-1]+ndimMalloced;
527  AulDest[i] = AulDest[i-1]+ndimMalloced;
528  AulEscp[i] = AulEscp[i-1]+ndimMalloced;
529  }
530 
531  /* the statistical weights of the levels
532  * and excitation potentials of each level relative to ground */
533  for( long j=0; j < ndimMalloced; j++ )
534  {
535  /* statistical weights for each level */
536  stat_levn[j] = states[j].g();
537  /* excitation energy of each level relative to ground, in K */
538  excit[j] = states[j].energy().K();
539  }
540  }
541  /* end malloc space and creating constant terms */
542 
543  /* this is test for call with too many rotation levels to handle -
544  * logic needs for largest model atom to be called first */
545  if( nXLevelsMatrix > ndimMalloced )
546  {
547  fprintf(ioQQQ," H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n");
548  cdEXIT(EXIT_FAILURE);
549  }
550 
551  /* all elements are used, and must be set to zero */
552  for( long i=0; i < nXLevelsMatrix; i++ )
553  {
554  pops[i] = 0.;
555  depart[i] = 0;
556  }
557 
558  /* do we need to reevaluate radiative quantities? only do this one time per zone */
559  if( nzone!=nzoneAsEval || iteration!=iterationAsEval || nXLevelsMatrix!=levelAsEval)
560  {
561  lgDoAs = true;
562  nzoneAsEval = nzone;
563  iterationAsEval = iteration;
564  levelAsEval = nXLevelsMatrix;
565  ASSERT( levelAsEval <= ndimMalloced );
566  }
567  else
568  {
569  lgDoAs = false;
570  }
571 
572  /* all elements are used, and must be set to zero */
573  if( lgDoAs )
574  {
575  for( long i=0; i < nXLevelsMatrix; i++ )
576  {
577  pops[i] = 0.;
578  depart[i] = 0;
579  for( long j=0; j < nXLevelsMatrix; j++ )
580  {
581  AulEscp[j][i] = 0.;
582  AulDest[j][i] = 0.;
583  AulPump[j][i] = 0.;
584  CollRate_levn[j][i] = 0.;
585  }
586  }
587  }
588 
589  /* find all radiative interactions within matrix, and between
590  * matrix and upper X and excited electronic states */
591  for( long ilo=0; ilo < nXLevelsMatrix; ilo++ )
592  {
593  long iRot = ipRot_H2_energy_sort[ilo];
594  long iVib = ipVib_H2_energy_sort[ilo];
595 
596  /* H2_X_sink[ilo] includes all processes that destroy H2 in one step,
597  * these include cosmic ray ionization and dissociation, photodissociation,
598  * BUT NOT THE SOLOMON process, which, directly, only goes to excited
599  * electronic states */
600  destroy[ilo] = H2_X_sink[ilo];
601 
602  /* rates H2 is created from grains and H- units cm-3 s-1, evaluated in mole_H2_form */
603  create[ilo] = H2_X_source[ilo];
604 
605  /* this loop does radiative decays from upper states inside matrix,
606  * and upward pumps within matrix region into this lower level */
607  if( lgDoAs )
608  {
609  for( long ihi=ilo+1; ihi<nXLevelsMatrix; ++ihi )
610  {
611  ASSERT( states[ihi].energy().WN() <= states[nXLevelsMatrix-1].energy().WN() );
612  /* general case - but line may not actually exist */
613  if( lgH2_radiative[ihi][ilo] )
614  {
615  const TransitionList::iterator&tr = trans.begin()+ ipTransitionSort[ihi][ilo] ;
616  ASSERT( (*tr).ipCont() > 0 );
617 
618  /* NB - the destruction probability is included in
619  * the total and the destruction is set to zero
620  * since we want to only count one ots rate, in
621  * main calling routine, and do not want matrix
622  * solver below to include it */
623  AulEscp[ihi][ilo] = (*tr).Emis().Aul()*(
624  (*tr).Emis().Pesc_total() );
625  AulDest[ihi][ilo] = (*tr).Emis().Aul()*(*tr).Emis().Pdest();
626  AulPump[ilo][ihi] = (*tr).Emis().pump();
627  AulPump[ihi][ilo] = AulPump[ilo][ihi]*stat_levn[ilo]/stat_levn[ihi];
628  }
629  }
630  }
631 
632  double rateout = 0.;
633  double ratein = 0.;
634  /* now do all levels within X, which are above nXLevelsMatrix,
635  * the highest level inside the matrix */
636  for( long ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
637  {
638  if( lgH2_radiative[ihi][ilo] )
639  {
640  const TransitionList::iterator&tr = trans.begin()+ ipTransitionSort[ihi][ilo] ;
641  ASSERT( (*tr).ipCont() > 0 );
642 
643  /* these will enter as net creation terms in creation vector, with
644  * units cm-3 s-1
645  * radiative transitions from above the matrix within X */
646  ratein +=
647  (*(*tr).Hi()).Pop() * (
648  (*tr).Emis().Aul()*( (*tr).Emis().Ploss() ) +
649  (*tr).Emis().pump() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g() );
650  /* rate out has units s-1 - destroys current lower level */
651  rateout += (*tr).Emis().pump();
652  }
653  }
654 
655  /* all states above the matrix but within X */
656  create[ilo] += ratein;
657 
658  /* rates out of matrix into levels in X but above matrix */
659  destroy[ilo] += rateout;
660 
661  /* Solomon process, this sum dos all pump and decays from all electronic excited states */
662  /* radiative rates [cm-3 s-1] from electronic excited states into X only vibration and rot */
663  create[ilo] += H2_X_rate_from_elec_excited[iVib][iRot];
664 
665  /* radiative & cosmic ray rates [s-1] to electronic excited states from X only vibration and rot */
666  destroy[ilo] += H2_X_rate_to_elec_excited[iVib][iRot];
667  }
668 
669  /* this flag set with atom H2 trace matrix */
670  if( nTRACE >= n_trace_matrix )
671  lgDeBug = true;
672  else
673  lgDeBug = false;
674 
675  /* now evaluate the rates for all transitions within matrix */
676  for( long ilo=0; ilo < nXLevelsMatrix; ilo++ )
677  {
678  double EloK = states[ilo].energy().K();
679 
680  if(lgDeBug)fprintf(ioQQQ,"DEBUG H2_Level_low_matrix, ilo=%li",ilo);
681  for( long ihi=ilo+1; ihi < nXLevelsMatrix; ihi++ )
682  {
683  /* >>chng 05 may 31, replace with simple expresion */
684  CollRate_levn[ihi][ilo] = H2_X_coll_rate[ihi][ilo];
685 
686  if(lgDeBug)fprintf(ioQQQ,"\t%.1e",CollRate_levn[ihi][ilo]);
687 
688  /* now get upward excitation rate - units s-1 */
689  // use dsexp() explicitly rather than divide Boltzmann() factors
690  // to avoid problems with underflow at low Te (see ticket #284)
691  CollRate_levn[ilo][ihi] = CollRate_levn[ihi][ilo]*
692  dsexp( (states[ihi].energy().K() - EloK) / phycon.te )*
693  states[ihi].g() / states[ilo].g();
694  }
695  if(lgDeBug)fprintf(ioQQQ,"\n");
696 
697  /* now do all collisions for levels within X, which are above nXLevelsMatrix,
698  * the highest level inside the matrix */
699  for( long ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
700  {
701  /* first do downward deexcitation rate */
702  /* >>chng 04 sep 14, do all levels */
703  /* >>chng 05 may 31, use summed rate */
704  double ratein = H2_X_coll_rate[ihi][ilo];
705  if(lgDeBug)fprintf(ioQQQ,"\t%.1e",ratein);
706 
707  /* now get upward excitation rate */
708  // use dsexp() explicitly rather than divide Boltzmann() factors
709  // to avoid problems with underflow at low Te (see ticket #284)
710  double rateout = ratein *
711  dsexp( (states[ihi].energy().K() - EloK) / phycon.te )*
712  states[ihi].g()/states[ilo].g();
713 
714  /* these are general entries and exits going into vector */
715  create[ilo] += ratein * states[ihi].Pop();
716  destroy[ilo] += rateout;
717  }
718  if(lgDeBug)fprintf(ioQQQ,"\n");
719  }
720 
721  /* H2 grain interactions */
722  {
723  for( long ihi=2; ihi < nXLevelsMatrix; ihi++ )
724  {
725  long iVibHi = ipVib_H2_energy_sort[ihi];
726  long iRotHi = ipRot_H2_energy_sort[ihi];
727 
728  /* collisions with grains goes to either J=1 or J=0 depending on
729  * spin of upper level - this conserves op ratio - following
730  * var is 1 if ortho, 0 if para, so this conserves op ratio
731  * units are s-1 */
732  CollRate_levn[ihi][H2_lgOrtho[0][iVibHi][iRotHi]] += rate_grain_op_conserve;
733  }
734 
735  /* H2 ortho - para conversion on grain surface,
736  * rate (s-1) all v,J levels go to 0 or 1 */
737  CollRate_levn[1][0] +=
738  (realnum)(rate_grain_J1_to_J0);
739  }
740 
741  /* now all levels in X above the matrix */
742  for( long ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
743  {
744  long iVibHi = ipVib_H2_energy_sort[ihi];
745  long iRotHi = ipRot_H2_energy_sort[ihi];
746 
747  /* these collisions all go into 0 or 1 depending on whether upper level was ortho or para
748  * units are cm-3 s-1 - rate new molecules appear in matrix */
749  create[H2_lgOrtho[0][iVibHi][iRotHi]] += states[ihi].Pop() * rate_grain_op_conserve;
750  }
751 
752  /* debug print individual contributors to matrix elements */
753  {
754  enum {DEBUG_LOC=false};
755  if( DEBUG_LOC || lgDeBug)
756  {
757  fprintf(ioQQQ,"DEBUG H2 matexcit");
758  for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
759  {
760  fprintf(ioQQQ,"\t%li",ilo );
761  }
762  fprintf(ioQQQ,"\n");
763  for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
764  {
765  fprintf(ioQQQ,"\t%.2e",excit[ihi] );
766  }
767  fprintf(ioQQQ,"\n");
768  for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
769  {
770  fprintf(ioQQQ,"\t%.2e",stat_levn[ihi] );
771  }
772  fprintf(ioQQQ,"\n");
773 
774  fprintf(ioQQQ,"AulEscp[n][]\\[][n] = Aul*Pesc\n");
775  for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
776  {
777  fprintf(ioQQQ,"\t%li",ilo );
778  }
779  fprintf(ioQQQ,"\n");
780  for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
781  {
782  fprintf(ioQQQ,"%li", ihi);
783  for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
784  {
785  fprintf(ioQQQ,"\t%.2e",AulEscp[ilo][ihi] );
786  }
787  fprintf(ioQQQ,"\n");
788  }
789 
790  fprintf(ioQQQ,"AulPump [n][]\\[][n]\n");
791  for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
792  {
793  fprintf(ioQQQ,"\t%li",ilo );
794  }
795  fprintf(ioQQQ,"\n");
796  for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
797  {
798  fprintf(ioQQQ,"%li", ihi);
799  for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
800  {
801  fprintf(ioQQQ,"\t%.2e",AulPump[ihi][ilo] );
802  }
803  fprintf(ioQQQ,"\n");
804  }
805 
806  fprintf(ioQQQ,"CollRate_levn [n][]\\[][n]\n");
807  for( long ilo=0; ilo<nXLevelsMatrix; ++ilo )
808  {
809  fprintf(ioQQQ,"\t%li",ilo );
810  }
811  fprintf(ioQQQ,"\n");
812  for( long ihi=0; ihi<nXLevelsMatrix;++ihi)
813  {
814  fprintf(ioQQQ,"%li", ihi);
815  for( long ilo=0; ilo<nXLevelsMatrix; ++ilo )
816  {
817  fprintf(ioQQQ,"\t%.2e",CollRate_levn[ihi][ilo] );
818  }
819  fprintf(ioQQQ,"\n");
820  }
821  fprintf(ioQQQ,"SOURCE");
822  for( long ihi=0; ihi<nXLevelsMatrix;++ihi)
823  {
824  fprintf(ioQQQ,"\t%.2e",create[ihi]);
825  }
826  fprintf(ioQQQ,"\nSINK");
827  for( long ihi=0; ihi<nXLevelsMatrix;++ihi)
828  {
829  fprintf(ioQQQ,"\t%.2e",destroy[ihi]);
830  }
831  fprintf(ioQQQ,"\n");
832  }
833  }
834 
835  static Atom_LevelN atom_levelN;
836  atom_levelN(
837  nXLevelsMatrix,
838  abundance,
839  &stat_levn[0],
840  &excit[0],
841  'K',
842  &pops[0],
843  &depart[0],
844  /* net transition rate, A * escape prob, s-1, indices are [upper][lower] */
845  AulEscp,
846  AulDest,
847  AulPump,
848  CollRate_levn,
849  &create[0],
850  &destroy[0],
851  &rot_cooling,
852  &dCoolDT,
853  " H2 ",
854  lgPrtMatrix,
855  /* nNegPop positive if negative pops occurred, negative if too cold */
856  &nNegPop,
857  &lgZeroPop,
858  lgDeBug ); /* option to print stuff - set to true for debug printout */
859 
860  if( nNegPop > 0 )
861  {
862  // Recovery procedure -- assume that sources have overshot
863  fprintf(ioQQQ," H2_Level_low_matrix called atom_levelN which returned negative populations.\n");
864  //ConvFail( "pops" , "H2" );
865  double totpop = 0.;
866  for( long i=0; i< nXLevelsMatrix; ++i )
867  {
868  if (pops[i] < 0.0)
869  pops[i] = 0.0;
870  totpop += pops[i];
871  }
872  totpop = abundance/totpop;
873  for( long i=0; i< nXLevelsMatrix; ++i )
874  {
875  pops[i] *= totpop;
876  }
877  }
878 
879  for( long i=0; i< nXLevelsMatrix; ++i )
880  {
881  states[i].Pop() = pops[i];
882  fixit("need to also store DepartCoef above X");
883  states[i].DepartCoef() = depart[i];
884  }
885 
886  if( 0 && nTRACE >= n_trace_full)
887  {
888  /* print pops that came out of matrix */
889  fprintf(ioQQQ,"\n DEBUG H2_Level_lowJ dense_total: %.3e matrix rel pops\n", *dense_total);
890  fprintf(ioQQQ,"v\tJ\tpop\n");
891  for( long i=0; i<nXLevelsMatrix; ++i )
892  {
893  long iRot = ipRot_H2_energy_sort[i];
894  long iVib = ipVib_H2_energy_sort[i];
895  fprintf(ioQQQ,"%3li\t%3li\t%.3e\t%.3e\t%.3e\n",
896  iVib , iRot , states[i].Pop(), create[i] , destroy[i]);
897  }
898  }
899 
900  /* nNegPop positive if negative pops occurred, negative if too cold */
901  return;
902 }
903 /* do level populations for H2, called by Hydrogenic after ionization and H chemistry
904  * has been recomputed */
905 void diatomics::H2_LevelPops( bool &lgPopsConverged, double &old_val, double &new_val )
906 {
907  DEBUG_ENTRY( "H2_LevelPops()" );
908 
909  string convLabel;
910  if( this==&h2 )
911  convLabel = "H2_LOOPS";
912  else if( this==&hd )
913  convLabel = "HD_LOOPS";
914  else
915  TotalInsanity();
916  static ConvergenceCounter cctr_diatom_l =conv.register_(convLabel);
917 
918  /* H2 not on, so space not allocated and return,
919  * also return if calculation has been declared a failure */
920  if( !lgEnabled || lgAbort )
921  {
922  // need to do this even if not doing big model
923  CalcPhotoionizationRate();
924  return;
925  }
926 
927  double old_solomon_rate=-1.;
928  long int n_pop_oscil = 0;
929  int kase=0;
930  bool lgConv_h2_soln,
931  lgPopsConv_total,
932  lgPopsConv_relative,
933  lgHeatConv,
934  lgSolomonConv,
935  lgOrthoParaRatioConv;
936  double quant_old=-1.,
937  quant_new=-1.;
938 
939  bool lgH2_pops_oscil=false,
940  lgH2_pops_ever_oscil=false;
941 
942  /* keep track of changes in population */
943  double PopChgMax_relative=0. , PopChgMaxOld_relative=0., PopChgMax_total=0., PopChgMaxOld_total=0.;
944  long int iRotMaxChng_relative , iVibMaxChng_relative,
945  iRotMaxChng_total , iVibMaxChng_total,
946  nXLevelsMatrix_save;
947  double popold_relative , popnew_relative , popold_total , popnew_total;
948  /* reason not converged */
949  char chReason[100];
950 
951  /* these are convergence criteria - will be increased during search phase */
952  double converge_pops_relative=conv.IonizErrorAllowed/3.,
953  converge_pops_total=1e-3,
954  converge_ortho_para=1e-2;
955 
956  double dens_rel_to_lim_react = mole.species[sp->index].xFracLim;
957 
958  if(nTRACE >= n_trace_full )
959  {
960  fprintf(ioQQQ,
961  "\n***************H2_LevelPops %s call %li this iteration, zone is %.2f, H2/H:%.e Te:%e ne:%e\n",
962  label.c_str(),
963  nCall_this_iteration,
964  fnzone,
965  dens_rel_to_lim_react,
966  phycon.te,
967  dense.eden
968  );
969  }
970  else if( nTRACE >= n_trace_final )
971  {
972  static long int nzone_prt=-1;
973  if( nzone!=nzone_prt )
974  {
975  nzone_prt = nzone;
976  fprintf(ioQQQ,"DEBUG zone %li species %s rel_to_lim:%.3e Te:%.3e *ne:%.3e n(%s):%.3e\n",
977  nzone,
978  label.c_str(),
979  dens_rel_to_lim_react,
980  phycon.te,
981  dense.eden,
982  label.c_str(),
983  *dense_total );
984  }
985  }
986 
987  CalcPhotoionizationRate();
988 
989  /* evaluate Boltzmann factors and LTE unit population - for trivial abundances
990  * LTE populations are used in place of full solution */
991  mole_H2_LTE();
992 
993  /* zero out populations and cooling, and return, if H2 fraction is small
994  * but, if H2 has ever been done, redo irregardless of abundance -
995  * if large H2 is ever evaluated then mole.H2_to_H_limit is ignored */
996  if( (!lgEvaluated && dens_rel_to_lim_react < H2_to_H_limit )
997  || *dense_total < 1e-20 )
998  {
999  /* will not compute the molecule */
1000  if( nTRACE >= n_trace_full )
1001  fprintf(ioQQQ,
1002  " H2_LevelPops %s pops too small, not computing, set to LTE and return, H2/H is %.2e and H2_to_H_limit is %.2e.\n",
1003  label.c_str(),
1004  dens_rel_to_lim_react,
1005  H2_to_H_limit);
1006  H2_zero_pops_too_low();
1007  fixit("set lgEvaluated = false here?");
1008  /* end of zero abundance branch */
1009  return;
1010  }
1011 
1012  /* check whether we need to update the H2_Boltzmann factors, LTE level populations,
1013  * and partition function. LTE level pops normalized by partition function,
1014  * so sum of pops is unity */
1015 
1016  /* say that H2 has been computed, ignore previous limit to abund
1017  * in future - this is to prevent oscillations as model is engaged */
1018  lgEvaluated = true;
1019  /* end loop setting H2_Boltzmann factors, partition function, and LTE populations */
1020 
1021  /* >>chng 05 jun 21,
1022  * during search phase we want to use full matrix - save number of levels so that
1023  * we can restore it */
1024  nXLevelsMatrix_save = nXLevelsMatrix;
1025  fixit("this does not appear to be necessary and may be counterproductive");
1026  if( conv.lgSearch )
1027  {
1028  nXLevelsMatrix = nLevels_per_elec[0];
1029  }
1030 
1031  /* 05 oct 27, had only reevaluated collision rates when 5% change in temperature
1032  * caused temp failures in large G0 sims -
1033  * do not check whether we need to update the collision rates but
1034  * reevaluate them always
1035  * >>chng 05 nov 04, above caused a 25% increase in the exec time for constant-T sims
1036  * in test suite- original code had reevaluated if > 0.05 change in T - was too much
1037  * change to 10x smaller, change > 0.005 */
1038  if( !fp_equal(phycon.te,TeUsedColl) )
1039  {
1040  H2_CollidRateEvalAll();
1041  TeUsedColl = phycon.te;
1042  }
1043 
1044  /* set the populations when this is the first call to this routine on
1045  * current iteration- will use LTE populations - populations were set by
1046  * call to mole_H2_LTE before above block */
1047  if( nCall_this_iteration==0 || lgLTE )
1048  {
1049  /* very first call so use LTE populations */
1050  if(nTRACE >= n_trace_full )
1051  fprintf(ioQQQ,"%s 1st call - using LTE level pops\n", label.c_str() );
1052 
1053  H2_den_s = 0.;
1054  H2_den_g = 0.;
1055  for( qList::iterator st = states.begin(); st != states.end(); ++st )
1056  {
1057  long iElec = (*st).n();
1058  long iVib = (*st).v();
1059  long iRot = (*st).J();
1060  /* LTE populations are for unit H2 density, so need to multiply
1061  * by total H2 density */
1062  double pop = H2_populations_LTE[iElec][iVib][iRot] * (*dense_total);
1063  H2_old_populations[iElec][iVib][iRot] = pop;
1064  (*st).Pop() = pop;
1065  /* find current population in H2s and H2g */
1066  if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
1067  {
1068  H2_den_s += pop;
1069  }
1070  else
1071  {
1072  H2_den_g += pop;
1073  }
1074  }
1075 
1076  /* first guess at ortho and para densities */
1077  ortho_density = 0.75* (*dense_total);
1078  para_density = 0.25* (*dense_total);
1079  {
1080  ortho_density_f = (realnum)ortho_density;
1081  para_density_f = (realnum)para_density;
1082  }
1083  ortho_para_current = ortho_density / SDIV( para_density );
1084  ortho_para_older = ortho_para_current;
1085  ortho_para_old = ortho_para_current;
1086  /* this is the fraction of the H2 pops that are within the levels done with a matrix */
1087  frac_matrix = 1.;
1088  }
1089 
1090  // make some population sums, normalize total to value handed from chemistry
1091  {
1092  pops_per_vib.zero();
1093  fill_n( pops_per_elec, N_ELEC, 0. );
1094  double pop_total = 0.;
1095  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1096  {
1097  long iElec = (*st).n();
1098  long iVib = (*st).v();
1099 
1100  pop_total += (*st).Pop();
1101  pops_per_elec[iElec] += (*st).Pop();
1102  pops_per_vib[iElec][iVib] += (*st).Pop();
1103  }
1104  ASSERT( pops_per_elec[0]>SMALLFLOAT );
1105  // Now renorm the old populations to the correct current H2 density.
1106  H2_renorm_chemistry = *dense_total/ SDIV(pop_total);
1107  }
1108 
1109  if(nTRACE >= n_trace_full)
1110  fprintf(ioQQQ,
1111  "%s H2_renorm_chemistry is %.4e, *dense_total is %.4e pops_per_elec[0] is %.4e\n",
1112  label.c_str(),
1113  H2_renorm_chemistry ,
1114  *dense_total,
1115  pops_per_elec[0]);
1116 
1117  /* renormalize all level populations for the current chemical solution */
1118  for( qList::iterator st = states.begin(); st != states.end(); ++st )
1119  {
1120  long iElec = (*st).n();
1121  long iVib = (*st).v();
1122  long iRot = (*st).J();
1123 
1124  (*st).Pop() *= H2_renorm_chemistry;
1125  H2_old_populations[iElec][iVib][iRot] = (*st).Pop();
1126  }
1127 
1128  if(nTRACE >= n_trace_full )
1129  fprintf(ioQQQ,
1130  " H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n",
1131  pops_per_elec[0],
1132  *dense_total);
1133 
1134  /* >>chng 05 feb 10, reset convergence criteria if we are in search phase */
1135  fixit("I suspect this is counterproductive. Test without it -- Ryan");
1136  if( conv.lgSearch )
1137  {
1138  converge_pops_relative *= 2.; /*def is 0.1 */
1139  converge_pops_total *= 3.; /*def is 1e-3*/
1140  converge_ortho_para *= 3.; /*def is 1e-2*/
1141  }
1142 
1143  if( !conv.nTotalIoniz )
1144  mole_update_species_cache();
1145 
1146  /* update state specific rates in H2_X_formation (cm-3 s-1) that H2 forms from grains and H- */
1147  mole_H2_form();
1148 
1149  /* evaluate total collision rates */
1150  H2_X_coll_rate_evaluate();
1151 
1152  /* this flag will say whether H2 populations have converged,
1153  * by comparing old and new values */
1154  lgConv_h2_soln = false;
1155  /* this will count number of passes around following loop */
1156  long loop_h2_pops = 0;
1157  {
1158  if( nzone != nzoneEval )
1159  {
1160  nzoneEval = nzone;
1161  /* this is number of zones with H2 solution in this iteration */
1162  ++nH2_zone;
1163  }
1164  }
1165 
1166  if( lgLTE )
1167  lgConv_h2_soln = true;
1168 
1169  /* begin - start level population solution
1170  * first do electronic excited states, Lyman, Werner, etc
1171  * using old solution for X
1172  * then do matrix if used, then solve for pops of rest of X
1173  * >>chng 04 apr 06, subtract number of oscillations from limit - don't waste loops
1174  * if solution is unstable */
1175  while( loop_h2_pops < LIM_H2_POP_LOOP-n_pop_oscil && !lgConv_h2_soln && !lgLTE )
1176  {
1177  ++cctr_diatom_l;
1178 
1179  /* this is number of trips around loop this time */
1180  ++loop_h2_pops;
1181  /* this is number of times through this loop in entire iteration */
1182  ++nH2_pops;
1183 
1184  /* radiative rates [cm-3 s-1] from electronic excited states into X vibration and rot */
1185  H2_X_rate_from_elec_excited.zero();
1186  /* radiative & cosmic ray rates [s-1] to electronic excited states from X */
1187  H2_X_rate_to_elec_excited.zero();
1188  H2_rad_rate_out.zero();
1189  H2_rad_rate_in.zero();
1190  pops_per_vib.zero();
1191  fill_n( pops_per_elec, N_ELEC, 0. );
1192 
1193  SolveExcitedElectronicLevels();
1194 
1195  /* evaluate rates that destroy or create ground electronic state */
1196  H2_X_sink_and_source();
1197 
1198  /* above set pops of excited electronic levels and found rates between them and X -
1199  * now solve highly excited levels within the X state by back-substitution */
1200  SolveSomeGroundElectronicLevels();
1201 
1202  /* now do lowest levels populations with matrix,
1203  * these should be collisionally dominated */
1204  if( nXLevelsMatrix )
1205  {
1206  H2_Level_low_matrix(
1207  /* the total abundance - frac_matrix is fraction of pop that was in these
1208  * levels the last time this was done */
1209  (*dense_total) * (realnum)frac_matrix );
1210  }
1211  if(nTRACE >= n_trace_full)
1212  {
1213  long iElecHi = 0;
1214  fprintf(ioQQQ," Rel pop(e=%li)" ,iElecHi);
1215  }
1216 
1217  /* find ortho and para densites, sum of pops in each vibration */
1218  /* this will become total pop is X, which will be renormed to equal *dense_total */
1219  {
1220  pops_per_elec[0] = 0.;
1221  for( md2i it = pops_per_vib.begin(0); it != pops_per_vib.end(0); ++it )
1222  *it = 0.;
1223 
1224  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1225  {
1226  long iElec = (*st).n();
1227  if( iElec > 0 ) continue;
1228  long iVib = (*st).v();
1229  pops_per_elec[iElec] += (*st).Pop();
1230  pops_per_vib[iElec][iVib] += (*st).Pop();
1231  }
1232 
1233  /* print sum of populations in each vibration if trace on */
1234  if(nTRACE >= n_trace_full)
1235  for( md2ci it = pops_per_vib.begin(0); it != pops_per_vib.end(0); ++it )
1236  fprintf(ioQQQ,"\t%.2e", *it/(*dense_total));
1237 
1238  ASSERT( pops_per_elec[0] > SMALLFLOAT );
1239  }
1240  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
1241  if(nTRACE >= n_trace_full)
1242  {
1243  fprintf(ioQQQ,"\n");
1244  /* print the ground vibration state */
1245  fprintf(ioQQQ," Rel pop(0,J)");
1246  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1247  {
1248  long iElec = (*st).n();
1249  if( iElec > 0 ) continue;
1250  long iVib = (*st).v();
1251  if( iVib > 0 ) continue;
1252  fprintf(ioQQQ,"\t%.2e", (*st).Pop()/(*dense_total) );
1253  }
1254  fprintf(ioQQQ,"\n");
1255  }
1256 
1257  {
1258  double pop_total = 0.;
1259  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1260  pop_total += (*st).Pop();
1261 
1262  // The ratio of H2 that came out of the chemistry network to what we just obtained.
1263  double H2_renorm_conserve = *dense_total/ SDIV(pop_total);
1264 
1265  if (0)
1266  fprintf(ioQQQ,"DEBUG H2 %ld %ld %g %g %g %g %g\n",
1267  nzone, loop_h2_pops, H2_renorm_conserve,frac_matrix,pop_total,states[0].Pop(),states[1].Pop());
1268 
1269  /* renormalize populations - were updated by renorm when routine entered,
1270  * before pops determined - but population determinations above do not have a sum rule on total
1271  * population - this renorm is to preserve total population */
1272  pops_per_vib.zero();
1273  fill_n( pops_per_elec, N_ELEC, 0. );
1274  for( qList::iterator st = states.begin(); st != states.end(); ++st )
1275  {
1276  (*st).Pop() *= H2_renorm_conserve;
1277  long iElec = (*st).n();
1278  long iVib = (*st).v();
1279  pops_per_elec[iElec] += (*st).Pop();
1280  pops_per_vib[iElec][iVib] += (*st).Pop();
1281  }
1282  }
1283 
1284  /* now find population in states done with matrix - this is only used to pass
1285  * to matrix solver */
1286  {
1287  double sum_pops_matrix = 0.;
1288  for( long i=0; i<nXLevelsMatrix; ++i )
1289  {
1290  sum_pops_matrix += states[i].Pop();
1291  }
1292  /* this is self consistent since pops_per_elec[0] came from current soln,
1293  * as did the matrix. pops will be renormalized by results from the chemistry
1294  * a few lines down */
1295  frac_matrix = sum_pops_matrix / SDIV(*dense_total);
1296  }
1297 
1298  /* these will do convergence check */
1299  PopChgMaxOld_relative = PopChgMax_relative;
1300  PopChgMaxOld_total = PopChgMax_total;
1301  PopChgMax_relative = 0.;
1302  PopChgMax_total = 0.;
1303  iRotMaxChng_relative =-1;
1304  iVibMaxChng_relative = -1;
1305  iRotMaxChng_total =-1;
1306  iVibMaxChng_total = -1;
1307  popold_relative = 0.;
1308  popnew_relative = 0.;
1309  popold_total = 0.;
1310  popnew_total = 0.;
1311 
1312  // *****************************************
1313  // *****************************************
1314  // *****************************************
1315  // *****************************************
1316  // Should be able to extract this loop!!!
1317  // *****************************************
1318  // *****************************************
1319  // *****************************************
1320  // *****************************************
1321  {
1322  /* this loop first checks for largest changes in populations, to determine whether
1323  * we have converged, then updates the population array with a new value,
1324  * which may be a mean of old and new
1325  * update populations check convergence converged */
1326  double sumold = 0.;
1327 
1328  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1329  {
1330  long iElec = (*st).n();
1331  long iVib = (*st).v();
1332  long iRot = (*st).J();
1333  double pop = states[ ipEnergySort[iElec][iVib][iRot] ].Pop();
1334  /* keep track of largest relative change in populations to
1335  * determines convergence */
1336  if(
1337  // >> chng 13 sep 21 -- pop test makes more sense as pre-condition
1338  pop>1e-6 * (*dense_total) &&
1339  /* >>chng 03 jul 19, this had simply been pop > SMALLFLOAT,
1340  * change to relative pops > 1e-15, spent too much time converging
1341  * levels at pops = 1e-37 */
1342  /* >>chng 03 dec 27, from rel pop 1e-15 to 1e-6 since converging heating will
1343  * be main convergence criteria check convergence */
1344  fabs( pop - H2_old_populations[iElec][iVib][iRot])
1345  /* on first call some very high J states can have zero pop ,
1346  * hence the SDIV, will retain sign for checks on oscilations,
1347  * hence the fabs */
1348  > fabs(PopChgMax_relative)*pop
1349  )
1350  {
1351  PopChgMax_relative =
1352  (pop - H2_old_populations[iElec][iVib][iRot])/SDIV(pop);
1353  iRotMaxChng_relative = iRot;
1354  iVibMaxChng_relative = iVib;
1355  popold_relative = H2_old_populations[iElec][iVib][iRot];
1356  popnew_relative = pop;
1357  }
1358  /* >>chng 05 feb 08, add largest rel change in total, this will be converged
1359  * down to higher accuracy than above
1360  * keep track of largest change in populations relative to total H2 to
1361  * determine convergence check convergence */
1362  const double rel_change = (pop - H2_old_populations[iElec][iVib][iRot])/SDIV(*dense_total);
1363  /* retain sign for checks on oscillations hence the fabs */
1364  if( fabs(rel_change) > fabs(PopChgMax_total) )
1365  {
1366  PopChgMax_total = rel_change;
1367  iRotMaxChng_total = iRot;
1368  iVibMaxChng_total = iVib;
1369  popold_total = H2_old_populations[iElec][iVib][iRot];
1370  popnew_total = pop;
1371  }
1372 
1373  kase = -1;
1374  /* update populations - we used the old populations to update the
1375  * current new populations - will do another iteration if they changed
1376  * by much. here old populations are updated for next sweep through molecule */
1377  /* pop oscillations have occurred - use small changes */
1378  /* >>chng 04 may 10, turn this back on - now with min on how small frac new
1379  * can become */
1380  const double abs_change = fabs( H2_old_populations[iElec][iVib][iRot] - pop );
1381 
1382  /* this branch very large changes, use mean of logs but onlly if both are positive*/
1383  if( abs_change > 3.*pop && H2_old_populations[iElec][iVib][iRot] * pop > 0 )
1384  {
1385  /* large changes or oscillations - take average in the log */
1386  H2_old_populations[iElec][iVib][iRot] = exp10(
1387  log10(H2_old_populations[iElec][iVib][iRot])/2. +
1388  log10(pop)/2. );
1389  kase = 2;
1390  }
1391 
1392  /* modest change, use means of old and new */
1393  else if( abs_change> 0.1*pop )
1394  {
1395  realnum frac_old=0.25f;
1396  /* large changes or oscillations - take average */
1397  H2_old_populations[iElec][iVib][iRot] =
1398  frac_old*H2_old_populations[iElec][iVib][iRot] +
1399  (1.f-frac_old)*pop;
1400  kase = 3;
1401  }
1402  else
1403  {
1404  /* small changes, use new value */
1405  H2_old_populations[iElec][iVib][iRot] = pop;
1406  kase = 4;
1407  }
1408  sumold += H2_old_populations[iElec][iVib][iRot];
1409  }
1410 
1411  /* will renormalize so that total population is correct */
1412  double H2_renorm_conserve_init = *dense_total/sumold;
1413 
1414  /* renormalize populations - were updated by renorm when routine entered,
1415  * before pops determined - but population determinations above do not have a sum rule on total
1416  * population - this renorm is to preserve total population */
1417  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1418  {
1419  long iElec = (*st).n();
1420  long iVib = (*st).v();
1421  long iRot = (*st).J();
1422  H2_old_populations[iElec][iVib][iRot] *= H2_renorm_conserve_init;
1423  }
1424  }
1425 
1426  /* get current ortho-para ratio, will be used as test on convergence */
1427  {
1428  ortho_density = 0.;
1429  para_density = 0.;
1430  H2_den_s = 0.;
1431  H2_den_g = 0.;
1432 
1433  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1434  {
1435  long iElec = (*st).n();
1436  long iVib = (*st).v();
1437  long iRot = (*st).J();
1438  const double& pop = (*st).Pop();
1439  /* find current population in H2s and H2g */
1440  if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
1441  {
1442  H2_den_s += pop;
1443  }
1444  else
1445  {
1446  H2_den_g += pop;
1447  }
1448  if( H2_lgOrtho[iElec][iVib][iRot] )
1449  {
1450  ortho_density += pop;
1451  }
1452  else
1453  {
1454  para_density += pop;
1455  }
1456  }
1457  ASSERT( fp_equal_tol( H2_den_s + H2_den_g, *dense_total, 1e-5 * (*dense_total) ) );
1458  }
1459 
1460  /* these will be used to determine whether solution has converged */
1461  ortho_para_older = ortho_para_old;
1462  ortho_para_old = ortho_para_current;
1463  ortho_para_current = ortho_density / SDIV( para_density );
1464 
1465  /* this will be evaluated in call to routine that follows - will check
1466  * whether this has converged */
1467  old_solomon_rate = Solomon_dissoc_rate_g;
1468 
1469  /* >>chng 05 jul 24, break code out into separate routine for clarify
1470  * located in mole_h2_etc.c - true says to only do Solomon rate */
1471  H2_Solomon_rate();
1472 
1473  /* are changes too large? must decide whether population shave converged,
1474  * will check whether populations themselves have changed by much,
1475  * but also change in heating by collisional deexcitation is stable */
1476  HeatChangeOld = HeatChange;
1477  HeatChange = HeatDexc_old - HeatDexc;
1478  {
1479  /* check whether pops are oscillating, as evidenced by change in
1480  * heating changing sign */
1481  if( loop_h2_pops>2 && (
1482  (HeatChangeOld*HeatChange<0. ) ||
1483  (PopChgMax_relative*PopChgMaxOld_relative<0. ) ) )
1484  {
1485  lgH2_pops_oscil = true;
1486  if( loop_h2_pops > 6 )
1487  {
1488  loop_h2_oscil = loop_h2_pops;
1489  lgH2_pops_ever_oscil = true;
1490  ++n_pop_oscil;
1491  }
1492  }
1493  else
1494  {
1495  lgH2_pops_oscil = false;
1496  /* turn off flag if no oscillations for a while */
1497  if( loop_h2_pops - loop_h2_oscil > 4 )
1498  {
1499  lgH2_pops_ever_oscil = false;
1500  }
1501  }
1502  }
1503 
1504  /* reevaluate heating - cooling */
1505  HeatDexc_old = HeatDexc;
1506  H2_Cooling();
1507 
1508  /* begin check on whether solution is converged */
1509  lgConv_h2_soln = true;
1510  lgPopsConv_total = true;
1511  lgPopsConv_relative = true;
1512  lgHeatConv = true;
1513  lgSolomonConv = true;
1514  lgOrthoParaRatioConv = true;
1515 
1516  /* these are all the convergence tests
1517  * check convergence converged */
1518  if( fabs(PopChgMax_relative)>converge_pops_relative )
1519  {
1520  /*lgPopsConv = (fabs(PopChgMax_relative)<=0.1);*/
1521  lgConv_h2_soln = false;
1522  lgPopsConv_relative = false;
1523  /* >>chng 04 sep 08, set quant_new to new chng max gs */
1524  /*quant_old = PopChgMax_relative;*/
1525  quant_old = PopChgMaxOld_relative;
1526  /*quant_new = 0.;*/
1527  quant_new = PopChgMax_relative;
1528 
1529  strcpy( chReason , "rel pops changed" );
1530  }
1531 
1532  /* check largest change in a level population relative to total h2
1533  * population convergence converged check */
1534  else if( fabs(PopChgMax_total)>converge_pops_total)
1535  {
1536  lgConv_h2_soln = false;
1537  lgPopsConv_total = false;
1538  /* >>chng 04 sep 08, set quant_new to new chng max gs */
1539  /*quant_old = PopChgMax_relative;*/
1540  quant_old = PopChgMaxOld_total;
1541  /*quant_new = 0.;*/
1542  quant_new = PopChgMax_total;
1543 
1544  strcpy( chReason , "tot pops changed" );
1545  }
1546 
1547  /* >>chng 04 apr 30, look at change in ortho-para ratio, also that is not
1548  * oscillating */
1549  /* >>chng 04 dec 15, only look at change, and don't make allowed change so tiny -
1550  * these were attempts at fixing problems that were due to shielding not thin*/
1551  else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> converge_ortho_para )
1552  /* else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> 1e-3
1553  && (ortho_para_current-ortho_para_old)*(ortho_para_old-ortho_para_older)>0. )*/
1554  {
1555  lgConv_h2_soln = false;
1556  lgOrthoParaRatioConv = false;
1557  quant_old = ortho_para_old;
1558  quant_new = ortho_para_current;
1559  strcpy( chReason , "ortho/para ratio changed" );
1560  }
1561  /* >>chng 04 dec 16, reduce error allowed fm /5 to /2, to be similar to
1562  * logic in conv_base */
1563  else if( !thermal.lgTemperatureConstant &&
1564  fabs(HeatDexc-HeatDexc_old)/MAX2(thermal.ctot,thermal.htot) >
1565  conv.HeatCoolRelErrorAllowed/10.
1566  /* >>chng 04 may 09, do not check on error in heating if constant temperature */
1567  /*&& !(thermal.lgTemperatureConstant || phycon.te <= phycon.TEMP_LIMIT_LOW )*/ )
1568  {
1569  /* default on HeatCoolRelErrorAllowed is 0.02 */
1570  /*lgHeatConv = (fabs(HeatDexc-HeatDexc_old)/thermal.ctot <=
1571  * conv.HeatCoolRelErrorAllowed/5.);*/
1572  lgConv_h2_soln = false;
1573  lgHeatConv = false;
1574  quant_old = HeatDexc_old/MAX2(thermal.ctot,thermal.htot);
1575  quant_new = HeatDexc/MAX2(thermal.ctot,thermal.htot);
1576  strcpy( chReason , "heating changed" );
1577  /*fprintf(ioQQQ,"DEBUG old new trip \t%.4e \t %.4e\n",
1578  HeatDexc_old,
1579  HeatDexc);*/
1580  }
1581 
1582  /* check on Solomon rate,
1583  * >>chng 04 aug 28, do not do this check if induced processes are disabled,
1584  * since Solomon process is then irrelevant */
1585  /* >>chng 04 sep 21, GS*/
1586  else if( rfield.lgInducProcess &&
1587  /* this is check that H2 abundance has not been set - if it has been
1588  * then we don't care what the Solomon rate is doing */
1589  hmi.H2_frac_abund_set==0 &&
1590  /*>>chng 05 feb 10, rather than checking change in Solomon relative to Solomon,
1591  * check it relative to total h2 destruction rate */
1592  fabs( Solomon_dissoc_rate_g - old_solomon_rate)/SDIV(hmi.H2_rate_destroy) >
1593  conv.EdenErrorAllowed/5.)
1594  {
1595  lgConv_h2_soln = false;
1596  lgSolomonConv = false;
1597  quant_old = old_solomon_rate;
1598  quant_new = Solomon_dissoc_rate_g;
1599  strcpy( chReason , "Solomon rate changed" );
1600  }
1601 
1602  /* did we pass all the convergence test */
1603  if( !lgConv_h2_soln )
1604  {
1605  /* this branch H2 populations within X are not converged,
1606  * print diagnostic */
1607 
1608  if( PRT_POPS || nTRACE >=n_trace_iterations )
1609  {
1610  /*fprintf(ioQQQ,"temppp\tnew\t%.4e\tnew\t%.4e\t%.4e\n",
1611  HeatDexc,
1612  HeatDexc_old,
1613  fabs(HeatDexc-HeatDexc_old)/thermal.ctot );*/
1614  fprintf(ioQQQ," %s loop %3li no conv oscl?%c why:%s ",
1615  label.c_str(),
1616  loop_h2_pops,
1617  TorF(lgH2_pops_ever_oscil),
1618  chReason );
1619  if( !lgPopsConv_relative )
1620  fprintf(ioQQQ," PopChgMax_relative:%.4e v:%li J:%li old:%.4e new:%.4e",
1621  PopChgMax_relative,
1622  iVibMaxChng_relative,
1623  iRotMaxChng_relative ,
1624  popold_relative ,
1625  popnew_relative );
1626  else if( !lgPopsConv_total )
1627  fprintf(ioQQQ," PopChgMax_total:%.4e v:%li J:%li old:%.4e new:%.4e",
1628  PopChgMax_total,
1629  iVibMaxChng_total,
1630  iRotMaxChng_total ,
1631  popold_total ,
1632  popnew_total );
1633  else if( !lgHeatConv )
1634  fprintf(ioQQQ," heat:%.4e old:%.4e new:%.4e",
1635  (HeatDexc-HeatDexc_old)/MAX2(thermal.ctot,thermal.htot),
1636  quant_old ,
1637  quant_new);
1638  /* Solomon rate changed */
1639  else if( !lgSolomonConv )
1640  fprintf(ioQQQ," d(sol rate)/tot dest\t%2e",(old_solomon_rate - Solomon_dissoc_rate_g)/SDIV(hmi.H2_rate_destroy));
1641  else if( !lgOrthoParaRatioConv )
1642  fprintf(ioQQQ," current, old, older ratios are %.4e %.4e %.4e",
1643  ortho_para_current , ortho_para_old, ortho_para_older );
1644  else
1645  TotalInsanity();
1646  fprintf(ioQQQ,"\n");
1647  }
1648  }
1649  /* end convergence criteria */
1650 
1651  if( trace.nTrConvg >= 5 )
1652  {
1653  fprintf( ioQQQ,
1654  " H2 5lev %li Conv?%c",
1655  loop_h2_pops ,
1656  TorF(lgConv_h2_soln) );
1657 
1658  if( fabs(PopChgMax_relative)>0.1 )
1659  fprintf(ioQQQ," pops, rel chng %.3e",PopChgMax_relative);
1660  else
1661  fprintf(ioQQQ," rel heat %.3e rel chng %.3e H2 heat/cool %.2e",
1662  HeatDexc/thermal.ctot ,
1663  fabs(HeatDexc-HeatDexc_old)/thermal.ctot ,
1664  HeatDexc/thermal.ctot);
1665 
1666  fprintf( ioQQQ,
1667  " Oscil?%c Ever Oscil?%c",
1668  TorF(lgH2_pops_oscil) ,
1669  TorF(lgH2_pops_ever_oscil) );
1670  fprintf(ioQQQ,"\n");
1671  }
1672 
1673  if( nTRACE >= n_trace_full )
1674  {
1675  fprintf(ioQQQ,
1676  "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",
1677  loop_h2_pops,
1678  kase,
1679  H2_renorm_chemistry,
1680  ortho_density / para_density ,
1681  frac_matrix);
1682 
1683  /* =======================INSIDE POPULATIONS CONVERGE LOOP =====================*/
1684  if( iVibMaxChng_relative>=0 && iRotMaxChng_relative>=0 && PopChgMax_relative>1e-10 )
1685  fprintf(ioQQQ,
1686  "end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n",
1687  loop_h2_pops,
1688  PopChgMax_relative ,
1689  H2_old_populations[0][iVibMaxChng_relative][iRotMaxChng_relative],
1690  states[ ipEnergySort[0][iVibMaxChng_relative][iRotMaxChng_relative] ].Pop(),
1691  iVibMaxChng_relative , iRotMaxChng_relative
1692  );
1693  }
1694  }
1695  /* =======================END POPULATIONS CONVERGE LOOP =====================*/
1696 
1697  /* >>chng 05 feb 08, do not print if we are in search phase */
1698  if( !lgConv_h2_soln && !conv.lgSearch )
1699  {
1700  conv.lgConvPops = false;
1701  lgPopsConverged = false;
1702  old_val = quant_old;
1703  new_val = quant_new;
1704  }
1705 
1706  for( qList::iterator st = states.begin(); st != states.end(); ++st )
1707  {
1708  ASSERT( (*st).Pop() >= 0. );
1709  }
1710 
1711  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1712  {
1713  /* following two heat exchange excitation, deexcitation */
1714  (*tr).Coll().cool() = 0.;
1715  (*tr).Coll().heat() = 0.;
1716 
1717  (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
1718 
1719  /* number of photons in the line
1720  * and line intensity */
1721  set_xIntensity( *tr );
1722  }
1723 
1724  average_energy_g = 0.;
1725  average_energy_s = 0.;
1726  /* determine average energy in ground and star */
1727  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1728  {
1729  double popTimesE = (*st).Pop() * (*st).energy().WN();
1730  if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
1731  average_energy_s += popTimesE;
1732  else
1733  average_energy_g += popTimesE;
1734  }
1735  /* average energy in ground and star */
1736  average_energy_g /= H2_den_g;
1737  if( H2_den_s > 1e-30 * (*dense_total) )
1738  average_energy_s /= H2_den_s;
1739  else
1740  average_energy_s = 0.;
1741 
1742  /* add up H2 + hnu => 2H, continuum photodissociation,
1743  * this is not the Solomon process, true continuum */
1744  /* >>chng 05 jun 16, GS, add dissociation to triplet states*/
1745  photodissoc_BigH2_H2s = 0.;
1746  photodissoc_BigH2_H2g = 0.;
1747  /* >>chng 05 jul 20, GS, add dissociation by H2 g and H2s*/
1748  Average_collH_dissoc_g = 0.;
1749  Average_collH_dissoc_s = 0.;
1750  Average_collH2_dissoc_g = 0.;
1751  Average_collH2_dissoc_s = 0.;
1752 
1753  rel_pop_LTE_g =0.;
1754  rel_pop_LTE_s = 0.;
1755 
1756  double exp_disoc = sexp(H2_DissocEnergies[0]/phycon.te_wn);
1757 
1758  /* >>chng 05 sep 12, TE, define a cutoff wavelength of 800 Angstrom
1759  * this is chosen as the cross sections given by
1760  *>>refer H2 photo cs Allison, A.C. & Dalgarno, A. 1969, Atomic Data, 1, 91
1761  * show a sharp decline in the cross section*/
1762  {
1763  static long ip_cut_off = -1;
1764  if( ip_cut_off < 0 )
1765  {
1766  /* one-time initialization of this pointer */
1767  ip_cut_off = ipoint( 1.14 );
1768  }
1769 
1770  /* >>chng 05 sep 12, TE, assume all H2s is at 2.5 eV
1771  * the dissociation threshold is at 1.07896 Rydberg*/
1772  double flux_accum_photodissoc_BigH2_H2s = 0;
1773  fixit("this 2.5 seems like a pretty bad (and unnecessary) approximation. Needs to be generalized at any rate.");
1774  long ip_H2_level = ipoint( 1.07896 - 2.5 / EVRYD);
1775  for( long i= ip_H2_level; i < ip_cut_off; ++i )
1776  {
1777  flux_accum_photodissoc_BigH2_H2s += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
1778  rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1] );
1779  }
1780 
1781  /* sum over all levels to obtain s and g populations and dissociation rates */
1782  for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
1783  {
1784  long iElec = (*st).n();
1785  if( iElec > 0 ) continue;
1786  long iVib = (*st).v();
1787  long iRot = (*st).J();
1788  const double &pop = (*st).Pop();
1789  fixit("generalize this factor (present value is (2m_e/m_H)^1.5/(2*2). See Robin's Feb 7, 2009 email.");
1790  const double mass_stat_factor = 3.634e-5/(2*2);
1791 
1792  /* >>chng 05 mar 22, TE, this should be for H2* rather than total */
1793  /* this is the total rate of direct photo-dissociation of excited electronic states into
1794  * the X continuum - this is continuum photodissociation, not the Solomon process */
1795  /* >>chng 03 sep 03, make sum of pops of excited states */
1796  if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
1797  {
1798  double arg_ratio;
1799  photodissoc_BigH2_H2s += pop * flux_accum_photodissoc_BigH2_H2s;
1800 
1801  /* >>chng 05 july 20, GS, collisional dissociation, unit s-1*/
1802  Average_collH_dissoc_s += pop * H2_coll_dissoc_rate_coef[iVib][iRot];
1803  Average_collH2_dissoc_s += pop * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
1804 
1805  /* >>chng 05 oct 17, GS, LTE populations of H2s*/
1806  arg_ratio = safe_div(exp_disoc,(*st).Boltzmann(),0.0);
1807  if( arg_ratio > 0. )
1808  {
1809  /* >>chng 05 oct 21, GS, only add ratio if Boltzmann factor > 0 */
1810  rel_pop_LTE_s += SAHA/SDIV(phycon.te32*arg_ratio)*
1811  (*st).g() * mass_stat_factor;
1812  }
1813  }
1814  else
1815  {
1816  double arg_ratio;
1817  /* >>chng 05 sep 12, TE, for H2g do the sum explicitly for every level*/
1818  double flux_accum_photodissoc_BigH2_H2g = 0;
1819  /* this is the dissociation energy needed for the level*/
1820  ip_H2_level = ipoint( 1.07896 - (*st).energy().Ryd() );
1821 
1822  for( long i= ip_H2_level; i < ip_cut_off; ++i )
1823  {
1824  flux_accum_photodissoc_BigH2_H2g += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
1825  rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1] );
1826  }
1827 
1828  photodissoc_BigH2_H2g += pop * flux_accum_photodissoc_BigH2_H2g;
1829 
1830  /* >>chng 05 jun 28, TE, determine average energy level in H2g */
1831  average_energy_g += (pop * (*st).energy().WN() );
1832 
1833  /* >>chng 05 july 20, GS, collisional dissociation, unit s-1*/
1834  Average_collH_dissoc_g += pop * H2_coll_dissoc_rate_coef[iVib][iRot];
1835  Average_collH2_dissoc_g += pop * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
1836 
1837  /* >>chng 05 oct 17, GS, LTE populations of H2g*/
1838  arg_ratio = safe_div(exp_disoc,(*st).Boltzmann(),0.0);
1839  if( arg_ratio > 0. )
1840  {
1841  rel_pop_LTE_g += SAHA/SDIV(phycon.te32*arg_ratio)*
1842  (*st).g() * mass_stat_factor;
1843  }
1844  }
1845  }
1846  }
1847 
1848  /* above sum was rate per unit vol since mult by H2 density, now div by H2* density to get rate s-1 */
1849  /* 0.25e-18 is wild guess of typical photodissociation cross section, from
1850  * >>refer H2 dissoc Allison, A.C. & Dalgarno, A. 1969, Atomic Data, 1, 91
1851  * this is based on an average of the highest v values they gave. unfortunately, we want
1852  * the highest J values -
1853  * final units are s-1*/
1854 
1855  if( H2_den_g > SMALLFLOAT )
1856  {
1857  Average_collH_dissoc_g /= SDIV(H2_den_g);/* unit cm3s-1*/
1858  Average_collH2_dissoc_g /= SDIV(H2_den_g);/* unit cm3s-1*/
1859  photodissoc_BigH2_H2g *= H2_DISS_ALLISON_DALGARNO / SDIV(H2_den_g);
1860  }
1861  else
1862  {
1863  Average_collH_dissoc_g = 0.;
1864  Average_collH2_dissoc_g = 0.;
1865  photodissoc_BigH2_H2g = 0.;
1866  }
1867  if( H2_den_s > SMALLFLOAT )
1868  {
1869  Average_collH_dissoc_s /= SDIV(H2_den_s);/* unit cm3s-1*/
1870  Average_collH2_dissoc_s /= SDIV(H2_den_s);/* unit cm3s-1*/
1871  photodissoc_BigH2_H2s *= H2_DISS_ALLISON_DALGARNO / SDIV(H2_den_s);
1872  }
1873  else
1874  {
1875  Average_collH_dissoc_s = 0.;
1876  Average_collH2_dissoc_s = 0.;
1877  photodissoc_BigH2_H2s = 0.;
1878  }
1879 
1880 
1881  // calculate some average rates from H2* to H2g
1882  H2_Calc_Average_Rates();
1883 
1884  if( nTRACE )
1885  {
1886  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",
1887  fnzone ,
1888  loop_h2_pops ,
1889  dens_rel_to_lim_react,
1890  old_solomon_rate,
1891  Solomon_dissoc_rate_g,
1892  Solomon_dissoc_rate_s,
1893  gs_rate(),
1894  photodissoc_BigH2_H2s );
1895  }
1896 
1897  /* >>chng 03 sep 01, add this population - before had just used H2star from chem network */
1898  /* if big H2 molecule is turned on and used for this zone, use its
1899  * value of H2* (pops of all states with v > 0 ) rather than simple network */
1900 
1901  /* update number of times we have been called */
1902  ++nCall_this_iteration;
1903 
1904  /* this will say how many times the large H2 molecule has been called in this zone -
1905  * if not called (due to low H2 abundance) then not need to update its line arrays */
1906  ++nCall_this_zone;
1907 
1908  /* >>chng 05 jun 21,
1909  * during search phase we want to use full matrix - save number of levels so that
1910  * we can restore it */
1911  nXLevelsMatrix = nXLevelsMatrix_save;
1912 
1913  /* >>chng 05 jan 19, check how many levels should be in the matrix if first call on
1914  * new zone, and we have a solution */
1915  /* end loop setting very first LTE populations */
1916  if( nCall_this_iteration && nzone != nzone_nlevel_set )
1917  {
1918  /* this is fraction of populations to include in matrix */
1919  const double FRAC = 0.99999;
1920  /* this loop is over increasing energy */
1921  double sum_pop = 0.;
1922  long nEner = 0;
1923  long iElec = 0;
1924  const bool PRT = false;
1925  if( PRT ) fprintf(ioQQQ,"DEBUG pops ");
1926  while( nEner < nLevels_per_elec[0] && sum_pop/(*dense_total) < FRAC )
1927  {
1928  /* array of energy sorted indices within X */
1929  ASSERT( iElec == ipElec_H2_energy_sort[nEner] );
1930  long iVib = ipVib_H2_energy_sort[nEner];
1931  long iRot = ipRot_H2_energy_sort[nEner];
1932  sum_pop += H2_old_populations[iElec][iVib][iRot];
1933  if( PRT ) fprintf(ioQQQ,"\t%.3e ", H2_old_populations[iElec][iVib][iRot]);
1934  ++nEner;
1935  }
1936  if( PRT ) fprintf(ioQQQ,"\n");
1937  nzone_nlevel_set = nzone;
1938  /*fprintf(ioQQQ,"DEBUG zone %.2f old nmatrix %li proposed nmatrix %li sum_pop %.4e H2_total %.4e\n",
1939  fnzone , nXLevelsMatrix ,nEner , sum_pop, *dense_total);
1940  nXLevelsMatrix = nEner;*/
1941  }
1942 
1943  return;
1944 }
1945 /*lint -e802 possible bad pointer */
1946 
1947 void diatomics::SolveExcitedElectronicLevels( void )
1948 {
1949  DEBUG_ENTRY( "diatomics::SolveExcitedElectronicLevels()" );
1950 
1951  multi_arr<double,3> rate_in;
1952  rate_in.alloc( H2_rad_rate_out.clone() );
1953  rate_in.zero();
1954  spon_diss_tot = 0.;
1955  double CosmicRayHILyaExcitationRate = ( hmi.lgLeidenCRHack ) ? secondaries.x12tot : 0.;
1956 
1957  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1958  {
1959  qList::iterator Lo = (*tr).Lo();
1960  long iElecLo = (*Lo).n();
1961  long iVibLo = (*Lo).v();
1962  long iRotLo = (*Lo).J();
1963  qList::iterator Hi = (*tr).Hi();
1964  long iElecHi = (*Hi).n();
1965  if( iElecHi < 1 ) continue;
1966  long iVibHi = (*Hi).v();
1967  long iRotHi = (*Hi).J();
1968  /* solve electronic excited state,
1969  * rate lower level in X goes to electronic excited state, s-1
1970  * first term is direct pump, second is cosmic ray excitation */
1971  /* collisional excitation of singlets by non-thermal electrons
1972  * this is stored ratio of electronic transition relative
1973  * cross section relative to the HI Lya cross section */
1974  double rate_up = (*tr).Emis().pump() + CosmicRayHILyaExcitationRate * (*tr).Coll().col_str();
1975  double rate_down =
1976  (*tr).Emis().Aul() * ( (*tr).Emis().Ploss() ) +
1977  rate_up * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
1978 
1979  /* this is a permitted electronic transition, must preserve nuclear spin */
1980  ASSERT( H2_lgOrtho[iElecHi][iVibHi][iRotHi] == H2_lgOrtho[iElecLo][iVibLo][iRotLo] );
1981 
1982  /* this is the rate [cm-3 s-1] electrons move into the upper level from X */
1983  rate_in[iElecHi][iVibHi][iRotHi] += H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_up;
1984 
1985  /* rate [s-1] from levels within X to electronic excited states,
1986  * includes photoexcitation and cosmic ray excitation */
1987  if( iElecLo==0 )
1988  H2_X_rate_to_elec_excited[iVibLo][iRotLo] += rate_up;
1989  H2_rad_rate_out[iElecLo][iVibLo][iRotLo] += rate_up;
1990 
1991  /* this is the rate [s-1] electrons leave the excited electronic upper level
1992  * and decay into X - will be used to get pops of electronic excited states */
1993  H2_rad_rate_out[iElecHi][iVibHi][iRotHi] += rate_down;
1994  ASSERT( rate_up >= 0. && rate_down >= 0. );
1995  }
1996 
1997  for( qList::iterator st = states.begin(); st != states.end(); ++st )
1998  {
1999  if( (*st).n() < 1 )
2000  continue;
2001 
2002  long iElec = (*st).n();
2003  long iVib = (*st).v();
2004  long iRot = (*st).J();
2005 
2006  H2_rad_rate_out[iElec][iVib][iRot] += H2_dissprob[iElec][iVib][iRot];
2007 
2008  /* update population [cm-3] of the electronic excited state this only includes
2009  * radiative processes between X and excited electronic states, and cosmic rays -
2010  * thermal collisions are neglected
2011  * X is done below and includes all processes */
2012  double pop = rate_in[iElec][iVib][iRot] / SDIV( H2_rad_rate_out[iElec][iVib][iRot] );
2013  (*st).Pop() = pop;
2014  spon_diss_tot += pop * H2_dissprob[iElec][iVib][iRot];
2015  if( H2_old_populations[iElec][iVib][iRot]==0. )
2016  H2_old_populations[iElec][iVib][iRot] = pop;
2017  /* this is total pop in this vibration state */
2018  pops_per_vib[iElec][iVib] += pop;
2019  /* total pop in each electronic state */
2020  pops_per_elec[iElec] += pop;
2021  }
2022 
2023  fixit("uncomment and test");
2024  if( H2_den_s > 1e-30 * (*dense_total) )
2025  spon_diss_tot /= H2_den_s;
2026  else
2027  spon_diss_tot = 0.;
2028 
2029  if(nTRACE >= n_trace_full)
2030  {
2031  for( long iElec=1; iElec<n_elec_states; ++iElec )
2032  {
2033  fprintf(ioQQQ," Pop(e=%li):",iElec);
2034  for( md2i it = pops_per_vib.begin(iElec); it != pops_per_vib.end(iElec); ++it )
2035  fprintf( ioQQQ,"\t%.2e", *it/(*dense_total) );
2036  fprintf(ioQQQ,"\n");
2037  }
2038  }
2039 
2040  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
2041  {
2042  qList::iterator Lo = (*tr).Lo();
2043  if( (*Lo).n() != 0 ) continue;
2044  qList::iterator Hi = (*tr).Hi();
2045  if( (*Hi).n() < 1 ) continue;
2046 
2047  /* radiative rates [cm-3 s-1] from electronic excited states to X */
2048  double rate = (*Hi).Pop() *
2049  ((*tr).Emis().Aul() * ( (*tr).Emis().Ploss() ) +
2050  (*tr).Emis().pump() * (*Lo).g() / (*Hi).g());
2051  H2_X_rate_from_elec_excited[(*Lo).v()][(*Lo).J()] += rate;
2052  H2_rad_rate_in[(*Lo).v()][(*Lo).J()] += rate;
2053  }
2054 
2055  return;
2056 }
2057 
2058 void diatomics::SolveSomeGroundElectronicLevels( void )
2059 {
2060  DEBUG_ENTRY( "diatomics::SolveSomeGroundElectronicLevels()" );
2061 
2062  /* these will be total rates into and out of the level */
2063  H2_col_rate_out.zero();
2064  H2_col_rate_in.zero();
2065 
2066  /* now evaluate total rates for all levels within X */
2067  for( long ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi)
2068  {
2069  /* array of energy sorted indices within X */
2070  ASSERT( ipElec_H2_energy_sort[ipHi] == 0 );
2071  long iVibHi = ipVib_H2_energy_sort[ipHi];
2072  long iRotHi = ipRot_H2_energy_sort[ipHi];
2073 
2074  realnum H2stat = states[ipHi].g();
2075  double EhiK = states[ipHi].energy().K();
2076 
2077  for( long ipLo=0; ipLo<ipHi; ++ipLo )
2078  {
2079  long iVibLo = ipVib_H2_energy_sort[ipLo];
2080  long iRotLo = ipRot_H2_energy_sort[ipLo];
2081 
2082  /* collision de-excitation [s-1] */
2083  realnum colldn = H2_X_coll_rate[ipHi][ipLo];
2084  /* inverse, rate up, [cm-3 s-1] */
2085  // use dsexp() explicitly rather than divide Boltzmann() factors
2086  // to avoid problems with underflow at low Te (see ticket #284)
2087  realnum collup = colldn *
2088  H2stat / states[ipLo].g() *
2089  dsexp( (EhiK - states[ipLo].energy().K()) / phycon.te );
2090 
2091  H2_col_rate_out[iVibHi][iRotHi] += colldn;
2092  H2_col_rate_in[iVibLo][iRotLo] += colldn * H2_old_populations[0][iVibHi][iRotHi];
2093 
2094  H2_col_rate_out[iVibLo][iRotLo] += collup;
2095  H2_col_rate_in[iVibHi][iRotHi] += collup * H2_old_populations[0][iVibLo][iRotLo];
2096  }
2097  }
2098 
2099  /* begin solving for X by back-substitution
2100  * this is the main loop that determines populations within X
2101  * units of all rates in are cm-3 s-1, all rates out are s-1
2102  * nLevels_per_elec is number of levels within electronic 0 - so nEner is one
2103  * beyond end of array here - but will be decremented at start of loop
2104  * this starts at the highest energy wihtin X and moves down to lower energies */
2105  long nEner = nLevels_per_elec[0];
2106  while( (--nEner) >= nXLevelsMatrix )
2107  {
2108  /* array of energy sorted indices within X - we are moving down
2109  * starting from highest level within X */
2110  long iElec = ipElec_H2_energy_sort[nEner];
2111  ASSERT( iElec == 0 );
2112  long iVib = ipVib_H2_energy_sort[nEner];
2113  long iRot = ipRot_H2_energy_sort[nEner];
2114 
2115  if( nEner+1 < nLevels_per_elec[0] )
2116  ASSERT( states[nEner].energy().WN() < states[nEner+1].energy().WN() ||
2117  fp_equal( states[nEner].energy().WN(), states[nEner+1].energy().WN() ) );
2118 
2119  /* >>chng 05 apr 30,GS, Instead of *dense_total, the specific populations are used because high levels have much less
2120  * populations than ground levels which consists most of the H2 population.
2121  * only do this if working level is not v=0, J=0, 1 */
2122  if( nEner >1 )
2123  {
2124  H2_col_rate_out[iVib][iRot] +=
2125  /* H2 grain interactions
2126  * rate (s-1) all v,J levels go to 0 or 1 preserving spin */
2127  (realnum)(rate_grain_op_conserve);
2128 
2129  /* this goes into v=0, and J=0 or 1 depending on whether initial
2130  * state is ortho or para */
2131  H2_col_rate_in[0][H2_lgOrtho[0][iVib][iRot]] +=
2132  /* H2 grain interactions
2133  * rate (cm-3 s-1) all v,J levels go to 0 or 1 preserving spin,
2134  * in above lgOrtho says whether should go to 0 or 1 */
2135  (realnum)(rate_grain_op_conserve*H2_old_populations[0][iVib][iRot]);
2136  }
2137  else if( nEner == 1 )
2138  {
2139  /* this is special J=1 to J=0 collision, which is only fast at
2140  * very low grain temperatures */
2141  H2_col_rate_out[0][1] +=
2142  /* H2 grain interactions
2143  * H2 ortho - para conversion on grain surface,
2144  * rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
2145  (realnum)(rate_grain_J1_to_J0);
2146 
2147  H2_col_rate_in[0][0] +=
2148  /* H2 grain interactions
2149  * H2 ortho - para conversion on grain surface,
2150  * rate (s-1) all v,J levels go to 0 or 1, preserving nuclear spin */
2151  (realnum)(rate_grain_J1_to_J0 *H2_old_populations[0][0][1]);
2152  }
2153 
2154  double pump_from_below = 0.;
2155  for( long ipLo = 0; ipLo<nEner; ++ipLo )
2156  {
2157  long iElecLo = ipElec_H2_energy_sort[ipLo];
2158  ASSERT( iElecLo == 0 );
2159  long iVibLo = ipVib_H2_energy_sort[ipLo];
2160  long iRotLo = ipRot_H2_energy_sort[ipLo];
2161  const TransitionList::iterator&tr = trans.begin() +ipTransitionSort[nEner][ipLo] ;
2162 
2163  /* the test on vibration is needed - the energies are ok but the space does not exist */
2164  if( ( abs(iRotLo-iRot) == 2 || iRotLo == iRot ) && (iVibLo <= iVib) && (*tr).ipCont() > 0 )
2165  {
2166  double rateone = (*tr).Emis().Aul() * ( (*tr).Emis().Ploss() );
2167  // Pumping and cosmic-ray excitation from these levels up to higher (than X) levels is already included in H2_rad_rate_out before reaching here.
2168  // The following lines take care of that process within X
2169  double CosmicRayHILyaExcitationRate = ( hmi.lgLeidenCRHack ) ? secondaries.x12tot : 0.;
2170  double pump_up = (*tr).Emis().pump() + CosmicRayHILyaExcitationRate * (*tr).Coll().col_str();
2171  pump_from_below += pump_up * H2_old_populations[iElecLo][iVibLo][iRotLo];
2172  rateone += pump_up * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
2173  ASSERT( rateone >=0 );
2174  H2_rad_rate_out[0][iVib][iRot] += rateone;
2175  H2_rad_rate_in[iVibLo][iRotLo] += rateone * H2_old_populations[iElec][iVib][iRot];
2176  }
2177  }
2178 
2179  /* we now have the total rates into and out of this level, get its population
2180  * units cm-3 */
2181  states[nEner].Pop() =
2182  (H2_col_rate_in[iVib][iRot]+ H2_rad_rate_in[iVib][iRot]+H2_X_source[nEner]+pump_from_below) /
2183  SDIV(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[0][iVib][iRot]+H2_X_sink[nEner]);
2184 
2185  ASSERT( states[nEner].Pop() >= 0. );
2186  }
2187 
2188  return;
2189 }
2190 
2191 /*H2_cooling evaluate cooling and heating due to H2 molecule */
2192 #if defined(__ICC) && defined(__i386)
2193 #pragma optimization_level 1
2194 #endif
2195 void diatomics::H2_Cooling( void )
2196 {
2197  DEBUG_ENTRY( "H2_Cooling()" );
2198 
2199  /* nCall_this_iteration is not incremented until after the level
2200  * populations have converged the first time. so for the first n calls
2201  * this will return zero, a good idea since populations will be wildly
2202  * incorrect during search for first valid pops */
2203  if( !lgEnabled || !nCall_this_iteration )
2204  {
2205  HeatDexc = 0.;
2206  HeatDiss = 0.;
2207  HeatDexc_deriv = 0.;
2208  return;
2209  }
2210 
2211  HeatDiss = 0.;
2212  /* heating due to dissociation of electronic excited states */
2213  for( qList::iterator st = states.begin(); st != states.end(); ++st )
2214  {
2215  long iElec = (*st).n();
2216  long iVib = (*st).v();
2217  long iRot = (*st).J();
2218  HeatDiss += (*st).Pop() * H2_dissprob[iElec][iVib][iRot] * H2_disske[iElec][iVib][iRot];
2219  }
2220 
2221  /* dissociation heating was in eV - convert to ergs */
2222  HeatDiss *= EN1EV;
2223 
2224  /* now work on collisional heating due to bound-bound
2225  * collisional transitions within X */
2226  HeatDexc = 0.;
2227  /* these are the colliders that will be considered as depopulating agents */
2228  /* the colliders are H, He, H2 ortho, H2 para, H+ */
2229  /* atomic hydrogen */
2230 
2231  /* this will be derivative */
2232  HeatDexc_deriv = 0.;
2233  long iElecHi = 0;
2234  for( long ipHi=1; ipHi<nLevels_per_elec[iElecHi]; ++ipHi )
2235  {
2236  realnum H2statHi = states[ipHi].g();
2237  double H2boltzHi = states[ipHi].Boltzmann();
2238  double H2popHi = states[ipHi].Pop();
2239  double ewnHi = states[ipHi].energy().WN();
2240 
2241  for( long ipLo=0; ipLo<ipHi; ++ipLo )
2242  {
2243  double rate_dn_heat = 0.;
2244 
2245  /* this sum is total downward heating summed over all colliders */
2246  mr3ci H2cr = CollRateCoeff.begin(ipHi, ipLo);
2247  for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
2248  /* downward collision rate */
2249  rate_dn_heat += H2cr[nColl]*collider_density[nColl];
2250 
2251  /* now get upward collisional cooling by detailed balance */
2252  double rate_up_cool = rate_dn_heat * states[ipLo].Pop() *
2253  /* rest converts into upward collision rate */
2254  H2statHi / states[ipLo].g() *
2255  safe_div(H2boltzHi, states[ipLo].Boltzmann(), 0.0 );
2256 
2257  rate_dn_heat *= H2popHi;
2258 
2259  /* net heating due to collisions within X -
2260  * positive if heating, negative is cooling
2261  * this will usually be heating if X is photo pumped
2262  * in printout and in save heating this is called "H2cX" */
2263  double conversion = (ewnHi - states[ipLo].energy().WN() ) * ERG1CM;
2264  double heatone = rate_dn_heat * conversion;
2265  double coolone = rate_up_cool * conversion;
2266  /* this is net heating, negative if cooling */
2267  double oneline = heatone - coolone;
2268  HeatDexc += oneline;
2269 
2270  /* derivative wrt temperature - assume exp wrt ground -
2271  * this needs to be divided by square of temperature in wn -
2272  * done at end of loop */
2273  HeatDexc_deriv += (realnum)(oneline * ewnHi);
2274 
2275  /* this would be a major logical error */
2276  ASSERT(
2277  (rate_up_cool==0 && rate_dn_heat==0) ||
2278  (states[ipHi].energy().WN() > states[ipLo].energy().WN()) );
2279  }/* end loop over lower levels, all collisions within X */
2280  }/* end loop over upper levels, all collisions within X */
2281 
2282  /* this is inside h2 cooling, and is called extra times when H2 heating is important */
2283  if( PRT_POPS )
2284  fprintf(ioQQQ,
2285  " DEBUG H2 heat fnzone\t%.2f\trenorm\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n",
2286  fnzone ,
2287  H2_renorm_chemistry ,
2288  phycon.te ,
2289  HeatDexc,
2290  HeatDexc/thermal.ctot);
2291 
2292  /* this is derivative of collisional heating wrt temperature - needs
2293  * to be divided by square of temperature in wn */
2294  HeatDexc_deriv /= (realnum)POW2(phycon.te_wn);
2295 
2296  if( nTRACE >= n_trace_full )
2297  fprintf(ioQQQ,
2298  " H2_Cooling Ctot\t%.4e\t HeatDiss \t%.4e\t HeatDexc \t%.4e\n" ,
2299  thermal.ctot ,
2300  HeatDiss ,
2301  HeatDexc );
2302 
2303  /* when we are very far from solution, during search phase, collisions within
2304  * X can be overwhelmingly large heating and cooling terms, which nearly
2305  * cancel out. Some dense cosmic ray heated clouds could not find correct
2306  * initial solution due to noise introduced by large net heating which was
2307  * the very noisy tiny difference between very large heating and cooling
2308  * terms. Do not include collisions with x as heat/cool during the
2309  * initial search phase */
2310  if( conv.lgSearch )
2311  {
2312  HeatDexc = 0.;
2313  HeatDexc_deriv = 0.;
2314  }
2315  return;
2316 }
2317 
2318 /*cdH2_colden return column density in H2, negative -1 if cannot find state,
2319  * header is cdDrive */
2320 double cdH2_colden( long iVib , long iRot )
2321 {
2322  diatomics& diatom = h2;
2323 
2324  /*if iVib is negative, return
2325  * total column density - iRot=0
2326  * ortho column density - iRot 1
2327  * para column density - iRot 2
2328  * else return column density in iVib, iRot */
2329  if( iVib < 0 )
2330  {
2331  if( iRot==0 )
2332  {
2333  /* return total column density */
2334  return( diatom.ortho_colden + diatom.para_colden );
2335  }
2336  else if( iRot==1 )
2337  {
2338  /* return ortho H2 column density */
2339  return diatom.ortho_colden;
2340  }
2341  else if( iRot==2 )
2342  {
2343  /* return para H2 column density */
2344  return diatom.para_colden;
2345  }
2346  else
2347  {
2348  fprintf(ioQQQ," iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n");
2349  return -1.;
2350  }
2351  }
2352  else if( diatom.lgEnabled )
2353  {
2354  return diatom.GetXColden( iVib, iRot );
2355  }
2356  /* error condition - no valid parameter */
2357  else
2358  return -1;
2359 }
2360 
2361 realnum diatomics::GetXColden( long iVib, long iRot )
2362 {
2363  DEBUG_ENTRY( "diatomics::GetXColden()" );
2364 
2365  /* this branch want state specific column density, which can only result from
2366  * evaluation of big molecule */
2367  int iElec = 0;
2368  if( iRot <0 || iVib >nVib_hi[iElec] || iRot > nRot_hi[iElec][iVib])
2369  {
2370  fprintf(ioQQQ," iVib and iRot must lie within X, returning -2.\n");
2371  fprintf(ioQQQ," iVib must be <= %li and iRot must be <= %li.\n",
2372  nVib_hi[iElec],nRot_hi[iElec][iVib]);
2373  return -2.;
2374  }
2375  else
2376  return H2_X_colden[iVib][iRot];
2377 }
2378 
2379 /*H2_Colden maintain H2 column densities within X */
2380 void diatomics::H2_Colden( const char *chLabel )
2381 {
2382  /* >>chng 05 jan 26, pops now set to LTE for small abundance case, so do this */
2383  if( !lgEnabled /*|| !nCall_this_zone*/ )
2384  return;
2385 
2386  DEBUG_ENTRY( "H2_Colden()" );
2387 
2388  if( strcmp(chLabel,"ZERO") == 0 )
2389  {
2390  /* zero out formation rates and column densites */
2391  H2_X_colden.zero();
2392  H2_X_colden_LTE.zero();
2393  }
2394 
2395  else if( strcmp(chLabel,"ADD ") == 0 )
2396  {
2397  /* add together column densities */
2398  for( qList::iterator st = states.begin(); st != states.end(); ++st )
2399  {
2400  long iElec = (*st).n();
2401  if( iElec > 0 ) continue;
2402  long iVib = (*st).v();
2403  long iRot = (*st).J();
2404  /* state specific H2 column density */
2405  H2_X_colden[iVib][iRot] += (realnum)( (*st).Pop() * radius.drad_x_fillfac);
2406  /* LTE state specific H2 column density - H2_populations_LTE is normed to unity
2407  * so must be multiplied by total H2 density */
2408  H2_X_colden_LTE[iVib][iRot] += (realnum)(H2_populations_LTE[0][iVib][iRot]*
2409  (*dense_total)*radius.drad_x_fillfac);
2410  }
2411  }
2412 
2413  /* we will not print column densities so skip that - if not print then we have a problem */
2414  else if( strcmp(chLabel,"PRIN") != 0 )
2415  {
2416  fprintf( ioQQQ, " H2_Colden does not understand the label %s\n",
2417  chLabel );
2418  cdEXIT(EXIT_FAILURE);
2419  }
2420 
2421  return;
2422 }
2423 
2424 /*H2_DR choose next zone thickness based on H2 big molecule */
2425 double diatomics::H2_DR(void)
2426 {
2427  return BIGFLOAT;
2428 }
2429 
2430 /*H2_RT_OTS - add H2 ots fields */
2431 void diatomics::H2_RT_OTS( void )
2432 {
2433  /* do not compute if H2 not turned on, or not computed for these conditions */
2434  if( !lgEnabled || !nCall_this_zone )
2435  return;
2436 
2437  DEBUG_ENTRY( "H2_RT_OTS()" );
2438 
2439  for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
2440  {
2441  qList::iterator Hi = (*tr).Hi();
2442  if( (*Hi).n() > 0 )
2443  continue;
2444 
2445  /* ots destruction rate */
2446  (*tr).Emis().ots() = (*(*tr).Hi()).Pop() * (*tr).Emis().Aul() * (*tr).Emis().Pdest();
2447 
2448  /* dump the ots rate into the stack - but only for ground electronic state*/
2449  RT_OTS_AddLine( (*tr).Emis().ots(), (*tr).ipCont() );
2450  }
2451 
2452  return;
2453 }
2454 
2455 void diatomics::H2_Calc_Average_Rates( void )
2456 {
2457  /* >>chng 05 jul 09, GS*/
2458  /* average Einstein value for H2* to H2g, GS*/
2459  double sumpop1 = 0.;
2460  double sumpopA1 = 0.;
2461  double sumpopcollH2O_deexcit = 0.;
2462  double sumpopcollH2p_deexcit = 0.;
2463  double sumpopcollH_deexcit = 0.;
2464  double popH2s = 0.;
2465  double sumpopcollH2O_excit = 0.;
2466  double sumpopcollH2p_excit = 0.;
2467  double sumpopcollH_excit = 0.;
2468  double popH2g = 0.;
2469 
2470  for( qList::const_iterator stHi = states.begin(); stHi != states.end(); ++stHi )
2471  {
2472  long iElecHi = (*stHi).n();
2473  if( iElecHi > 0 ) continue;
2474  long iVibHi = (*stHi).v();
2475  long iRotHi = (*stHi).J();
2476  double ewnHi = (*stHi).energy().WN();
2477  for( qList::const_iterator stLo = states.begin(); stLo != stHi; ++stLo )
2478  {
2479  long iVibLo = (*stLo).v();
2480  long iRotLo = (*stLo).J();
2481  double ewnLo2 = (*stLo).energy().WN();
2482  if( ewnHi > ENERGY_H2_STAR && ewnLo2 < ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
2483  {
2484  /* >>chng 05 jul 10, GS*/
2485  /* average collisional rate for H2* to H2g, GS*/
2486  if( H2_lgOrtho[0][iVibHi][iRotHi] == H2_lgOrtho[0][iVibLo][iRotLo] )
2487  {
2488  long ihi = ipEnergySort[0][iVibHi][iRotHi];
2489  long ilo = ipEnergySort[0][iVibLo][iRotLo];
2490  const TransitionList::iterator&tr =
2491  trans.begin()+ ipTransitionSort[ihi][ilo];
2492  double popHi = (*(*tr).Hi()).Pop();
2493  double popLo = (*(*tr).Lo()).Pop();
2494 
2495  /* sums of populations */
2496  popH2s += popHi;
2497  popH2g += popLo;
2498 
2499  /* sums of deexcitation rates - H2* to H2g */
2500  sumpopcollH_deexcit += popHi * CollRateCoeff[ihi][ilo][0];
2501  sumpopcollH2O_deexcit += popHi * CollRateCoeff[ihi][ilo][2];
2502  sumpopcollH2p_deexcit += popHi * CollRateCoeff[ihi][ilo][3];
2503 
2504  double temp = popLo *
2505  (*stHi).g() / (*stLo).g() *
2506  safe_div((*stHi).Boltzmann(), (*stLo).Boltzmann(), 0.0 );
2507 
2508  /* sums of excitation rates - H2g to H2* */
2509  sumpopcollH_excit += temp * CollRateCoeff[ihi][ilo][0];
2510  sumpopcollH2O_excit += temp * CollRateCoeff[ihi][ilo][2];
2511  sumpopcollH2p_excit += temp * CollRateCoeff[ihi][ilo][3];
2512 
2513  if( lgH2_radiative[ihi][ilo] )
2514  {
2515  sumpop1 += popHi;
2516  sumpopA1 += popHi * (*tr).Emis().Aul();
2517  }
2518  }
2519  }
2520  }
2521  }
2522  Average_A = sumpopA1/SDIV(sumpop1);
2523 
2524  /* collisional excitation and deexcitation of H2g and H2s */
2525  Average_collH2_deexcit = (sumpopcollH2O_deexcit+sumpopcollH2p_deexcit)/SDIV(popH2s);
2526  Average_collH2_excit = (sumpopcollH2O_excit+sumpopcollH2p_excit)/SDIV(popH2g);
2527  Average_collH_excit = sumpopcollH_excit/SDIV(popH2g);
2528  Average_collH_deexcit = sumpopcollH_deexcit/SDIV(popH2s);
2529 
2530  /*fprintf(ioQQQ,
2531  "DEBUG Average_collH_excit sumpop = %.2e %.2e %.2e %.2e %.2e %.2e \n",
2532  popH2g,popH2s,sumpopcollH_deexcit ,sumpopcollH_excit ,
2533  sumpopcollH_deexcit/SDIV(popH2s) ,sumpopcollH_excit/SDIV(popH2g));*/
2534  /*fprintf(ioQQQ,"sumpop = %le sumpopA = %le Av= %le\n",
2535  sumpop1,sumpopA1 , Average_A );*/
2536 
2537  return;
2538 }
2539 
2540 double diatomics::GetExcitedElecDensity( void )
2541 {
2542  double H2_sum_excit_elec_den = 0.;
2543  for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
2544  {
2545  if( iElecHi > 0 )
2546  H2_sum_excit_elec_den += pops_per_elec[iElecHi];
2547  }
2548 
2549  return H2_sum_excit_elec_den;
2550 }
2551 /*lint +e802 possible bad pointer */
2552 
diatomics::H2_rad_rate_in
multi_arr< double, 2 > H2_rad_rate_in
Definition: h2_priv.h:514
diatomics::nTRACE
int nTRACE
Definition: h2_priv.h:406
t_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:65
multi_arr::begin
iterator begin(size_type i1)
Definition: container_classes.h:1514
diatomics::H2_col_rate_out
multi_arr< double, 2 > H2_col_rate_out
Definition: h2_priv.h:513
Atom_LevelN
Definition: atoms.h:78
md2ci
multi_arr< double, 2 >::const_iterator md2ci
Definition: container_classes.h:1829
N_ELEC
const int N_ELEC
Definition: h2_priv.h:34
t_mole_local::sink_rate_tot
double sink_rate_tot(const char chSpecies[]) const
Definition: mole_reactions.cpp:4184
mole_global
t_mole_global mole_global
Definition: mole.cpp:7
diatomics::H2_dissprob
multi_arr< realnum, 3 > H2_dissprob
Definition: h2_priv.h:496
t_mole_global::lgStancil
bool lgStancil
Definition: mole.h:337
t_thermal::htot
double htot
Definition: thermal.h:169
diatomics::ENERGY_H2_STAR
const double ENERGY_H2_STAR
Definition: h2_priv.h:449
diatomics::Average_A
double Average_A
Definition: h2_priv.h:303
thermal
t_thermal thermal
Definition: thermal.cpp:6
diatomics::gs_rate
double gs_rate(void)
Definition: mole_h2_etc.cpp:111
diatomics::GetXColden
realnum GetXColden(long iVib, long iRot)
Definition: mole_h2.cpp:2361
diatomics::renorm_min
double renorm_min
Definition: h2_priv.h:343
diatomics::rel_pop_LTE_s
double rel_pop_LTE_s
Definition: h2_priv.h:290
diatomics::rel_pop_LTE_g
double rel_pop_LTE_g
Definition: h2_priv.h:289
diatomics::H2_DissocEnergies
double H2_DissocEnergies[N_ELEC]
Definition: h2_priv.h:473
diatomics::SolveExcitedElectronicLevels
void SolveExcitedElectronicLevels(void)
Definition: mole_h2.cpp:1947
diatomics::H2_LineZero
void H2_LineZero(void)
Definition: mole_h2.cpp:439
t_conv::EdenErrorAllowed
double EdenErrorAllowed
Definition: conv.h:263
exp10
double exp10(double x)
Definition: cddefines.h:1368
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
multi_arr::alloc
void alloc()
Definition: container_classes.h:1113
diatomics::Average_collH2_excit
double Average_collH2_excit
Definition: h2_priv.h:307
diatomics::n_elec_states
long int n_elec_states
Definition: h2_priv.h:416
t_rfield::flux
realnum ** flux
Definition: rfield.h:68
diatomics::HeatChangeOld
double HeatChangeOld
Definition: h2_priv.h:300
multi_arr::zero
void zero()
Definition: container_classes.h:1048
diatomics::mass_amu
realnum mass_amu
Definition: h2_priv.h:403
diatomics::spon_diss_tot
double spon_diss_tot
Definition: h2_priv.h:269
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:246
diatomics::H2_LevelPops
void H2_LevelPops(bool &lgPopsConverged, double &old_value, double &new_value)
Definition: mole_h2.cpp:905
t_hmi::lgLeiden_Keep_ipMH2s
bool lgLeiden_Keep_ipMH2s
Definition: hmi.h:219
diatomics::n_trace_full
int n_trace_full
Definition: h2_priv.h:409
diatomics::ortho_para_older
double ortho_para_older
Definition: h2_priv.h:339
t_conv::register_
ConvergenceCounter register_(const string name)
Definition: conv.cpp:87
diatomics::average_energy_s
double average_energy_s
Definition: h2_priv.h:294
t_rfield::outlin_noplot
realnum * outlin_noplot
Definition: rfield.h:189
diatomics::pops_per_vib
multi_arr< double, 2 > pops_per_vib
Definition: h2_priv.h:464
t_hmi::H2_rate_destroy
double H2_rate_destroy
Definition: hmi.h:30
diatomics::HeatDexc_deriv
double HeatDexc_deriv
Definition: h2_priv.h:299
set_xIntensity
void set_xIntensity(const TransitionProxy &t)
Definition: lines_service.cpp:811
diatomics::H2_coll_dissoc_rate_coef_H2
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef_H2
Definition: h2_priv.h:532
TorF
char TorF(bool l)
Definition: cddefines.h:749
diatomics::lgPrtMatrix
bool lgPrtMatrix
Definition: h2_priv.h:592
diatomics::sp_star
molecule * sp_star
Definition: h2_priv.h:428
diatomics::H2_RTMake
void H2_RTMake(linefunc line_one)
Definition: mole_h2.cpp:387
diatomics::ipVib_H2_energy_sort
valarray< long > ipVib_H2_energy_sort
Definition: h2_priv.h:548
diatomics::ortho_colden
double ortho_colden
Definition: h2_priv.h:335
conv
t_conv conv
Definition: conv.cpp:5
diatomics::H2_to_H_limit
double H2_to_H_limit
Definition: h2_priv.h:401
diatomics::AulEscp
double ** AulEscp
Definition: h2_priv.h:558
diatomics::lgEvaluated
bool lgEvaluated
Definition: h2_priv.h:317
phycon
t_phycon phycon
Definition: phycon.cpp:6
fp_equal_tol
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:904
diatomics::rate_grain_op_conserve
double rate_grain_op_conserve
Definition: h2_priv.h:280
mr3ci
multi_arr< realnum, 3 >::const_iterator mr3ci
Definition: container_classes.h:1820
t_conv::lgConvPops
bool lgConvPops
Definition: conv.h:136
TransitionList::begin
iterator begin(void)
Definition: transition.h:339
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:999
PRT_POPS
#define PRT_POPS
Definition: mole_h2.cpp:22
ipFineCont
long ipFineCont(double energy_ryd)
Definition: cont_ipoint.cpp:161
t_thermal::lgTemperatureConstant
bool lgTemperatureConstant
Definition: thermal.h:44
t_rfield::outlin
realnum ** outlin
Definition: rfield.h:189
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
diatomics::H2_den_g
double H2_den_g
Definition: h2_priv.h:543
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:831
diatomics::H2_coll_dissoc_rate_coef
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef
Definition: h2_priv.h:529
diatomics::HeatChange
double HeatChange
Definition: h2_priv.h:300
H2_DISS_ALLISON_DALGARNO
#define H2_DISS_ALLISON_DALGARNO
Definition: mole_h2.cpp:27
nzone
long int nzone
Definition: cddefines.cpp:14
t_mole_local::source_rate_tot
double source_rate_tot(const char chSpecies[]) const
Definition: mole_reactions.cpp:4320
diatomics::trans
TransitionList trans
Definition: h2_priv.h:430
MIN2
#define MIN2(a, b)
Definition: cddefines.h:803
collider_density
static realnum collider_density[N_X_COLLIDER]
Definition: mole_h2.cpp:48
diatomics::Average_collH2_dissoc_g
double Average_collH2_dissoc_g
Definition: h2_priv.h:312
dsexp
double dsexp(double x)
Definition: service.cpp:1038
t_rfield::flux_accum
realnum * flux_accum
Definition: rfield.h:77
diatomics::H2_X_sink_and_source
void H2_X_sink_and_source(void)
Definition: mole_h2.cpp:51
diatomics::Average_collH2_deexcit
double Average_collH2_deexcit
Definition: h2_priv.h:305
diatomics::n_trace_matrix
int n_trace_matrix
Definition: h2_priv.h:409
dense
t_dense dense
Definition: global.cpp:15
diatomics::H2_X_formation
multi_arr< realnum, 2 > H2_X_formation
Definition: h2_priv.h:517
diatomics::Solomon_dissoc_rate_g
double Solomon_dissoc_rate_g
Definition: h2_priv.h:271
diatomics::TeUsedColl
double TeUsedColl
Definition: h2_priv.h:423
diatomics::CollRateCoeff
multi_arr< realnum, 3 > CollRateCoeff
Definition: h2_priv.h:485
diatomics::dense_total
const double *const dense_total
Definition: h2_priv.h:453
collider_density_total_not_H2
static realnum collider_density_total_not_H2
Definition: mole_h2.cpp:49
diatomics::H2_X_sink
valarray< realnum > H2_X_sink
Definition: h2_priv.h:536
diatomics::levelAsEval
long int levelAsEval
Definition: h2_priv.h:564
diatomics::CollRate_levn
double ** CollRate_levn
Definition: h2_priv.h:558
diatomics::Cont_Dissoc_Rate
multi_arr< double, 3 > Cont_Dissoc_Rate
Definition: h2_priv.h:286
diatomics::ipEnergySort
multi_arr< long int, 3 > ipEnergySort
Definition: h2_priv.h:551
diatomics::H2_old_populations
multi_arr< double, 3 > H2_old_populations
Definition: h2_priv.h:501
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
diatomics::GetExcitedElecDensity
double GetExcitedElecDensity(void)
Definition: mole_h2.cpp:2540
diatomics::sp
molecule * sp
Definition: h2_priv.h:427
diatomics::Average_collH2_dissoc_s
double Average_collH2_dissoc_s
Definition: h2_priv.h:313
TransitionList::states
qList *& states()
Definition: transition.h:359
iteration
long int iteration
Definition: cddefines.cpp:16
diatomics::H2_CollidRateEvalAll
void H2_CollidRateEvalAll(void)
Definition: mole_h2_coll.cpp:13
diatomics::iterationAsEval
long int iterationAsEval
Definition: h2_priv.h:507
t_conv::lgSearch
bool lgSearch
Definition: conv.h:168
trace
t_trace trace
Definition: trace.cpp:5
multi_arr::clone
const multi_geom< d, ALLOC > & clone() const
Definition: container_classes.h:1779
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:554
diatomics::H2_disske
multi_arr< realnum, 3 > H2_disske
Definition: h2_priv.h:497
diatomics::ortho_density
double ortho_density
Definition: h2_priv.h:326
lines_service.h
diatomics::HeatDexc_old
double HeatDexc_old
Definition: h2_priv.h:298
diatomics::n_trace_iterations
int n_trace_iterations
Definition: h2_priv.h:409
diatomics::para_density_f
realnum para_density_f
Definition: h2_priv.h:331
diatomics::mole_H2_LTE
void mole_H2_LTE(void)
Definition: mole_h2_etc.cpp:229
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:15
cdH2_colden
double cdH2_colden(long iVib, long iRot)
Definition: mole_h2.cpp:2320
diatomics::H2_RT_diffuse
void H2_RT_diffuse(void)
Definition: mole_h2.cpp:368
diatomics::mole_H2_form
void mole_H2_form(void)
Definition: mole_h2_form.cpp:14
diatomics::para_density
double para_density
Definition: h2_priv.h:326
diatomics::H2_RT_tau_reset
void H2_RT_tau_reset(void)
Definition: mole_h2.cpp:455
t_hmi::H2_frac_abund_set
double H2_frac_abund_set
Definition: hmi.h:196
diatomics::ortho_para_old
double ortho_para_old
Definition: h2_priv.h:339
diatomics::lgLTE
bool lgLTE
Definition: h2_priv.h:376
qList::Boltzmann
double * Boltzmann()
Definition: quantumstate.h:149
diatomics::photodissoc_BigH2_H2s
double photodissoc_BigH2_H2s
Definition: h2_priv.h:264
POW2
#define POW2
Definition: cddefines.h:979
energy
double energy(const genericState &gs)
Definition: generic_state.cpp:49
multi_arr< double, 3 >
diatomics::nzoneAsEval
long int nzoneAsEval
Definition: h2_priv.h:507
diatomics::nLevels_per_elec
long int nLevels_per_elec[N_ELEC]
Definition: h2_priv.h:482
diatomics::destroy
vector< double > destroy
Definition: h2_priv.h:562
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:352
diatomics::H2_den_s
double H2_den_s
Definition: h2_priv.h:543
ipLineEnergy
long ipLineEnergy(double energy, const char *chLabel, long ipIonEnergy)
Definition: cont_ipoint.cpp:68
PressureRadiationLine
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
Definition: pressure.cpp:50
diatomics::nzoneEval
long int nzoneEval
Definition: h2_priv.h:395
diatomics::H2_populations_LTE
multi_arr< double, 3 > H2_populations_LTE
Definition: h2_priv.h:502
diatomics::label
string label
Definition: h2_priv.h:435
diatomics::nCall_this_zone
long int nCall_this_zone
Definition: h2_priv.h:348
diatomics::H2_RT_tau_inc
void H2_RT_tau_inc(void)
Definition: mole_h2.cpp:409
mole
t_mole_local mole
Definition: mole.cpp:8
linefunc
void(* linefunc)(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
Definition: h2_priv.h:10
rfield
t_rfield rfield
Definition: rfield.cpp:9
diatomics::Average_collH_deexcit
double Average_collH_deexcit
Definition: h2_priv.h:306
N_X_COLLIDER
const int N_X_COLLIDER
Definition: h2_priv.h:20
diatomics::lgFirst
bool lgFirst
Definition: h2_priv.h:565
diatomics::nCall_this_iteration
long int nCall_this_iteration
Definition: h2_priv.h:586
diatomics::H2_InterEnergy
double H2_InterEnergy(void)
diatomics::ndimMalloced
long int ndimMalloced
Definition: h2_priv.h:557
diatomics::H2_RT_OTS
void H2_RT_OTS(void)
Definition: mole_h2.cpp:2431
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:154
realnum
float realnum
Definition: cddefines.h:124
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:468
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
diatomics::H2_X_colden_LTE
multi_arr< realnum, 2 > H2_X_colden_LTE
Definition: h2_priv.h:523
diatomics::photodissoc_BigH2_H2g
double photodissoc_BigH2_H2g
Definition: h2_priv.h:265
t_conv::IonizErrorAllowed
realnum IonizErrorAllowed
Definition: conv.h:276
diatomics::lgH2_radiative
multi_arr< bool, 2 > lgH2_radiative
Definition: h2_priv.h:574
diatomics::nH2_pops
long int nH2_pops
Definition: h2_priv.h:577
BIGFLOAT
const realnum BIGFLOAT
Definition: cpu.h:244
diatomics::nXLevelsMatrix
long int nXLevelsMatrix
Definition: h2_priv.h:556
diatomics::H2_X_Hmin_back
multi_arr< realnum, 2 > H2_X_Hmin_back
Definition: h2_priv.h:519
diatomics::AulDest
double ** AulDest
Definition: h2_priv.h:558
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
diatomics::H2_Cooling
void H2_Cooling(void)
Definition: mole_h2.cpp:2195
molecule::index
int index
Definition: mole.h:194
diatomics::H2_rad_rate_out
multi_arr< double, 3 > H2_rad_rate_out
Definition: h2_priv.h:498
RT_OTS_AddLine
void RT_OTS_AddLine(double ots, long int ip)
Definition: rt_ots.cpp:390
h2
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
diatomics::lgColl_deexec_Calc
bool lgColl_deexec_Calc
Definition: h2_priv.h:366
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1011
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:376
diatomics::average_energy_g
double average_energy_g
Definition: h2_priv.h:293
diatomics::H2_DR
double H2_DR(void)
Definition: mole_h2.cpp:2425
diatomics::ipElec_H2_energy_sort
valarray< long > ipElec_H2_energy_sort
Definition: h2_priv.h:549
t_trace::nTrConvg
int nTrConvg
Definition: trace.h:27
t_conv::HeatCoolRelErrorAllowed
realnum HeatCoolRelErrorAllowed
Definition: conv.h:274
diatomics::H2_ContPoint
void H2_ContPoint(void)
Definition: mole_h2.cpp:276
RT_line_one_tauinc
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)
Definition: rt_line_one_tauinc.cpp:14
radius
t_radius radius
Definition: radius.cpp:5
diatomics::loop_h2_oscil
long int loop_h2_oscil
Definition: h2_priv.h:394
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:159
diatomics::ortho_density_f
realnum ortho_density_f
Definition: h2_priv.h:331
diatomics::rate_grain_J1_to_J0
double rate_grain_J1_to_J0
Definition: h2_priv.h:281
FRAC
#define FRAC
diatomics::HeatDexc
double HeatDexc
Definition: h2_priv.h:297
diatomics::H2_Calc_Average_Rates
void H2_Calc_Average_Rates(void)
Definition: mole_h2.cpp:2455
t_hmi::lgLeidenCRHack
bool lgLeidenCRHack
Definition: hmi.h:220
diatomics::ortho_para_current
double ortho_para_current
Definition: h2_priv.h:339
qList::end
iterator end()
Definition: quantumstate.h:373
diatomics::frac_matrix
double frac_matrix
Definition: h2_priv.h:419
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
diatomics::H2_Solomon_rate
void H2_Solomon_rate(void)
Definition: mole_h2_etc.cpp:24
diatomics::HeatDiss
double HeatDiss
Definition: h2_priv.h:296
diatomics::stat_levn
vector< double > stat_levn
Definition: h2_priv.h:562
diatomics::H2_X_rate_to_elec_excited
multi_arr< double, 2 > H2_X_rate_to_elec_excited
Definition: h2_priv.h:527
diatomics::H2_zero_pops_too_low
void H2_zero_pops_too_low(void)
Definition: mole_h2_etc.cpp:188
diatomics::H2_X_rate_from_elec_excited
multi_arr< double, 2 > H2_X_rate_from_elec_excited
Definition: h2_priv.h:525
diatomics::renorm_max
double renorm_max
Definition: h2_priv.h:343
t_radius::drad_x_fillfac
double drad_x_fillfac
Definition: radius.h:77
diatomics::states
qList states
Definition: h2_priv.h:429
diatomics::H2_Colden
void H2_Colden(const char *chLabel)
Definition: mole_h2.cpp:2380
molezone::den
double den
Definition: mole.h:421
diatomics::Average_collH_excit
double Average_collH_excit
Definition: h2_priv.h:308
diatomics::H2_RadPress
double H2_RadPress(void)
Definition: mole_h2.cpp:317
diatomics::H2_X_coll_rate_evaluate
void H2_X_coll_rate_evaluate(void)
Definition: mole_h2.cpp:198
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
diatomics::Solomon_dissoc_rate_s
double Solomon_dissoc_rate_s
Definition: h2_priv.h:272
qList::begin
iterator begin()
Definition: quantumstate.h:365
diatomics::H2_Level_low_matrix
void H2_Level_low_matrix(realnum abundance)
Definition: mole_h2.cpp:471
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:350
multi_arr::end
iterator end(size_type i1)
Definition: container_classes.h:1555
md2i
multi_arr< double, 2 >::iterator md2i
Definition: container_classes.h:1828
diatomics::ipTransitionSort
multi_arr< long int, 2 > ipTransitionSort
Definition: h2_priv.h:552
mole_update_species_cache
void mole_update_species_cache(void)
Definition: mole_species.cpp:866
t_hmi::H2_total
double H2_total
Definition: hmi.h:25
hd
diatomics hd("hd", 4100.,&hmi.HD_total, Yan_H2_CS)
t_dense::eden
double eden
Definition: dense.h:201
diatomics::pops
vector< double > pops
Definition: h2_priv.h:562
t_hmi::H2_total_f
realnum H2_total_f
Definition: hmi.h:26
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
diatomics::rad_end
TransitionList::iterator rad_end
Definition: h2_priv.h:431
t_rfield::lgInducProcess
bool lgInducProcess
Definition: rfield.h:233
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
diatomics::nRot_hi
valarray< long > nRot_hi[N_ELEC]
Definition: h2_priv.h:477
diatomics::pops_per_elec
double pops_per_elec[N_ELEC]
Definition: h2_priv.h:484
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
diatomics::para_colden
double para_colden
Definition: h2_priv.h:335
diatomics::H2_renorm_chemistry
double H2_renorm_chemistry
Definition: h2_priv.h:467
diatomics::create
vector< double > create
Definition: h2_priv.h:562
diatomics::Average_collH_dissoc_g
double Average_collH_dissoc_g
Definition: h2_priv.h:310
diatomics::ipRot_H2_energy_sort
valarray< long > ipRot_H2_energy_sort
Definition: h2_priv.h:550
diatomics::AulPump
double ** AulPump
Definition: h2_priv.h:558
diatomics::H2_X_source
valarray< realnum > H2_X_source
Definition: h2_priv.h:535
diatomics::H2_Accel
double H2_Accel(void)
Definition: mole_h2.cpp:294
diatomics::nzone_nlevel_set
long int nzone_nlevel_set
Definition: h2_priv.h:581
t_phycon::te_wn
double te_wn
Definition: phycon.h:30
diatomics::Average_collH_dissoc_s
double Average_collH_dissoc_s
Definition: h2_priv.h:311
diatomics::CalcPhotoionizationRate
void CalcPhotoionizationRate(void)
Definition: mole_h2_etc.cpp:411
fixit
#define fixit(a)
Definition: cddefines.h:417
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
hmi
t_hmi hmi
Definition: hmi.cpp:5
fnzone
double fnzone
Definition: cddefines.cpp:15
LIM_H2_POP_LOOP
#define LIM_H2_POP_LOOP
Definition: mole_h2.cpp:24
t_phycon::te
double te
Definition: phycon.h:21
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:349
RT_line_one_tau_reset
void RT_line_one_tau_reset(const TransitionProxy &t)
Definition: rt_line_one_tau_reset.cpp:13
diatomics::depart
vector< double > depart
Definition: h2_priv.h:562
diatomics::excit
vector< double > excit
Definition: h2_priv.h:562
diatomics::H2_X_coll_rate
multi_arr< realnum, 2 > H2_X_coll_rate
Definition: h2_priv.h:470
diatomics::H2_itrzn
double H2_itrzn(void)
Definition: mole_h2.cpp:263
diatomics::nVib_hi
long int nVib_hi[N_ELEC]
Definition: h2_priv.h:475
diatomics::SolveSomeGroundElectronicLevels
void SolveSomeGroundElectronicLevels(void)
Definition: mole_h2.cpp:2058
t_phycon::te32
double te32
Definition: phycon.h:58
diatomics::H2_col_rate_in
multi_arr< double, 2 > H2_col_rate_in
Definition: h2_priv.h:512
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
diatomics::n_trace_final
int n_trace_final
Definition: h2_priv.h:409
diatomics::H2_X_colden
multi_arr< realnum, 2 > H2_X_colden
Definition: h2_priv.h:521
t_thermal::ctot
double ctot
Definition: thermal.h:130
diatomics::nH2_zone
long int nH2_zone
Definition: h2_priv.h:578
diatomics::H2_ipPhoto
multi_arr< int, 2 > H2_ipPhoto
Definition: h2_priv.h:511
diatomics::H2_lgOrtho
multi_arr< bool, 3 > H2_lgOrtho
Definition: h2_priv.h:504
Generated by   doxygen 1.8.5