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