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