cloudy trunk
Loading...
Searching...
No Matches
mole_h2_io.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_ParseSave parse the save h2 command */
4/*H2_PunchDo save some properties of the large H2 molecule */
5/*chMolBranch returns a char with the spectroscopic branch of a transition */
6/*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
7/*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from SaveLineStuff */
8/*H2_Punch_line_data save line data for H2 molecule */
9/*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */
10/*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
11/*H2_ReadEnergies read energies for all electronic levels */
12/*H2_ReadTransprob read transition probabilities */
13/*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
14/*H2_ParseSave parse the save h2 command */
15/*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */
16/*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
17 /*cdH2_Line returns 1 if we found the line,
18 * or false==0 if we did not find the line because ohoto-para transition
19 * or upper level has lower energy than lower level */
20#include "cddefines.h"
21#include "physconst.h"
22#include "save.h"
23#include "hmi.h"
24#include "prt.h"
25#include "secondaries.h"
26#include "grainvar.h"
27#include "input.h"
28#include "phycon.h"
29#include "rfield.h"
30#include "hyperfine.h"
31#include "thermal.h"
32#include "lines.h"
33#include "lines_service.h"
34#include "mole.h"
35#include "dense.h"
36#include "radius.h"
37#include "colden.h"
38#include "taulines.h"
39#include "h2.h"
40#include "h2_priv.h"
41#include "cddrive.h"
42#include "doppvel.h"
43#include "doppvel.h"
44#include "parser.h"
45
47
48/*H2_LinesAdd add in explicit lines from the large H2 molecule, called by lines_molecules */
50{
51 /* H2 not on, so space not allocated */
52 if( !lgEnabled )
53 return;
54
55 DEBUG_ENTRY( "H2_LinesAdd()" );
56
57 if( strcmp( "H2 ", label.c_str() ) == 0 )
58 {
59 /* >>chng 05 nov 04, make info copies of these lines up here
60 * these are among the strongest lines in the 2 micron window and some have nearly the same
61 * wavelength as far weaker lines that may come before them in the line stack. in that case
62 * cdLine would find the much weaker line with the same wavelength.
63 * put strong H2 lines first so that line search will find these, and not far weaker
64 * lines with nearly the same wavelength - these will be duplicated in the output but
65 * these are here for into (the 'i') so does no harm
66 */
67
68 /* >>chng 05 dec 22, had hand entered wavelength in A as second parameter. This gave
69 * rounded off result when set line precision 5 was used. now uses same logic that
70 * PutLine will eventually use - simply enter same wl in Ang */
71 /* 1-0 S(4) - 18910 */
72 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][6] ][ ipEnergySort[0][0][4] ] ], "H2 ", 'i', false, "H2 line");
73 /* 1-0 S(3) - 19570 */
74 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][5] ][ ipEnergySort[0][0][3] ] ], "H2 ", 'i', false, "H2 line");
75 /* 1-0 S(2) - 20330 */
76 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][4] ][ ipEnergySort[0][0][2] ] ], "H2 ", 'i', false, "H2 line");
77 /* 1-0 S1) - 21210 */
78 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][3] ][ ipEnergySort[0][0][1] ] ], "H2 ", 'i', false, "H2 line");
79 /* 1-0 S(0) - 22230 */
80 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][2] ][ ipEnergySort[0][0][0] ] ], "H2 ", 'i', false, "H2 line");
81 /* start Q branch - selection rule requires that J be non-zero, so no Q(0) */
82 /* 1-0 Q(2) - 24130 */
83 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][2] ][ ipEnergySort[0][0][2] ] ], "H2 ", 'i', false, "H2 line");
84 /* 1-0 Q(1) - 24060 */
85 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][1] ][ ipEnergySort[0][0][1] ] ], "H2 ", 'i', false, "H2 line");
86 }
87
88
89 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
90 {
91 qList::iterator Hi = ( (*tr).Hi() );
92 if( (*Hi).n() >= nElecLevelOutput ) continue;
93 qList::iterator Lo = ( (*tr).Lo() );
94 /* all ground vib state rotation lines - first is J to J-2 */
95 PutLine( *tr, "diatoms lines", label.c_str() );
96 if( LineSave.ipass == 0 )
97 {
98 H2_SaveLine[(*Hi).n()][(*Hi).v()][(*Hi).J()][(*Lo).n()][(*Lo).v()][(*Lo).J()] = 0.;
99 }
100 else if( LineSave.ipass == 1 )
101 {
102 H2_SaveLine[(*Hi).n()][(*Hi).v()][(*Hi).J()][(*Lo).n()][(*Lo).v()][(*Lo).J()] +=
103 (realnum)( radius.dVeffAper * (*tr).Emis().xIntensity() );
104 }
105 }
106
107 return;
108}
109
110/*H2_ParseSave parse the save h2 command */
112 char *chHeader)
113{
114 DEBUG_ENTRY( "H2_ParseSave()" );
115
116 save.whichDiatomToPrint[save.nsave] = &(*this);
117
118 /* this provides info on the large H2 molecule */
119 if( p.nMatch("COLU") )
120 {
121 /* save column density */
122 strcpy( save.chSave[save.nsave], "H2cl" );
123
124 /* this is an option to scan off highest vib and rot states
125 * to save pops - first is limit to vibration, then rotation
126 * if no number is entered then 0 is set and all levels punched */
127 /* now get vib limit */
128 save.punarg[save.nsave][0] = (realnum)p.getNumberDefault(
129 "H2 vibration state",0.0);
130
131 /* highest rotation */
132 save.punarg[save.nsave][1] = (realnum)p.getNumberDefault(
133 "H2 rotation state",0.0);
134 /* this says whether to save triplets or a matrix for output -
135 * default is triplets, so only check for matrix */
136 if( p.nMatch( "MATR" ) )
137 {
138 /* matrix */
139 save.punarg[save.nsave][2] = 1;
140 sprintf( chHeader, "#vib\trot\tcolumn density\n" );
141 }
142 else
143 {
144 /* triplets */
145 save.punarg[save.nsave][2] = -1;
146 sprintf( chHeader, "#vib\trot\tEner(K)\tcolden\tcolden/stat wght\tLTE colden\tLTE colden/stat wght\n" );
147 }
148 }
149 else if( p.nMatch("COOL") )
150 {
151 /* heating and cooling rates */
152 strcpy( save.chSave[save.nsave], "H2co" );
153 sprintf( chHeader,
154 "#H2 depth\ttot cool\tTH Sol\tBig Sol\tTH pht dis\tpht dis\tTH Xcool\tXcool \n" );
155 }
156
157 else if( p.nMatch("CREA") )
158 {
159 /* H2 creation rates */
160 fprintf( ioQQQ, " This command has been superseded by the \"creation\" option of the \"save chemistry rates\" command.\n" );
161 fprintf( ioQQQ, " Sorry.\n" );
163 }
164 else if( p.nMatch("DEST") )
165 {
166 /* save H2 destruction - output destruction rates */
167 fprintf( ioQQQ, " This command has been superseded by the \"destruction\" option of the \"save chemistry rates\" command.\n" );
168 fprintf( ioQQQ, " Sorry.\n" );
170 }
171
172 else if( p.nMatch("HEAT") )
173 {
174 /* heating and cooling rates */
175 strcpy( save.chSave[save.nsave], "H2he" );
176 sprintf( chHeader,
177 "#H2 depth\ttot Heat\tHeat(big)\tHeat(TH85)\tDissoc(Big)\tDissoc(TH85) \n" );
178 }
179
180 else if( p.nMatch("LEVE") )
181 {
182 /* save H2 level energies */
183 strcpy( save.chSave[save.nsave], "H2le" );
184 sprintf( chHeader,
185 "#H2 v\tJ\tenergy(wn)\tstat wght\tSum As" );
186 char chHoldit[chN_X_COLLIDER+12];
187 for( int nColl=0; nColl<N_X_COLLIDER; ++nColl )
188 {
189 /* labels for all colliders */
190 sprintf(chHoldit,"\tCritDen %s",chH2ColliderLabels[nColl]);
191 strcat( chHeader , chHoldit );
192 }
193 strcat( chHeader , "\n" );
194 }
195
196 else if( p.nMatch("LINE") )
197 {
198 /* save H2 lines - all in X */
199 strcpy( save.chSave[save.nsave], "H2ln" );
200 sprintf( chHeader,
201 "#H2 line\tEhi\tVhi\tJhi\tElo\tVlo\tJlo\twl(mic)\twl(lab)\tlog L or I\tI/Inorm\tExcit(hi, K)\tg_u h nu * Aul\n" );
202 /* first optional number changes the threshold of weakest line to print*/
203 /* fe2thresh is intensity relative to normalization line,
204 * normally Hbeta, and is set to zero in zero.c */
205
206 /* threshold for faintest line to save, default is 1e-4 of norm line */
208 "faintest line to save",1e-4);
209
210 /* lines from how many electronic states? default is one, just X, and is
211 * obtained with GROUND keyword. ALL will produce all lines from all levels.
212 * else, if a number is present, will be the number. if no number, no keyword,
213 * appear then just ground */
214 if( p.nMatch( "ELEC" ) )
215 {
216 if( p.nMatch(" ALL") )
217 {
218 /* all electronic levels - when done, will set upper limit, the
219 * number of electronic levels actually computed, don't know this yet,
220 * so signify with negative number */
221 nElecLevelOutput = -1;
222 }
223 else if( p.nMatch("GROU") )
224 {
225 /* just the ground electronic state */
227 }
228 else
229 {
231 "electronic levels for output",1.0);
232 }
233 }
234 }
235
236 else if( p.nMatch(" PDR") )
237 {
238 /* creation and destruction processes */
239 strcpy( save.chSave[save.nsave], "H2pd" );
240 sprintf( chHeader, "#H2 creation, destruction. \n" );
241 }
242 else if( p.nMatch("POPU") )
243 {
244 /* save populations */
245 strcpy( save.chSave[save.nsave], "H2po" );
246
247 /* this is an option to scan off highest vib and rot states
248 * to save pops - first is limit to vibration, then rotation
249 * if no number is entered then 0 is set and all levels punched */
250 /* now get vib lim */
251 save.punarg[save.nsave][0] = (realnum)p.getNumberDefault(
252 "highest H2 save vibration state",0.0);
253
254 /* this is limit to rotation quantum index */
255 save.punarg[save.nsave][1] = (realnum)p.getNumberDefault(
256 "highest H2 save rotation state",0.0);
257
258 if( p.nMatch( "ZONE" ) )
259 {
260 /* save v=0 pops for each zone, all along one line */
261 save.punarg[save.nsave][2] = 0;
262 sprintf( chHeader, "#depth\torth\tpar\te=1 rel pop\te=2 rel pop\tv,J rel pops\n" );
263 }
264 else
265 {
266 /* will not do zone output, only output at the end of the calculation
267 * now check whether to save triplets or a matrix for output -
268 * default is triplets, so only check for matrix */
269 if( p.nMatch( "MATR" ) )
270 {
271 /* matrix */
272 save.punarg[save.nsave][2] = 1;
273 sprintf( chHeader, "#vib\trot\tpops\n" );
274 }
275 else
276 {
277 /* triplets */
278 save.punarg[save.nsave][2] = -1;
279 sprintf( chHeader, "#vib\trot\ts\tenergy(wn)\tpops/H2\told/H2\tpops/g/H2\tdep coef\tFin(Col)\tFout(col)\tRCout\tRRout\tRCin\tRRin\n" );
280 }
281 }
282 }
283
284 else if( p.nMatch("RATE") )
285 {
286 /* save h2 rates - creation and destruction rates */
287 strcpy( save.chSave[save.nsave], "H2ra" );
288 sprintf( chHeader,
289 "#depth\tN(H2)\tN(H2)/u(H2)\tA_V(star)\tn(Eval)"
290 "\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tELWERT\tBigH2\telec->H2g\telec->H2s"
291 "\tG(TH85)\tG(DB96)\tCR\tEleclife\tShield(BD96)\tShield(H2)\tBigh2/G0(spc)\ttot dest"
292 "\tHeatH2Dish_TH85\tHeatH2Dexc_TH85\tHeatDish_BigH2\tHeatDexc_BigH2\thtot\n" );
293 }
294 else if( p.nMatch("SOLO") )
295 {
296 /* rate of Solomon process then fracs of exits from each v, J level */
297 strcpy( save.chSave[save.nsave], "H2so" );
298 sprintf( chHeader,
299 "#depth\tSol tot\tpump/dissoc\tpump/dissoc BigH2\tavH2g\tavH2s\tH2g chem/big H2\tH2s chem/big H2\tfrac H2g BigH2\tfrac H2s BigH2\teHi\tvHi\tJHi\tvLo\tJLo\tfrac\twl(A)\n" );
300 }
301 else if( p.nMatch("SPEC") )
302 {
303 /* special save command*/
304 strcpy( save.chSave[save.nsave], "H2sp" );
305 sprintf( chHeader,
306 "#depth\tspecial\n" );
307 }
308 else if( p.nMatch("TEMP") )
309 {
310 /* various temperatures for neutral/molecular gas */
311 strcpy( save.chSave[save.nsave], "H2te" );
312 sprintf( chHeader,
313 "#depth\tH2/H\tn(1/0)\tn(ortho/para)\tT(1/0)\tT(2/0)\tT(3/0)\tT(3/1)\tT(4/0)\tT(kin)\tT(21cm)\tT_sum(1/0)\tT_sum(2/0)\tT_sum(3/0)\tT_sum(3/1)\tT_sum(4/0) \n");
314 }
315 else if( p.nMatch("THER") )
316 {
317 /* thermal heating cooling processes involving H2 */
318 strcpy( save.chSave[save.nsave], "H2th" );
319 sprintf( chHeader,
320 "#depth\tH2/H\tn(1/0)\tn(ortho/para)\tT(1/0)\tT(2/0)\tT(3/0)\tT(3/1)\tT(4/0)\tT(kin)\tT(21cm)\tT_sum(1/0)\tT_sum(2/0)\tT_sum(3/0)\tT_sum(3/1)\tT_sum(4/0) \n");
321 }
322 else
323 {
324 fprintf( ioQQQ,
325 " There must be a second key; they are RATE, LINE, COOL, COLUMN, _PDR, SOLOmon, TEMP, and POPUlations\n" );
327 }
328 return;
329}
330
331
332/*H2_Prt_Zone print H2 info into zone results, called from prtzone for each printed zone */
334{
335 /* no print if H2 not turned on, or not computed for these conditions */
336 if( !lgEnabled || !nCall_this_zone )
337 return;
338
339 DEBUG_ENTRY( "H2_Prt_Zone()" );
340
341 fprintf( ioQQQ, " %s density ", label.c_str() );
342 fprintf(ioQQQ,PrintEfmt("%9.2e", (*dense_total)));
343
344 fprintf( ioQQQ, " orth/par");
345 fprintf(ioQQQ,PrintEfmt("%9.2e", ortho_density / SDIV( para_density )));
346
347 fprintf( ioQQQ, " v0 J=0,3");
348 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][0] ].Pop() / (*dense_total)));
349 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][1] ].Pop() / (*dense_total)));
350 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][2] ].Pop() / (*dense_total)));
351 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][3] ].Pop() / (*dense_total)));
352
353 fprintf( ioQQQ, " TOTv=0,3");
354 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][0] / (*dense_total)));
355 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][1] / (*dense_total)));
356 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][2] / (*dense_total)));
357 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][3] / (*dense_total)));
358 fprintf( ioQQQ, "\n");
359 return;
360}
361
363{
364 /* no print if H2 not turned on, or not computed for these conditions */
365 if( !lgEnabled || !nCall_this_zone )
366 return;
367
368 DEBUG_ENTRY( "H2_PrtDepartCoef()" );
369
370 // print departure coefficients
371 fprintf( ioQQQ, " %s departure coefficients\n", label.c_str() );
372 for( long iElec=0; iElec<n_elec_states; ++iElec )
373 {
374 fprintf( ioQQQ, "%li electronic\n", iElec );
375 for( long iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
376 {
377 for( long iRot=0; iRot<Jlowest[iElec]; ++iRot )
378 fprintf( ioQQQ, " -----" );
379 for( long iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
380 {
381 long i = ipEnergySort[iElec][iVib][iRot];
382 fprintf( ioQQQ, " %5.3f", depart[i] );
383 }
384 fprintf( ioQQQ, "\n" );
385 }
386 fprintf( ioQQQ, "\n" );
387 if( iElec==0 )
388 break;
389 }
390
391 return;
392}
393
394/*H2_Prt_column_density print H2 info into zone results, called from prtzone for each printed zone */
396 /* this is stream used for io, is stdout when called by final,
397 * is save unit when save output generated */
398 FILE *ioMEAN )
399
400{
401 int iVibHi;
402
403 /* no print if H2 not turned on, or not computed for these conditions */
404 if( !lgEnabled || !nCall_this_zone )
405 return;
406
407 DEBUG_ENTRY( "H2_Prt_column_density()" );
408
409 fprintf( ioMEAN, " H2 total ");
410 fprintf(ioMEAN,"%7.3f", log10(SDIV(ortho_colden + para_colden)));
411
412 fprintf( ioMEAN, " H2 ortho ");
413 fprintf(ioMEAN,"%7.3f", log10(SDIV(ortho_colden)));
414
415 fprintf( ioMEAN, " para");
416 fprintf(ioMEAN,"%7.3f", log10(SDIV(para_colden)));
417
418 iVibHi = 0;
419 fprintf( ioMEAN, " v0 J=0,3");
420 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][0])));
421 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][1])));
422 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][2])));
423 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][3])));
424
425 return;
426}
427
428
429/*H2_ReadTransprob read transition probabilities */
430void diatomics::H2_ReadTransprob( long int nelec , TransitionList &trns)
431{
432 const char* cdDATAFILE[N_ELEC] =
433 {
434 "transprob_X.dat",
435 "transprob_B.dat",
436 "transprob_C_plus.dat",
437 "transprob_C_minus.dat",
438 "transprob_B_primed.dat",
439 "transprob_D_plus.dat",
440 "transprob_D_minus.dat"
441 };
442 FILE *ioDATA;
443 char chLine[FILENAME_PATH_LENGTH_2];
444 long int i, n1, n2, n3;
445 long int iVibHi , iVibLo , iRotHi , iRotLo , iElecHi , iElecLo;
446 bool lgEOL;
447
448 DEBUG_ENTRY( "H2_ReadTransprob()" );
449
450 /* now open the data file */
451 char chPath[FILENAME_PATH_LENGTH_2];
452 strcpy( chPath, path.c_str() );
453 strcat( chPath, input.chDelimiter );
454 strcat( chPath, cdDATAFILE[nelec] );
455 ioDATA = open_data( chPath , "r" );
456
457 /* read the first line and check that magic number is ok */
458 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
459 {
460 fprintf( ioQQQ, " H2_ReadTransprob could not read first line of %s\n", cdDATAFILE[nelec]);
462 }
463 i = 1;
464 /* magic number */
465 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
466 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
467 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
468
469 /* magic number
470 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
471 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
472 {
473 fprintf( ioQQQ,
474 " H2_ReadTransprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
475 fprintf( ioQQQ,
476 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
477 n1 , n2 , n3 );
478 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
480 }
481
482 long nlines = 0;
483 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
484 {
485 /* skip comment */
486 if( chLine[0]=='#' )
487 continue;
488 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
489 break;
490
491 double Aul;
492 int n = sscanf(chLine,"%li\t%li\t%li\t%li\t%li\t%li\t%le",
493 &iElecHi , &iVibHi ,&iRotHi , &iElecLo , &iVibLo , &iRotLo , &Aul );
494 ASSERT( n == 7 );
495 ASSERT( iElecHi == nelec );
496 ASSERT( iElecHi < N_ELEC );
497 ASSERT( iElecLo < N_ELEC );
498
499 /* check that we actually included the levels in the model representation */
500 if( iVibHi <= nVib_hi[iElecHi] &&
501 iVibLo <= nVib_hi[iElecLo] &&
502 iRotHi <= nRot_hi[iElecHi][iVibHi] &&
503 iRotLo <= nRot_hi[iElecLo][iVibLo])
504 {
505 long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
506 long ipLo = ipEnergySort[iElecLo][iVibLo][iRotLo];
507 double ener = states[ipHi].energy().WN() - states[ipLo].energy().WN();
508 long lineIndex = ipTransitionSort[ipHi][ipLo];
509
510 /* only lines that have real Aul are added to stack. */
511 trns[lineIndex].AddLine2Stack();
512 trns[lineIndex].Emis().Aul() = (realnum)Aul;
513
514 /* say that this line exists */
515 lgH2_radiative[ipHi][ipLo] = true;
516 ++nlines;
517
518 /* prints transitions with negative energies - should not happen */
519 if( ener <= 0. )
520 {
521 fprintf(ioQQQ,"negative energy H2 transition\t%li\t%li\t%li\t%li\t%.2e\t%.2e\n",
522 iVibHi,iVibLo,iRotHi,iRotLo,Aul,ener);
523 ShowMe();
525 }
526 }
527 }
528
529 fclose( ioDATA );
530 return;
531}
532
533#if 0
534/*H2_Read_Cosmicray_distribution read distribution function for H2 population following cosmic ray collisional excitation */
535void H2_Read_Cosmicray_distribution(void)
536{
537 //CR_PRINT (false), CR_X (1), CR_VIB(15), CR_J(10), CR_EXIT(3)
538
539 /*>>refer H2 cr excit Tine, S., Lepp, S., Gredel, R., & Dalgarno, A. 1997, ApJ, 481, 282 */
540 FILE *ioDATA;
541 char chLine[FILENAME_PATH_LENGTH_2];
542 long int i, n1, n2, n3, iVib , iRot;
543 long neut_frac;
544 bool lgEOL;
545
546 DEBUG_ENTRY( "H2_Read_Cosmicray_distribution()" );
547
548 /* now open the data file */
549 char chPath[FILENAME_PATH_LENGTH_2];
550 strcpy( chPath, path.c_str() );
551 strcat( chPath, input.chDelimiter );
552 strcat( chPath, "H2_CosmicRay_collision.dat" );
553 ioDATA = open_data( chPath, "r" );
554
555 /* read the first line and check that magic number is ok */
556 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
557 {
558 fprintf( ioQQQ, " H2_Read_Cosmicray_distribution could not read first line of %s\n", "H2_Cosmic_collision.dat");
560 }
561
562 i = 1;
563 /* magic number */
564 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
565 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
566 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
567
568 /* magic number
569 * the following is the set of numbers that appear at the start of H2_Cosmic_collision.dat 01 21 03 */
570 if( ( n1 != 1 ) || ( n2 != 21 ) || ( n3 != 3 ) )
571 {
572 fprintf( ioQQQ,
573 " H2_Read_Cosmicray_distribution: the version of %s is not the current version.\n", "H2_Cosmic_collision.dat" );
574 fprintf( ioQQQ,
575 " I expected to find the number 1 21 3 and got %li %li %li instead.\n" ,
576 n1 , n2 , n3 );
577 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
579 }
580
581 /* read until not a comment */
582 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
583 BadRead();
584
585 while( chLine[0]=='#' )
586 {
587 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
588 BadRead();
589 }
590
591 iRot = 1;
592 iVib = 1;
593 neut_frac = 0;
594 while( iVib >= 0 )
595 {
596 long int j_minus_ji;
597 double a[10];
598
599 sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf",
600 &iVib ,&j_minus_ji , &a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9]
601 );
602 /* negative iVib says end of data */
603 if( iVib < 0 )
604 continue;
605
606 /* cr_rate[CR_X][CR_VIB][CR_J][CR_EXIT];*/
607 /* check that we actually included the levels in the model representation */
608 ASSERT( iVib < CR_VIB );
609 ASSERT( j_minus_ji == -2 || j_minus_ji == +2 || j_minus_ji == 0 );
610 ASSERT( neut_frac < CR_X );
611
612 /* now make i_minus_ji an array index */
613 j_minus_ji = 1 + j_minus_ji/2;
614 ASSERT( j_minus_ji>=0 && j_minus_ji<=2 );
615
616 /* option to add Gaussian random mole */
617 for( iRot=0; iRot<CR_J; ++iRot )
618 {
619 cr_rate[neut_frac][iVib][iRot][j_minus_ji] = (realnum)a[iRot];
620 }
621 if( lgH2_NOISECOSMIC )
622 {
623 realnum r;
624 r = (realnum)RandGauss( xMeanNoise , xSTDNoise );
625
626 for( iRot=0; iRot<CR_J; ++iRot )
627 {
628 cr_rate[neut_frac][iVib][iRot][j_minus_ji] *= (realnum)pow(10.,(double)r);
629 }
630 }
631
632 if( CR_PRINT )
633 {
634 fprintf(ioQQQ,"cr rate\t%li\t%li", iVib , j_minus_ji );
635 for( iRot=0; iRot<CR_J; ++iRot )
636 {
637 fprintf(ioQQQ,"\t%.3e", cr_rate[neut_frac][iVib][iRot][j_minus_ji] );
638 }
639 fprintf(ioQQQ,"\n" );
640 }
641
642 /* now get next line */
643 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
644 BadRead();
645 while( chLine[0]=='#' )
646 {
647 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
648 BadRead();
649 }
650 }
651 fclose( ioDATA );
652
653 return;
654}
655#endif
656
658{
659public:
660 bool operator<( const level_tmp& second ) const
661 {
662 if( eWN < second.eWN )
663 return true;
664 else
665 return false;
666 }
667 long n, v, J;
668 double eWN;
669};
670
671/*H2_ReadEnergies read energies for all electronic levels */
673{
674 DEBUG_ENTRY( "H2_ReadEnergies()" );
675
676 vector<int> n, v, J;
677 vector<double>eWN;
678
679 for( long nelec=0; nelec<n_elec_states; ++nelec )
680 {
681 /* get energies out of files */
682 H2_ReadEnergies(nelec,n,v,J,eWN);
683 }
684
685 vector<level_tmp> levels;
686 levels.resize( n.size() );
687 ASSERT( levels.size() > 0 );
688 for( unsigned i = 0; i < n.size(); ++i )
689 {
690 levels[i].n = n[i];
691 levels[i].v = v[i];
692 levels[i].J = J[i];
693 levels[i].eWN = eWN[i];
694 }
695
696 // now get levels in energy order (comparison done by operator< of level_tmp class)
697 sort( levels.begin(), levels.end() );
698
699 // populate states
700 for( vector<level_tmp>::iterator lev = levels.begin(); lev != levels.end(); ++lev )
701 {
702 states.resize( states.size() + 1 );
703 long i = states.size() - 1;
704 states[i].n() = lev->n;
705 states[i].v() = lev->v;
706 states[i].J() = lev->J;
707 states[i].energy().set( lev->eWN, "cm^-1" );
708 /* NB this must be kept parallel with nelem and ionstag in transition struc,
709 * since that struc expects to find the abundances here - abund set in hmole.c */
710 states[i].nelem() = -1;
711 /* this does not mean anything for a molecule */
712 states[i].IonStg() = -1;
713 strcpy( states[i].chLabel(), label.c_str() );
714 }
715
716 ASSERT( states.size() > 0 );
717 ASSERT( states.size() == levels.size() );
718
719 for( long nelec=0; nelec<n_elec_states; ++nelec )
720 {
721 ASSERT( nLevels_per_elec[nelec] > 0 );
722 ASSERT( nVib_hi[nelec] > 0 );
723 ASSERT( nVib_hi[nelec] > Jlowest[nelec] );
724
725 nRot_hi[nelec].resize( nVib_hi[nelec]+1 );
726 nRot_hi[nelec] = 0;
727 }
728
729 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
730 {
731 nRot_hi[ (*st).n() ][ (*st).v() ] = MAX2( nRot_hi[ (*st).n() ][ (*st).v() ], (*st).J() );
732 }
733
734 return;
735}
736
737void diatomics::H2_ReadEnergies( long int nelec, vector<int>& n, vector<int>& v, vector<int>&J, vector<double>& eWN )
738{
739 DEBUG_ENTRY("diatomics::H2_ReadEnergies()");
740 const char* cdDATAFILE[N_ELEC] =
741 {
742 "energy_X.dat",
743 "energy_B.dat",
744 "energy_C_plus.dat",
745 "energy_C_minus.dat",
746 "energy_B_primed.dat",
747 "energy_D_plus.dat",
748 "energy_D_minus.dat"
749 };
750 /* now open the data file */
751 char chPath[FILENAME_PATH_LENGTH_2];
752 strcpy( chPath, path.c_str() );
753 strcat( chPath, input.chDelimiter );
754 strcat( chPath, cdDATAFILE[nelec] );
755 FILE *ioDATA = open_data( chPath, "r" );
756
757 char chLine[FILENAME_PATH_LENGTH_2];
758 long int i, n1, n2, n3;
759 bool lgEOL;
760
761 /* read the first line and check that magic number is ok */
762 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
763 {
764 fprintf( ioQQQ, " H2_ReadEnergies could not read first line of %s\n", cdDATAFILE[nelec]);
766 }
767 i = 1;
768 /* magic number */
769 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
770 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
771 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
772
773 /* magic number
774 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
775 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
776 {
777 fprintf( ioQQQ,
778 " H2_ReadEnergies: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
779 fprintf( ioQQQ,
780 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
781 n1 , n2 , n3 );
782 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
784 }
785
786 /* this will count the number of levels within each electronic state */
787 nLevels_per_elec[nelec] = 0;
788 nVib_hi[nelec] = 0;
789 Jlowest[nelec] = LONG_MAX;
790
791 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
792 {
793 /* skip comment */
794 if( chLine[0]=='#' )
795 continue;
796 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
797 break;
798 long iVib, iRot;
799 double energyWN;
800 int nReads = sscanf(chLine,"%li\t%li\t%le", &iVib, &iRot, &energyWN );
801 ASSERT( nReads == 3 );
802 ASSERT( iVib >= 0 );
803 ASSERT( iRot >= 0 );
804 ASSERT( energyWN > 0. || (nelec==0 && iVib==0 && iRot==0 ) );
805
806 n.push_back( nelec );
807 v.push_back( iVib );
808 J.push_back( iRot );
809 eWN.push_back( energyWN );
810
811 // update limits
812 nVib_hi[nelec] = MAX2( nVib_hi[nelec], iVib );
813 Jlowest[nelec] = MIN2( Jlowest[nelec], iRot );
814 /* increment number of levels within this electronic state */
815 ++nLevels_per_elec[nelec];
816 }
817
818 ASSERT( n.size() > 0 );
819 ASSERT( nLevels_per_elec[nelec] > 0 );
820 ASSERT( nVib_hi[nelec] > 0 );
821 ASSERT( nVib_hi[nelec] > Jlowest[nelec] );
822
823 fclose( ioDATA );
824
825 return;
826}
827
828/*H2_ReadDissocEnergies read energies for all electronic levels */
830{
831 const char* cdDATAFILE = "energy_dissoc.dat";
832 FILE *ioDATA;
833 char chLine[FILENAME_PATH_LENGTH_2];
834 long int i, n1, n2, n3;
835 bool lgEOL;
836
837 DEBUG_ENTRY( "H2_ReadDissocEnergies()" );
838
839 /* now open the data file */
840 char chPath[FILENAME_PATH_LENGTH_2];
841 strcpy( chPath, path.c_str() );
842 strcat( chPath, input.chDelimiter );
843 strcat( chPath, cdDATAFILE );
844 ioDATA = open_data( chPath, "r" );
845
846 /* read the first line and check that magic number is ok */
847 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
848 {
849 fprintf( ioQQQ, " H2_ReadDissocEnergies could not read first line of %s\n", cdDATAFILE );
851 }
852 i = 1;
853 /* magic number */
854 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
855 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
856 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
857
858 /* magic number
859 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
860 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
861 {
862 fprintf( ioQQQ,
863 " H2_ReadDissocEnergies: the version of %s is not the current version.\n", cdDATAFILE );
864 fprintf( ioQQQ,
865 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
866 n1 , n2 , n3 );
867 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
869 }
870
871 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
872 {
873 /* skip comment */
874 if( chLine[0]=='#' )
875 continue;
876 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
877 break;
878 long iElec;
879 double energyWN;
880 int n = sscanf(chLine,"%li\t%le", &iElec, &energyWN );
881 ASSERT( n == 2 );
882 ASSERT( iElec >= 0 );
883 ASSERT( iElec < N_ELEC );
884 ASSERT( energyWN > 0. );
885 H2_DissocEnergies[iElec] = energyWN;
886 }
887 fclose( ioDATA );
888
889 return;
890}
891
892/*H2_ReadDissprob read dissociation probabilities and kinetic energies for all electronic levels */
893void diatomics::H2_ReadDissprob( long int nelec )
894{
895 const char* cdDATAFILE[N_ELEC] =
896 {
897 "dissprob_X.dat",/* this does not exist and nelec == 0 is not valid */
898 "dissprob_B.dat",
899 "dissprob_C_plus.dat",
900 "dissprob_C_minus.dat",
901 "dissprob_B_primed.dat",
902 "dissprob_D_plus.dat",
903 "dissprob_D_minus.dat"
904 };
905 char chLine[FILENAME_PATH_LENGTH_2];
906 bool lgEOL;
907
908 DEBUG_ENTRY( "H2_ReadDissprob()" );
909
910 ASSERT( nelec > 0 );
911
912 /* now open the data file */
913 char chPath[FILENAME_PATH_LENGTH_2];
914 strcpy( chPath, path.c_str() );
915 strcat( chPath, input.chDelimiter );
916 strcat( chPath, cdDATAFILE[nelec] );
917 FILE *ioDATA = open_data( chPath, "r" );
918
919 /* read the first line and check that magic number is ok */
920 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
921 {
922 fprintf( ioQQQ, " H2_ReadDissprob could not read first line of %s\n", cdDATAFILE[nelec]);
924 }
925 long i = 1;
926 /* magic number */
927 long n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
928 long n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
929 long n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
930
931 /* magic number
932 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
933 if( ( n1 != 3 ) || ( n2 != 2 ) || ( n3 != 11 ) )
934 {
935 fprintf( ioQQQ,
936 " H2_ReadDissprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
937 fprintf( ioQQQ,
938 " I expected to find the number 3 2 11 and got %li %li %li instead.\n" ,
939 n1 , n2 , n3 );
940 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
942 }
943
944 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
945 {
946 /* skip comment */
947 if( chLine[0]=='#' )
948 continue;
949 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
950 break;
951
952 long iVib, iRot;
953 double a, b;
954 i = 1;
955 sscanf(chLine,"%li\t%li\t%le\t%le",
956 &iVib, &iRot,
957 /* dissociation probability */
958 &a ,
959 /* dissociation kinetic energy - eV not ergs */
960 &b);
961
962 /* these have to agree if data file is valid */
963 //ASSERT( iVib >= 0 );
964 //ASSERT( iVib <= nVib_hi[nelec] );
965 //ASSERT( iRot >= Jlowest[nelec] );
966 //ASSERT( iRot <= nRot_hi[nelec][iVib] );
967 if( ( iVib < 0 ) ||
968 ( iVib > nVib_hi[nelec] ) ||
969 ( iRot < Jlowest[nelec] ) ||
970 ( iRot > nRot_hi[nelec][iVib] ) )
971 continue;
972
973 /* dissociation probability */
974 H2_dissprob[nelec][iVib][iRot] = (realnum)a;
975 /* dissociation kinetic energy - eV not ergs */
976 H2_disske[nelec][iVib][iRot] = (realnum)b;
977 }
978 fclose( ioDATA );
979 return;
980}
981
982
983/*H2_Read_hminus_distribution read distribution function for H2 population following formation from H minus */
985{
986 FILE *ioDATA;
987 char chLine[FILENAME_PATH_LENGTH_2];
988 long int i, n1, n2, n3, iVib , iRot;
989 bool lgEOL;
990 double sumrate[nTE_HMINUS] = {0};
991 /* set true for lots of printout */
992 const bool lgH2HMINUS_PRT = false;
993
994 DEBUG_ENTRY( "H2_Read_hminus_distribution()" );
995
996 /* now open the data file */
997 char chPath[FILENAME_PATH_LENGTH_2];
998 strcpy( chPath, path.c_str() );
999 strcat( chPath, input.chDelimiter );
1000 strcat( chPath, "hminus_deposit.dat" );
1001 ioDATA = open_data( chPath, "r" );
1002
1003 /* read the first line and check that magic number is ok */
1004 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1005 {
1006 fprintf( ioQQQ, " H2_Read_hminus_distribution could not read first line of %s\n", chPath );
1008 }
1009
1010 i = 1;
1011 /* magic number */
1012 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1013 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1014 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
1015
1016 /* magic number
1017 * the following is the set of numbers that appear at the start of H2_hminus_deposit.dat 01 08 10 */
1018 if( ( n1 != 2 ) || ( n2 != 10 ) || ( n3 != 17 ) )
1019 {
1020 fprintf( ioQQQ,
1021 " H2_Read_hminus_distribution: the version of %s is not the current version.\n", chPath );
1022 fprintf( ioQQQ,
1023 " I expected to find the number 2 10 17 and got %li %li %li instead.\n" ,
1024 n1 , n2 , n3 );
1025 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
1027 }
1028
1029 /* read until not a comment */
1030 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1031 BadRead();
1032
1033 while( chLine[0]=='#' )
1034 {
1035 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1036 BadRead();
1037 }
1038
1039 iRot = 1;
1040 iVib = 1;
1041 while( iVib >= 0 )
1042 {
1043 /* set true to print rates */
1044
1045 double a[nTE_HMINUS] , ener;
1046 sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf",
1047 &iVib ,&iRot , &ener, &a[0],&a[1],&a[2] , &a[3],&a[4],&a[5] ,&a[6]
1048 );
1049 /* negative iVib says end of data */
1050 if( iVib < 0 )
1051 continue;
1052
1053 /* check that we actually included the levels in the model representation */
1054 ASSERT( iVib <= nVib_hi[0] &&
1055 iRot <= nRot_hi[0][iVib] );
1056
1057 if( lgH2HMINUS_PRT )
1058 fprintf(ioQQQ,"hminusss\t%li\t%li", iVib , iRot );
1059 for( i=0; i<nTE_HMINUS; ++i )
1060 {
1061 H2_X_hminus_formation_distribution[i][iVib][iRot] = (realnum)pow(10.,-a[i]);
1062 sumrate[i] += H2_X_hminus_formation_distribution[i][iVib][iRot];
1063 if( lgH2HMINUS_PRT )
1064 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
1065 }
1066 if( lgH2HMINUS_PRT )
1067 fprintf(ioQQQ,"\n" );
1068
1069 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1070 BadRead();
1071 while( chLine[0]=='#' )
1072 {
1073 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
1074 BadRead();
1075 }
1076 }
1077 fclose( ioDATA );
1078
1079 if( lgH2HMINUS_PRT )
1080 {
1081 /* print total rate */
1082 fprintf(ioQQQ," total H- formation rate ");
1083 /* convert temps to log */
1084 for(i=0; i<nTE_HMINUS; ++i )
1085 {
1086 fprintf(ioQQQ,"\t%.3e" , sumrate[i]);
1087 }
1088 fprintf(ioQQQ,"\n" );
1089 }
1090
1091 /* convert to dimensionless factors that add to unity */
1092 for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
1093 {
1094 for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
1095 {
1096 for(i=0; i<nTE_HMINUS; ++i )
1097 {
1098 H2_X_hminus_formation_distribution[i][iVib][iRot] /= (realnum)sumrate[i];
1099 }
1100 }
1101 }
1102
1103 if( lgH2HMINUS_PRT )
1104 {
1105 /* print total rate */
1106 fprintf(ioQQQ," H- distribution function ");
1107 for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
1108 {
1109 for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
1110 {
1111 fprintf(ioQQQ,"%li\t%li", iVib , iRot );
1112 for(i=0; i<nTE_HMINUS; ++i )
1113 {
1114 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
1115 }
1116 fprintf(ioQQQ,"\n" );
1117 }
1118 }
1119 }
1120 return;
1121}
1122
1123/* ===================================================================== */
1124/*H2_Punch_line_data save line data for H2 molecule */
1126 /* io unit for save */
1127 FILE* ioPUN ,
1128 /* save all levels if true, only subset if false */
1129 bool lgDoAll )
1130{
1131 if( !lgEnabled )
1132 return;
1133
1134 DEBUG_ENTRY( "H2_Punch_line_data()" );
1135
1136 if( lgDoAll )
1137 {
1138 fprintf( ioQQQ,
1139 " H2_Punch_line_data ALL option not implemented in H2_Punch_line_data yet 1\n" );
1141 }
1142 else
1143 {
1144 bool lgPrint = false;
1145 fprintf( ioPUN, "#Eu\tVu\tJu\tEl\tVl\tJl\tWL\tgl\tgu\tgf\tA\tCS\tn(crt)\n" );
1146 /* save line date, looping over all possible lines */
1147 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
1148 {
1149 if( (*tr).ipCont() <= 0 )
1150 continue;
1151 (*tr).Coll().col_str() = 0.;
1152 qList::iterator Hi = ( (*tr).Hi() );
1153 qList::iterator Lo = ( (*tr).Lo() );
1154 /* print quantum indices */
1155 fprintf(ioPUN,"%2li\t%2li\t%2li\t%2li\t%2li\t%2li\t",
1156 (*Hi).n(), (*Hi).v(), (*Hi).J(),
1157 (*Lo).n(), (*Lo).v(), (*Lo).J() );
1158 Save1LineData( *tr, ioPUN, false , lgPrint);
1159 }
1160
1161 fprintf( ioPUN , "\n");
1162 }
1163 return;
1164}
1165
1166/*H2_PunchLineStuff include H2 lines in punched optical depths, etc, called from SaveLineStuff */
1167void diatomics::H2_PunchLineStuff( FILE * io , realnum xLimit , long index)
1168{
1169 if( !lgEnabled )
1170 return;
1171
1172 DEBUG_ENTRY( "H2_PunchLineStuff()" );
1173
1174 /* loop over all possible lines */
1175 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
1176 {
1177 if( (*tr).ipCont() <= 0 )
1178 continue;
1179 Save1Line( *tr, io, xLimit, index, GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN]));
1180 }
1181
1182 return;
1183}
1184
1185
1186/*H2_Prt_line_tau print line optical depths, called from premet in response to print line optical depths command*/
1188{
1189 if( !lgEnabled )
1190 return;
1191
1192 DEBUG_ENTRY( "H2_Prt_line_tau()" );
1193
1194 /* loop over all possible lines */
1195 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
1196 {
1197 if( (*tr).ipCont() <= 0 )
1198 continue;
1199 prme( false, *tr );
1200 }
1201
1202 return;
1203}
1204
1205
1206/*chMolBranch returns a char with the spectroscopic branch of a transition */
1207STATIC char chMolBranch( long iRotHi , long int iRotLo )
1208{
1209 /* these are the spectroscopic branches */
1210 char chBranch[5] = {'O','P','Q','R','S'};
1211 /* this is the index within the chBranch array */
1212 int ip = 2 + (iRotHi - iRotLo);
1213 if( ip<0 || ip>=5 )
1214 {
1215 fprintf(ioQQQ," chMolBranch called with insane iRotHi=%li iRotLo=%li ip=%i\n",
1216 iRotHi , iRotLo , ip );
1217 ip = 0;
1218 }
1219
1220 return( chBranch[ip] );
1221}
1222
1223/*H2_PunchDo save some properties of the large H2 molecule */
1224void diatomics::H2_PunchDo( FILE* io , char chJOB[] , const char chTime[] , long int ipPun )
1225{
1226 DEBUG_ENTRY( "H2_PunchDo()" );
1227
1228 /* which job are we supposed to do? This routine is active even when H2 is not turned on
1229 * so do not test on lgEnabled initially */
1230
1231 /* H2 populations computed in last zone -
1232 * give all of molecule in either matrix or triplet format */
1233 if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") == 0) &&
1234 (save.punarg[ipPun][2] != 0) )
1235 {
1236 /* >>chng 04 feb 19, do not save if H2 not yet evaluated */
1237 if( lgEnabled && lgEvaluated )
1238 {
1239 long iVibHi= 0;
1240 long iRotHi = 0;
1241 long iElecHi=0;
1242 long LimVib, LimRot;
1243 /* the limit to the number of vibration levels punched -
1244 * default is all, but first two numbers on save h2 pops command
1245 * reset limit */
1246 /* this is limit to vibration */
1247 if( save.punarg[ipPun][0] > 0 )
1248 {
1249 LimVib = (long)save.punarg[ipPun][0];
1250 }
1251 else
1252 {
1253 LimVib = nVib_hi[iElecHi];
1254 }
1255
1256 /* first save the current ortho, para, and total H2 density */
1257 fprintf(io,"%i\t%i\t%.3e\tortho\n",
1258 103 ,
1259 103 ,
1260 ortho_density );
1261 fprintf(io,"%i\t%i\t%.3e\tpara\n",
1262 101 ,
1263 101 ,
1264 para_density );
1265 fprintf(io,"%i\t%i\t%.3e\ttotal\n",
1266 0 ,
1267 0 ,
1268 (*dense_total) );
1269
1270 /* now save the actual populations, first part both matrix and triplets */
1271 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
1272 {
1273 /* this is limit to rotation quantum index */
1274 if( save.punarg[ipPun][1] > 0 )
1275 {
1276 LimRot = (long)MIN2(
1277 save.punarg[ipPun][1] , (realnum)nRot_hi[iElecHi][iVibHi]);
1278 }
1279 else
1280 {
1281 LimRot = nRot_hi[iElecHi][iVibHi];
1282 }
1283 if( save.punarg[ipPun][2] > 0 )
1284 {
1285 long int i;
1286 /* this option save matrix */
1287 if( iVibHi == 0 )
1288 {
1289 fprintf(io,"vib\\rot");
1290 /* this is first vib, so make row of rot numbs */
1291 for( i=0; i<=LimRot; ++i )
1292 {
1293 fprintf(io,"\t%li",i);
1294 }
1295 fprintf(io,"\n");
1296 }
1297 fprintf(io,"%li",iVibHi );
1298 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1299 {
1300 fprintf(io,"\t%.3e",
1301 states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].Pop()/(*dense_total) );
1302 }
1303 fprintf(io,"\n" );
1304 }
1305 else if( save.punarg[ipPun][2] < 0 )
1306 {
1307 /* this option save triplets - the default */
1308 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1309 {
1310 /* this will say whether ortho or para,
1311 * H2_lgOrtho is 0 or 1 depending on whether or not ortho,
1312 * so chlgPara[H2_lgOrtho] gives P or O for printing */
1313 const char chlgPara[2]={'P','O'};
1314 const long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
1315
1316 /* intensity, relative to normalization line, for faintest line to save */
1317 fprintf(io,"%li\t%li\t%c\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
1318 /* upper vibration and rotation quantum numbers */
1319 iVibHi , iRotHi ,
1320 /* an 'O' or 'P' for ortho or para */
1321 chlgPara[H2_lgOrtho[iElecHi][iVibHi][iRotHi]],
1322 /* the level excitation energy in wavenumbers */
1323 states[ipHi].energy().WN(),
1324 /* actual population relative to total H2 */
1325 states[ipHi].Pop()/(*dense_total) ,
1326 /* old level populations for comparison */
1327 H2_old_populations[iElecHi][iVibHi][iRotHi]/(*dense_total) ,
1328 /* populations per h2 and per statistical weight */
1329 states[ipHi].Pop()/(*dense_total)/H2_stat[iElecHi][iVibHi][iRotHi] ,
1330 /* LTE departure coefficient */
1331 /* >>chng 05 jan 26, missing factor of H2 abundance LTE is norm to unity, not tot abund */
1332 states[ipHi].Pop()/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi]*(*dense_total) ) ,
1333 /* fraction of exits that were collisional */
1334 H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[0][iVibHi][iRotHi]) ,
1335 /* fraction of entries that were collisional */
1336 H2_col_rate_in[iVibHi][iRotHi]/SDIV(H2_col_rate_in[iVibHi][iRotHi]+H2_rad_rate_in[iVibHi][iRotHi]),
1337 /* collisions out */
1338 H2_col_rate_out[iVibHi][iRotHi],
1339 /* radiation out */
1340 H2_rad_rate_out[0][iVibHi][iRotHi] ,
1341 /* radiation out */
1342 H2_col_rate_in[iVibHi][iRotHi],
1343 /* radiation in */
1344 H2_rad_rate_in[iVibHi][iRotHi]
1345 );
1346 }
1347 }
1348 }
1349 }
1350 }
1351 /* save H2 populations for each zone
1352 * populations of v=0 for each zone */
1353 else if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") != 0) &&
1354 (save.punarg[ipPun][2] == 0) )
1355 {
1356 /* >>chng 04 feb 19, do not save if h2 not yet evaluated */
1357 if( lgEnabled && lgEvaluated )
1358 {
1359 fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth_mid_zone ,
1361 /* rel pops of first two excited electronic states */
1362 fprintf(io,"\t%.3e\t%.3e",
1364 long iElecHi = 0;
1365 long iVibHi = 0;
1366 long LimVib, LimRot;
1367 /* this is limit to vibration quantum index */
1368 if( save.punarg[ipPun][0] > 0 )
1369 {
1370 LimVib = (long)save.punarg[ipPun][1];
1371 }
1372 else
1373 {
1374 LimVib = nRot_hi[iElecHi][iVibHi];
1375 }
1376 LimVib = MIN2( LimVib , nVib_hi[iElecHi] );
1377 /* this is limit to rotation quantum index */
1378 if( save.punarg[ipPun][1] > 0 )
1379 {
1380 LimRot = (long)save.punarg[ipPun][1];
1381 }
1382 else
1383 {
1384 LimRot = nRot_hi[iElecHi][iVibHi];
1385 }
1386 for( iVibHi = 0; iVibHi<=LimVib; ++iVibHi )
1387 {
1388 fprintf(io,"\tv=%li",iVibHi);
1389 long int LimRotVib = MIN2( LimRot , nRot_hi[iElecHi][iVibHi] );
1390 for( long iRotHi=Jlowest[iElecHi]; iRotHi<=LimRotVib; ++iRotHi )
1391 {
1392 fprintf(io,"\t%.3e",
1393 states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].Pop()/(*dense_total) );
1394 }
1395 }
1396 fprintf(io,"\n");
1397 }
1398 }
1399
1400 /* save column densities */
1401 else if( (strcmp( chJOB , "H2cl" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1402 {
1403 long iVibHi= 0;
1404 long iRotHi = 0;
1405 long iElecHi=0;
1406 long LimVib, LimRot;
1407 /* the limit to the number of vibration levels punched -
1408 * default is all, but first two numbers on save h2 pops command
1409 * reset limit */
1410 /* this is limit to vibration */
1411 if( save.punarg[ipPun][0] > 0 )
1412 {
1413 LimVib = (long)save.punarg[ipPun][0];
1414 }
1415 else
1416 {
1417 LimVib = nVib_hi[iElecHi];
1418 }
1419
1420 /* first save ortho and para populations */
1421 fprintf(io,"%i\t%i\t%.3e\tortho\n",
1422 103 ,
1423 103 ,
1424 ortho_colden );
1425 fprintf(io,"%i\t%i\t%.3e\tpara\n",
1426 101 ,
1427 101 ,
1428 para_colden );
1429 /* total H2 column density */
1430 fprintf(io,"%i\t%i\t%.3e\ttotal\n",
1431 0 ,
1432 0 ,
1433 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]);
1434
1435 /* save level column densities */
1436 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
1437 {
1438 if( lgEnabled )
1439 {
1440 /* this is limit to rotation quantum index */
1441 if( save.punarg[ipPun][1] > 0 )
1442 {
1443 LimRot = (long)save.punarg[ipPun][1];
1444 }
1445 else
1446 {
1447 LimRot = nRot_hi[iElecHi][iVibHi];
1448 }
1449 if( save.punarg[ipPun][2] > 0 )
1450 {
1451 long int i;
1452 /* save matrix */
1453 if( iVibHi == 0 )
1454 {
1455 fprintf(io,"vib\\rot");
1456 /* this is first vib, so make row of rot numbs */
1457 for( i=0; i<=LimRot; ++i )
1458 {
1459 fprintf(io,"\t%li",i);
1460 }
1461 fprintf(io,"\n");
1462 }
1463 fprintf(io,"%li",iVibHi );
1464 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1465 {
1466 fprintf(io,"\t%.3e",
1467 H2_X_colden[iVibHi][iRotHi]/(*dense_total) );
1468 }
1469 fprintf(io,"\n" );
1470 }
1471 else
1472 {
1473 /* save triplets - the default */
1474 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
1475 {
1476 fprintf(io,"%li\t%li\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\n",
1477 iVibHi ,
1478 iRotHi ,
1479 /* energy relative to 0,0, T1CM converts wavenumber to K */
1480 states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].energy().K(),
1481 /* these are column densities for actual molecule */
1482 H2_X_colden[iVibHi][iRotHi] ,
1483 H2_X_colden[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi] ,
1484 /* these are same column densities but for LTE populations */
1485 H2_X_colden_LTE[iVibHi][iRotHi] ,
1486 H2_X_colden_LTE[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi]);
1487 }
1488 }
1489 }
1490 }
1491 }
1492 else if( (strcmp(chJOB , "H2pd" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1493 {
1494 /* save PDR
1495 * output some PDR information (densities, rates) for each zone */
1496 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1497 /* depth in cm */
1498 radius.depth_mid_zone ,
1499 /* the computed ortho and para densities */
1500 ortho_density ,
1501 para_density ,
1502 /* the Lyman Werner band dissociation, Tielens & Hollenbach */
1503 hmi.H2_Solomon_dissoc_rate_TH85_H2g ,
1504 /* the Lyman Werner band dissociation, Bertoldi & Draine */
1505 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
1506 /* the Lyman Werner band dissociation, big H2 mole */
1508 }
1509 else if( (strcmp(chJOB , "H2co" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1510 {
1511 /* save H2 cooling - do heating cooling for each zone old new H2 */
1512 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1513 /* depth in cm */
1514 radius.depth_mid_zone ,
1515 /* total cooling, equal to total heating */
1516 thermal.ctot ,
1517 /* H2 destruction by Solomon process, TH85 rate */
1518 hmi.H2_Solomon_dissoc_rate_TH85_H2g,
1519 /* H2 destruction by Solomon process, big H2 model rate */
1522 /* H2 photodissociation heating, eqn A9 of Tielens & Hollenbach 1985a */
1523 hmi.HeatH2Dish_TH85,
1524 /* heating due to dissociation of electronic excited states */
1525 HeatDiss ,
1526 /* cooling (usually neg and so heating) due to collisions within X */
1527 hmi.HeatH2Dexc_TH85,
1528 HeatDexc
1529 );
1530
1531 }
1532
1533 else if( (strcmp(chJOB , "H2le" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1534 {
1535 /* save H2 levels */
1536 for( long int ipHi=0; ipHi < nLevels_per_elec[0]; ipHi++ )
1537 {
1538 long iRotHi = ipRot_H2_energy_sort[ipHi];
1539 long iVibHi = ipVib_H2_energy_sort[ipHi];
1540 long int nColl;
1541 double Asum , Csum[N_X_COLLIDER];
1542 Asum = 0;
1543 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1544 Csum[nColl] = 0.;
1545 for( long int ipLo=0; ipLo<ipHi; ++ipLo )
1546 {
1547 /* all lower levels */
1548 long iRotLo = ipRot_H2_energy_sort[ipLo];
1549 long iVibLo = ipVib_H2_energy_sort[ipLo];
1550 EmissionProxy em = trans[ ipTransitionSort[ipHi][ipLo] ].Emis();
1551
1552 /* radiative decays down */
1553 if( ( abs(iRotHi-iRotLo) == 2 || (iRotHi-iRotLo) == 0 ) && iVibLo <= iVibHi &&
1554 lgH2_radiative[ipHi][ipLo] )
1555 {
1556 Asum += em.Aul() * ( em.Pesc() + em.Pdest() + em.Pelec_esc() );
1557 }
1558 /* all collisions down */
1559 mr3ci H2cr = CollRateCoeff.begin(ipHi,ipLo);
1560 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1561 Csum[nColl] += H2cr[nColl];
1562 }
1563
1564 /* save H2 level energies */
1565 fprintf(io,"%li\t%li\t%.2f\t%li\t%.3e",
1566 iVibHi , iRotHi,
1567 states[ipHi].energy().WN(),
1568 (long)states[ipHi].g(),
1569 Asum );
1570 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
1571 /* sum over all lower levels */
1572 fprintf(io,"\t%.3e",Csum[nColl]);
1573 fprintf(io,"\n");
1574 }
1575 }
1576
1577 else if( (strcmp(chJOB , "H2ra" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1578 {
1579 /* save h2 rates - some rates and lifetimes */
1580 double sumpop = 0. , sumlife = 0.;
1581
1582 /* this block, find lifetime against photo excitation into excited electronic states */
1583 if( lgEnabled && lgEvaluated )
1584 {
1585 /* only do if radiative transition exists */
1586 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1587 {
1588 qList::iterator Lo = ( (*tr).Lo() );
1589 if( (*Lo).n() > 0 || (*Lo).v() > 0 )
1590 continue;
1591 sumlife += (*tr).Emis().pump() * (*(*tr).Lo()).Pop();
1592 sumpop += (*(*tr).Lo()).Pop();
1593 }
1594 }
1595
1596 /* continue output from save h2 rates command */
1597 /* find photoexcitation rates from v=0 */
1598 /* PDR information for each zone */
1599 fprintf(io,
1600 "%.5e\t%.3e\t%.3e\t%.3e\t%li",
1601 /* depth in cm */
1602 radius.depth_mid_zone ,
1603 /* the column density (cm^-2) in H2 */
1604 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
1605 /* this is a special form of column density - should be proportional
1606 * to total shielding */
1607 colden.coldenH2_ov_vel ,
1608 /* visual extinction due to dust alone, of point source (star)*/
1609 rfield.extin_mag_V_point,
1610 /* number of large molecule evaluations in this zone */
1612 fprintf(io,
1613 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1614 /* total H2 fraction */
1615 (*dense_total)/dense.gas_phase[ipHYDROGEN] ,
1616 /* chemistry renorm factor */
1618 /* rate H2 forms on grains */
1619 gv.rate_h2_form_grains_used_total ,
1620 /* rate H2 forms by H minus route */
1621 findspecieslocal("H-")->den*1.35e-9,
1622 /* H2 destruction by Solomon process, TH85 rate */
1623 hmi.H2_Solomon_dissoc_rate_TH85_H2g + hmi.H2_Solomon_dissoc_rate_TH85_H2s,
1624 /* H2 destruction by Solomon process, Bertoldi & Draine rate */
1625 hmi.H2_Solomon_dissoc_rate_BD96_H2g + hmi.H2_Solomon_dissoc_rate_BD96_H2s,
1626 /* H2 destruction by Solomon process, Elwert et al. in preparation */
1627 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + hmi.H2_Solomon_dissoc_rate_ELWERT_H2g,
1628 /* H2 destruction by Solomon process, big H2 model rate */
1630 /* rate s-1 H2 electronic excit states decay into H2g */
1632 /* rate s-1 H2 electronic excit states decay into H2s */
1634 );
1635 fprintf(io,
1636 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
1637 /* The TH85 estimate of the radiation field relative to the Habing value */
1638 hmi.UV_Cont_rel2_Habing_TH85_depth,
1639 /* The DB96 estimate of the radiation field relative to the Habing value */
1640 hmi.UV_Cont_rel2_Draine_DB96_depth,
1641 /* cosmic ray ionization rate */
1642 secondaries.csupra[ipHYDROGEN][0]*0.93,
1643 sumlife/SDIV( sumpop ) ,
1644 hmi.H2_Solomon_dissoc_rate_BD96_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth) ,
1645 Solomon_dissoc_rate_g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth),
1646 Solomon_dissoc_rate_g/SDIV(hmi.UV_Cont_rel2_Habing_spec_depth),
1647 hmi.H2_rate_destroy);
1648 fprintf(io,
1649 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1650 hmi.HeatH2Dish_TH85,
1651 hmi.HeatH2Dexc_TH85,
1652 HeatDiss,
1653 HeatDexc,
1654 thermal.htot);
1655 }
1656 /* save h2 solomon */
1657 else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1658 {
1659 /* remember as many as nSOL lines contributing to total Solomon process */
1660 const int nSOL = 100;
1661 double sum, one;
1662 long int jlosave[nSOL] , ivlosave[nSOL],
1663 iehisave[nSOL] ,jhisave[nSOL] , ivhisave[nSOL],
1664 nsave,
1665 ipOrdered[nSOL];
1666 int nFail;
1667 realnum fsave[nSOL], wlsave[nSOL];
1668 /* Solomon process, and where it came from */
1669 fprintf(io,"%.5e\t%.3e",
1670 /* depth in cm */
1671 radius.depth_mid_zone ,
1672 /* H2 destruction by Solomon process, big H2 model rate */
1675 sum = 0.;
1676 /* find sum of all radiative exits from X into excited electronic states */
1677 if( lgEnabled && lgEvaluated )
1678 {
1679 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1680 {
1681 qList::iterator Lo = ( (*tr).Lo() );
1682 if( (*Lo).n() > 0 )
1683 continue;
1684 sum += (*(*tr).Lo()).Pop() * (*tr).Emis().pump();
1685 }
1686
1687 /* make sure it is safe to div by sum */
1688 sum = SDIV( sum );
1689 nsave = 0;
1690 /* now loop back over X and print all those which contribute more than frac of the total */
1691 const double frac = 0.01;
1692 /* only do if radiative transition exists */
1693 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1694 {
1695 qList::iterator Lo = ( (*tr).Lo() );
1696 if( (*Lo).n() > 0 )
1697 continue;
1698 one = (*(*tr).Lo()).Pop() * (*tr).Emis().pump();
1699 if( one/sum > frac && nsave<nSOL)
1700 {
1701 qList::iterator Hi = ( (*tr).Hi() );
1702 fsave[nsave] = (realnum)(one/sum);
1703 jlosave[nsave] = (*Lo).J();
1704 ivlosave[nsave] = (*Lo).v();
1705 jhisave[nsave] = (*Hi).J();
1706 ivhisave[nsave] = (*Hi).v();
1707 iehisave[nsave] = (*Hi).n();
1708 wlsave[nsave] = (*tr).WLAng();
1709 ++nsave;
1710 }
1711 }
1712
1713 /* now sort by decreasing importance */
1714 /*spsort netlib routine to sort array returning sorted indices */
1715 spsort(
1716 /* input array to be sorted */
1717 fsave,
1718 /* number of values in x */
1719 nsave,
1720 /* permutation output array */
1721 ipOrdered,
1722 /* flag saying what to do - 1 sorts into increasing order, not changing
1723 * the original routine */
1724 -1,
1725 /* error condition, should be 0 */
1726 &nFail);
1727
1728 /* print ratio of pumps to dissociations - this is 9:1 in TH85 */
1729 /*>>chng 05 jul 20, TE, save average energy in H2s and renormalization factors for H2g and H2s */
1730 /* >>chng 05 sep 16, TE, chng denominator to do g and s with proper dissoc rates */
1731 fprintf(io,"\t%.3f\t%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e",
1732 /* this is sum of photons and CRs */
1733 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f)/SDIV((Solomon_dissoc_rate_g * findspecieslocal("H2")->den +
1734 Solomon_dissoc_rate_s * findspecieslocal("H2*")->den) ),
1735 /* this is sum of photons and CRs */
1736 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f) /SDIV((Solomon_dissoc_rate_g * H2_den_g+
1739 findspecieslocal("H2")->den/SDIV(H2_den_g), findspecieslocal("H2*")->den/SDIV(H2_den_s),
1741 );
1742 for( long i=0; i<nsave; ++i )
1743 {
1744 long ip = ipOrdered[i];
1745 /*lint -e644 not init */
1746 fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f\t%.3f",
1747 iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[ip] , wlsave[ip] );
1748 /*lint +e644 not init */
1749 }
1750 fprintf(io,"\n");
1751 }
1752 }
1753
1754 else if( (strcmp(chJOB , "H2te" ) == 0) && (strcmp(chTime,"LAST") != 0) )
1755 {
1756 /* save h2 temperatures */
1757 double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40;
1758 double T10,T20,T30,T31,T40;
1759 /* subscript"sum" denotes integrated quantities */
1760 double T10_sum,T20_sum,T30_sum,T31_sum,T40_sum;
1761 double pop_ratio10_sum,pop_ratio20_sum,pop_ratio30_sum,pop_ratio31_sum,pop_ratio40_sum;
1762 if( lgEnabled && nCall_this_zone )
1763 {
1764 double pop0 = states[0].Pop();
1765 double pop1 = states[1].Pop();
1766 double pop2 = states[2].Pop();
1767 double pop3 = states[3].Pop();
1768 double pop4 = states[4].Pop();
1769
1770 double energyK = states[1].energy().K() - states[0].energy().K();
1771 /* the ratio of populations of J=1 to 0 */
1772 pop_ratio10 = pop1/SDIV(pop0);
1773 pop_ratio10_sum = H2_X_colden[0][1]/SDIV(H2_X_colden[0][0]);
1774 /* the corresponding temperature */
1775 T10 = -170.5/log(SDIV(pop_ratio10) * H2_stat[0][0][0]/H2_stat[0][0][1]);
1776 T10_sum = -170.5/log(SDIV(pop_ratio10_sum) * H2_stat[0][0][0]/H2_stat[0][0][1]);
1777
1778 energyK = states[2].energy().K() - states[0].energy().K();
1779 pop_ratio20 = pop2/SDIV(pop0);
1780 T20 = -energyK/log(SDIV(pop_ratio20) * H2_stat[0][0][0]/H2_stat[0][0][2]);
1781
1782 pop_ratio20_sum = H2_X_colden[0][2]/SDIV(H2_X_colden[0][0]);
1783 T20_sum = -energyK/log(SDIV(pop_ratio20_sum) * H2_stat[0][0][0]/H2_stat[0][0][2]);
1784
1785 energyK = states[3].energy().K() - states[0].energy().K();
1786 pop_ratio30 = pop3/SDIV(pop0);
1787 T30 = -energyK/log(SDIV(pop_ratio30) * H2_stat[0][0][0]/H2_stat[0][0][3]);
1788
1789 pop_ratio30_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][0]);
1790 T30_sum = -energyK/log(SDIV(pop_ratio30_sum) * H2_stat[0][0][0]/H2_stat[0][0][3]);
1791
1792 energyK = states[3].energy().K() - states[1].energy().K();
1793 pop_ratio31 = pop3/SDIV(pop1);
1794 T31 = -energyK/log(SDIV(pop_ratio31) * H2_stat[0][0][1]/H2_stat[0][0][3]);
1795
1796 pop_ratio31_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][1]);
1797 T31_sum = -energyK/log(SDIV(pop_ratio31_sum) * H2_stat[0][0][1]/H2_stat[0][0][3]);
1798
1799 energyK = states[4].energy().K() - states[0].energy().K();
1800 pop_ratio40 = pop4/SDIV(pop0);
1801 T40 = -energyK/log(SDIV(pop_ratio40) * H2_stat[0][0][0]/H2_stat[0][0][4]);
1802
1803 pop_ratio40_sum = H2_X_colden[0][4]/SDIV(H2_X_colden[0][0]);
1804 T40_sum = -energyK/log(SDIV(pop_ratio40_sum) * H2_stat[0][0][0]/H2_stat[0][0][4]);
1805 }
1806 else
1807 {
1808 pop_ratio10 = 0.;
1809 pop_ratio10_sum = 0.;
1810 T10 = 0.;
1811 T20 = 0.;
1812 T30 = 0.;
1813 T31 = 0.;
1814 T40 = 0.;
1815 T10_sum = 0.;
1816 T20_sum = 0.;
1817 T30_sum = 0.;
1818 T31_sum = 0.;
1819 T40_sum = 0.;
1820 }
1821
1822 /* various temperatures for neutral/molecular gas */
1823 fprintf( io,
1824 "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n" ,
1825 /* depth in cm */
1826 radius.depth_mid_zone ,
1827 /* total H2 fraction */
1828 (*dense_total)/dense.gas_phase[ipHYDROGEN] ,
1829 /* ratio of populations of 1 to 0 only */
1830 pop_ratio10 ,
1831 /* sum of all ortho and para */
1833 T10,T20,T30,T31,T40,
1834 phycon.te ,
1835 hyperfine.Tspin21cm,T10_sum,T20_sum,T30_sum,T31_sum,T40_sum );
1836 }
1837 else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) )
1838 {
1839 /* save H2 lines - output the full emission-line spectrum */
1840 double thresh;
1841 double renorm;
1842 /* first test, is H2 turned on? Second test, have lines arrays
1843 * been set up - nsum is negative if abort occurs before lines
1844 * are set up */
1845 if( lgEnabled && LineSave.nsum > 0)
1846 {
1847 ASSERT( LineSave.ipNormWavL >= 0 );
1848 /* get the normalization line */
1849 if( LineSv[LineSave.ipNormWavL].SumLine[0] > SMALLFLOAT )
1850 renorm = LineSave.ScaleNormLine/
1851 LineSv[LineSave.ipNormWavL].SumLine[0];
1852 else
1853 renorm = 1.;
1854
1855 if( renorm > SMALLFLOAT )
1856 {
1857 /* this is threshold for faintest line, normally 0, set with
1858 * number on save H2 command */
1859 thresh = thresh_punline_h2/(realnum)renorm;
1860 }
1861 else
1862 thresh = 0.f;
1863
1864 /* save H2 line intensities at end of iteration
1865 * nElecLevelOutput is electronic level with 1 for ground, so this loop is < nElecLevelOutput */
1866 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
1867 {
1868 qList::iterator Hi = ( (*tr).Hi() );
1869 qList::iterator Lo = ( (*tr).Lo() );
1870 long iElecHi = (*Hi).n();
1871 long iVibHi = (*Hi).v();
1872 long iRotHi = (*Hi).J();
1873 long iElecLo = (*Lo).n();
1874 long iVibLo = (*Lo).v();
1875 long iRotLo = (*Lo).J();
1876 if( iElecHi >= nElecLevelOutput )
1877 continue;
1878 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > thresh )
1879 {
1880 /* air wavelength in microns */
1881 /* WLAng contains correction for index of refraction of air */
1882 double wl = (*tr).WLAng()/1e4;
1883 fprintf(io, "%li-%li %c(%li)", iVibHi, iVibLo, chMolBranch( iRotHi, iRotLo ), iRotLo );
1884 fprintf( io, "\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld", iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo);
1885 /* WLAng contains correction for index of refraction of air */
1886 fprintf( io, "\t%.7f\t", wl );
1887 /*prt_wl print floating wavelength in Angstroms, in output format */
1888 prt_wl( io , (*tr).WLAng() );
1889 /* the log of the line intensity or luminosity */
1890 fprintf( io, "\t%.3f\t%.3e",
1891 log10(MAX2(1e-37,H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo])) + radius.Conv2PrtInten,
1892 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]*renorm );
1893 /* excitation energy of upper level in K */
1894 fprintf( io, "\t%.3f", (*Hi).energy().K() );
1895 /* the product g_hi h nu * Aul */
1896 fprintf( io, "\t%.3e", (*tr).Emis().Aul() * (*tr).EnergyErg() * (*(*tr).Hi()).g() );
1897 fprintf( io, "\n");
1898 }
1899 }
1900 }
1901 }
1902 else if( (strcmp(chJOB , "H2sp" ) == 0) )
1903 {
1904 fprintf( io, "PUT SOMETHING HERE!\n" );
1905 }
1906 return;
1907}
1908 /*cdH2_Line determines intensity and luminosity of and H2 line. The first
1909 * six arguments give the upper and lower quantum designation of the levels.
1910 * The function returns 1 if we found the line,
1911 * and false==0 if we did not find the line because ohoto-para transition
1912 * or upper level has lower energy than lower level */
1913long int cdH2_Line(
1914 /* indices for the upper level */
1915 long int iElecHi,
1916 long int iVibHi ,
1917 long int iRotHi ,
1918 /* indices for lower level */
1919 long int iElecLo,
1920 long int iVibLo ,
1921 long int iRotLo ,
1922 /* linear intensity relative to normalization line*/
1923 double *relint,
1924 /* log of luminosity or intensity of line */
1925 double *absint )
1926{
1927 DEBUG_ENTRY( "cdH2_Line()" );
1928 return h2.getLine( iElecHi, iVibHi, iRotHi, iElecLo, iVibLo, iRotLo, relint, absint );
1929};
1930
1931long int diatomics::getLine( long iElecHi, long iVibHi, long iRotHi, long iElecLo, long iVibLo, long iRotLo, double *relint, double *absint )
1932{
1933
1934 DEBUG_ENTRY( "diatomics::getline()" );
1935
1936 /* these will be return values if we can't find the line */
1937 *relint = 0.;
1938 *absint = 0.;
1939
1940 /* for now both electronic levels must be zero */
1941 if( iElecHi!=0 || iElecLo!=0 )
1942 {
1943 return 0;
1944 }
1945
1946 long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
1947 long ipLo = ipEnergySort[iElecLo][iVibLo][iRotLo];
1948
1949 /* check that energy of first level is higher than energy of second level */
1950 if( states[ipHi].energy().WN() < states[ipLo].energy().WN() )
1951 {
1952 return 0;
1953 }
1954
1955 /* check that ortho-para does not change */
1956 if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] - H2_lgOrtho[iElecLo][iVibLo][iRotLo] != 0 )
1957 {
1958 return 0;
1959 }
1960
1961 /* exit if lines does not exist */
1962 if( !lgH2_radiative[ipHi][ipLo] )
1963 {
1964 return 0;
1965 }
1966
1967 ASSERT( LineSave.ipNormWavL >= 0 );
1968 /* does the normalization line have a positive intensity*/
1969 if( LineSv[LineSave.ipNormWavL].SumLine[0] > 0. )
1970 {
1971 *relint = H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]/
1972 LineSv[LineSave.ipNormWavL].SumLine[0] * LineSave.ScaleNormLine;
1973 }
1974 else
1975 {
1976 *relint = 0.;
1977 }
1978
1979 /* return log of line intensity if it is positive */
1980 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > 0. )
1981 {
1982 *absint = log10(H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]) +
1983 radius.Conv2PrtInten;
1984 }
1985 else
1986 {
1987 /* line intensity is actually zero, return small number */
1988 *absint = -37.;
1989 }
1990 /* this indicates success */
1991 return 1;
1992}
1993
1995{
1996 if( !lgREAD_DATA )
1997 nXLevelsMatrix = numLevels;
1998}
1999
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define MIN2
Definition cddefines.h:761
NORETURN void BadRead(void)
Definition service.cpp:901
#define EXIT_FAILURE
Definition cddefines.h:140
double RandGauss(double xMean, double s)
Definition service.cpp:1643
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
void spsort(realnum x[], long int n, long int iperm[], int kflag, int *ier)
Definition service.cpp:1100
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
LinSv * LineSv
Definition cdinit.cpp:70
realnum & Pesc() const
Definition emission.h:523
realnum & Pelec_esc() const
Definition emission.h:533
realnum & Aul() const
Definition emission.h:613
realnum & Pdest() const
Definition emission.h:543
bool nMatch(const char *chKey) const
Definition parser.h:135
double getNumberDefault(const char *chDesc, double fdef)
Definition parser.cpp:282
double getNumberDefaultNegImplLog(const char *chDesc, double fdef)
Definition parser.cpp:336
TransitionProxy::iterator iterator
Definition transition.h:280
EmissionList & Emis()
Definition transition.h:329
void H2_PunchDo(FILE *io, char chJOB[], const char chTime[], long int ipPun)
long int n_elec_states
Definition h2_priv.h:409
char chH2ColliderLabels[N_X_COLLIDER][chN_X_COLLIDER]
Definition h2_priv.h:591
double para_density
Definition h2_priv.h:321
multi_arr< double, 2 > H2_col_rate_out
Definition h2_priv.h:652
double Solomon_elec_decay_s
Definition h2_priv.h:269
multi_arr< bool, 3 > H2_lgOrtho
Definition h2_priv.h:643
double H2_renorm_chemistry
Definition h2_priv.h:603
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
double pops_per_elec[N_ELEC]
Definition h2_priv.h:620
long int getLine(long iElecHi, long iVibHi, long iRotHi, long iElecLo, long iVibLo, long iRotLo, double *relint, double *absint)
bool lgEnabled
Definition h2_priv.h:345
vector< double > depart
Definition h2_priv.h:702
double Solomon_elec_decay_g
Definition h2_priv.h:268
multi_arr< realnum, 3 > H2_dissprob
Definition h2_priv.h:632
string path
Definition h2_priv.h:573
void set_numLevelsMatrix(long numLevels)
long int nXLevelsMatrix
Definition h2_priv.h:695
void H2_ReadDissprob(long int nelec)
double ortho_colden
Definition h2_priv.h:328
qList states
Definition h2_priv.h:565
string label
Definition h2_priv.h:571
void H2_Read_hminus_distribution(void)
void H2_ParseSave(Parser &p, char *chHeader)
long int Jlowest[N_ELEC]
Definition h2_priv.h:616
multi_arr< long int, 2 > ipTransitionSort
Definition h2_priv.h:691
int nElecLevelOutput
Definition h2_priv.h:349
double average_energy_g
Definition h2_priv.h:286
void H2_Prt_Zone(void)
valarray< long > ipVib_H2_energy_sort
Definition h2_priv.h:687
double HeatDiss
Definition h2_priv.h:289
const double *const dense_total
Definition h2_priv.h:589
void H2_ReadEnergies()
void H2_Prt_line_tau(void)
multi_arr< realnum, 3 > H2_stat
Definition h2_priv.h:641
double Solomon_dissoc_rate_g
Definition h2_priv.h:264
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
bool lgEvaluated
Definition h2_priv.h:310
multi_arr< double, 3 > H2_populations_LTE
Definition h2_priv.h:639
double H2_den_s
Definition h2_priv.h:682
valarray< long > ipRot_H2_energy_sort
Definition h2_priv.h:689
multi_arr< realnum, 3 > CollRateCoeff
Definition h2_priv.h:621
void H2_Prt_column_density(FILE *ioMEAN)
void H2_LinesAdd(void)
void H2_ReadDissocEnergies(void)
void H2_PrtDepartCoef(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
double HeatDexc
Definition h2_priv.h:290
double H2_den_g
Definition h2_priv.h:682
double ortho_density
Definition h2_priv.h:319
multi_arr< realnum, 2 > H2_X_colden_LTE
Definition h2_priv.h:662
void H2_PunchLineStuff(FILE *io, realnum xLimit, long index)
multi_arr< realnum, 2 > H2_X_colden
Definition h2_priv.h:660
double para_colden
Definition h2_priv.h:329
TransitionList trans
Definition h2_priv.h:566
double average_energy_s
Definition h2_priv.h:287
double Solomon_dissoc_rate_s
Definition h2_priv.h:265
long int nCall_this_zone
Definition h2_priv.h:341
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
long int nLevels_per_elec[N_ELEC]
Definition h2_priv.h:618
void H2_Punch_line_data(FILE *ioPUN, bool lgDoAll)
multi_arr< long int, 3 > ipEnergySort
Definition h2_priv.h:690
bool operator<(const level_tmp &second) const
ProxyIterator< qStateConstProxy, qStateConstProxy > const_iterator
ProxyIterator< qStateProxy, qStateConstProxy > iterator
t_colden colden
Definition colden.cpp:5
#define ipCOL_H2s
Definition colden.h:18
#define ipCOL_H2g
Definition colden.h:16
multi_arr< realnum, 3 >::const_iterator mr3ci
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
realnum GetDopplerWidth(realnum massAMU)
GrainVar gv
Definition grainvar.cpp:5
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
const int N_ELEC
Definition h2_priv.h:27
const int nTE_HMINUS
Definition h2_priv.h:24
const int chN_X_COLLIDER
Definition h2_priv.h:15
const int N_X_COLLIDER
Definition h2_priv.h:13
t_hmi hmi
Definition hmi.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_input input
Definition input.cpp:12
t_LineSave LineSave
Definition lines.cpp:5
void lindst(double xInten, realnum wavelength, const char *chLab, long int ipnt, char chInfo, bool lgOutToo, const char *chComment)
molezone * findspecieslocal(const char buf[])
STATIC char chMolBranch(long iRotHi, long int iRotLo)
long int cdH2_Line(long int iElecHi, long int iVibHi, long int iRotHi, long int iElecLo, long int iVibLo, long int iRotLo, double *relint, double *absint)
static realnum thresh_punline_h2
t_phycon phycon
Definition phycon.cpp:6
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
void prme(const bool lgReset, const TransitionProxy &t)
Definition prt_met.cpp:97
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_save save
Definition save.cpp:5
void Save1Line(const TransitionProxy &t, FILE *io, realnum xLimit, long index, realnum DopplerWidth)
Definition save_do.cpp:4347
void Save1LineData(const TransitionProxy &t, FILE *io, bool lgCS_2, bool &lgPrint)
t_secondaries secondaries
static double * g
Definition species2.cpp:28
t_thermal thermal
Definition thermal.cpp:5
void PutLine(const TransitionProxy &t, const char *chComment, const char *chLabelTemp)