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