00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "cddefines.h"
00021 #include "physconst.h"
00022 #include "save.h"
00023 #include "hmi.h"
00024 #include "prt.h"
00025 #include "secondaries.h"
00026 #include "grainvar.h"
00027 #include "input.h"
00028 #include "phycon.h"
00029 #include "rfield.h"
00030 #include "hyperfine.h"
00031 #include "thermal.h"
00032 #include "lines.h"
00033 #include "lines_service.h"
00034 #include "mole.h"
00035 #include "dense.h"
00036 #include "radius.h"
00037 #include "colden.h"
00038 #include "taulines.h"
00039 #include "h2.h"
00040 #include "h2_priv.h"
00041 #include "cddrive.h"
00042 #include "doppvel.h"
00043 #include "doppvel.h"
00044 #include "parser.h"
00045
00046 static realnum thresh_punline_h2;
00047
00048
00049 void diatomics::H2_LinesAdd(void)
00050 {
00051
00052 if( !lgEnabled )
00053 return;
00054
00055 DEBUG_ENTRY( "H2_LinesAdd()" );
00056
00057 if( strcmp( "H2 ", label.c_str() ) == 0 )
00058 {
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][6] ][ ipEnergySort[0][0][4] ] ], "H2 ", 'i', false, "H2 line");
00073
00074 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][5] ][ ipEnergySort[0][0][3] ] ], "H2 ", 'i', false, "H2 line");
00075
00076 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][4] ][ ipEnergySort[0][0][2] ] ], "H2 ", 'i', false, "H2 line");
00077
00078 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][3] ][ ipEnergySort[0][0][1] ] ], "H2 ", 'i', false, "H2 line");
00079
00080 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][2] ][ ipEnergySort[0][0][0] ] ], "H2 ", 'i', false, "H2 line");
00081
00082
00083 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][2] ][ ipEnergySort[0][0][2] ] ], "H2 ", 'i', false, "H2 line");
00084
00085 lindst( trans[ ipTransitionSort[ ipEnergySort[0][1][1] ][ ipEnergySort[0][0][1] ] ], "H2 ", 'i', false, "H2 line");
00086 }
00087
00088
00089 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00090 {
00091 qList::iterator Hi = ( (*tr).Hi() );
00092 if( (*Hi).n() >= nElecLevelOutput ) continue;
00093 qList::iterator Lo = ( (*tr).Lo() );
00094
00095 PutLine( *tr, "diatoms lines", label.c_str() );
00096 if( LineSave.ipass == 0 )
00097 {
00098 H2_SaveLine[(*Hi).n()][(*Hi).v()][(*Hi).J()][(*Lo).n()][(*Lo).v()][(*Lo).J()] = 0.;
00099 }
00100 else if( LineSave.ipass == 1 )
00101 {
00102 H2_SaveLine[(*Hi).n()][(*Hi).v()][(*Hi).J()][(*Lo).n()][(*Lo).v()][(*Lo).J()] +=
00103 (realnum)( radius.dVeffAper * (*tr).Emis().xIntensity() );
00104 }
00105 }
00106
00107 return;
00108 }
00109
00110
00111 void diatomics::H2_ParseSave( Parser &p ,
00112 char *chHeader)
00113 {
00114 DEBUG_ENTRY( "H2_ParseSave()" );
00115
00116 save.whichDiatomToPrint[save.nsave] = &(*this);
00117
00118
00119 if( p.nMatch("COLU") )
00120 {
00121
00122 strcpy( save.chSave[save.nsave], "H2cl" );
00123
00124
00125
00126
00127
00128 save.punarg[save.nsave][0] = (realnum)p.getNumberDefault(
00129 "H2 vibration state",0.0);
00130
00131
00132 save.punarg[save.nsave][1] = (realnum)p.getNumberDefault(
00133 "H2 rotation state",0.0);
00134
00135
00136 if( p.nMatch( "MATR" ) )
00137 {
00138
00139 save.punarg[save.nsave][2] = 1;
00140 sprintf( chHeader, "#vib\trot\tcolumn density\n" );
00141 }
00142 else
00143 {
00144
00145 save.punarg[save.nsave][2] = -1;
00146 sprintf( chHeader, "#vib\trot\tEner(K)\tcolden\tcolden/stat wght\tLTE colden\tLTE colden/stat wght\n" );
00147 }
00148 }
00149 else if( p.nMatch("COOL") )
00150 {
00151
00152 strcpy( save.chSave[save.nsave], "H2co" );
00153 sprintf( chHeader,
00154 "#H2 depth\ttot cool\tTH Sol\tBig Sol\tTH pht dis\tpht dis\tTH Xcool\tXcool \n" );
00155 }
00156
00157 else if( p.nMatch("CREA") )
00158 {
00159
00160 fprintf( ioQQQ, " This command has been superseded by the \"creation\" option of the \"save chemistry rates\" command.\n" );
00161 fprintf( ioQQQ, " Sorry.\n" );
00162 cdEXIT(EXIT_FAILURE);
00163 }
00164 else if( p.nMatch("DEST") )
00165 {
00166
00167 fprintf( ioQQQ, " This command has been superseded by the \"destruction\" option of the \"save chemistry rates\" command.\n" );
00168 fprintf( ioQQQ, " Sorry.\n" );
00169 cdEXIT(EXIT_FAILURE);
00170 }
00171
00172 else if( p.nMatch("HEAT") )
00173 {
00174
00175 strcpy( save.chSave[save.nsave], "H2he" );
00176 sprintf( chHeader,
00177 "#H2 depth\ttot Heat\tHeat(big)\tHeat(TH85)\tDissoc(Big)\tDissoc(TH85) \n" );
00178 }
00179
00180 else if( p.nMatch("LEVE") )
00181 {
00182
00183 strcpy( save.chSave[save.nsave], "H2le" );
00184 sprintf( chHeader,
00185 "#H2 v\tJ\tenergy(wn)\tstat wght\tSum As" );
00186 char chHoldit[chN_X_COLLIDER+12];
00187 for( int nColl=0; nColl<N_X_COLLIDER; ++nColl )
00188 {
00189
00190 sprintf(chHoldit,"\tCritDen %s",chH2ColliderLabels[nColl]);
00191 strcat( chHeader , chHoldit );
00192 }
00193 strcat( chHeader , "\n" );
00194 }
00195
00196 else if( p.nMatch("LINE") )
00197 {
00198
00199 strcpy( save.chSave[save.nsave], "H2ln" );
00200 sprintf( chHeader,
00201 "#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" );
00202
00203
00204
00205
00206
00207 thresh_punline_h2 = (realnum)p.getNumberDefaultNegImplLog(
00208 "faintest line to save",1e-4);
00209
00210
00211
00212
00213
00214 if( p.nMatch( "ELEC" ) )
00215 {
00216 if( p.nMatch(" ALL") )
00217 {
00218
00219
00220
00221 nElecLevelOutput = -1;
00222 }
00223 else if( p.nMatch("GROU") )
00224 {
00225
00226 nElecLevelOutput = 1;
00227 }
00228 else
00229 {
00230 nElecLevelOutput = (int)p.getNumberDefault(
00231 "electronic levels for output",1.0);
00232 }
00233 }
00234 }
00235
00236 else if( p.nMatch(" PDR") )
00237 {
00238
00239 strcpy( save.chSave[save.nsave], "H2pd" );
00240 sprintf( chHeader, "#H2 creation, destruction. \n" );
00241 }
00242 else if( p.nMatch("POPU") )
00243 {
00244
00245 strcpy( save.chSave[save.nsave], "H2po" );
00246
00247
00248
00249
00250
00251 save.punarg[save.nsave][0] = (realnum)p.getNumberDefault(
00252 "highest H2 save vibration state",0.0);
00253
00254
00255 save.punarg[save.nsave][1] = (realnum)p.getNumberDefault(
00256 "highest H2 save rotation state",0.0);
00257
00258 if( p.nMatch( "ZONE" ) )
00259 {
00260
00261 save.punarg[save.nsave][2] = 0;
00262 sprintf( chHeader, "#depth\torth\tpar\te=1 rel pop\te=2 rel pop\tv,J rel pops\n" );
00263 }
00264 else
00265 {
00266
00267
00268
00269 if( p.nMatch( "MATR" ) )
00270 {
00271
00272 save.punarg[save.nsave][2] = 1;
00273 sprintf( chHeader, "#vib\trot\tpops\n" );
00274 }
00275 else
00276 {
00277
00278 save.punarg[save.nsave][2] = -1;
00279 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" );
00280 }
00281 }
00282 }
00283
00284 else if( p.nMatch("RATE") )
00285 {
00286
00287 strcpy( save.chSave[save.nsave], "H2ra" );
00288 sprintf( chHeader,
00289 "#depth\tN(H2)\tN(H2)/u(H2)\tA_V(star)\tn(Eval)"
00290 "\tH2/Htot\trenorm\tfrm grn\tfrmH-\tdstTH85\tBD96\tELWERT\tBigH2\telec->H2g\telec->H2s"
00291 "\tG(TH85)\tG(DB96)\tCR\tEleclife\tShield(BD96)\tShield(H2)\tBigh2/G0(spc)\ttot dest"
00292 "\tHeatH2Dish_TH85\tHeatH2Dexc_TH85\tHeatDish_BigH2\tHeatDexc_BigH2\thtot\n" );
00293 }
00294 else if( p.nMatch("SOLO") )
00295 {
00296
00297 strcpy( save.chSave[save.nsave], "H2so" );
00298 sprintf( chHeader,
00299 "#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" );
00300 }
00301 else if( p.nMatch("SPEC") )
00302 {
00303
00304 strcpy( save.chSave[save.nsave], "H2sp" );
00305 sprintf( chHeader,
00306 "#depth\tspecial\n" );
00307 }
00308 else if( p.nMatch("TEMP") )
00309 {
00310
00311 strcpy( save.chSave[save.nsave], "H2te" );
00312 sprintf( chHeader,
00313 "#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");
00314 }
00315 else if( p.nMatch("THER") )
00316 {
00317
00318 strcpy( save.chSave[save.nsave], "H2th" );
00319 sprintf( chHeader,
00320 "#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");
00321 }
00322 else
00323 {
00324 fprintf( ioQQQ,
00325 " There must be a second key; they are RATE, LINE, COOL, COLUMN, _PDR, SOLOmon, TEMP, and POPUlations\n" );
00326 cdEXIT(EXIT_FAILURE);
00327 }
00328 return;
00329 }
00330
00331
00332
00333 void diatomics::H2_Prt_Zone(void)
00334 {
00335
00336 if( !lgEnabled || !nCall_this_zone )
00337 return;
00338
00339 DEBUG_ENTRY( "H2_Prt_Zone()" );
00340
00341 fprintf( ioQQQ, " %s density ", label.c_str() );
00342 fprintf(ioQQQ,PrintEfmt("%9.2e", (*dense_total)));
00343
00344 fprintf( ioQQQ, " orth/par");
00345 fprintf(ioQQQ,PrintEfmt("%9.2e", ortho_density / SDIV( para_density )));
00346
00347 fprintf( ioQQQ, " v0 J=0,3");
00348 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][0] ].Pop() / (*dense_total)));
00349 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][1] ].Pop() / (*dense_total)));
00350 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][2] ].Pop() / (*dense_total)));
00351 fprintf(ioQQQ,PrintEfmt("%9.2e", states[ ipEnergySort[0][0][3] ].Pop() / (*dense_total)));
00352
00353 fprintf( ioQQQ, " TOTv=0,3");
00354 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][0] / (*dense_total)));
00355 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][1] / (*dense_total)));
00356 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][2] / (*dense_total)));
00357 fprintf(ioQQQ,PrintEfmt("%9.2e", pops_per_vib[0][3] / (*dense_total)));
00358 fprintf( ioQQQ, "\n");
00359 return;
00360 }
00361
00362 void diatomics::H2_PrtDepartCoef(void)
00363 {
00364
00365 if( !lgEnabled || !nCall_this_zone )
00366 return;
00367
00368 DEBUG_ENTRY( "H2_PrtDepartCoef()" );
00369
00370
00371 fprintf( ioQQQ, " %s departure coefficients\n", label.c_str() );
00372 for( long iElec=0; iElec<n_elec_states; ++iElec )
00373 {
00374 fprintf( ioQQQ, "%li electronic\n", iElec );
00375 for( long iVib=0; iVib<=nVib_hi[iElec]; ++iVib )
00376 {
00377 for( long iRot=0; iRot<Jlowest[iElec]; ++iRot )
00378 fprintf( ioQQQ, " -----" );
00379 for( long iRot=Jlowest[iElec]; iRot<=nRot_hi[iElec][iVib]; ++iRot )
00380 {
00381 long i = ipEnergySort[iElec][iVib][iRot];
00382 fprintf( ioQQQ, " %5.3f", depart[i] );
00383 }
00384 fprintf( ioQQQ, "\n" );
00385 }
00386 fprintf( ioQQQ, "\n" );
00387 if( iElec==0 )
00388 break;
00389 }
00390
00391 return;
00392 }
00393
00394
00395 void diatomics::H2_Prt_column_density(
00396
00397
00398 FILE *ioMEAN )
00399
00400 {
00401 int iVibHi;
00402
00403
00404 if( !lgEnabled || !nCall_this_zone )
00405 return;
00406
00407 DEBUG_ENTRY( "H2_Prt_column_density()" );
00408
00409 fprintf( ioMEAN, " H2 total ");
00410 fprintf(ioMEAN,"%7.3f", log10(SDIV(ortho_colden + para_colden)));
00411
00412 fprintf( ioMEAN, " H2 ortho ");
00413 fprintf(ioMEAN,"%7.3f", log10(SDIV(ortho_colden)));
00414
00415 fprintf( ioMEAN, " para");
00416 fprintf(ioMEAN,"%7.3f", log10(SDIV(para_colden)));
00417
00418 iVibHi = 0;
00419 fprintf( ioMEAN, " v0 J=0,3");
00420 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][0])));
00421 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][1])));
00422 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][2])));
00423 fprintf(ioMEAN,"%7.3f", log10(SDIV(H2_X_colden[iVibHi][3])));
00424
00425 return;
00426 }
00427
00428
00429
00430 void diatomics::H2_ReadTransprob( long int nelec , TransitionList &trns)
00431 {
00432 const char* cdDATAFILE[N_ELEC] =
00433 {
00434 "transprob_X.dat",
00435 "transprob_B.dat",
00436 "transprob_C_plus.dat",
00437 "transprob_C_minus.dat",
00438 "transprob_B_primed.dat",
00439 "transprob_D_plus.dat",
00440 "transprob_D_minus.dat"
00441 };
00442 FILE *ioDATA;
00443 char chLine[FILENAME_PATH_LENGTH_2];
00444 long int i, n1, n2, n3;
00445 long int iVibHi , iVibLo , iRotHi , iRotLo , iElecHi , iElecLo;
00446 bool lgEOL;
00447
00448 DEBUG_ENTRY( "H2_ReadTransprob()" );
00449
00450
00451 char chPath[FILENAME_PATH_LENGTH_2];
00452 strcpy( chPath, path.c_str() );
00453 strcat( chPath, input.chDelimiter );
00454 strcat( chPath, cdDATAFILE[nelec] );
00455 ioDATA = open_data( chPath , "r" );
00456
00457
00458 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00459 {
00460 fprintf( ioQQQ, " H2_ReadTransprob could not read first line of %s\n", cdDATAFILE[nelec]);
00461 cdEXIT(EXIT_FAILURE);
00462 }
00463 i = 1;
00464
00465 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00466 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00467 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00468
00469
00470
00471 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
00472 {
00473 fprintf( ioQQQ,
00474 " H2_ReadTransprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
00475 fprintf( ioQQQ,
00476 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
00477 n1 , n2 , n3 );
00478 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00479 cdEXIT(EXIT_FAILURE);
00480 }
00481
00482 long nlines = 0;
00483 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00484 {
00485
00486 if( chLine[0]=='#' )
00487 continue;
00488 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
00489 break;
00490
00491 double Aul;
00492 int n = sscanf(chLine,"%li\t%li\t%li\t%li\t%li\t%li\t%le",
00493 &iElecHi , &iVibHi ,&iRotHi , &iElecLo , &iVibLo , &iRotLo , &Aul );
00494 ASSERT( n == 7 );
00495 ASSERT( iElecHi == nelec );
00496 ASSERT( iElecHi < N_ELEC );
00497 ASSERT( iElecLo < N_ELEC );
00498
00499
00500 if( iVibHi <= nVib_hi[iElecHi] &&
00501 iVibLo <= nVib_hi[iElecLo] &&
00502 iRotHi <= nRot_hi[iElecHi][iVibHi] &&
00503 iRotLo <= nRot_hi[iElecLo][iVibLo])
00504 {
00505 long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
00506 long ipLo = ipEnergySort[iElecLo][iVibLo][iRotLo];
00507 double ener = states[ipHi].energy().WN() - states[ipLo].energy().WN();
00508 long lineIndex = ipTransitionSort[ipHi][ipLo];
00509
00510
00511 trns[lineIndex].AddLine2Stack();
00512 trns[lineIndex].Emis().Aul() = (realnum)Aul;
00513
00514
00515 lgH2_radiative[ipHi][ipLo] = true;
00516 ++nlines;
00517
00518
00519 if( ener <= 0. )
00520 {
00521 fprintf(ioQQQ,"negative energy H2 transition\t%li\t%li\t%li\t%li\t%.2e\t%.2e\n",
00522 iVibHi,iVibLo,iRotHi,iRotLo,Aul,ener);
00523 ShowMe();
00524 cdEXIT(EXIT_FAILURE);
00525 }
00526 }
00527 }
00528 if( nTRACE >= n_trace_full )
00529 fprintf(ioQQQ," There are a total of %li lines in the entire H2 molecule.\n", nlines );
00530
00531 fclose( ioDATA );
00532 return;
00533 }
00534
00535 #if 0
00536
00537 void H2_Read_Cosmicray_distribution(void)
00538 {
00539
00540
00541
00542 FILE *ioDATA;
00543 char chLine[FILENAME_PATH_LENGTH_2];
00544 long int i, n1, n2, n3, iVib , iRot;
00545 long neut_frac;
00546 bool lgEOL;
00547
00548 DEBUG_ENTRY( "H2_Read_Cosmicray_distribution()" );
00549
00550
00551 char chPath[FILENAME_PATH_LENGTH_2];
00552 strcpy( chPath, path.c_str() );
00553 strcat( chPath, input.chDelimiter );
00554 strcat( chPath, "H2_CosmicRay_collision.dat" );
00555 ioDATA = open_data( chPath, "r" );
00556
00557
00558 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00559 {
00560 fprintf( ioQQQ, " H2_Read_Cosmicray_distribution could not read first line of %s\n", "H2_Cosmic_collision.dat");
00561 cdEXIT(EXIT_FAILURE);
00562 }
00563
00564 i = 1;
00565
00566 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00567 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00568 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00569
00570
00571
00572 if( ( n1 != 1 ) || ( n2 != 21 ) || ( n3 != 3 ) )
00573 {
00574 fprintf( ioQQQ,
00575 " H2_Read_Cosmicray_distribution: the version of %s is not the current version.\n", "H2_Cosmic_collision.dat" );
00576 fprintf( ioQQQ,
00577 " I expected to find the number 1 21 3 and got %li %li %li instead.\n" ,
00578 n1 , n2 , n3 );
00579 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00580 cdEXIT(EXIT_FAILURE);
00581 }
00582
00583
00584 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00585 BadRead();
00586
00587 while( chLine[0]=='#' )
00588 {
00589 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00590 BadRead();
00591 }
00592
00593 iRot = 1;
00594 iVib = 1;
00595 neut_frac = 0;
00596 while( iVib >= 0 )
00597 {
00598 long int j_minus_ji;
00599 double a[10];
00600
00601 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",
00602 &iVib ,&j_minus_ji , &a[0],&a[1],&a[2],&a[3],&a[4],&a[5],&a[6],&a[7],&a[8],&a[9]
00603 );
00604
00605 if( iVib < 0 )
00606 continue;
00607
00608
00609
00610 ASSERT( iVib < CR_VIB );
00611 ASSERT( j_minus_ji == -2 || j_minus_ji == +2 || j_minus_ji == 0 );
00612 ASSERT( neut_frac < CR_X );
00613
00614
00615 j_minus_ji = 1 + j_minus_ji/2;
00616 ASSERT( j_minus_ji>=0 && j_minus_ji<=2 );
00617
00618
00619 for( iRot=0; iRot<CR_J; ++iRot )
00620 {
00621 cr_rate[neut_frac][iVib][iRot][j_minus_ji] = (realnum)a[iRot];
00622 }
00623 if( lgH2_NOISECOSMIC )
00624 {
00625 realnum r;
00626 r = (realnum)RandGauss( xMeanNoise , xSTDNoise );
00627
00628 for( iRot=0; iRot<CR_J; ++iRot )
00629 {
00630 cr_rate[neut_frac][iVib][iRot][j_minus_ji] *= (realnum)pow(10.,(double)r);
00631 }
00632 }
00633
00634 if( CR_PRINT )
00635 {
00636 fprintf(ioQQQ,"cr rate\t%li\t%li", iVib , j_minus_ji );
00637 for( iRot=0; iRot<CR_J; ++iRot )
00638 {
00639 fprintf(ioQQQ,"\t%.3e", cr_rate[neut_frac][iVib][iRot][j_minus_ji] );
00640 }
00641 fprintf(ioQQQ,"\n" );
00642 }
00643
00644
00645 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00646 BadRead();
00647 while( chLine[0]=='#' )
00648 {
00649 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00650 BadRead();
00651 }
00652 }
00653 fclose( ioDATA );
00654
00655 return;
00656 }
00657 #endif
00658
00659 class level_tmp
00660 {
00661 public:
00662 bool operator<( const level_tmp& second ) const
00663 {
00664 if( eWN < second.eWN )
00665 return true;
00666 else
00667 return false;
00668 }
00669 long n, v, J;
00670 double eWN;
00671 };
00672
00673
00674 void diatomics::H2_ReadEnergies( )
00675 {
00676 DEBUG_ENTRY( "H2_ReadEnergies()" );
00677
00678 vector<int> n, v, J;
00679 vector<double>eWN;
00680
00681 for( long nelec=0; nelec<n_elec_states; ++nelec )
00682 {
00683
00684 H2_ReadEnergies(nelec,n,v,J,eWN);
00685 }
00686
00687 vector<level_tmp> levels;
00688 levels.resize( n.size() );
00689 ASSERT( levels.size() > 0 );
00690 for( unsigned i = 0; i < n.size(); ++i )
00691 {
00692 levels[i].n = n[i];
00693 levels[i].v = v[i];
00694 levels[i].J = J[i];
00695 levels[i].eWN = eWN[i];
00696 }
00697
00698
00699 sort( levels.begin(), levels.end() );
00700
00701
00702 for( vector<level_tmp>::iterator lev = levels.begin(); lev != levels.end(); ++lev )
00703 {
00704 states.resize( states.size() + 1 );
00705 long i = states.size() - 1;
00706 states[i].n() = lev->n;
00707 states[i].v() = lev->v;
00708 states[i].J() = lev->J;
00709 states[i].energy().set( lev->eWN, "cm^-1" );
00710
00711
00712 states[i].nelem() = -1;
00713
00714 states[i].IonStg() = -1;
00715 strcpy( states[i].chLabel(), label.c_str() );
00716 }
00717
00718 ASSERT( states.size() > 0 );
00719 ASSERT( states.size() == levels.size() );
00720
00721 for( long nelec=0; nelec<n_elec_states; ++nelec )
00722 {
00723 ASSERT( nLevels_per_elec[nelec] > 0 );
00724 ASSERT( nVib_hi[nelec] > 0 );
00725 ASSERT( nVib_hi[nelec] > Jlowest[nelec] );
00726
00727 nRot_hi[nelec].resize( nVib_hi[nelec]+1 );
00728 nRot_hi[nelec] = 0;
00729 }
00730
00731 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
00732 {
00733 nRot_hi[ (*st).n() ][ (*st).v() ] = MAX2( nRot_hi[ (*st).n() ][ (*st).v() ], (*st).J() );
00734 }
00735
00736 return;
00737 }
00738
00739 void diatomics::H2_ReadEnergies( long int nelec, vector<int>& n, vector<int>& v, vector<int>&J, vector<double>& eWN )
00740 {
00741 DEBUG_ENTRY("diatomics::H2_ReadEnergies()");
00742 const char* cdDATAFILE[N_ELEC] =
00743 {
00744 "energy_X.dat",
00745 "energy_B.dat",
00746 "energy_C_plus.dat",
00747 "energy_C_minus.dat",
00748 "energy_B_primed.dat",
00749 "energy_D_plus.dat",
00750 "energy_D_minus.dat"
00751 };
00752
00753 char chPath[FILENAME_PATH_LENGTH_2];
00754 strcpy( chPath, path.c_str() );
00755 strcat( chPath, input.chDelimiter );
00756 strcat( chPath, cdDATAFILE[nelec] );
00757 FILE *ioDATA = open_data( chPath, "r" );
00758
00759 char chLine[FILENAME_PATH_LENGTH_2];
00760 long int i, n1, n2, n3;
00761 bool lgEOL;
00762
00763
00764 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00765 {
00766 fprintf( ioQQQ, " H2_ReadEnergies could not read first line of %s\n", cdDATAFILE[nelec]);
00767 cdEXIT(EXIT_FAILURE);
00768 }
00769 i = 1;
00770
00771 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00772 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00773 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00774
00775
00776
00777 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
00778 {
00779 fprintf( ioQQQ,
00780 " H2_ReadEnergies: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
00781 fprintf( ioQQQ,
00782 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
00783 n1 , n2 , n3 );
00784 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00785 cdEXIT(EXIT_FAILURE);
00786 }
00787
00788
00789 nLevels_per_elec[nelec] = 0;
00790 nVib_hi[nelec] = 0;
00791 Jlowest[nelec] = LONG_MAX;
00792
00793 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00794 {
00795
00796 if( chLine[0]=='#' )
00797 continue;
00798 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
00799 break;
00800 long iVib, iRot;
00801 double energyWN;
00802 int nReads = sscanf(chLine,"%li\t%li\t%le", &iVib, &iRot, &energyWN );
00803 ASSERT( nReads == 3 );
00804 ASSERT( iVib >= 0 );
00805 ASSERT( iRot >= 0 );
00806 ASSERT( energyWN > 0. || (nelec==0 && iVib==0 && iRot==0 ) );
00807
00808 n.push_back( nelec );
00809 v.push_back( iVib );
00810 J.push_back( iRot );
00811 eWN.push_back( energyWN );
00812
00813
00814 nVib_hi[nelec] = MAX2( nVib_hi[nelec], iVib );
00815 Jlowest[nelec] = MIN2( Jlowest[nelec], iRot );
00816
00817 ++nLevels_per_elec[nelec];
00818 }
00819
00820 ASSERT( n.size() > 0 );
00821 ASSERT( nLevels_per_elec[nelec] > 0 );
00822 ASSERT( nVib_hi[nelec] > 0 );
00823 ASSERT( nVib_hi[nelec] > Jlowest[nelec] );
00824
00825 fclose( ioDATA );
00826
00827 return;
00828 }
00829
00830
00831 void diatomics::H2_ReadDissocEnergies( void )
00832 {
00833 const char* cdDATAFILE = "energy_dissoc.dat";
00834 FILE *ioDATA;
00835 char chLine[FILENAME_PATH_LENGTH_2];
00836 long int i, n1, n2, n3;
00837 bool lgEOL;
00838
00839 DEBUG_ENTRY( "H2_ReadDissocEnergies()" );
00840
00841
00842 char chPath[FILENAME_PATH_LENGTH_2];
00843 strcpy( chPath, path.c_str() );
00844 strcat( chPath, input.chDelimiter );
00845 strcat( chPath, cdDATAFILE );
00846 ioDATA = open_data( chPath, "r" );
00847
00848
00849 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00850 {
00851 fprintf( ioQQQ, " H2_ReadDissocEnergies could not read first line of %s\n", cdDATAFILE );
00852 cdEXIT(EXIT_FAILURE);
00853 }
00854 i = 1;
00855
00856 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00857 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00858 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00859
00860
00861
00862 if( ( n1 != 2 ) || ( n2 != 4 ) || ( n3 != 29 ) )
00863 {
00864 fprintf( ioQQQ,
00865 " H2_ReadDissocEnergies: the version of %s is not the current version.\n", cdDATAFILE );
00866 fprintf( ioQQQ,
00867 " I expected to find the number 2 4 29 and got %li %li %li instead.\n" ,
00868 n1 , n2 , n3 );
00869 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00870 cdEXIT(EXIT_FAILURE);
00871 }
00872
00873 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00874 {
00875
00876 if( chLine[0]=='#' )
00877 continue;
00878 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
00879 break;
00880 long iElec;
00881 double energyWN;
00882 int n = sscanf(chLine,"%li\t%le", &iElec, &energyWN );
00883 ASSERT( n == 2 );
00884 ASSERT( iElec >= 0 );
00885 ASSERT( iElec < N_ELEC );
00886 ASSERT( energyWN > 0. );
00887 H2_DissocEnergies[iElec] = energyWN;
00888 }
00889 fclose( ioDATA );
00890
00891 return;
00892 }
00893
00894
00895 void diatomics::H2_ReadDissprob( long int nelec )
00896 {
00897 const char* cdDATAFILE[N_ELEC] =
00898 {
00899 "dissprob_X.dat",
00900 "dissprob_B.dat",
00901 "dissprob_C_plus.dat",
00902 "dissprob_C_minus.dat",
00903 "dissprob_B_primed.dat",
00904 "dissprob_D_plus.dat",
00905 "dissprob_D_minus.dat"
00906 };
00907 char chLine[FILENAME_PATH_LENGTH_2];
00908 bool lgEOL;
00909
00910 DEBUG_ENTRY( "H2_ReadDissprob()" );
00911
00912 ASSERT( nelec > 0 );
00913
00914
00915 char chPath[FILENAME_PATH_LENGTH_2];
00916 strcpy( chPath, path.c_str() );
00917 strcat( chPath, input.chDelimiter );
00918 strcat( chPath, cdDATAFILE[nelec] );
00919 FILE *ioDATA = open_data( chPath, "r" );
00920
00921
00922 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00923 {
00924 fprintf( ioQQQ, " H2_ReadDissprob could not read first line of %s\n", cdDATAFILE[nelec]);
00925 cdEXIT(EXIT_FAILURE);
00926 }
00927 long i = 1;
00928
00929 long n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00930 long n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00931 long n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00932
00933
00934
00935 if( ( n1 != 3 ) || ( n2 != 2 ) || ( n3 != 11 ) )
00936 {
00937 fprintf( ioQQQ,
00938 " H2_ReadDissprob: the version of %s is not the current version.\n", cdDATAFILE[nelec] );
00939 fprintf( ioQQQ,
00940 " I expected to find the number 3 2 11 and got %li %li %li instead.\n" ,
00941 n1 , n2 , n3 );
00942 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00943 cdEXIT(EXIT_FAILURE);
00944 }
00945
00946 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00947 {
00948
00949 if( chLine[0]=='#' )
00950 continue;
00951 if( chLine[0]=='\n' || chLine[0]=='\0' || chLine[0]==' ' )
00952 break;
00953
00954 long iVib, iRot;
00955 double a, b;
00956 i = 1;
00957 sscanf(chLine,"%li\t%li\t%le\t%le",
00958 &iVib, &iRot,
00959
00960 &a ,
00961
00962 &b);
00963
00964
00965
00966
00967
00968
00969 if( ( iVib < 0 ) ||
00970 ( iVib > nVib_hi[nelec] ) ||
00971 ( iRot < Jlowest[nelec] ) ||
00972 ( iRot > nRot_hi[nelec][iVib] ) )
00973 continue;
00974
00975
00976 H2_dissprob[nelec][iVib][iRot] = (realnum)a;
00977
00978 H2_disske[nelec][iVib][iRot] = (realnum)b;
00979 }
00980 fclose( ioDATA );
00981 return;
00982 }
00983
00984
00985
00986 void diatomics::H2_Read_hminus_distribution(void)
00987 {
00988 FILE *ioDATA;
00989 char chLine[FILENAME_PATH_LENGTH_2];
00990 long int i, n1, n2, n3, iVib , iRot;
00991 bool lgEOL;
00992 double sumrate[nTE_HMINUS] = {0};
00993
00994 const bool lgH2HMINUS_PRT = false;
00995
00996 DEBUG_ENTRY( "H2_Read_hminus_distribution()" );
00997
00998
00999 char chPath[FILENAME_PATH_LENGTH_2];
01000 strcpy( chPath, path.c_str() );
01001 strcat( chPath, input.chDelimiter );
01002 strcat( chPath, "hminus_deposit.dat" );
01003 ioDATA = open_data( chPath, "r" );
01004
01005
01006 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01007 {
01008 fprintf( ioQQQ, " H2_Read_hminus_distribution could not read first line of %s\n", chPath );
01009 cdEXIT(EXIT_FAILURE);
01010 }
01011
01012 i = 1;
01013
01014 n1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01015 n2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01016 n3 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
01017
01018
01019
01020 if( ( n1 != 2 ) || ( n2 != 10 ) || ( n3 != 17 ) )
01021 {
01022 fprintf( ioQQQ,
01023 " H2_Read_hminus_distribution: the version of %s is not the current version.\n", chPath );
01024 fprintf( ioQQQ,
01025 " I expected to find the number 2 10 17 and got %li %li %li instead.\n" ,
01026 n1 , n2 , n3 );
01027 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
01028 cdEXIT(EXIT_FAILURE);
01029 }
01030
01031
01032 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01033 BadRead();
01034
01035 while( chLine[0]=='#' )
01036 {
01037 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01038 BadRead();
01039 }
01040
01041 iRot = 1;
01042 iVib = 1;
01043 while( iVib >= 0 )
01044 {
01045
01046
01047 double a[nTE_HMINUS] , ener;
01048 sscanf(chLine,"%li\t%li\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf",
01049 &iVib ,&iRot , &ener, &a[0],&a[1],&a[2] , &a[3],&a[4],&a[5] ,&a[6]
01050 );
01051
01052 if( iVib < 0 )
01053 continue;
01054
01055
01056 ASSERT( iVib <= nVib_hi[0] &&
01057 iRot <= nRot_hi[0][iVib] );
01058
01059 if( lgH2HMINUS_PRT )
01060 fprintf(ioQQQ,"hminusss\t%li\t%li", iVib , iRot );
01061 for( i=0; i<nTE_HMINUS; ++i )
01062 {
01063 H2_X_hminus_formation_distribution[i][iVib][iRot] = (realnum)pow(10.,-a[i]);
01064 sumrate[i] += H2_X_hminus_formation_distribution[i][iVib][iRot];
01065 if( lgH2HMINUS_PRT )
01066 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
01067 }
01068 if( lgH2HMINUS_PRT )
01069 fprintf(ioQQQ,"\n" );
01070
01071 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01072 BadRead();
01073 while( chLine[0]=='#' )
01074 {
01075 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01076 BadRead();
01077 }
01078 }
01079 fclose( ioDATA );
01080
01081 if( lgH2HMINUS_PRT )
01082 {
01083
01084 fprintf(ioQQQ," total H- formation rate ");
01085
01086 for(i=0; i<nTE_HMINUS; ++i )
01087 {
01088 fprintf(ioQQQ,"\t%.3e" , sumrate[i]);
01089 }
01090 fprintf(ioQQQ,"\n" );
01091 }
01092
01093
01094 for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
01095 {
01096 for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
01097 {
01098 for(i=0; i<nTE_HMINUS; ++i )
01099 {
01100 H2_X_hminus_formation_distribution[i][iVib][iRot] /= (realnum)sumrate[i];
01101 }
01102 }
01103 }
01104
01105 if( lgH2HMINUS_PRT )
01106 {
01107
01108 fprintf(ioQQQ," H- distribution function ");
01109 for( iVib=0; iVib<=nVib_hi[0]; ++iVib )
01110 {
01111 for( iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
01112 {
01113 fprintf(ioQQQ,"%li\t%li", iVib , iRot );
01114 for(i=0; i<nTE_HMINUS; ++i )
01115 {
01116 fprintf(ioQQQ,"\t%.3e", H2_X_hminus_formation_distribution[i][iVib][iRot] );
01117 }
01118 fprintf(ioQQQ,"\n" );
01119 }
01120 }
01121 }
01122 return;
01123 }
01124
01125
01126
01127 void diatomics::H2_Punch_line_data(
01128
01129 FILE* ioPUN ,
01130
01131 bool lgDoAll )
01132 {
01133 if( !lgEnabled )
01134 return;
01135
01136 DEBUG_ENTRY( "H2_Punch_line_data()" );
01137
01138 if( lgDoAll )
01139 {
01140 fprintf( ioQQQ,
01141 " H2_Punch_line_data ALL option not implemented in H2_Punch_line_data yet 1\n" );
01142 cdEXIT(EXIT_FAILURE);
01143 }
01144 else
01145 {
01146 bool lgPrint = false;
01147 fprintf( ioPUN, "#Eu\tVu\tJu\tEl\tVl\tJl\tWL\tgl\tgu\tgf\tA\tCS\tn(crt)\n" );
01148
01149 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
01150 {
01151 if( (*tr).ipCont() <= 0 )
01152 continue;
01153 (*tr).Coll().col_str() = 0.;
01154 qList::iterator Hi = ( (*tr).Hi() );
01155 qList::iterator Lo = ( (*tr).Lo() );
01156
01157 fprintf(ioPUN,"%2li\t%2li\t%2li\t%2li\t%2li\t%2li\t",
01158 (*Hi).n(), (*Hi).v(), (*Hi).J(),
01159 (*Lo).n(), (*Lo).v(), (*Lo).J() );
01160 Save1LineData( *tr, ioPUN, false , lgPrint);
01161 }
01162
01163 fprintf( ioPUN , "\n");
01164 }
01165 return;
01166 }
01167
01168
01169 void diatomics::H2_PunchLineStuff( FILE * io , realnum xLimit , long index)
01170 {
01171 if( !lgEnabled )
01172 return;
01173
01174 DEBUG_ENTRY( "H2_PunchLineStuff()" );
01175
01176
01177 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
01178 {
01179 if( (*tr).ipCont() <= 0 )
01180 continue;
01181 Save1Line( *tr, io, xLimit, index, GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN]));
01182 }
01183
01184 return;
01185 }
01186
01187
01188
01189 void diatomics::H2_Prt_line_tau(void)
01190 {
01191 if( !lgEnabled )
01192 return;
01193
01194 DEBUG_ENTRY( "H2_Prt_line_tau()" );
01195
01196
01197 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
01198 {
01199 if( (*tr).ipCont() <= 0 )
01200 continue;
01201 prme( false, *tr );
01202 }
01203
01204 return;
01205 }
01206
01207
01208
01209 STATIC char chMolBranch( long iRotHi , long int iRotLo )
01210 {
01211
01212 char chBranch[5] = {'O','P','Q','R','S'};
01213
01214 int ip = 2 + (iRotHi - iRotLo);
01215 if( ip<0 || ip>=5 )
01216 {
01217 fprintf(ioQQQ," chMolBranch called with insane iRotHi=%li iRotLo=%li ip=%i\n",
01218 iRotHi , iRotLo , ip );
01219 ip = 0;
01220 }
01221
01222 return( chBranch[ip] );
01223 }
01224
01225
01226 void diatomics::H2_PunchDo( FILE* io , char chJOB[] , const char chTime[] , long int ipPun )
01227 {
01228 DEBUG_ENTRY( "H2_PunchDo()" );
01229
01230
01231
01232
01233
01234
01235 if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") == 0) &&
01236 (save.punarg[ipPun][2] != 0) )
01237 {
01238
01239 if( lgEnabled && lgEvaluated )
01240 {
01241 long iVibHi= 0;
01242 long iRotHi = 0;
01243 long iElecHi=0;
01244 long LimVib, LimRot;
01245
01246
01247
01248
01249 if( save.punarg[ipPun][0] > 0 )
01250 {
01251 LimVib = (long)save.punarg[ipPun][0];
01252 }
01253 else
01254 {
01255 LimVib = nVib_hi[iElecHi];
01256 }
01257
01258
01259 fprintf(io,"%i\t%i\t%.3e\tortho\n",
01260 103 ,
01261 103 ,
01262 ortho_density );
01263 fprintf(io,"%i\t%i\t%.3e\tpara\n",
01264 101 ,
01265 101 ,
01266 para_density );
01267 fprintf(io,"%i\t%i\t%.3e\ttotal\n",
01268 0 ,
01269 0 ,
01270 (*dense_total) );
01271
01272
01273 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
01274 {
01275
01276 if( save.punarg[ipPun][1] > 0 )
01277 {
01278 LimRot = (long)MIN2(
01279 save.punarg[ipPun][1] , (realnum)nRot_hi[iElecHi][iVibHi]);
01280 }
01281 else
01282 {
01283 LimRot = nRot_hi[iElecHi][iVibHi];
01284 }
01285 if( save.punarg[ipPun][2] > 0 )
01286 {
01287 long int i;
01288
01289 if( iVibHi == 0 )
01290 {
01291 fprintf(io,"vib\\rot");
01292
01293 for( i=0; i<=LimRot; ++i )
01294 {
01295 fprintf(io,"\t%li",i);
01296 }
01297 fprintf(io,"\n");
01298 }
01299 fprintf(io,"%li",iVibHi );
01300 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01301 {
01302 fprintf(io,"\t%.3e",
01303 states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].Pop()/(*dense_total) );
01304 }
01305 fprintf(io,"\n" );
01306 }
01307 else if( save.punarg[ipPun][2] < 0 )
01308 {
01309
01310 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01311 {
01312
01313
01314
01315 const char chlgPara[2]={'P','O'};
01316 const long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
01317
01318
01319 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",
01320
01321 iVibHi , iRotHi ,
01322
01323 chlgPara[H2_lgOrtho[iElecHi][iVibHi][iRotHi]],
01324
01325 states[ipHi].energy().WN(),
01326
01327 states[ipHi].Pop()/(*dense_total) ,
01328
01329 H2_old_populations[iElecHi][iVibHi][iRotHi]/(*dense_total) ,
01330
01331 states[ipHi].Pop()/(*dense_total)/H2_stat[iElecHi][iVibHi][iRotHi] ,
01332
01333
01334 states[ipHi].Pop()/SDIV(H2_populations_LTE[iElecHi][iVibHi][iRotHi]*(*dense_total) ) ,
01335
01336 H2_col_rate_out[iVibHi][iRotHi]/SDIV(H2_col_rate_out[iVibHi][iRotHi]+H2_rad_rate_out[0][iVibHi][iRotHi]) ,
01337
01338 H2_col_rate_in[iVibHi][iRotHi]/SDIV(H2_col_rate_in[iVibHi][iRotHi]+H2_rad_rate_in[iVibHi][iRotHi]),
01339
01340 H2_col_rate_out[iVibHi][iRotHi],
01341
01342 H2_rad_rate_out[0][iVibHi][iRotHi] ,
01343
01344 H2_col_rate_in[iVibHi][iRotHi],
01345
01346 H2_rad_rate_in[iVibHi][iRotHi]
01347 );
01348 }
01349 }
01350 }
01351 }
01352 }
01353
01354
01355 else if( (strcmp( chJOB , "H2po" ) == 0) && (strcmp(chTime,"LAST") != 0) &&
01356 (save.punarg[ipPun][2] == 0) )
01357 {
01358
01359 if( lgEnabled && lgEvaluated )
01360 {
01361 fprintf(io,"%.5e\t%.3e\t%.3e", radius.depth_mid_zone ,
01362 ortho_density , para_density);
01363
01364 fprintf(io,"\t%.3e\t%.3e",
01365 pops_per_elec[1] , pops_per_elec[2]);
01366 long iElecHi = 0;
01367 long iVibHi = 0;
01368 long LimVib, LimRot;
01369
01370 if( save.punarg[ipPun][0] > 0 )
01371 {
01372 LimVib = (long)save.punarg[ipPun][1];
01373 }
01374 else
01375 {
01376 LimVib = nRot_hi[iElecHi][iVibHi];
01377 }
01378 LimVib = MIN2( LimVib , nVib_hi[iElecHi] );
01379
01380 if( save.punarg[ipPun][1] > 0 )
01381 {
01382 LimRot = (long)save.punarg[ipPun][1];
01383 }
01384 else
01385 {
01386 LimRot = nRot_hi[iElecHi][iVibHi];
01387 }
01388 for( iVibHi = 0; iVibHi<=LimVib; ++iVibHi )
01389 {
01390 fprintf(io,"\tv=%li",iVibHi);
01391 long int LimRotVib = MIN2( LimRot , nRot_hi[iElecHi][iVibHi] );
01392 for( long iRotHi=Jlowest[iElecHi]; iRotHi<=LimRotVib; ++iRotHi )
01393 {
01394 fprintf(io,"\t%.3e",
01395 states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].Pop()/(*dense_total) );
01396 }
01397 }
01398 fprintf(io,"\n");
01399 }
01400 }
01401
01402
01403 else if( (strcmp( chJOB , "H2cl" ) == 0) && (strcmp(chTime,"LAST") == 0) )
01404 {
01405 long iVibHi= 0;
01406 long iRotHi = 0;
01407 long iElecHi=0;
01408 long LimVib, LimRot;
01409
01410
01411
01412
01413 if( save.punarg[ipPun][0] > 0 )
01414 {
01415 LimVib = (long)save.punarg[ipPun][0];
01416 }
01417 else
01418 {
01419 LimVib = nVib_hi[iElecHi];
01420 }
01421
01422
01423 fprintf(io,"%i\t%i\t%.3e\tortho\n",
01424 103 ,
01425 103 ,
01426 ortho_colden );
01427 fprintf(io,"%i\t%i\t%.3e\tpara\n",
01428 101 ,
01429 101 ,
01430 para_colden );
01431
01432 fprintf(io,"%i\t%i\t%.3e\ttotal\n",
01433 0 ,
01434 0 ,
01435 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s]);
01436
01437
01438 for( iVibHi=0; iVibHi<=LimVib; ++iVibHi )
01439 {
01440 if( lgEnabled )
01441 {
01442
01443 if( save.punarg[ipPun][1] > 0 )
01444 {
01445 LimRot = (long)save.punarg[ipPun][1];
01446 }
01447 else
01448 {
01449 LimRot = nRot_hi[iElecHi][iVibHi];
01450 }
01451 if( save.punarg[ipPun][2] > 0 )
01452 {
01453 long int i;
01454
01455 if( iVibHi == 0 )
01456 {
01457 fprintf(io,"vib\\rot");
01458
01459 for( i=0; i<=LimRot; ++i )
01460 {
01461 fprintf(io,"\t%li",i);
01462 }
01463 fprintf(io,"\n");
01464 }
01465 fprintf(io,"%li",iVibHi );
01466 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01467 {
01468 fprintf(io,"\t%.3e",
01469 H2_X_colden[iVibHi][iRotHi]/(*dense_total) );
01470 }
01471 fprintf(io,"\n" );
01472 }
01473 else
01474 {
01475
01476 for( iRotHi=Jlowest[iElecHi]; iRotHi<=LimRot; ++iRotHi )
01477 {
01478 fprintf(io,"%li\t%li\t%.1f\t%.3e\t%.3e\t%.3e\t%.3e\n",
01479 iVibHi ,
01480 iRotHi ,
01481
01482 states[ ipEnergySort[iElecHi][iVibHi][iRotHi] ].energy().K(),
01483
01484 H2_X_colden[iVibHi][iRotHi] ,
01485 H2_X_colden[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi] ,
01486
01487 H2_X_colden_LTE[iVibHi][iRotHi] ,
01488 H2_X_colden_LTE[iVibHi][iRotHi]/H2_stat[iElecHi][iVibHi][iRotHi]);
01489 }
01490 }
01491 }
01492 }
01493 }
01494 else if( (strcmp(chJOB , "H2pd" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01495 {
01496
01497
01498 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01499
01500 radius.depth_mid_zone ,
01501
01502 ortho_density ,
01503 para_density ,
01504
01505 hmi.H2_Solomon_dissoc_rate_TH85_H2g ,
01506
01507 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
01508
01509 Solomon_dissoc_rate_g);
01510 }
01511 else if( (strcmp(chJOB , "H2co" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01512 {
01513
01514 fprintf(io,"%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01515
01516 radius.depth_mid_zone ,
01517
01518 thermal.ctot ,
01519
01520 hmi.H2_Solomon_dissoc_rate_TH85_H2g,
01521
01522 Solomon_dissoc_rate_g +
01523 Solomon_dissoc_rate_s,
01524
01525 hmi.HeatH2Dish_TH85,
01526
01527 HeatDiss ,
01528
01529 hmi.HeatH2Dexc_TH85,
01530 HeatDexc
01531 );
01532
01533 }
01534
01535 else if( (strcmp(chJOB , "H2le" ) == 0) && (strcmp(chTime,"LAST") == 0) )
01536 {
01537
01538 for( long int ipHi=0; ipHi < nLevels_per_elec[0]; ipHi++ )
01539 {
01540 long iRotHi = ipRot_H2_energy_sort[ipHi];
01541 long iVibHi = ipVib_H2_energy_sort[ipHi];
01542 long int nColl;
01543 double Asum , Csum[N_X_COLLIDER];
01544 Asum = 0;
01545 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
01546 Csum[nColl] = 0.;
01547 for( long int ipLo=0; ipLo<ipHi; ++ipLo )
01548 {
01549
01550 long iRotLo = ipRot_H2_energy_sort[ipLo];
01551 long iVibLo = ipVib_H2_energy_sort[ipLo];
01552 EmissionProxy em = trans[ ipTransitionSort[ipHi][ipLo] ].Emis();
01553
01554
01555 if( ( abs(iRotHi-iRotLo) == 2 || (iRotHi-iRotLo) == 0 ) && iVibLo <= iVibHi &&
01556 lgH2_radiative[ipHi][ipLo] )
01557 {
01558 Asum += em.Aul() * ( em.Pesc() + em.Pdest() + em.Pelec_esc() );
01559 }
01560
01561 mr3ci H2cr = CollRateCoeff.begin(ipHi,ipLo);
01562 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
01563 Csum[nColl] += H2cr[nColl];
01564 }
01565
01566
01567 fprintf(io,"%li\t%li\t%.2f\t%li\t%.3e",
01568 iVibHi , iRotHi,
01569 states[ipHi].energy().WN(),
01570 (long)states[ipHi].g(),
01571 Asum );
01572 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
01573
01574 fprintf(io,"\t%.3e",Csum[nColl]);
01575 fprintf(io,"\n");
01576 }
01577 }
01578
01579 else if( (strcmp(chJOB , "H2ra" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01580 {
01581
01582 double sumpop = 0. , sumlife = 0.;
01583
01584
01585 if( lgEnabled && lgEvaluated )
01586 {
01587
01588 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
01589 {
01590 qList::iterator Lo = ( (*tr).Lo() );
01591 if( (*Lo).n() > 0 || (*Lo).v() > 0 )
01592 continue;
01593 sumlife += (*tr).Emis().pump() * (*(*tr).Lo()).Pop();
01594 sumpop += (*(*tr).Lo()).Pop();
01595 }
01596 }
01597
01598
01599
01600
01601 fprintf(io,
01602 "%.5e\t%.3e\t%.3e\t%.3e\t%li",
01603
01604 radius.depth_mid_zone ,
01605
01606 colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
01607
01608
01609 colden.coldenH2_ov_vel ,
01610
01611 rfield.extin_mag_V_point,
01612
01613 nCall_this_zone );
01614 fprintf(io,
01615 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
01616
01617 (*dense_total)/dense.gas_phase[ipHYDROGEN] ,
01618
01619 H2_renorm_chemistry,
01620
01621 gv.rate_h2_form_grains_used_total ,
01622
01623 findspecieslocal("H-")->den*1.35e-9,
01624
01625 hmi.H2_Solomon_dissoc_rate_TH85_H2g + hmi.H2_Solomon_dissoc_rate_TH85_H2s,
01626
01627 hmi.H2_Solomon_dissoc_rate_BD96_H2g + hmi.H2_Solomon_dissoc_rate_BD96_H2s,
01628
01629 hmi.H2_Solomon_dissoc_rate_ELWERT_H2g + hmi.H2_Solomon_dissoc_rate_ELWERT_H2g,
01630
01631 Solomon_dissoc_rate_g + Solomon_dissoc_rate_s,
01632
01633 Solomon_elec_decay_g ,
01634
01635 Solomon_elec_decay_s
01636 );
01637 fprintf(io,
01638 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e",
01639
01640 hmi.UV_Cont_rel2_Habing_TH85_depth,
01641
01642 hmi.UV_Cont_rel2_Draine_DB96_depth,
01643
01644 secondaries.csupra[ipHYDROGEN][0]*0.93,
01645 sumlife/SDIV( sumpop ) ,
01646 hmi.H2_Solomon_dissoc_rate_BD96_H2g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth) ,
01647 Solomon_dissoc_rate_g/SDIV(hmi.UV_Cont_rel2_Habing_TH85_depth),
01648 Solomon_dissoc_rate_g/SDIV(hmi.UV_Cont_rel2_Habing_spec_depth),
01649 hmi.H2_rate_destroy);
01650 fprintf(io,
01651 "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01652 hmi.HeatH2Dish_TH85,
01653 hmi.HeatH2Dexc_TH85,
01654 HeatDiss,
01655 HeatDexc,
01656 thermal.htot);
01657 }
01658
01659 else if( (strcmp(chJOB , "H2so" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01660 {
01661
01662 const int nSOL = 100;
01663 double sum, one;
01664 long int jlosave[nSOL] , ivlosave[nSOL],
01665 iehisave[nSOL] ,jhisave[nSOL] , ivhisave[nSOL],
01666 nsave,
01667 ipOrdered[nSOL];
01668 int nFail;
01669 realnum fsave[nSOL], wlsave[nSOL];
01670
01671 fprintf(io,"%.5e\t%.3e",
01672
01673 radius.depth_mid_zone ,
01674
01675 Solomon_dissoc_rate_g +
01676 Solomon_dissoc_rate_s);
01677 sum = 0.;
01678
01679 if( lgEnabled && lgEvaluated )
01680 {
01681 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
01682 {
01683 qList::iterator Lo = ( (*tr).Lo() );
01684 if( (*Lo).n() > 0 )
01685 continue;
01686 sum += (*(*tr).Lo()).Pop() * (*tr).Emis().pump();
01687 }
01688
01689
01690 sum = SDIV( sum );
01691 nsave = 0;
01692
01693 const double frac = 0.01;
01694
01695 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
01696 {
01697 qList::iterator Lo = ( (*tr).Lo() );
01698 if( (*Lo).n() > 0 )
01699 continue;
01700 one = (*(*tr).Lo()).Pop() * (*tr).Emis().pump();
01701 if( one/sum > frac && nsave<nSOL)
01702 {
01703 qList::iterator Hi = ( (*tr).Hi() );
01704 fsave[nsave] = (realnum)(one/sum);
01705 jlosave[nsave] = (*Lo).J();
01706 ivlosave[nsave] = (*Lo).v();
01707 jhisave[nsave] = (*Hi).J();
01708 ivhisave[nsave] = (*Hi).v();
01709 iehisave[nsave] = (*Hi).n();
01710 wlsave[nsave] = (*tr).WLAng();
01711 ++nsave;
01712 }
01713 }
01714
01715
01716
01717 spsort(
01718
01719 fsave,
01720
01721 nsave,
01722
01723 ipOrdered,
01724
01725
01726 -1,
01727
01728 &nFail);
01729
01730
01731
01732
01733 fprintf(io,"\t%.3f\t%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e",
01734
01735 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f)/SDIV((Solomon_dissoc_rate_g * findspecieslocal("H2")->den +
01736 Solomon_dissoc_rate_s * findspecieslocal("H2*")->den) ),
01737
01738 (sum + secondaries.csupra[ipHYDROGEN][0]*2.02f) /SDIV((Solomon_dissoc_rate_g * H2_den_g+
01739 Solomon_dissoc_rate_s * H2_den_s) ),
01740 average_energy_g, average_energy_s,
01741 findspecieslocal("H2")->den/SDIV(H2_den_g), findspecieslocal("H2*")->den/SDIV(H2_den_s),
01742 H2_den_g/SDIV((*dense_total)), H2_den_s/SDIV((*dense_total))
01743 );
01744 for( long i=0; i<nsave; ++i )
01745 {
01746 long ip = ipOrdered[i];
01747
01748 fprintf(io,"\t%li\t%li\t%li\t%li\t%li\t%.3f\t%.3f",
01749 iehisave[ip],ivhisave[ip],jhisave[ip],ivlosave[ip] , jlosave[ip] , fsave[ip] , wlsave[ip] );
01750
01751 }
01752 fprintf(io,"\n");
01753 }
01754 }
01755
01756 else if( (strcmp(chJOB , "H2te" ) == 0) && (strcmp(chTime,"LAST") != 0) )
01757 {
01758
01759 double pop_ratio10,pop_ratio20,pop_ratio30,pop_ratio31,pop_ratio40;
01760 double T10,T20,T30,T31,T40;
01761
01762 double T10_sum,T20_sum,T30_sum,T31_sum,T40_sum;
01763 double pop_ratio10_sum,pop_ratio20_sum,pop_ratio30_sum,pop_ratio31_sum,pop_ratio40_sum;
01764 if( lgEnabled && nCall_this_zone )
01765 {
01766 double pop0 = states[0].Pop();
01767 double pop1 = states[1].Pop();
01768 double pop2 = states[2].Pop();
01769 double pop3 = states[3].Pop();
01770 double pop4 = states[4].Pop();
01771
01772 double energyK = states[1].energy().K() - states[0].energy().K();
01773
01774 pop_ratio10 = pop1/SDIV(pop0);
01775 pop_ratio10_sum = H2_X_colden[0][1]/SDIV(H2_X_colden[0][0]);
01776
01777 T10 = -170.5/log(SDIV(pop_ratio10) * H2_stat[0][0][0]/H2_stat[0][0][1]);
01778 T10_sum = -170.5/log(SDIV(pop_ratio10_sum) * H2_stat[0][0][0]/H2_stat[0][0][1]);
01779
01780 energyK = states[2].energy().K() - states[0].energy().K();
01781 pop_ratio20 = pop2/SDIV(pop0);
01782 T20 = -energyK/log(SDIV(pop_ratio20) * H2_stat[0][0][0]/H2_stat[0][0][2]);
01783
01784 pop_ratio20_sum = H2_X_colden[0][2]/SDIV(H2_X_colden[0][0]);
01785 T20_sum = -energyK/log(SDIV(pop_ratio20_sum) * H2_stat[0][0][0]/H2_stat[0][0][2]);
01786
01787 energyK = states[3].energy().K() - states[0].energy().K();
01788 pop_ratio30 = pop3/SDIV(pop0);
01789 T30 = -energyK/log(SDIV(pop_ratio30) * H2_stat[0][0][0]/H2_stat[0][0][3]);
01790
01791 pop_ratio30_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][0]);
01792 T30_sum = -energyK/log(SDIV(pop_ratio30_sum) * H2_stat[0][0][0]/H2_stat[0][0][3]);
01793
01794 energyK = states[3].energy().K() - states[1].energy().K();
01795 pop_ratio31 = pop3/SDIV(pop1);
01796 T31 = -energyK/log(SDIV(pop_ratio31) * H2_stat[0][0][1]/H2_stat[0][0][3]);
01797
01798 pop_ratio31_sum = H2_X_colden[0][3]/SDIV(H2_X_colden[0][1]);
01799 T31_sum = -energyK/log(SDIV(pop_ratio31_sum) * H2_stat[0][0][1]/H2_stat[0][0][3]);
01800
01801 energyK = states[4].energy().K() - states[0].energy().K();
01802 pop_ratio40 = pop4/SDIV(pop0);
01803 T40 = -energyK/log(SDIV(pop_ratio40) * H2_stat[0][0][0]/H2_stat[0][0][4]);
01804
01805 pop_ratio40_sum = H2_X_colden[0][4]/SDIV(H2_X_colden[0][0]);
01806 T40_sum = -energyK/log(SDIV(pop_ratio40_sum) * H2_stat[0][0][0]/H2_stat[0][0][4]);
01807 }
01808 else
01809 {
01810 pop_ratio10 = 0.;
01811 pop_ratio10_sum = 0.;
01812 T10 = 0.;
01813 T20 = 0.;
01814 T30 = 0.;
01815 T31 = 0.;
01816 T40 = 0.;
01817 T10_sum = 0.;
01818 T20_sum = 0.;
01819 T30_sum = 0.;
01820 T31_sum = 0.;
01821 T40_sum = 0.;
01822 }
01823
01824
01825 fprintf( io,
01826 "%.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" ,
01827
01828 radius.depth_mid_zone ,
01829
01830 (*dense_total)/dense.gas_phase[ipHYDROGEN] ,
01831
01832 pop_ratio10 ,
01833
01834 ortho_density / SDIV(para_density),
01835 T10,T20,T30,T31,T40,
01836 phycon.te ,
01837 hyperfine.Tspin21cm,T10_sum,T20_sum,T30_sum,T31_sum,T40_sum );
01838 }
01839 else if( (strcmp(chJOB , "H2ln" ) == 0) && (strcmp(chTime,"LAST") == 0) )
01840 {
01841
01842 double thresh;
01843 double renorm;
01844
01845
01846
01847 if( lgEnabled && LineSave.nsum > 0)
01848 {
01849 ASSERT( LineSave.ipNormWavL >= 0 );
01850
01851 if( LineSv[LineSave.ipNormWavL].SumLine[0] > SMALLFLOAT )
01852 renorm = LineSave.ScaleNormLine/
01853 LineSv[LineSave.ipNormWavL].SumLine[0];
01854 else
01855 renorm = 1.;
01856
01857 if( renorm > SMALLFLOAT )
01858 {
01859
01860
01861 thresh = thresh_punline_h2/(realnum)renorm;
01862 }
01863 else
01864 thresh = 0.f;
01865
01866
01867
01868 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
01869 {
01870 qList::iterator Hi = ( (*tr).Hi() );
01871 qList::iterator Lo = ( (*tr).Lo() );
01872 long iElecHi = (*Hi).n();
01873 long iVibHi = (*Hi).v();
01874 long iRotHi = (*Hi).J();
01875 long iElecLo = (*Lo).n();
01876 long iVibLo = (*Lo).v();
01877 long iRotLo = (*Lo).J();
01878 if( iElecHi >= nElecLevelOutput )
01879 continue;
01880 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > thresh )
01881 {
01882
01883
01884 double wl = (*tr).WLAng()/1e4;
01885 fprintf(io, "%li-%li %c(%li)", iVibHi, iVibLo, chMolBranch( iRotHi, iRotLo ), iRotLo );
01886 fprintf( io, "\t%ld\t%ld\t%ld\t%ld\t%ld\t%ld", iElecHi , iVibHi , iRotHi , iElecLo , iVibLo , iRotLo);
01887
01888 fprintf( io, "\t%.7f\t", wl );
01889
01890 prt_wl( io , (*tr).WLAng() );
01891
01892 fprintf( io, "\t%.3f\t%.3e",
01893 log10(MAX2(1e-37,H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo])) + radius.Conv2PrtInten,
01894 H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]*renorm );
01895
01896 fprintf( io, "\t%.3f", (*Hi).energy().K() );
01897
01898 fprintf( io, "\t%.3e", (*tr).Emis().Aul() * (*tr).EnergyErg() * (*(*tr).Hi()).g() );
01899 fprintf( io, "\n");
01900 }
01901 }
01902 }
01903 }
01904 else if( (strcmp(chJOB , "H2sp" ) == 0) )
01905 {
01906 fprintf( io, "PUT SOMETHING HERE!\n" );
01907 }
01908 return;
01909 }
01910
01911
01912
01913
01914
01915 long int cdH2_Line(
01916
01917 long int iElecHi,
01918 long int iVibHi ,
01919 long int iRotHi ,
01920
01921 long int iElecLo,
01922 long int iVibLo ,
01923 long int iRotLo ,
01924
01925 double *relint,
01926
01927 double *absint )
01928 {
01929 DEBUG_ENTRY( "cdH2_Line()" );
01930 return h2.getLine( iElecHi, iVibHi, iRotHi, iElecLo, iVibLo, iRotLo, relint, absint );
01931 };
01932
01933 long int diatomics::getLine( long iElecHi, long iVibHi, long iRotHi, long iElecLo, long iVibLo, long iRotLo, double *relint, double *absint )
01934 {
01935
01936 DEBUG_ENTRY( "diatomics::getline()" );
01937
01938
01939 *relint = 0.;
01940 *absint = 0.;
01941
01942
01943 if( iElecHi!=0 || iElecLo!=0 )
01944 {
01945 return 0;
01946 }
01947
01948 long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
01949 long ipLo = ipEnergySort[iElecLo][iVibLo][iRotLo];
01950
01951
01952 if( states[ipHi].energy().WN() < states[ipLo].energy().WN() )
01953 {
01954 return 0;
01955 }
01956
01957
01958 if( H2_lgOrtho[iElecHi][iVibHi][iRotHi] - H2_lgOrtho[iElecLo][iVibLo][iRotLo] != 0 )
01959 {
01960 return 0;
01961 }
01962
01963
01964 if( !lgH2_radiative[ipHi][ipLo] )
01965 {
01966 return 0;
01967 }
01968
01969 ASSERT( LineSave.ipNormWavL >= 0 );
01970
01971 if( LineSv[LineSave.ipNormWavL].SumLine[0] > 0. )
01972 {
01973 *relint = H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]/
01974 LineSv[LineSave.ipNormWavL].SumLine[0] * LineSave.ScaleNormLine;
01975 }
01976 else
01977 {
01978 *relint = 0.;
01979 }
01980
01981
01982 if( H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] > 0. )
01983 {
01984 *absint = log10(H2_SaveLine[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo]) +
01985 radius.Conv2PrtInten;
01986 }
01987 else
01988 {
01989
01990 *absint = -37.;
01991 }
01992
01993 return 1;
01994 }
01995
01996 void diatomics::set_numLevelsMatrix( long numLevels )
01997 {
01998 if( !lgREAD_DATA )
01999 nXLevelsMatrix = numLevels;
02000 }
02001