cloudy  trunk
 All Data Structures 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-2008 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 "physconst.h"
7 #include "mole.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 
17 /* if this is set true then code will print energies and stop */
18 /*@-redef@*/
19 enum {DEBUG_ENER=false};
20 /*@+redef@*/
21 
22 /* this is equation 8 of Takahashi 2001, clearer definition is given in
23  * equation 5 and following discussion of
24  * >>refer H2 formation Takahashi, J., & Uehara, H., 2001, ApJ, 561, 843-857
25  * 0.27eV, convert into wavenumbers */
26 static double XVIB[H2_TOP] = { 0.70 , 0.60 , 0.20 };
27 static double Xdust[H2_TOP] = { 0.04 , 0.10 , 0.40 };
28 
29 /* this is energy difference between bottom of potential well and 0,0
30  * the Takahashi energy scale is from the bottom,
31  * 2201.9 wavenumbers */
32 static const double energy_off = 0.273*FREQ_1EV/SPEEDLIGHT;
33 
34 STATIC double EH2_eval( long int iVib , int ipH2 )
35 {
36  double EH2_here;
37  double Evm = H2_DissocEnergies[0]* XVIB[ipH2] + energy_off;
38 
39  double Ev = (energy_wn[0][iVib][0]+energy_off);
40  /* equation 9 of Takahashi 2001 which is only an approximation
41  double EH2 = H2_DissocEnergies[0] * (1. - Xdust[ipH2] ); */
42  /* equation 1, 2 of
43  * Takahashi, Junko, & Uehara, Hideya, 2001, ApJ, 561, 843-857,
44  * this is heat deposited on grain by H2 formation in this state */
45  double Edust = H2_DissocEnergies[0] * Xdust[ipH2] *
46  ( 1. - ( (Ev - Evm) / (H2_DissocEnergies[0]+energy_off-Evm)) *
47  ( (1.-Xdust[ipH2])/2.) );
48  ASSERT( Edust >= 0. );
49 
50  /* energy is total binding energy less energy lost on grain surface
51  * and energy offset */
52  EH2_here = H2_DissocEnergies[0]+energy_off - Edust;
53  ASSERT( EH2_here >= 0.);
54 
55  return EH2_here;
56 }
57 
58 /*H2_vib_dist evaluates the vibration distribution for H2 formed on grains */
59 STATIC double H2_vib_dist( long int iVib , int ipH2 , double EH2)
60 {
61  double G1[H2_TOP] = { 0.3 , 0.4 , 0.9 };
62  double G2[H2_TOP] = { 0.6 , 0.6 , 0.4 };
63  double Evm = H2_DissocEnergies[0]* XVIB[ipH2] + energy_off;
64  double Fv;
65  if( (energy_wn[0][iVib][0]+energy_off) <= Evm )
66  {
67  /* equation 4 of Takahashi 2001 */
68  Fv = sexp( POW2( (energy_wn[0][iVib][0]+energy_off - Evm)/(G1[ipH2]* Evm ) ) );
69  }
70  else
71  {
72  /* equation 5 of Takahashi 2001 */
73  Fv = sexp( POW2( (energy_wn[0][iVib][0]+energy_off - Evm)/(G2[ipH2]*(EH2 - Evm ) ) ) );
74  }
75  return Fv;
76 }
77 
78 
79 /*H2_Create create variables for the H2 molecule, called by
80  * ContCreatePointers after continuum mesh has been set up */
81 void H2_Create(void)
82 {
83  long int i , iElecHi , iElecLo;
84  long int iVibHi , iVibLo;
85  long int iRotHi , iRotLo;
86  long int iElec, iVib , iRot;
87  long int nColl,
88  nlines;
89  int ier;
90  int nEner;
91  /* >>chng 03 nov 26, from enum H2_type to int */
92  int ipH2;
93  realnum sum , sumj , sumv , sumo , sump;
94 
95  /* this is flag set above - when true h2 code is not executed - this is way to
96  * avoid this code when it is not working */
97  /* only malloc vectors one time per core load */
98  if( lgH2_READ_DATA || !h2.lgH2ON )
99  return;
100 
101  DEBUG_ENTRY( "H2_Create()" );
102 
103  /* print string if H2 debugging is enabled */
104  if( mole.nH2_TRACE )
105  fprintf(ioQQQ," H2_Create called in DEBUG mode.\n");
106 
107  /* this was option to print all electronic states in the main printout - but
108  * number of electronic states was not known at initialization so set to -1,
109  * will set properly now */
110  if( h2.nElecLevelOutput < 1 )
112 
113  /* this var is in h2.h and prevents h2 from being change once committed here */
114  lgH2_READ_DATA = true;
115 
116  /* create special vector that saves collision rates within ground */
117  /* this will contain a vector for collisions within the X ground electronic state,
118  * CollRateFit[vib_up][rot_up][vib_lo][rot_lo][coll_type][3] */
119  /* N_X_COLLIDER is number of different species that collide */
120  iElecHi = 0;
121  /* the current data set is limited to vib hi <= 3 */
122  /* will define collision rates for all possible transitions within X */
123  CollRateFit.reserve(h2.nVib_hi[iElecHi]+1);
124  H2_CollRate.reserve(h2.nVib_hi[iElecHi]+1);
125  for( iVibHi = 0; iVibHi <= h2.nVib_hi[iElecHi]; ++iVibHi )
126  {
127  CollRateFit.reserve(iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
128  H2_CollRate.reserve(iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
129  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
130  {
131  CollRateFit.reserve(iVibHi,iRotHi,h2.nVib_hi[iElecHi]+1);
132  H2_CollRate.reserve(iVibHi,iRotHi,h2.nVib_hi[iElecHi]+1);
133  for( iVibLo=0; iVibLo<(h2.nVib_hi[iElecHi]+1); ++iVibLo )
134  {
135  CollRateFit.reserve(iVibHi,iRotHi,iVibLo,h2.nRot_hi[iElecHi][iVibLo]+1);
136  H2_CollRate.reserve(iVibHi,iRotHi,iVibLo,h2.nRot_hi[iElecHi][iVibLo]+1);
137  for( iRotLo=0; iRotLo<=h2.nRot_hi[iElecHi][iVibLo]; ++iRotLo )
138  {
139  H2_CollRate.reserve(iVibHi,iRotHi,iVibLo,iRotLo,N_X_COLLIDER);
140  CollRateFit.reserve(iVibHi,iRotHi,iVibLo,iRotLo,N_X_COLLIDER);
141  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
142  {
143  /* the last one - the three coefficients */
144  CollRateFit.reserve(iVibHi,iRotHi,iVibLo,iRotLo,nColl,3);
145  }
146  }
147  }
148  }
149  }
150 
151  CollRateFit.alloc();
152  H2_CollRate.alloc();
153 
154  /* zero out the collisional rates since only a minority of them are known*/
155  CollRateFit.zero();
156  H2_CollRate.zero();
157 
158  /* create space for the electronic levels */
162 
163  for( iElec = 0; iElec<mole.n_h2_elec_states; ++iElec )
164  {
165 
167  fprintf(ioQQQ,"elec %li highest vib= %li\n", iElec , h2.nVib_hi[iElec] );
168 
169  ASSERT( h2.nVib_hi[iElec] > 0 );
170 
171  /* h2.nVib_hi is now the highest vibrational level before dissociation,
172  * now allocate space to hold the number of rotation levels */
173  H2_populations.reserve(iElec,h2.nVib_hi[iElec]+1);
174  pops_per_vib.reserve(iElec,h2.nVib_hi[iElec]+1);
175  if( iElec > 0 )
176  H2_dissprob.reserve(iElec,h2.nVib_hi[iElec]+1);
177 
178  /* now loop over all vibrational levels, and find out how many rotation levels there are */
179  /* ground is special since use tabulated data - there are 14 vib states,
180  * ivib=14 is highest */
181  for( iVib = 0; iVib <= h2.nVib_hi[iElec]; ++iVib )
182  {
183  /* lastly create the space for the rotation quantum number */
184  H2_populations.reserve(iElec,iVib,h2.nRot_hi[iElec][iVib]+1);
185  if( iElec > 0 )
186  H2_dissprob.reserve(iElec,iVib,h2.nRot_hi[iElec][iVib]+1);
187  }
188  }
189 
197  H2_lgOrtho.alloc( H2_populations.clone() );
198 
199  pops_per_vib.alloc();
200 
201  H2_dissprob.alloc();
203 
204  /* set this one time, will never be set again, but might be printed */
206 
207  /* these do not have electronic levels - all within X */
208  H2_ipPhoto.reserve(h2.nVib_hi[0]+1);
209 
210  /* space for the vibration levels */
211  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
212  {
213  /* space for the rotation quantum number */
214  H2_ipPhoto.reserve(iVib,h2.nRot_hi[0][iVib]+1);
215  }
216 
217  H2_ipPhoto.alloc();
218  H2_col_rate_in.alloc( H2_ipPhoto.clone() );
219  H2_col_rate_out.alloc( H2_ipPhoto.clone() );
220  H2_rad_rate_in.alloc( H2_ipPhoto.clone() );
221  H2_coll_dissoc_rate_coef.alloc( H2_ipPhoto.clone() );
222  H2_coll_dissoc_rate_coef_H2.alloc( H2_ipPhoto.clone() );
223  H2_X_colden.alloc( H2_ipPhoto.clone() );
224  H2_X_rate_from_elec_excited.alloc( H2_ipPhoto.clone() );
225  H2_X_rate_to_elec_excited.alloc( H2_ipPhoto.clone() );
226  H2_X_colden_LTE.alloc( H2_ipPhoto.clone() );
227  H2_X_formation.alloc( H2_ipPhoto.clone() );
228  H2_X_Hmin_back.alloc( H2_ipPhoto.clone() );
229 
230  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
231  {
232  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
233  {
234  /* >>chng 04 jun 14, set these to bad numbers */
235  H2_rad_rate_in[iVib][iRot] = -BIGFLOAT;
236  H2_coll_dissoc_rate_coef[iVib][iRot] = -BIGFLOAT;
237  H2_coll_dissoc_rate_coef_H2[iVib][iRot] = -BIGFLOAT;
238  }
239  }
240  /* zero out the matrices */
241  H2_X_colden.zero();
242  H2_X_colden_LTE.zero();
243  H2_X_formation.zero();
244  H2_X_Hmin_back.zero();
245  /* rates [cm-3 s-1] from elec excited states into X only vib and rot */
247  /* rates [s-1] to elec excited states from X only vib and rot */
249 
250  /* distribution function for populations following formation from H minus H- */
252  for( i=0; i<nTE_HMINUS; ++i )
253  {
255  /* space for the vibration levels */
256  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
257  {
259  }
260  }
264 
265  /* >>chng 05 jun 20, do not use this, which is highly processed - use ab initio
266  * rates of excitation to electronic levels instead */
267  /* read in cosmic ray distribution information
268  H2_Read_Cosmicray_distribution(); */
269 
270  /* grain formation matrix */
272  for( ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
273  {
275 
276  /* space for the vibration levels */
277  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
278  {
279  H2_X_grain_formation_distribution.reserve(ipH2,iVib,h2.nRot_hi[0][iVib]+1);
280  }
281  }
284 
285  /* space for the energy vector is now malloced, must fill it in,
286  * defines array energy_wn[nelec][iVib][iRot] */
287  for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
288  {
289  /* get energies out of files into array energy_wn[nelec][iVib][iRot] */
290  H2_ReadEnergies(iElec);
291 
292  /* get dissociation probabilities and energies - ground state is stable */
293  if( iElec > 0 )
294  H2_ReadDissprob(iElec);
295  }
296 
297  /* >>02 oct 18, add photodissociation, H2 + hnu => 2H + KE */
298  /* we now have ro-vib energies, now set up threshold array offsets
299  * for photodissociation */
300  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
301  {
302  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
303  {
304  /* this is energy needed to get up to n=3 electronic continuum
305  * H2 cannot dissociate following absorption of a continuum photon into the
306  * continuum above X, which would require little energy, and would be given
307  * by H2_DissocEnergies[0], because that process violates momentum conservation
308  * these would be the triplet states - permitted are into singlets
309  * the effective full wavelength range of this process is from Lya to
310  * Lyman limit in shielded regions
311  * tests show limits are between 850A and 1220A - so Lya is included */
313  /*>>KEYWORD Allison & Dalgarno; continuum dissociation; */
314  double thresh = (H2_DissocEnergies[1] - energy_wn[0][iVib][iRot])*WAVNRYD;
315  /*fprintf(ioQQQ,"DEBUG\t%.2f\t%f\n", RYDLAM/thresh , thresh);*/
316  /* in theory we should be able to assert that thesh just barely reaches
317  * lya, but actual numbers reach down to 0.749 ryd */
318  ASSERT( thresh > 0.74 );
319  H2_ipPhoto[iVib][iRot] = ipoint(thresh);
320  }
321  }
322 
323  nH2_energies = 0;
324  for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec)
325  {
326  /* the number of levels within the molecule */
327  nH2_energies += nLevels_per_elec[iElec];
328  }
329 
331  {
332  for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec)
333  {
334  /* print the number of levels within iElec */
335  fprintf(ioQQQ,"\t(%li %li)", iElec , nLevels_per_elec[iElec] );
336  }
337  fprintf(ioQQQ,
338  " H2_Create: there are %li electronic levels, in each level there are",
340  fprintf(ioQQQ,
341  " for a total of %li levels.\n", nH2_energies );
342  }
343 
344  /* now read in the various sets of collision data */
345  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
346  {
347  /* ground state has tabulated data */
348  H2_CollidRateRead(nColl);
349  }
350  /* option to add gaussian random mole */
351  if( mole.lgH2_NOISE )
352  {
353  iElecHi = 0;
354  /* loop over all transitions */
355  for( iVibHi = 0; iVibHi <= VIB_COLLID; ++iVibHi )
356  {
357  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
358  {
359  for( iVibLo=0; iVibLo<(VIB_COLLID+1); ++iVibLo )
360  {
361  for( iRotLo=0; iRotLo<=h2.nRot_hi[iElecHi][iVibLo]; ++iRotLo )
362  {
363  /* first set of expressions are series that adds to log of rate,
364  * so we will add the gaussian noise */
365  /* >>chng 05 dec 13, GS, last two fits are different,
366  * loop had been to N_X_COLLIDER-1 and so included Stancil data
367  * with noise became negative */
368  for( nColl=0; nColl<N_X_COLLIDER-2; ++nColl )
369  {
370  /* the gaussian random number, many possible collision rates
371  * have no data, and CollRateFit[][][][][][0] is zero - do not
372  * scramble these, only scramble the non-zero rates */
373  if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] != 0. )
374  {
375  /* this returns the log of the random noise */
377  /* check that coefficient 0 is the one we want to hit with the mole */
378  /* these are used at line 2990 below, */
379  CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] += r;
380  }
381  }
382  /* >>chng 04 feb 19, break out last one which is linear and must be treated separately */
383  /* for late one is linear so use pow */
384  if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][N_X_COLLIDER-2][0] != 0. )
385  {
386  /* this returns the log of the random noise */
388  /* check that coefficient 0 is the one we want to hit with the mole */
389  /* these are used at line 2990 below, */
390  CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][N_X_COLLIDER-2][0] *= pow((realnum)10.f,r);
391  }
392  }
393  }
394  }
395  }
396  }
397 
398 
399  /* create arrays for energy sorted referencing of e, v, J */
400  H2_energies = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nH2_energies );
401  H2_ipX_ener_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
402  ipElec_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
403  ipVib_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
404  ipRot_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
405 
406  /* this will be total collision rate from an upper to a lower level within X */
407  H2_X_source = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nLevels_per_elec[0] );
408  H2_X_sink = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nLevels_per_elec[0] );
409 
410  H2_X_coll_rate.reserve(nLevels_per_elec[0]);
411  /* now expand out to include all lower levels as lower state */
412  for( i=1; i<nLevels_per_elec[0]; ++i )
413  {
414  H2_X_coll_rate.reserve(i,i);
415  }
416  H2_X_coll_rate.alloc();
417 
418  /* create a vector of sorted energies for X */
420  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
421  {
422  ipEnergySort.reserve(iElecHi,h2.nVib_hi[iElecHi]+1);
423  for( iVibHi = 0; iVibHi <= h2.nVib_hi[iElecHi]; ++iVibHi )
424  {
425  ipEnergySort.reserve(iElecHi,iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
426  }
427  }
429 
430  nEner = 0;
431  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
432  {
433  /* get set of energies */
434  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
435  {
436  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
437  {
438  H2_energies[nEner] = (realnum)energy_wn[iElecHi][iVibHi][iRotHi];
439  ipElec_H2_energy_sort[nEner] = iElecHi;
440  ipVib_H2_energy_sort[nEner] = iVibHi;
441  ipRot_H2_energy_sort[nEner] = iRotHi;
442  ipEnergySort[iElecHi][iVibHi][iRotHi] = -1;
443  ++nEner;
444  }
445  }
446  }
447 
448  ASSERT( nH2_energies == nEner );
449 
450  /* sort the energy levels so that we can do top-down trickle of states */
451  /*spsort netlib routine to sort array returning sorted indices */
452  spsort(
453  /* input array to be sorted */
454  H2_energies,
455  /* number of values in the molecule */
456  nH2_energies,
457  /* permutation output array */
459  /* flag saying what to do - 1 sorts into increasing order, not changing
460  * the original vector, -1 sorts into decreasing order. 2, -2 change vector */
461  1,
462  /* error condition, should be 0 */
463  &ier);
464 
465  /* now loop over the energies confirming the order */
466  for( nEner=0; nEner<nH2_energies; ++nEner )
467  {
468  if( nEner+1 < nLevels_per_elec[0] )
470  H2_energies[H2_ipX_ener_sort[nEner+1]] );
471  /* following will print quantum indices and energies */
472  /*fprintf(ioQQQ,"%li\t%li\t%.3e\n",
473  ipVib_H2_energy_sort[H2_ipX_ener_sort[nEner]],
474  ipRot_H2_energy_sort[H2_ipX_ener_sort[nEner]],
475  H2_energies[H2_ipX_ener_sort[nEner]]);*/
476  i = H2_ipX_ener_sort[nEner];
477  iElec = ipElec_H2_energy_sort[i];
478  iRot = ipRot_H2_energy_sort[i];
479  iVib = ipVib_H2_energy_sort[i];
480  /* this allows v,J to map into energy sorted array */
481  ipEnergySort[iElec][iVib][iRot] = nEner;
482  }
483 
484  /* now find number of levels in H2g */
485  for( nEner=0; nEner<nLevels_per_elec[0]; ++nEner )
486  {
487  i = H2_ipX_ener_sort[nEner];
488  iRot = ipRot_H2_energy_sort[i];
489  iVib = ipVib_H2_energy_sort[i];
490  if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
491  break;
492  nEner_H2_ground = nEner;
493  }
494  /* need to increment it so that this is the number of levels, not the index
495  * of the highest level */
496  ++nEner_H2_ground;
497 
498  /* this is the number of levels to do with the matrix - set with the
499  * atom h2 matrix command, keyword ALL means to do all of X in the matrix
500  * but number of levels within X was not known when the command was parsed,
501  * so this was set to -1 to defer setting to all until now */
502  if( nXLevelsMatrix<0 )
503  {
504  nXLevelsMatrix = nLevels_per_elec[0];
505  }
506  else if( nXLevelsMatrix > nLevels_per_elec[0] )
507  {
508  fprintf( ioQQQ,
509  " 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",
511  nLevels_per_elec[0]);
512  cdEXIT(EXIT_FAILURE);
513  }
514 
515  /* create the main array of lines */
516  H2Lines.reserve(mole.n_h2_elec_states);
517 
518  nlines = 0;
519  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
520  {
521  H2Lines.reserve(iElecHi,h2.nVib_hi[iElecHi]+1);
522  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
523  {
524  H2Lines.reserve(iElecHi,iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
525  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
526  {
527  /* now the lower levels */
528  /* NB - X is the only lower level considered here, since we are only
529  * concerned with excited electronic levels as a photodissociation process
530  * code exists to relax this assumption - simply change following to iElecHi */
531  long int lim_elec_lo = 0;
532  H2Lines.reserve(iElecHi,iVibHi,iRotHi,1);
533  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
534  {
535  /* want to include all vib states in lower level if
536  * different elec level, but only lower vib levels if
537  * same elec level */
538  long int nv = h2.nVib_hi[iElecLo];
539  /* within X, no transitions v_hi < v_lo transitions exist */
540  if( iElecLo==iElecHi )
541  nv = iVibHi;
542  H2Lines.reserve(iElecHi,iVibHi,iRotHi,iElecLo,nv+1);
543  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
544  {
545  long nr = h2.nRot_hi[iElecLo][iVibLo];
546  if( iElecLo==iElecHi && iVibHi==iVibLo )
547  /* max because cannot malloc 0 bytes */
548  nr = MAX2(1,iRotHi-1);
549  H2Lines.reserve(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,nr+1);
550  nlines += nr+1;
551  }
552  }
553  }
554  }
555  }
556 
557  H2Lines.alloc();
558  H2_SaveLine.alloc( H2Lines.clone() );
559  lgH2_line_exists.alloc( H2Lines.clone() );
560 
561  /* set flag saying that space exists */
562  H2_SaveLine.zero();
563  lgH2_line_exists.zero();
564 
566  fprintf(ioQQQ," There are a total of %li lines in the entire H2 molecule.\n", nlines );
567 
568  /* junk the transitions */
569  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
570  {
571  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
572  {
573  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
574  {
575  /* NB - X is the only lower level considered here, since we are only
576  * concerned with excited electronic levels as a photodissociation process
577  * code exists to relax this assumption - simply change following to iElecHi */
578  long int lim_elec_lo = 0;
579  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
580  {
581  /* want to include all vib states in lower level if different elec level,
582  * but only lower vib levels if same elec level */
583  long int nv = h2.nVib_hi[iElecLo];
584  if( iElecLo==iElecHi )
585  nv = iVibHi;
586  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
587  {
588  long nr = h2.nRot_hi[iElecLo][iVibLo];
589  if( iElecLo==iElecHi && iVibHi==iVibLo )
590  nr = iRotHi-1;
591  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
592  {
593  TransitionJunk( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
594  }
595  }
596  }
597  }
598  }
599  }
600 
601  /* now set up state pointers and zero out the transitions */
602  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
603  {
604  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
605  {
606  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
607  {
608  H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi = AddState2Stack();
609 
610  /* NB - X is the only lower level considered here, since we are only
611  * concerned with excited electronic levels as a photodissociation process
612  * code exists to relax this assumption - simply change following to iElecHi */
613  long int lim_elec_lo = 0;
614  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
615  {
616  /* want to include all vib states in lower level if different elec level,
617  * but only lower vib levels if same elec level */
618  long int nv = h2.nVib_hi[iElecLo];
619  if( iElecLo==iElecHi )
620  nv = iVibHi;
621  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
622  {
623  long nr = h2.nRot_hi[iElecLo][iVibLo];
624  if( iElecLo==iElecHi && iVibHi==iVibLo )
625  nr = iRotHi-1;
626  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
627  {
628  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi =
629  H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi;
630  /* lower level is higher level of a previous transition. */
631  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo =
632  H2Lines[iElecLo][iVibLo][iRotLo][0][0][0].Hi;
633 
634  /* set initial values for each line structure */
635  TransitionZero( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] );
636  }
637  }
638  }
639  }
640  }
641  }
642 
643  /* space for the energy vector is now malloced, must read trans probs from table */
644  for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
645  {
646  /* ground state has tabulated data */
647  H2_ReadTransprob(iElec);
648  }
649 
650  /* set all statistical weights - ours is total statistical weight -
651  * including nuclear spin */
652  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
653  {
654  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
655  {
656  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
657  {
658  /* unlike atoms, for H2 nuclear spin is taken into account - so the
659  * statistical weight of even and odd J states differ by factor of 3 - see page 166, sec par
660  * >>>refer H2 H2_stat wght Shull, J.M., & Beckwith, S., 1982, ARAA, 20, 163-188 */
661  if( is_odd(iRotHi+H2_nRot_add_ortho_para[iElecHi]) )
662  {
663  /* ortho */
664  H2_lgOrtho[iElecHi][iVibHi][iRotHi] = true;
665  H2_stat[iElecHi][iVibHi][iRotHi] = 3.f*(2.f*iRotHi+1.f);
666  }
667  else
668  {
669  /* para */
670  H2_lgOrtho[iElecHi][iVibHi][iRotHi] = false;
671  H2_stat[iElecHi][iVibHi][iRotHi] = (2.f*iRotHi+1.f);
672  }
673  }
674  }
675  }
676 
677  /* set up transition parameters */
678  for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
679  {
680  for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
681  {
682  for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
683  {
684  /* now the lower levels */
685  /* NB - X is the only lower level considered here, since we are only
686  * concerned with excited electronic levels as a photodissociation process
687  * code exists to relax this assumption - simply change following to iElecHi */
688  long int lim_elec_lo = 0;
689  for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
690  {
691  /* want to include all vib states in lower level if different elec level,
692  * but only lower vib levels if same elec level */
693  long int nv = h2.nVib_hi[iElecLo];
694  if( iElecLo==iElecHi )
695  nv = iVibHi;
696  for( iVibLo=0; iVibLo<=nv; ++iVibLo )
697  {
698  long nr = h2.nRot_hi[iElecLo][iVibLo];
699 
700  if( iElecLo==iElecHi && iVibHi==iVibLo )
701  nr = iRotHi-1;
702  for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
703  {
704  /* this is special flag for H2 - these are used in velset (in tfidle.c) to
705  * set doppler velocities for species */
706  /* NB this must be kept parallel with nelem and ionstag in H2Lines transition struc,
707  * since that struc expects to find the abundances here - abund set in hmole.c */
708  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->nelem = LIMELM+3;
709  /* this does not mean anything for a molecule */
710  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->IonStg = 1;
711 
712  /* statistical weights of lower and upper levels */
713  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g = H2_stat[iElecLo][iVibLo][iRotLo];
714  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g = H2_stat[iElecHi][iVibHi][iRotHi];
715 
716  /* energy of the transition in wavenumbers */
717  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN =
718  (realnum)(energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]);
719 
720  /*wavelength of transition in Angstroms */
721  if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN > SMALLFLOAT)
722  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng =
723  (realnum)(1.e8f/H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN /
724  RefIndex( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN));
725 
726  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK =
727  (realnum)(T1CM)*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN;
728 
729  /* energy of photon in ergs */
730  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg =
731  (realnum)(ERG1CM)*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN;
732 
733  /* only do this if radiative transition exists */
734  if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
735  {
736  /* line redistribution function - will use complete redistribution */
737  /* >>chng 04 mar 26, should include damping wings, especially for electronic
738  * transitions, had used doppler core only */
739  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->iRedisFun = ipCRDW;
740 
741  /* line optical depths in direction towards source of ionizing radiation */
742  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauIn = opac.taumin;
743  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauCon = opac.taumin;
744  /* outward optical depth */
745  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauTot = 1e20f;
746 
747 
748  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->dampXvel =
749  (realnum)(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul/
750  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN/PI4);
751 
752  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->gf =
753  (realnum)(GetGF(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul,
754  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN,
755  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g ) );
756 
757  /* derive the absorption coefficient, call to function is gf, wl (A), g_low */
758  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->opacity = (realnum)(
759  abscf(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->gf,
760  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN,
761  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g));
762 
763  if( iElecHi > 0 )
764  {
765  /* cosmic ray and non-thermal suprathermal excitation
766  * to singlet state of H2 (B,C,B',D) from
767  *>>refer H2 cr exc Dalgarno, Yan, & Liu 1999 ApJs;
768  * cross section is scaled from transition probability
769  * and equation 5 of
770  *>>refer H2 cr excit Dalgarno, A., Yan, Min, & Liu, Weihong 1999, ApJS, 125, 237
771  * in terms of hydrogen ionization cross-section
772  *>>refer H2 cr exc Shemansky et al. 1985
773  * this is used in mole_h2.cpp at place where Solomon
774  * rate is added
775  * cosmic ray excitation of triplets is
776  * done
777  */
778  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = (realnum)(
779  pow(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng*1e-8,3)*
780  (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g/
781  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g)*
782  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul*
783  log(3.)*HPLANCK/(160.f*pow(PI,3)*0.5*1e-8*EN1EV));
784  ASSERT(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs>0.);
785  /*fprintf(ioQQQ,"DEBUG %3li %3li %3li %3li %3li %3li %.2e\n",
786  iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,iRotLo,
787  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cs
788  );*/
789  }
790  else
791  {
792  /* excitation within X - not treated this way */
793  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = 0.;
794  }
795  }
796  else
797  {
798  /* Aul is zero but cosmic ray collisions are not
799  * zero in this case - this is the Aul = 0 branch */
800  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = 0.;
801  }
802  }
803  }
804  }
805  }
806  }
807  }
808 
809  /* define branching ratios for deposition of H2 formed on grain surfaces,
810  * set true to use Takahashi distribution, false to use Draine & Bertoldi */
811 
812  /* loop over all types of grain surfaces */
813  /* >>chng 02 oct 08, resolved grain types */
814  /* number of different grain types H2_TOP is set in grainvar.h,
815  * types are ice, silicate, graphite */
816  for( ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
817  {
818  sum = 0.;
819  sumj = 0.;
820  sumv = 0.;
821  sumo = 0.;
822  sump = 0.;
823  iElec = 0;
824  /* first is Draine distribution */
825  if( hmi.chGrainFormPump == 'D' )
826  {
827  /* H2 formation temperature, for equation 19, page 271, of
828  * >>refer H2 formation distribution Draine, B.T., & Bertoldi, F., 1996, ApJ, 468, 269-289
829  */
830 # define T_H2_FORM 50000.
831  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
832  {
833  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
834  {
835  /* no distinction between grain surface composition */
836  H2_X_grain_formation_distribution[ipH2][iVib][iRot] =
837  /* first term is nuclear H2_stat weight */
838  (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * (1.f+iVib) *
839  (realnum)sexp( energy_wn[iElec][iVib][iRot]*T1CM/T_H2_FORM );
840  sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
841  sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
842  sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
843  if( H2_lgOrtho[iElec][iVib][iRot] )
844  {
845  sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
846  }
847  else
848  {
849  /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */
850  sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
851  }
852  }
853  }
854  }
855  else if( hmi.chGrainFormPump == 'T' )
856  {
857  /* Takahashi 2001 distribution */
858  double Xrot[H2_TOP] = { 0.14 , 0.15 , 0.15 };
859  double Xtrans[H2_TOP] = { 0.12 , 0.15 , 0.25 };
860  /* first normalize the vibration distribution function */
861  double sumvib = 0.;
862  double EH2;
863 
864  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
865  {
866  double vibdist;
867  EH2 = EH2_eval( iVib , ipH2 );
868  vibdist = H2_vib_dist( iVib , ipH2 , EH2);
869  sumvib += vibdist;
870  }
871  /* this branch, use distribution function from
872  * >>refer grain physics Takahashi, Junko, 2001, ApJ, 561, 254-263 */
873  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
874  {
875  double Ev = (energy_wn[iElec][iVib][0]+energy_off);
876  double Fv;
877  /* equation 10 of Takahashi 2001, extra term is energy offset between bottom of potential
878  * the 0,0 level */
879  double Erot;
880  /*fprintf(ioQQQ," Evvvv\t%i\t%li\t%.3e\n", ipH2 ,iVib , Ev*WAVNRYD*EVRYD);*/
881 
882  EH2 = EH2_eval( iVib , ipH2 );
883 
884  /* equation 3 of Taktahashi & Uehara */
885  Erot = (EH2 - Ev) * Xrot[ipH2] / (Xrot[ipH2] + Xtrans[ipH2]);
886 
887  /* email exchange with Junko Takahashi -
888  Thank you for your E-mail.
889  I did not intend to generate negative Erot.
890  I cut off the populations if their energy levels are negative, and made the total
891  population be unity by using normalization factors (see, e.g., Eq. 12).
892 
893  I hope that my answer is of help to you and your work is going well.
894  With best wishes,
895  Junko
896 
897  >Thanks for the reply. By cutting off the population, should we set the
898  >population to zero when Erot becomes negative, or should we set Erot to
899  >a small positive number?
900 
901  I just set the population to zero when Erot becomes negative.
902  Our model is still a rough one for the vibration-rotation distribution function
903  of H2 newly formed on dust, because we have not yet had any exact
904  experimental or theoretical data about it.
905  With best wishes,
906  Junko
907 
908  */
909 
910  if( Erot > 0. )
911  {
912  /* the vibrational distribution */
913  Fv = H2_vib_dist( iVib , ipH2 , EH2) / sumvib;
914  /*fprintf(ioQQQ," vibbb\t%li\t%.3e\n", iVib , Fv );*/
915 
916  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
917  {
918  /* equation 6 of Takahashi 2001 */
919  double gaussian =
920  sexp( POW2( (energy_wn[iElec][iVib][iRot] - energy_wn[iElec][iVib][0] - Erot) /
921  (0.5 * Erot) ) );
922  /* equation 7 of Takahashi 2001 */
923  double thermal_dist =
924  sexp( (energy_wn[iElec][iVib][iRot] - energy_wn[iElec][iVib][0]) /
925  Erot );
926 
927  /* take the mean of the two */
928  double aver = ( gaussian + thermal_dist ) / 2.;
929  /*fprintf(ioQQQ,"rottt\t%i\t%li\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
930  ipH2,iVib,iRot,
931  (energy_wn[iElec][iVib][iRot]+energy_off)*WAVNRYD*EVRYD,
932  gaussian, thermal_dist , aver );*/
933 
934  /* thermal_dist does become > 1 since Erot can become negative */
935  ASSERT( gaussian <= 1. /*&& thermal_dist <= 10.*/ );
936 
937  H2_X_grain_formation_distribution[ipH2][iVib][iRot] = (realnum)(
938  /* first term is nuclear H2_stat weight */
939  (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * Fv * (2.*iRot+1.) * aver );
940 
941  sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
942  sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
943  sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
944  if( H2_lgOrtho[iElec][iVib][iRot] )
945  {
946  sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
947  }
948  else
949  {
950  sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
951  }
952 
953  }
954  }
955  else
956  {
957  /* this branch Erot is non-positive, so no distribution */
958  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
959  {
960  H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 0.;
961  }
962  }
963  }
964  }
965 # undef T_H2_FORM
966  else if( hmi.chGrainFormPump == 't' )
967  {
968  /* thermal distribution at 1.5 eV, as suggested by Amiel & Jaques */
969  /* thermal distribution, upper right column of page 239 of
970  *>>refer H2 formation Le Bourlot, J, 1991, A&A, 242, 235
971  * set with command
972  * set h2 grain formation pumping thermal */
973 # define T_H2_FORM 17329.
974  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
975  {
976  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
977  {
978  /* no distinction between grain surface composition */
979  H2_X_grain_formation_distribution[ipH2][iVib][iRot] =
980  /* first term is nuclear H2_stat weight */
981  H2_stat[0][iVib][iRot] *
982  (realnum)sexp( energy_wn[0][iVib][iRot]*T1CM/T_H2_FORM );
983  sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
984  sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
985  sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
986  if( H2_lgOrtho[iElec][iVib][iRot] )
987  {
988  sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
989  }
990  else
991  {
992  /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */
993  sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
994  }
995  }
996  }
997  }
998  else
999  TotalInsanity();
1000 
1002  fprintf(ioQQQ, "H2 form grains mean J= %.3f mean v = %.3f ortho/para= %.3f\n",
1003  sumj/sum , sumv/sum , sumo/sump );
1004 
1005  iElec = 0;
1006  /* now rescale so that integral is unity */
1007  for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
1008  {
1009  for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
1010  {
1011  H2_X_grain_formation_distribution[ipH2][iVib][iRot] /= sum;
1012  /* print the distribution function */
1013  /*if( energy_wn[iElec][iVib][iRot] < 5200. )
1014  fprintf(ioQQQ,"disttt\t%i\t%li\t%li\t%li\t%.4e\t%.4e\t%.4e\t%.4e\n",
1015  ipH2, iVib , iRot, (long)H2_stat[0][iVib][iRot] ,
1016  energy_wn[iElec][iVib][iRot] ,
1017  energy_wn[iElec][iVib][iRot]*T1CM ,
1018  H2_X_grain_formation_distribution[ipH2][iVib][iRot],
1019  H2_X_grain_formation_distribution[ipH2][iVib][iRot]/H2_stat[0][iVib][iRot]
1020  );*/
1021  }
1022  }
1023  }
1024 
1025  /* at this stage the full electronic, vibration, and rotation energies have been defined,
1026  * this is an option to print the energies */
1027  {
1028  /* set following true to get printout, false to not print energies */
1029  if( DEBUG_ENER )
1030  {
1031  /* print title for quantum numbers and energies */
1032  /*fprintf(ioQQQ,"elec\tvib\trot\tenergy\n");*/
1033  for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
1034  {
1035  /* now must specify the number of rotation levels within the vib levels */
1036  for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
1037  {
1038  for( iRot=0; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
1039  {
1040  fprintf(ioQQQ,"%li\t%li\t%li\t%.5e\n",
1041  iElec, iVib, iRot ,
1042  energy_wn[iElec][iVib][iRot]);
1043  }
1044  }
1045  }
1046  /* this will exit the program after printing the level energies */
1047  cdEXIT(EXIT_SUCCESS);
1048  }
1049  }
1050 
1051  return;
1052 }

Generated for cloudy by doxygen 1.8.3.1