cloudy trunk
Loading...
Searching...
No Matches
mole_h2_create.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_Create create variables for the H2 molecule, called by ContCreatePointers after continuum
4 * mesh has been set up */
5#include "cddefines.h"
6#include "collision.h"
7#include "physconst.h"
8#include "taulines.h"
9#include "lines_service.h"
10#include "opacity.h"
11#include "hmi.h"
12#include "ipoint.h"
13#include "grainvar.h"
14#include "h2.h"
15#include "h2_priv.h"
16#include "mole.h"
17#include "dense.h"
18
19/* this integer is added to rotation quantum number J for the test of whether
20 * a particular J state is ortho or para - the state is ortho if J+below is odd,
21 * and para if J+below is even */
22int H2_nRot_add_ortho_para[N_ELEC] = {0 , 1 , 1 , 0, 1, 1 , 0};
23
24/* if this is set true then code will print energies and stop */
25enum {DEBUG_ENER=false};
26
27/* this is equation 8 of Takahashi 2001, clearer definition is given in
28 * equation 5 and following discussion of
29 * >>refer H2 formation Takahashi, J., & Uehara, H., 2001, ApJ, 561, 843-857
30 * 0.27eV, convert into wavenumbers */
31static double XVIB[H2_TOP] = { 0.70 , 0.60 , 0.20 };
32static double Xdust[H2_TOP] = { 0.04 , 0.10 , 0.40 };
33
34/* this is energy difference between bottom of potential well and 0,0
35 * the Takahashi energy scale is from the bottom,
36 * 2201.9 wavenumbers */
37static const double energy_off = 0.273*FREQ_1EV/SPEEDLIGHT;
38
39STATIC double EH2_eval( int ipH2, double DissocEnergy, double energy_wn )
40{
41 double EH2_here;
42 double Evm = DissocEnergy * XVIB[ipH2] + energy_off;
43
44 double Ev = (energy_wn+energy_off);
45 /* equation 9 of Takahashi 2001 which is only an approximation
46 * equation 1, 2 of
47 * Takahashi, Junko, & Uehara, Hideya, 2001, ApJ, 561, 843-857,
48 * this is heat deposited on grain by H2 formation in this state */
49 double Edust = DissocEnergy * Xdust[ipH2] *
50 ( 1. - ( (Ev - Evm) / (DissocEnergy+energy_off-Evm)) *
51 ( (1.-Xdust[ipH2])/2.) );
52 ASSERT( Edust >= 0. );
53
54 /* energy is total binding energy less energy lost on grain surface
55 * and energy offset */
56 EH2_here = DissocEnergy +energy_off - Edust;
57 ASSERT( EH2_here >= 0.);
58
59 return EH2_here;
60}
61
62/*H2_vib_dist evaluates the vibration distribution for H2 formed on grains */
63STATIC double H2_vib_dist( int ipH2 , double EH2, double DissocEnergy, double energy_wn)
64{
65 double G1[H2_TOP] = { 0.3 , 0.4 , 0.9 };
66 double G2[H2_TOP] = { 0.6 , 0.6 , 0.4 };
67 double Evm = DissocEnergy * XVIB[ipH2] + energy_off;
68 double Fv;
69 if( (energy_wn+energy_off) <= Evm )
70 {
71 /* equation 4 of Takahashi 2001 */
72 Fv = sexp( POW2( (energy_wn+energy_off - Evm)/(G1[ipH2]* Evm ) ) );
73 }
74 else
75 {
76 /* equation 5 of Takahashi 2001 */
77 Fv = sexp( POW2( (energy_wn+energy_off - Evm)/(G2[ipH2]*(EH2 - Evm ) ) ) );
78 }
79 return Fv;
80}
81
83{
84 //quantumState_diatoms *Hi = static_cast<quantumState_diatoms*>( tr.Hi );
85 qList::iterator Lo = (*tr).Lo() ;
86 //if( Lo->n==0 && lgH2_radiative[ ipEnergySort[(*Hi).n][(*Hi).v][(*Hi).J] ][ ipEnergySort[Lo->n][Lo->v][Lo->J] ] )
87 if( (*Lo).n()==0 && (*tr).Emis().Aul() > 1.01e-30 )
88 return true;
89 else
90 return false;
91}
92
94{
95 if( lgRadiative( tr1 ) && !lgRadiative( tr2 ) )
96 return true;
97 else
98 return false;
99}
100
102{
103 if( st1.energy().Ryd() <= st2.energy().Ryd() )
104 return true;
105 else
106 return false;
107}
108
109/* create variables for the H2 molecule, called by
110 * ContCreatePointers after continuum mesh has been set up */
112{
113 /* this is flag set above - when true h2 code is not executed - this is way to
114 * avoid this code when it is not working */
115 /* only malloc vectors one time per core load */
116 if( lgREAD_DATA || !lgEnabled )
117 return;
118
119 DEBUG_ENTRY( "H2_Create()" );
120
121 bool lgDebugPrt = false;
122
123 /* print string if H2 debugging is enabled */
124 if( lgDebugPrt )
125 fprintf(ioQQQ," H2_Create called in DEBUG mode.\n");
126
127 // find the species in the chemistry network
128 sp = findspecies( label.c_str() );
129 ASSERT( sp != null_mole );
130
131 // find the same species with the excitation marker
132 string strSpStar = shortlabel + "*";
133 sp_star = findspecies( strSpStar.c_str() );
134 ASSERT( sp != null_mole );
135
136 mass_amu = sp->mole_mass / ATOMIC_MASS_UNIT;
137 ASSERT( mass_amu > 0. );
138
140
141 /* This does more than just read energies... it actually defines the states */
143
144 // now sort states into energy order
145 fixit(); // this doesn't work!
146 //sort( states.begin(), states.end(), compareEnergies );
147
148 ASSERT( n_elec_states > 0 );
149 /* create a vector of sorted energies */
151 for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
152 {
153 ipEnergySort.reserve(iElecHi,nVib_hi[iElecHi]+1);
154 for( long iVibHi = 0; iVibHi <= nVib_hi[iElecHi]; ++iVibHi )
155 {
156 ipEnergySort.reserve(iElecHi,iVibHi,nRot_hi[iElecHi][iVibHi]+1);
157 }
158 }
159 ipEnergySort.alloc();
160
161 /* create arrays for energy sorted referencing of e, v, J */
162 ipElec_H2_energy_sort.resize( states.size() );
163 ipVib_H2_energy_sort.resize( states.size() );
164 ipRot_H2_energy_sort.resize( states.size() );
165 for( unsigned nEner = 0; nEner < states.size(); ++nEner )
166 {
167 long iElec = states[nEner].n();
168 long iVib = states[nEner].v();
169 long iRot = states[nEner].J();
170 ipElec_H2_energy_sort[nEner] = iElec;
171 ipVib_H2_energy_sort[nEner] = iVib;
172 ipRot_H2_energy_sort[nEner] = iRot;
173 ipEnergySort[iElec][iVib][iRot] = nEner;
174 /* following will print quantum indices and energies */
175 /*fprintf(ioQQQ,"%li\t%li\t%li\t%.3e\n", iElec, iVib, iRot,
176 states[nEner].energy().WN() );*/
177 }
178
179 H2_stat.alloc( ipEnergySort.clone() );
180 H2_lgOrtho.alloc( ipEnergySort.clone() );
181
182 /* set all statistical weights - ours is total statistical weight -
183 * including nuclear spin */
184 for( qList::iterator st = states.begin(); st != states.end(); ++st )
185 {
186 /* unlike atoms, for H2 nuclear spin is taken into account - so the
187 * statistical weight of even and odd J states differ by factor of 3 - see page 166, sec par
188 * >>>refer H2 H2_stat wght Shull, J.M., & Beckwith, S., 1982, ARAA, 20, 163-188 */
189 if( is_odd( (*st).J() + H2_nRot_add_ortho_para[(*st).n()]) )
190 {
191 /* ortho */
192 H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] = true;
193 H2_stat[(*st).n()][(*st).v()][(*st).J()] = 3.f*(2.f*(*st).J()+1.f);
194 }
195 else
196 {
197 /* para */
198 H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] = false;
199 H2_stat[(*st).n()][(*st).v()][(*st).J()] = (2.f*(*st).J()+1.f);
200 }
201 (*st).g() = H2_stat[(*st).n()][(*st).v()][(*st).J()];
202 }
203
204 if( lgDebugPrt )
205 {
206 for( long iElec=0; iElec<n_elec_states; ++iElec)
207 {
208 /* print the number of levels within iElec */
209 fprintf(ioQQQ,"\t(%li %li)", iElec , nLevels_per_elec[iElec] );
210 }
211 fprintf(ioQQQ,
212 " H2_Create: there are %li electronic levels, in each level there are",
214 fprintf(ioQQQ,
215 " for a total of %li levels.\n", (long int) states.size() );
216 }
217
218 /* now find number of levels in H2g */
219 for( long nEner=0; nEner<nLevels_per_elec[0]; ++nEner )
220 {
221 if( states[nEner].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
222 break;
223 nEner_H2_ground = nEner;
224 }
225 /* need to increment it so that this is the number of levels, not the index
226 * of the highest level */
228
229 /* this is the number of levels to do with the matrix - set with the
230 * atom h2 matrix command, keyword ALL means to do all of X in the matrix
231 * but number of levels within X was not known when the command was parsed,
232 * so this was set to -1 to defer setting to all until now */
233 if( nXLevelsMatrix<0 )
234 {
236 }
237 else if( nXLevelsMatrix > nLevels_per_elec[0] )
238 {
239 fprintf( ioQQQ,
240 " The total number of levels used in the matrix solver was set to %li but there are only %li levels in X.\n Sorry.\n",
244 }
245
246 /* at this stage the full electronic, vibration, and rotation energies have been defined,
247 * this is an option to print the energies */
248 {
249 /* set following true to get printout, false to not print energies */
250 if( DEBUG_ENER )
251 {
252 /* print title for quantum numbers and energies */
253 /*fprintf(ioQQQ,"elec\tvib\trot\tenergy\n");*/
254 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
255 {
256 fprintf(ioQQQ,"%li\t%li\t%li\t%.5e\n", (*st).n(), (*st).v(), (*st).J(), (*st).energy().WN() );
257 }
258 /* this will exit the program after printing the level energies */
260 }
261 }
262
263 /* this will prevent data from being read twice */
264 lgREAD_DATA = true;
265
266 /* create space for the electronic levels */
269 H2_dissprob.reserve(n_elec_states);
270
271 for( long iElec = 0; iElec<n_elec_states; ++iElec )
272 {
273
274 if( lgDebugPrt )
275 fprintf(ioQQQ,"elec %li highest vib= %li\n", iElec , nVib_hi[iElec] );
276
277 ASSERT( nVib_hi[iElec] > 0 );
278
279 /* nVib_hi is now the highest vibrational level before dissociation,
280 * now allocate space to hold the number of rotation levels */
281 H2_populations_LTE.reserve(iElec,nVib_hi[iElec]+1);
282 pops_per_vib.reserve(iElec,nVib_hi[iElec]+1);
283
284 /* now loop over all vibrational levels, and find out how many rotation levels there are */
285 /* ground is special since use tabulated data - there are 14 vib states,
286 * ivib=14 is highest */
287 for( long iVib = 0; iVib <= nVib_hi[iElec]; ++iVib )
288 {
289 /* lastly create the space for the rotation quantum number */
290 H2_populations_LTE.reserve(iElec,iVib,nRot_hi[iElec][iVib]+1);
291 }
292 }
293
294 H2_populations_LTE.alloc();
295 H2_old_populations.alloc( H2_populations_LTE.clone() );
296 H2_Boltzmann.alloc( H2_populations_LTE.clone() );
297 H2_rad_rate_out.alloc( H2_populations_LTE.clone() );
298 H2_dissprob.alloc( H2_populations_LTE.clone() );
299 H2_disske.alloc( H2_populations_LTE.clone() );
300 // this has a different geometry than the above
301 pops_per_vib.alloc();
302
303 H2_dissprob.zero();
304 H2_disske.zero();
305
306 /* set this one time, will never be set again, but might be printed */
307 H2_rad_rate_out.zero();
308
309 /* these do not have electronic levels - all within X */
310 H2_ipPhoto.reserve(nVib_hi[0]+1);
311
312 /* space for the vibration levels */
313 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
314 {
315 /* space for the rotation quantum number */
316 H2_ipPhoto.reserve(iVib,nRot_hi[0][iVib]+1);
317 }
318
319 H2_ipPhoto.alloc();
320 H2_col_rate_in.alloc( H2_ipPhoto.clone() );
321 H2_col_rate_out.alloc( H2_ipPhoto.clone() );
322 H2_rad_rate_in.alloc( H2_ipPhoto.clone() );
323 H2_coll_dissoc_rate_coef.alloc( H2_ipPhoto.clone() );
325 H2_X_colden.alloc( H2_ipPhoto.clone() );
327 H2_X_rate_to_elec_excited.alloc( H2_ipPhoto.clone() );
328 H2_X_colden_LTE.alloc( H2_ipPhoto.clone() );
329 H2_X_formation.alloc( H2_ipPhoto.clone() );
330 H2_X_Hmin_back.alloc( H2_ipPhoto.clone() );
331
332 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
333 {
334 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
335 {
336 /* >>chng 04 jun 14, set these to bad numbers */
337 H2_rad_rate_in[iVib][iRot] = -BIGFLOAT;
338 H2_coll_dissoc_rate_coef[iVib][iRot] = -BIGFLOAT;
340 }
341 }
342 /* zero out the matrices */
343 H2_X_colden.zero();
344 H2_X_colden_LTE.zero();
345 H2_X_formation.zero();
346 H2_X_Hmin_back.zero();
347 /* rates [cm-3 s-1] from elec excited states into X only vib and rot */
349 /* rates [s-1] to elec excited states from X only vib and rot */
351
352 /* distribution function for populations following formation from H minus H- */
354 for( long i=0; i<nTE_HMINUS; ++i )
355 {
357 /* space for the vibration levels */
358 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
359 {
360 H2_X_hminus_formation_distribution.reserve(i,iVib,nRot_hi[0][iVib]+1);
361 }
362 }
366
367 /* >>chng 05 jun 20, do not use this, which is highly processed - use ab initio
368 * rates of excitation to electronic levels instead */
369 /* read in cosmic ray distribution information
370 H2_Read_Cosmicray_distribution(); */
371
372 /* grain formation matrix */
374 for( long ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
375 {
377
378 /* space for the vibration levels */
379 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
380 {
381 H2_X_grain_formation_distribution.reserve(ipH2,iVib,nRot_hi[0][iVib]+1);
382 }
383 }
386
387 for( long iElec=0; iElec<n_elec_states; ++iElec )
388 {
389 /* get dissociation probabilities and energies - ground state is stable */
390 if( iElec > 0 )
391 H2_ReadDissprob(iElec);
392 }
393
394 /* >>02 oct 18, add photodissociation, H2 + hnu => 2H + KE */
395 /* we now have ro-vib energies, now set up threshold array offsets
396 * for photodissociation */
397 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
398 {
399 long iElec = (*st).n();
400 if( iElec > 0 ) continue;
401 long iVib = (*st).v();
402 long iRot = (*st).J();
403 /* this is energy needed to get up to n=3 electronic continuum
404 * H2 cannot dissociate following absorption of a continuum photon into the
405 * continuum above X (which would require little energy)
406 * because that process violates momentum conservation
407 * these would be the triplet states - permitted are into singlets
408 * the effective full wavelength range of this process is from Lya to
409 * Lyman limit in shielded regions
410 * tests show limits are between 850A and 1220A - so Lya is included */
412 /*>>KEYWORD Allison & Dalgarno; continuum dissociation; */
413 double thresh = (H2_DissocEnergies[1] - (*st).energy().WN())*WAVNRYD;
414 /*fprintf(ioQQQ,"DEBUG\t%.2f\t%f\n", RYDLAM/thresh , thresh);*/
415 /* in theory we should be able to assert that thesh just barely reaches
416 * lya, but actual numbers reach down to 0.749 ryd */
417 ASSERT( thresh > 0.74 );
418 H2_ipPhoto[iVib][iRot] = ipoint(thresh);
419 fixit(); // this needs to be generalized
420 }
421
422 CollRateCoeff.reserve( nLevels_per_elec[0] );
423 for( long j = 1; j < nLevels_per_elec[0]; ++j )
424 {
425 CollRateCoeff.reserve( j, j );
426 for( long k = 0; k < j; ++k )
427 {
428 CollRateCoeff.reserve( j, k, N_X_COLLIDER );
429 }
430 }
431 CollRateCoeff.alloc();
432 CollRateCoeff.zero();
433
434 fixit(); // Does RateCoefTable only need to be N_X_COLLIDER long?
436 /* now read in the various sets of collision data */
437 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
438 {
439 /* ground state has tabulated data */
440 H2_CollidRateRead(nColl);
441 }
442
443 CollRateErrFac.alloc( CollRateCoeff.clone() );
444 CollRateErrFac.zero();
445
446 /* option to add gaussian random mole */
447 if( lgH2_NOISE )
448 {
449 for( long ipHi = 1; ipHi < nLevels_per_elec[0]; ++ipHi )
450 {
451 for( long ipLo = 0; ipLo < ipHi; ++ipLo )
452 {
453 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
454 {
455 /* this returns the log of the random noise */
457 CollRateErrFac[ipHi][ipLo][nColl] = pow((realnum)10.f,r);
458 }
459 }
460 }
461 }
462
463 /* this will be total collision rate from an upper to a lower level within X */
464 H2_X_source.resize( nLevels_per_elec[0] );
465 H2_X_sink.resize( nLevels_per_elec[0] );
466
468 /* now expand out to include all lower levels as lower state */
469 for( long i=1; i<nLevels_per_elec[0]; ++i )
470 {
471 H2_X_coll_rate.reserve(i,i);
472 }
473 H2_X_coll_rate.alloc();
474
475 ipTransitionSort.reserve( states.size() );
476 for( unsigned nEner = 1; nEner < states.size(); ++nEner )
477 ipTransitionSort.reserve( nEner, nEner );
478 ipTransitionSort.alloc();
479 ipTransitionSort.zero();
480 lgH2_radiative.alloc( ipTransitionSort.clone() );
481 lgH2_radiative.zero();
482
483 trans.resize( (states.size() * (states.size()-1))/2 );
484 AllTransitions.push_back(trans);
485 qList initStates(1);
486 TransitionList initlist("H2InitList",&initStates);
487 vector<TransitionList::iterator> initptrs;
488 initlist.resize(trans.size());
489 initlist.states() = &states;
490 initptrs.resize(trans.size());
491
492 {
493 long lineIndex = 0;
494 TransitionList::iterator tr = initlist.begin();
495 for( unsigned ipHi=1; ipHi< states.size(); ++ipHi )
496 {
497 for( unsigned ipLo=0; ipLo<ipHi; ++ipLo )
498 {
499 (*tr).Junk();
500 (*tr).setHi(ipHi);
501 (*tr).setLo(ipLo);
502 (*tr).Zero();
503 initptrs[lineIndex] = tr;
504 ipTransitionSort[ipHi][ipLo] = lineIndex;
505 lineIndex++;
506 ++tr;
507 }
508 }
509 }
510
511 /* create the main array of lines */
512 H2_SaveLine.reserve(n_elec_states);
513 for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
514 {
515 H2_SaveLine.reserve(iElecHi,nVib_hi[iElecHi]+1);
516 for( long iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
517 {
518 H2_SaveLine.reserve(iElecHi,iVibHi,nRot_hi[iElecHi][iVibHi]+1);
519 for( long iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
520 {
521 /* now the lower levels */
522 /* NB - X is the only lower level considered here, since we are only
523 * concerned with excited electronic levels as a photodissociation process
524 * code exists to relax this assumption - simply change following to iElecHi */
525 long int lim_elec_lo = 0;
526 H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,lim_elec_lo+1);
527 for( long iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
528 {
529 H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,iElecLo,nVib_hi[iElecLo]+1);
530 for( long iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
531 {
532 H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,nRot_hi[iElecLo][iVibLo]+1);
533 }
534 }
535 }
536 }
537 }
538
539 H2_SaveLine.alloc();
540
541 /* zero out array used to save emission line intensities */
542 H2_SaveLine.zero();
543
544 /* space for the energy vector is now malloced, must read trans probs from table */
545 for( long iElec=0; iElec<n_elec_states; ++iElec )
546 {
547 /* ground state has tabulated data */
548 H2_ReadTransprob(iElec,initlist);
549 }
550
551 // sort sys.trans so that radiative lines are at the beginning
552 stable_sort( initptrs.begin(), initptrs.end(), compareEmis );
553 rad_end = trans.begin();
554 // rad_end will be used for the end of the range (non-inclusive) for operations on radiative lines
555 {
556 TransitionList::iterator tr = trans.begin();
557 vector<TransitionList::iterator>::iterator ptr = initptrs.begin();
558 for (size_t i=0; i < initptrs.size(); ++i, ++tr, ++ptr)
559 {
560 (*tr).copy(*(*ptr));
561 if (lgRadiative(tr))
562 {
563 rad_end = tr;
564 }
565 }
566 }
567 ++rad_end;
568 ASSERT( rad_end != trans.end() );
569 // after above sorting, ipTransitionSort is now invalid. Fix here
570 for( unsigned i = 0; i < trans.size(); ++i )
571 {
572 qList::iterator Hi = trans[i].Hi();
573 qList::iterator Lo = trans[i].Lo();
574 ipTransitionSort[ ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()] ][ ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()] ] = i;
575 trans[i].ipHi() = ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()];
576 trans[i].ipLo() = ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()];
577 }
578
579 // link level and line stacks to species in the chemistry network
580 mole.species[ sp->index ].levels = &states;
581 mole.species[ sp->index ].lines = &trans;
582
583 // after above sorting, ipTransitionSort is now invalid. Fix here
584 for( unsigned i = 0; i < trans.size(); ++i )
585 {
586 qList::iterator Hi = trans[i].Hi();
587 qList::iterator Lo = trans[i].Lo();
588 ipTransitionSort[ ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()] ][ ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()] ] = i;
589 //trans[i].ipHi() = ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()];
590 //trans[i].ipLo() = ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()];
591 }
592
593 // this loop is over all transitions
594 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
595 {
596 (*tr).EnergyWN() = (realnum)((*(*tr).Hi()).energy().WN() - (*(*tr).Lo()).energy().WN());
597 /*wavelength of transition in Angstroms */
598 if( (*tr).EnergyWN() > SMALLFLOAT)
599 (*tr).WLAng() = (realnum)(1.e8f/(*tr).EnergyWN() / RefIndex( (*tr).EnergyWN() ) );
600
601 (*tr).Coll().col_str() = 0.;
602 }
603
604 // this loop is over only radiative transitions. Notice the end iterator.
605 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
606 {
607 /* line redistribution function - will use complete redistribution */
608 /* >>chng 04 mar 26, should include damping wings, especially for electronic
609 * transitions, had used doppler core only */
610 (*tr).resetEmis();
611 (*tr).Emis().iRedisFun() = ipCRDW;
612 /* line optical depths in direction towards source of ionizing radiation */
613 (*tr).Emis().TauIn() = opac.taumin;
614 (*tr).Emis().TauCon() = opac.taumin;
615 /* outward optical depth */
616 (*tr).Emis().TauTot() = 1e20f;
617
618 (*tr).Emis().dampXvel() = (realnum)( (*tr).Emis().Aul()/(*tr).EnergyWN()/PI4);
619 (*tr).Emis().gf() = (realnum)(GetGF( (*tr).Emis().Aul(),(*tr).EnergyWN(), (*(*tr).Hi()).g() ) );
620
621 /* derive the absorption coefficient, call to function is gf, wl (A), g_low */
622 (*tr).Emis().opacity() = (realnum)( abscf( (*tr).Emis().gf(), (*tr).EnergyWN(), (*(*tr).Lo()).g()) );
623
624 qList::iterator Hi = (*tr).Hi() ;
625
626 if( (*Hi).n() == 0 )
627 {
628 /* the ground electronic state, most excitations are not direct pumping
629 * (rather indirect, which does not count for ColOvTot) */
630 (*tr).Emis().ColOvTot() = 1.;
631 }
632 else
633 {
634 /* these are excited electronic states, mostly pumped, except for supras */
636 (*tr).Emis().ColOvTot() = 0.;
637 }
638
639 fixit(); // why not include this for excitations within X as well?
640 if( (*Hi).n() > 0 )
641 {
642 /* cosmic ray and non-thermal suprathermal excitation
643 * to singlet state of H2 (B,C,B',D)
644 * cross section is equation 5 of
645 *>>refer H2 cs Liu, W. & Dalgarno, A. 1994, ApJ, 428, 769
646 * relative to H I Lya cross section
647 * this is used to derive H2 electronic excitations
648 * from the H I Lya rate
649 * the following is dimensionless scale factor for excitation
650 * relative to H I Lya
651 */
652 (*tr).Coll().col_str() = (realnum)(
653 pow3( (*tr).WLAng()*1e-8 ) *
654 ( (*(*tr).Hi()).g()/(*(*tr).Lo()).g()) *
655 (*tr).Emis().Aul() *
656 log(3.)*HPLANCK/(160.f*pow3(PI)*0.5*1e-8*EN1EV)/6.e-17);
657 ASSERT( (*tr).Coll().col_str()>0.);
658 }
659 }
660
661 /* Read continuum photodissociation cross section files */
662 if( mole_global.lgStancil )
664
665 /* define branching ratios for deposition of H2 formed on grain surfaces,
666 * set true to use Takahashi distribution, false to use Draine & Bertoldi */
667
668 /* loop over all types of grain surfaces */
669 /* >>chng 02 oct 08, resolved grain types */
670 /* number of different grain types H2_TOP is set in grainvar.h,
671 * types are ice, silicate, graphite */
672 for( long ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
673 {
674 realnum sum = 0., sumj = 0., sumv = 0., sumo = 0., sump = 0.;
675
676 /* first is Draine distribution */
677 if( hmi.chGrainFormPump == 'D' )
678 {
679 long iElec = 0;
680 /* H2 formation temperature, for equation 19, page 271, of
681 * >>refer H2 formation distribution Draine, B.T., & Bertoldi, F., 1996, ApJ, 468, 269-289
682 */
683 double T_H2_FORM = 50000.;
684 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
685 {
686 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
687 {
688 /* no distinction between grain surface composition */
690 /* first term is nuclear H2_stat weight */
691 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * (1.f+iVib) *
692 (realnum)sexp( states[ ipEnergySort[iElec][iVib][iRot] ].energy().K()/T_H2_FORM );
693 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
694 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
695 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
696 if( H2_lgOrtho[iElec][iVib][iRot] )
697 {
698 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
699 }
700 else
701 {
702 /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */
703 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
704 }
705 }
706 }
707 }
708 else if( hmi.chGrainFormPump == 'T' )
709 {
710 /* Takahashi 2001 distribution */
711 double Xrot[H2_TOP] = { 0.14 , 0.15 , 0.15 };
712 double Xtrans[H2_TOP] = { 0.12 , 0.15 , 0.25 };
713 /* first normalize the vibration distribution function */
714 double sumvib = 0.;
715 double EH2;
716 long iElec = 0;
717
718 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
719 {
720 double vibdist;
721 EH2 = EH2_eval( ipH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
722 vibdist = H2_vib_dist( ipH2 , EH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
723 sumvib += vibdist;
724 }
725 /* this branch, use distribution function from
726 * >>refer grain physics Takahashi, Junko, 2001, ApJ, 561, 254-263 */
727 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
728 {
729 double Ev = states[ ipEnergySort[iElec][iVib][0] ].energy().WN()+energy_off;
730 double Fv;
731 /* equation 10 of Takahashi 2001, extra term is energy offset between bottom of potential
732 * the 0,0 level */
733 double Erot;
734 /*fprintf(ioQQQ," Evvvv\t%i\t%li\t%.3e\n", ipH2 ,iVib , Ev*WAVNRYD*EVRYD);*/
735
736 EH2 = EH2_eval( ipH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
737
738 /* equation 3 of Taktahashi & Uehara */
739 Erot = (EH2 - Ev) * Xrot[ipH2] / (Xrot[ipH2] + Xtrans[ipH2]);
740
741 /* email exchange with Junko Takahashi -
742 Thank you for your E-mail.
743 I did not intend to generate negative Erot.
744 I cut off the populations if their energy levels are negative, and made the total
745 population be unity by using normalization factors (see, e.g., Eq. 12).
746
747 I hope that my answer is of help to you and your work is going well.
748 With best wishes,
749 Junko
750
751 >Thanks for the reply. By cutting off the population, should we set the
752 >population to zero when Erot becomes negative, or should we set Erot to
753 >a small positive number?
754
755 I just set the population to zero when Erot becomes negative.
756 Our model is still a rough one for the vibration-rotation distribution function
757 of H2 newly formed on dust, because we have not yet had any exact
758 experimental or theoretical data about it.
759 With best wishes,
760 Junko
761
762 */
763
764 if( Erot > 0. )
765 {
766 /* the vibrational distribution */
767 Fv = H2_vib_dist( ipH2 , EH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() ) / sumvib;
768 /*fprintf(ioQQQ," vibbb\t%li\t%.3e\n", iVib , Fv );*/
769
770 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
771 {
772 double deltaE = states[ ipEnergySort[iElec][iVib][iRot] ].energy().WN() -
773 states[ ipEnergySort[iElec][iVib][0] ].energy().WN();
774 /* equation 6 of Takahashi 2001 */
775 double gaussian = sexp( POW2( (deltaE - Erot) / (0.5 * Erot) ) );
776 /* equation 7 of Takahashi 2001 */
777 double thermal_dist = sexp( deltaE / Erot );
778
779 /* take the mean of the two */
780 double aver = ( gaussian + thermal_dist ) / 2.;
781 /*fprintf(ioQQQ,"rottt\t%i\t%li\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
782 ipH2,iVib,iRot,
783 deltaE*WAVNRYD*EVRYD,
784 gaussian, thermal_dist , aver );*/
785
786 /* thermal_dist does become > 1 since Erot can become negative */
787 ASSERT( gaussian <= 1. /*&& thermal_dist <= 10.*/ );
788
790 /* first term is nuclear H2_stat weight */
791 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * Fv * (2.*iRot+1.) * aver );
792
793 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
794 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
795 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
796 if( H2_lgOrtho[iElec][iVib][iRot] )
797 {
798 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
799 }
800 else
801 {
802 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
803 }
804
805 }
806 }
807 else
808 {
809 /* this branch Erot is non-positive, so no distribution */
810 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
811 {
812 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 0.;
813 }
814 }
815 }
816 }
817 else if( hmi.chGrainFormPump == 't' )
818 {
819 /* thermal distribution at 1.5 eV, as suggested by Amiel & Jaques */
820 /* thermal distribution, upper right column of page 239 of
821 *>>refer H2 formation Le Bourlot, J, 1991, A&A, 242, 235
822 * set with command
823 * set h2 grain formation pumping thermal */
824 double T_H2_FORM = 17329.;
825 long iElec = 0;
826 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
827 {
828 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
829 {
830 /* no distinction between grain surface composition */
832 /* first term is nuclear H2_stat weight */
833 H2_stat[0][iVib][iRot] *
834 (realnum)sexp( states[ ipEnergySort[0][iVib][iRot] ].energy().K()/T_H2_FORM );
835 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
836 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
837 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
838 if( H2_lgOrtho[iElec][iVib][iRot] )
839 {
840 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
841 }
842 else
843 {
844 /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */
845 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
846 }
847 }
848 }
849 }
850 else
852
853 if( lgDebugPrt )
854 fprintf(ioQQQ, "H2 form grains mean J= %.3f mean v = %.3f ortho/para= %.3f\n",
855 sumj/sum , sumv/sum , sumo/sump );
856
857 //iElec = 0;
858 /* now rescale so that integral is unity */
859 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
860 {
861 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
862 {
863 H2_X_grain_formation_distribution[ipH2][iVib][iRot] /= sum;
864 /* print the distribution function */
865 /*if( states[ ipEnergySort[iElec][iVib][iRot] ].energy().WN() < 5200. )
866 fprintf(ioQQQ,"disttt\t%i\t%li\t%li\t%li\t%.4e\t%.4e\t%.4e\t%.4e\n",
867 ipH2, iVib , iRot, (long)H2_stat[0][iVib][iRot] ,
868 states[ ipEnergySort[iElec][iVib][iRot] ].energy().WN(),
869 states[ ipEnergySort[iElec][iVib][iRot] ].energy().K(),
870 H2_X_grain_formation_distribution[ipH2][iVib][iRot],
871 H2_X_grain_formation_distribution[ipH2][iVib][iRot]/H2_stat[0][iVib][iRot]
872 );*/
873 }
874 }
875 }
876
877 return;
878}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
bool is_odd(int j)
Definition cddefines.h:714
#define STATIC
Definition cddefines.h:97
#define POW2
Definition cddefines.h:929
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
double RandGauss(double xMean, double s)
Definition service.cpp:1643
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
const int ipCRDW
Definition cddefines.h:294
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
T pow3(T a)
Definition cddefines.h:938
void fixit(void)
Definition service.cpp:991
double Ryd() const
Definition energy.h:26
TransitionProxy::iterator iterator
Definition transition.h:280
void resize(size_t newsize)
Definition transition.h:285
iterator begin(void)
Definition transition.h:305
qList *& states()
Definition transition.h:325
long int n_elec_states
Definition h2_priv.h:409
multi_arr< realnum, 2 > H2_X_coll_rate
Definition h2_priv.h:606
multi_arr< double, 2 > H2_col_rate_out
Definition h2_priv.h:652
multi_arr< bool, 3 > H2_lgOrtho
Definition h2_priv.h:643
double xSTDNoise
Definition h2_priv.h:391
bool lgREAD_DATA
Definition h2_priv.h:252
long int nVib_hi[N_ELEC]
Definition h2_priv.h:611
multi_arr< realnum, 3 > H2_X_hminus_formation_distribution
Definition h2_priv.h:685
molecule * sp_star
Definition h2_priv.h:564
realnum mass_amu
Definition h2_priv.h:396
bool lgEnabled
Definition h2_priv.h:345
multi_arr< realnum, 3 > H2_dissprob
Definition h2_priv.h:632
multi_arr< realnum, 2 > H2_X_formation
Definition h2_priv.h:656
valarray< long > ipElec_H2_energy_sort
Definition h2_priv.h:688
long int nXLevelsMatrix
Definition h2_priv.h:695
void H2_ReadDissprob(long int nelec)
qList states
Definition h2_priv.h:565
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef
Definition h2_priv.h:668
string shortlabel
Definition h2_priv.h:572
string label
Definition h2_priv.h:571
void H2_Read_hminus_distribution(void)
long int Jlowest[N_ELEC]
Definition h2_priv.h:616
multi_arr< long int, 2 > ipTransitionSort
Definition h2_priv.h:691
void init(void)
valarray< long > ipVib_H2_energy_sort
Definition h2_priv.h:687
void H2_ReadEnergies()
multi_arr< realnum, 3 > CollRateErrFac
Definition h2_priv.h:622
multi_arr< realnum, 3 > H2_stat
Definition h2_priv.h:641
double H2_DissocEnergies[N_ELEC]
Definition h2_priv.h:609
multi_arr< double, 2 > H2_col_rate_in
Definition h2_priv.h:651
multi_arr< realnum, 6 > H2_SaveLine
Definition h2_priv.h:710
vector< CollRateCoeffArray > RateCoefTable
Definition h2_priv.h:623
multi_arr< double, 3 > H2_populations_LTE
Definition h2_priv.h:639
valarray< long > ipRot_H2_energy_sort
Definition h2_priv.h:689
multi_arr< realnum, 3 > CollRateCoeff
Definition h2_priv.h:621
multi_arr< double, 3 > H2_Boltzmann
Definition h2_priv.h:638
const double ENERGY_H2_STAR
Definition h2_priv.h:585
void Read_Mol_Diss_cross_sections(void)
bool lgH2_NOISE
Definition h2_priv.h:383
multi_arr< double, 2 > H2_X_rate_from_elec_excited
Definition h2_priv.h:664
double xMeanNoise
Definition h2_priv.h:391
multi_arr< int, 2 > H2_ipPhoto
Definition h2_priv.h:650
void H2_ReadDissocEnergies(void)
void H2_ReadTransprob(long int nelec, TransitionList &trans)
multi_arr< bool, 2 > lgH2_radiative
Definition h2_priv.h:714
multi_arr< double, 2 > pops_per_vib
Definition h2_priv.h:600
valarray< realnum > H2_X_source
Definition h2_priv.h:674
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
multi_arr< realnum, 2 > H2_X_colden
Definition h2_priv.h:660
TransitionList trans
Definition h2_priv.h:566
multi_arr< realnum, 2 > H2_X_Hmin_back
Definition h2_priv.h:658
multi_arr< realnum, 3 > H2_X_grain_formation_distribution
Definition h2_priv.h:679
valarray< long > nRot_hi[N_ELEC]
Definition h2_priv.h:613
multi_arr< double, 3 > H2_old_populations
Definition h2_priv.h:637
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
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
void H2_CollidRateRead(long int nColl)
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
long int nEner_H2_ground
Definition h2_priv.h:597
ProxyIterator< qStateConstProxy, qStateConstProxy > const_iterator
ProxyIterator< qStateProxy, qStateConstProxy > iterator
Energy & energy() const
@ ipH2
Definition collision.h:17
@ ipNCOLLIDER
Definition collision.h:18
long ipoint(double energy_ryd)
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
const realnum SMALLFLOAT
Definition cpu.h:191
@ H2_TOP
Definition grainvar.h:93
const int N_ELEC
Definition h2_priv.h:27
const int nTE_HMINUS
Definition h2_priv.h:24
const int N_X_COLLIDER
Definition h2_priv.h:13
t_hmi hmi
Definition hmi.cpp:5
double RefIndex(double EnergyWN)
double abscf(double gf, double enercm, double gl)
double GetGF(double trans_prob, double enercm, double gup)
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molecule * findspecies(const char buf[])
molecule * null_mole
static const double energy_off
STATIC double EH2_eval(int ipH2, double DissocEnergy, double energy_wn)
static double XVIB[H2_TOP]
int H2_nRot_add_ortho_para[N_ELEC]
STATIC bool compareEmis(const TransitionList::iterator &tr1, const TransitionList::iterator &tr2)
STATIC double H2_vib_dist(int ipH2, double EH2, double DissocEnergy, double energy_wn)
bool compareEnergies(qStateProxy st1, qStateProxy st2)
static double Xdust[H2_TOP]
STATIC bool lgRadiative(const TransitionList::iterator &tr)
@ DEBUG_ENER
t_opac opac
Definition opacity.cpp:5
UNUSED const double PI
Definition physconst.h:29
UNUSED const double HPLANCK
Definition physconst.h:103
UNUSED const double PI4
Definition physconst.h:35
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double FREQ_1EV
Definition physconst.h:213
UNUSED const double EN1EV
Definition physconst.h:192
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
UNUSED const double WAVNRYD
Definition physconst.h:173
vector< TransitionList > AllTransitions
Definition taulines.cpp:8