00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "collision.h"
00007 #include "physconst.h"
00008 #include "taulines.h"
00009 #include "lines_service.h"
00010 #include "opacity.h"
00011 #include "hmi.h"
00012 #include "ipoint.h"
00013 #include "grainvar.h"
00014 #include "h2.h"
00015 #include "h2_priv.h"
00016 #include "mole.h"
00017 #include "dense.h"
00018
00019
00020
00021
00022 int H2_nRot_add_ortho_para[N_ELEC] = {0 , 1 , 1 , 0, 1, 1 , 0};
00023
00024
00025 enum {DEBUG_ENER=false};
00026
00027
00028
00029
00030
00031 static double XVIB[H2_TOP] = { 0.70 , 0.60 , 0.20 };
00032 static double Xdust[H2_TOP] = { 0.04 , 0.10 , 0.40 };
00033
00034
00035
00036
00037 static const double energy_off = 0.273*FREQ_1EV/SPEEDLIGHT;
00038
00039 STATIC double EH2_eval( int ipH2, double DissocEnergy, double energy_wn )
00040 {
00041 double EH2_here;
00042 double Evm = DissocEnergy * XVIB[ipH2] + energy_off;
00043
00044 double Ev = (energy_wn+energy_off);
00045
00046
00047
00048
00049 double Edust = DissocEnergy * Xdust[ipH2] *
00050 ( 1. - ( (Ev - Evm) / (DissocEnergy+energy_off-Evm)) *
00051 ( (1.-Xdust[ipH2])/2.) );
00052 ASSERT( Edust >= 0. );
00053
00054
00055
00056 EH2_here = DissocEnergy +energy_off - Edust;
00057 ASSERT( EH2_here >= 0.);
00058
00059 return EH2_here;
00060 }
00061
00062
00063 STATIC double H2_vib_dist( int ipH2 , double EH2, double DissocEnergy, double energy_wn)
00064 {
00065 double G1[H2_TOP] = { 0.3 , 0.4 , 0.9 };
00066 double G2[H2_TOP] = { 0.6 , 0.6 , 0.4 };
00067 double Evm = DissocEnergy * XVIB[ipH2] + energy_off;
00068 double Fv;
00069 if( (energy_wn+energy_off) <= Evm )
00070 {
00071
00072 Fv = sexp( POW2( (energy_wn+energy_off - Evm)/(G1[ipH2]* Evm ) ) );
00073 }
00074 else
00075 {
00076
00077 Fv = sexp( POW2( (energy_wn+energy_off - Evm)/(G2[ipH2]*(EH2 - Evm ) ) ) );
00078 }
00079 return Fv;
00080 }
00081
00082 STATIC bool lgRadiative( const TransitionList::iterator &tr )
00083 {
00084
00085 qList::iterator Lo = (*tr).Lo() ;
00086
00087 if( (*Lo).n()==0 && (*tr).Emis().Aul() > 1.01e-30 )
00088 return true;
00089 else
00090 return false;
00091 }
00092
00093 STATIC bool compareEmis( const TransitionList::iterator& tr1, const TransitionList::iterator &tr2 )
00094 {
00095 if( lgRadiative( tr1 ) && !lgRadiative( tr2 ) )
00096 return true;
00097 else
00098 return false;
00099 }
00100
00101 bool compareEnergies( qStateProxy st1, qStateProxy st2 )
00102 {
00103 if( st1.energy().Ryd() <= st2.energy().Ryd() )
00104 return true;
00105 else
00106 return false;
00107 }
00108
00109
00110
00111 void diatomics::init(void)
00112 {
00113
00114
00115
00116 if( lgREAD_DATA || !lgEnabled )
00117 return;
00118
00119 DEBUG_ENTRY( "H2_Create()" );
00120
00121
00122 if( nTRACE )
00123 fprintf(ioQQQ," H2_Create called in DEBUG mode.\n");
00124
00125
00126 sp = findspecies( label.c_str() );
00127 ASSERT( sp != null_mole );
00128
00129
00130 string strSpStar = shortlabel + "*";
00131 sp_star = findspecies( strSpStar.c_str() );
00132 ASSERT( sp != null_mole );
00133
00134 mass_amu = sp->mole_mass / ATOMIC_MASS_UNIT;
00135 ASSERT( mass_amu > 0. );
00136
00137 H2_ReadDissocEnergies();
00138
00139
00140 H2_ReadEnergies();
00141
00142
00143 fixit();
00144
00145
00146 ASSERT( n_elec_states > 0 );
00147
00148 ipEnergySort.reserve(n_elec_states);
00149 for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
00150 {
00151 ipEnergySort.reserve(iElecHi,nVib_hi[iElecHi]+1);
00152 for( long iVibHi = 0; iVibHi <= nVib_hi[iElecHi]; ++iVibHi )
00153 {
00154 ipEnergySort.reserve(iElecHi,iVibHi,nRot_hi[iElecHi][iVibHi]+1);
00155 }
00156 }
00157 ipEnergySort.alloc();
00158
00159
00160 ipElec_H2_energy_sort.resize( states.size() );
00161 ipVib_H2_energy_sort.resize( states.size() );
00162 ipRot_H2_energy_sort.resize( states.size() );
00163 for( unsigned nEner = 0; nEner < states.size(); ++nEner )
00164 {
00165 long iElec = states[nEner].n();
00166 long iVib = states[nEner].v();
00167 long iRot = states[nEner].J();
00168 ipElec_H2_energy_sort[nEner] = iElec;
00169 ipVib_H2_energy_sort[nEner] = iVib;
00170 ipRot_H2_energy_sort[nEner] = iRot;
00171 ipEnergySort[iElec][iVib][iRot] = nEner;
00172
00173
00174
00175 }
00176
00177 H2_stat.alloc( ipEnergySort.clone() );
00178 H2_lgOrtho.alloc( ipEnergySort.clone() );
00179
00180
00181
00182 for( qList::iterator st = states.begin(); st != states.end(); ++st )
00183 {
00184
00185
00186
00187 if( is_odd( (*st).J() + H2_nRot_add_ortho_para[(*st).n()]) )
00188 {
00189
00190 H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] = true;
00191 H2_stat[(*st).n()][(*st).v()][(*st).J()] = 3.f*(2.f*(*st).J()+1.f);
00192 }
00193 else
00194 {
00195
00196 H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] = false;
00197 H2_stat[(*st).n()][(*st).v()][(*st).J()] = (2.f*(*st).J()+1.f);
00198 }
00199 (*st).g() = H2_stat[(*st).n()][(*st).v()][(*st).J()];
00200 }
00201
00202 if( nTRACE >= n_trace_full )
00203 {
00204 for( long iElec=0; iElec<n_elec_states; ++iElec)
00205 {
00206
00207 fprintf(ioQQQ,"\t(%li %li)", iElec , nLevels_per_elec[iElec] );
00208 }
00209 fprintf(ioQQQ,
00210 " H2_Create: there are %li electronic levels, in each level there are",
00211 n_elec_states);
00212 fprintf(ioQQQ,
00213 " for a total of %li levels.\n", (long int) states.size() );
00214 }
00215
00216
00217 for( long nEner=0; nEner<nLevels_per_elec[0]; ++nEner )
00218 {
00219 if( states[nEner].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
00220 break;
00221 nEner_H2_ground = nEner;
00222 }
00223
00224
00225 ++nEner_H2_ground;
00226
00227
00228
00229
00230
00231 if( nXLevelsMatrix<0 )
00232 {
00233 nXLevelsMatrix = nLevels_per_elec[0];
00234 }
00235 else if( nXLevelsMatrix > nLevels_per_elec[0] )
00236 {
00237 fprintf( ioQQQ,
00238 " The total number of levels used in the matrix solver was set to %li but there are only %li levels in X.\n Sorry.\n",
00239 nXLevelsMatrix ,
00240 nLevels_per_elec[0]);
00241 cdEXIT(EXIT_FAILURE);
00242 }
00243
00244
00245
00246 {
00247
00248 if( DEBUG_ENER )
00249 {
00250
00251
00252 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
00253 {
00254 fprintf(ioQQQ,"%li\t%li\t%li\t%.5e\n", (*st).n(), (*st).v(), (*st).J(), (*st).energy().WN() );
00255 }
00256
00257 cdEXIT(EXIT_SUCCESS);
00258 }
00259 }
00260
00261
00262 lgREAD_DATA = true;
00263
00264
00265 H2_populations_LTE.reserve(n_elec_states);
00266 pops_per_vib.reserve(n_elec_states);
00267 H2_dissprob.reserve(n_elec_states);
00268
00269 for( long iElec = 0; iElec<n_elec_states; ++iElec )
00270 {
00271
00272 if( nTRACE >= n_trace_full )
00273 fprintf(ioQQQ,"elec %li highest vib= %li\n", iElec , nVib_hi[iElec] );
00274
00275 ASSERT( nVib_hi[iElec] > 0 );
00276
00277
00278
00279 H2_populations_LTE.reserve(iElec,nVib_hi[iElec]+1);
00280 pops_per_vib.reserve(iElec,nVib_hi[iElec]+1);
00281
00282
00283
00284
00285 for( long iVib = 0; iVib <= nVib_hi[iElec]; ++iVib )
00286 {
00287
00288 H2_populations_LTE.reserve(iElec,iVib,nRot_hi[iElec][iVib]+1);
00289 }
00290 }
00291
00292 H2_populations_LTE.alloc();
00293 H2_old_populations.alloc( H2_populations_LTE.clone() );
00294 H2_Boltzmann.alloc( H2_populations_LTE.clone() );
00295 H2_rad_rate_out.alloc( H2_populations_LTE.clone() );
00296 H2_dissprob.alloc( H2_populations_LTE.clone() );
00297 H2_disske.alloc( H2_populations_LTE.clone() );
00298
00299 pops_per_vib.alloc();
00300
00301 H2_dissprob.zero();
00302 H2_disske.zero();
00303
00304
00305 H2_rad_rate_out.zero();
00306
00307
00308 H2_ipPhoto.reserve(nVib_hi[0]+1);
00309
00310
00311 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00312 {
00313
00314 H2_ipPhoto.reserve(iVib,nRot_hi[0][iVib]+1);
00315 }
00316
00317 H2_ipPhoto.alloc();
00318 H2_col_rate_in.alloc( H2_ipPhoto.clone() );
00319 H2_col_rate_out.alloc( H2_ipPhoto.clone() );
00320 H2_rad_rate_in.alloc( H2_ipPhoto.clone() );
00321 H2_coll_dissoc_rate_coef.alloc( H2_ipPhoto.clone() );
00322 H2_coll_dissoc_rate_coef_H2.alloc( H2_ipPhoto.clone() );
00323 H2_X_colden.alloc( H2_ipPhoto.clone() );
00324 H2_X_rate_from_elec_excited.alloc( H2_ipPhoto.clone() );
00325 H2_X_rate_to_elec_excited.alloc( H2_ipPhoto.clone() );
00326 H2_X_colden_LTE.alloc( H2_ipPhoto.clone() );
00327 H2_X_formation.alloc( H2_ipPhoto.clone() );
00328 H2_X_Hmin_back.alloc( H2_ipPhoto.clone() );
00329
00330 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00331 {
00332 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00333 {
00334
00335 H2_rad_rate_in[iVib][iRot] = -BIGFLOAT;
00336 H2_coll_dissoc_rate_coef[iVib][iRot] = -BIGFLOAT;
00337 H2_coll_dissoc_rate_coef_H2[iVib][iRot] = -BIGFLOAT;
00338 }
00339 }
00340
00341 H2_X_colden.zero();
00342 H2_X_colden_LTE.zero();
00343 H2_X_formation.zero();
00344 H2_X_Hmin_back.zero();
00345
00346 H2_X_rate_from_elec_excited.zero();
00347
00348 H2_X_rate_to_elec_excited.zero();
00349
00350
00351 H2_X_hminus_formation_distribution.reserve(nTE_HMINUS);
00352 for( long i=0; i<nTE_HMINUS; ++i )
00353 {
00354 H2_X_hminus_formation_distribution.reserve(i,nVib_hi[0]+1);
00355
00356 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00357 {
00358 H2_X_hminus_formation_distribution.reserve(i,iVib,nRot_hi[0][iVib]+1);
00359 }
00360 }
00361 H2_X_hminus_formation_distribution.alloc();
00362 H2_X_hminus_formation_distribution.zero();
00363 H2_Read_hminus_distribution();
00364
00365
00366
00367
00368
00369
00370
00371 H2_X_grain_formation_distribution.reserve(H2_TOP);
00372 for( long ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
00373 {
00374 H2_X_grain_formation_distribution.reserve(ipH2,nVib_hi[0]+1);
00375
00376
00377 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00378 {
00379 H2_X_grain_formation_distribution.reserve(ipH2,iVib,nRot_hi[0][iVib]+1);
00380 }
00381 }
00382 H2_X_grain_formation_distribution.alloc();
00383 H2_X_grain_formation_distribution.zero();
00384
00385 for( long iElec=0; iElec<n_elec_states; ++iElec )
00386 {
00387
00388 if( iElec > 0 )
00389 H2_ReadDissprob(iElec);
00390 }
00391
00392
00393
00394
00395 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
00396 {
00397 long iElec = (*st).n();
00398 if( iElec > 0 ) continue;
00399 long iVib = (*st).v();
00400 long iRot = (*st).J();
00401
00402
00403
00404
00405
00406
00407
00408
00410
00411 double thresh = (H2_DissocEnergies[1] - (*st).energy().WN())*WAVNRYD;
00412
00413
00414
00415 ASSERT( thresh > 0.74 );
00416 H2_ipPhoto[iVib][iRot] = ipoint(thresh);
00417 fixit();
00418 }
00419
00420 CollRateCoeff.reserve( nLevels_per_elec[0] );
00421 for( long j = 1; j < nLevels_per_elec[0]; ++j )
00422 {
00423 CollRateCoeff.reserve( j, j );
00424 for( long k = 0; k < j; ++k )
00425 {
00426 CollRateCoeff.reserve( j, k, N_X_COLLIDER );
00427 }
00428 }
00429 CollRateCoeff.alloc();
00430 CollRateCoeff.zero();
00431
00432 fixit();
00433 RateCoefTable.resize(ipNCOLLIDER);
00434
00435 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
00436 {
00437
00438 H2_CollidRateRead(nColl);
00439 }
00440
00441 CollRateErrFac.alloc( CollRateCoeff.clone() );
00442 CollRateErrFac.zero();
00443
00444
00445 if( lgH2_NOISE )
00446 {
00447 for( long ipHi = 1; ipHi < nLevels_per_elec[0]; ++ipHi )
00448 {
00449 for( long ipLo = 0; ipLo < ipHi; ++ipLo )
00450 {
00451 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
00452 {
00453
00454 realnum r = (realnum)RandGauss( xMeanNoise , xSTDNoise );
00455 CollRateErrFac[ipHi][ipLo][nColl] = pow((realnum)10.f,r);
00456 }
00457 }
00458 }
00459 }
00460
00461
00462 H2_X_source.resize( nLevels_per_elec[0] );
00463 H2_X_sink.resize( nLevels_per_elec[0] );
00464
00465 H2_X_coll_rate.reserve(nLevels_per_elec[0]);
00466
00467 for( long i=1; i<nLevels_per_elec[0]; ++i )
00468 {
00469 H2_X_coll_rate.reserve(i,i);
00470 }
00471 H2_X_coll_rate.alloc();
00472
00473 ipTransitionSort.reserve( states.size() );
00474 for( unsigned nEner = 1; nEner < states.size(); ++nEner )
00475 ipTransitionSort.reserve( nEner, nEner );
00476 ipTransitionSort.alloc();
00477 ipTransitionSort.zero();
00478 lgH2_radiative.alloc( ipTransitionSort.clone() );
00479 lgH2_radiative.zero();
00480
00481 trans.resize( (states.size() * (states.size()-1))/2 );
00482 AllTransitions.push_back(trans);
00483 qList initStates(1);
00484 TransitionList initlist("H2InitList",&initStates);
00485 vector<TransitionList::iterator> initptrs;
00486 initlist.resize(trans.size());
00487 initlist.states() = &states;
00488 initptrs.resize(trans.size());
00489
00490 {
00491 long lineIndex = 0;
00492 TransitionList::iterator tr = initlist.begin();
00493 for( unsigned ipHi=1; ipHi< states.size(); ++ipHi )
00494 {
00495 for( unsigned ipLo=0; ipLo<ipHi; ++ipLo )
00496 {
00497 (*tr).Junk();
00498 (*tr).setHi(ipHi);
00499 (*tr).setLo(ipLo);
00500 (*tr).Zero();
00501 initptrs[lineIndex] = tr;
00502 ipTransitionSort[ipHi][ipLo] = lineIndex;
00503 lineIndex++;
00504 ++tr;
00505 }
00506 }
00507 }
00508
00509
00510 H2_SaveLine.reserve(n_elec_states);
00511 for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
00512 {
00513 H2_SaveLine.reserve(iElecHi,nVib_hi[iElecHi]+1);
00514 for( long iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
00515 {
00516 H2_SaveLine.reserve(iElecHi,iVibHi,nRot_hi[iElecHi][iVibHi]+1);
00517 for( long iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00518 {
00519
00520
00521
00522
00523 long int lim_elec_lo = 0;
00524 H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,lim_elec_lo+1);
00525 for( long iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00526 {
00527 H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,iElecLo,nVib_hi[iElecLo]+1);
00528 for( long iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
00529 {
00530 H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,nRot_hi[iElecLo][iVibLo]+1);
00531 }
00532 }
00533 }
00534 }
00535 }
00536
00537 H2_SaveLine.alloc();
00538
00539
00540 H2_SaveLine.zero();
00541
00542
00543 for( long iElec=0; iElec<n_elec_states; ++iElec )
00544 {
00545
00546 H2_ReadTransprob(iElec,initlist);
00547 }
00548
00549
00550 stable_sort( initptrs.begin(), initptrs.end(), compareEmis );
00551 rad_end = trans.begin();
00552
00553 {
00554 TransitionList::iterator tr = trans.begin();
00555 vector<TransitionList::iterator>::iterator ptr = initptrs.begin();
00556 for (size_t i=0; i < initptrs.size(); ++i, ++tr, ++ptr)
00557 {
00558 (*tr).copy(*(*ptr));
00559 if (lgRadiative(tr))
00560 {
00561 rad_end = tr;
00562 }
00563 }
00564 }
00565 ++rad_end;
00566 ASSERT( rad_end != trans.end() );
00567
00568 for( unsigned i = 0; i < trans.size(); ++i )
00569 {
00570 qList::iterator Hi = trans[i].Hi();
00571 qList::iterator Lo = trans[i].Lo();
00572 ipTransitionSort[ ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()] ][ ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()] ] = i;
00573 trans[i].ipHi() = ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()];
00574 trans[i].ipLo() = ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()];
00575 }
00576
00577 findspecieslocal("H2")->levels = &states;
00578 findspecieslocal("H2")->lines = &trans;
00579
00580
00581 for( unsigned i = 0; i < trans.size(); ++i )
00582 {
00583 qList::iterator Hi = trans[i].Hi();
00584 qList::iterator Lo = trans[i].Lo();
00585 ipTransitionSort[ ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()] ][ ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()] ] = i;
00586
00587
00588 }
00589
00590
00591 for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
00592 {
00593 (*tr).EnergyWN() = (realnum)((*(*tr).Hi()).energy().WN() - (*(*tr).Lo()).energy().WN());
00594
00595 if( (*tr).EnergyWN() > SMALLFLOAT)
00596 (*tr).WLAng() = (realnum)(1.e8f/(*tr).EnergyWN() / RefIndex( (*tr).EnergyWN() ) );
00597
00598 (*tr).Coll().col_str() = 0.;
00599 }
00600
00601
00602 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00603 {
00604
00605
00606
00607 (*tr).resetEmis();
00608 (*tr).Emis().iRedisFun() = ipCRDW;
00609
00610 (*tr).Emis().TauIn() = opac.taumin;
00611 (*tr).Emis().TauCon() = opac.taumin;
00612
00613 (*tr).Emis().TauTot() = 1e20f;
00614
00615 (*tr).Emis().dampXvel() = (realnum)( (*tr).Emis().Aul()/(*tr).EnergyWN()/PI4);
00616 (*tr).Emis().gf() = (realnum)(GetGF( (*tr).Emis().Aul(),(*tr).EnergyWN(), (*(*tr).Hi()).g() ) );
00617
00618
00619 (*tr).Emis().opacity() = (realnum)( abscf( (*tr).Emis().gf(), (*tr).EnergyWN(), (*(*tr).Lo()).g()) );
00620
00621 qList::iterator Hi = (*tr).Hi() ;
00622
00623 if( (*Hi).n() == 0 )
00624 {
00625
00626
00627 (*tr).Emis().ColOvTot() = 1.;
00628 }
00629 else
00630 {
00631
00633 (*tr).Emis().ColOvTot() = 0.;
00634 }
00635
00636 fixit();
00637 if( (*Hi).n() > 0 )
00638 {
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649 (*tr).Coll().col_str() = (realnum)(
00650 pow3( (*tr).WLAng()*1e-8 ) *
00651 ( (*(*tr).Hi()).g()/(*(*tr).Lo()).g()) *
00652 (*tr).Emis().Aul() *
00653 log(3.)*HPLANCK/(160.f*pow3(PI)*0.5*1e-8*EN1EV)/6.e-17);
00654 ASSERT( (*tr).Coll().col_str()>0.);
00655 }
00656 }
00657
00658
00659 if( mole_global.lgStancil )
00660 Read_Mol_Diss_cross_sections();
00661
00662
00663
00664
00665
00666
00667
00668
00669 for( long ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
00670 {
00671 realnum sum = 0., sumj = 0., sumv = 0., sumo = 0., sump = 0.;
00672
00673
00674 if( hmi.chGrainFormPump == 'D' )
00675 {
00676 long iElec = 0;
00677
00678
00679
00680 double T_H2_FORM = 50000.;
00681 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00682 {
00683 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00684 {
00685
00686 H2_X_grain_formation_distribution[ipH2][iVib][iRot] =
00687
00688 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * (1.f+iVib) *
00689 (realnum)sexp( states[ ipEnergySort[iElec][iVib][iRot] ].energy().K()/T_H2_FORM );
00690 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00691 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00692 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00693 if( H2_lgOrtho[iElec][iVib][iRot] )
00694 {
00695 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00696 }
00697 else
00698 {
00699
00700 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00701 }
00702 }
00703 }
00704 }
00705 else if( hmi.chGrainFormPump == 'T' )
00706 {
00707
00708 double Xrot[H2_TOP] = { 0.14 , 0.15 , 0.15 };
00709 double Xtrans[H2_TOP] = { 0.12 , 0.15 , 0.25 };
00710
00711 double sumvib = 0.;
00712 double EH2;
00713 long iElec = 0;
00714
00715 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00716 {
00717 double vibdist;
00718 EH2 = EH2_eval( ipH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
00719 vibdist = H2_vib_dist( ipH2 , EH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
00720 sumvib += vibdist;
00721 }
00722
00723
00724 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00725 {
00726 double Ev = states[ ipEnergySort[iElec][iVib][0] ].energy().WN()+energy_off;
00727 double Fv;
00728
00729
00730 double Erot;
00731
00732
00733 EH2 = EH2_eval( ipH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
00734
00735
00736 Erot = (EH2 - Ev) * Xrot[ipH2] / (Xrot[ipH2] + Xtrans[ipH2]);
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761 if( Erot > 0. )
00762 {
00763
00764 Fv = H2_vib_dist( ipH2 , EH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() ) / sumvib;
00765
00766
00767 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00768 {
00769 double deltaE = states[ ipEnergySort[iElec][iVib][iRot] ].energy().WN() -
00770 states[ ipEnergySort[iElec][iVib][0] ].energy().WN();
00771
00772 double gaussian = sexp( POW2( (deltaE - Erot) / (0.5 * Erot) ) );
00773
00774 double thermal_dist = sexp( deltaE / Erot );
00775
00776
00777 double aver = ( gaussian + thermal_dist ) / 2.;
00778
00779
00780
00781
00782
00783
00784 ASSERT( gaussian <= 1. );
00785
00786 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = (realnum)(
00787
00788 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * Fv * (2.*iRot+1.) * aver );
00789
00790 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00791 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00792 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00793 if( H2_lgOrtho[iElec][iVib][iRot] )
00794 {
00795 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00796 }
00797 else
00798 {
00799 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00800 }
00801
00802 }
00803 }
00804 else
00805 {
00806
00807 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00808 {
00809 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 0.;
00810 }
00811 }
00812 }
00813 }
00814 else if( hmi.chGrainFormPump == 't' )
00815 {
00816
00817
00818
00819
00820
00821 double T_H2_FORM = 17329.;
00822 long iElec = 0;
00823 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00824 {
00825 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00826 {
00827
00828 H2_X_grain_formation_distribution[ipH2][iVib][iRot] =
00829
00830 H2_stat[0][iVib][iRot] *
00831 (realnum)sexp( states[ ipEnergySort[0][iVib][iRot] ].energy().K()/T_H2_FORM );
00832 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00833 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00834 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00835 if( H2_lgOrtho[iElec][iVib][iRot] )
00836 {
00837 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00838 }
00839 else
00840 {
00841
00842 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00843 }
00844 }
00845 }
00846 }
00847 else
00848 TotalInsanity();
00849
00850 if( nTRACE >= n_trace_full )
00851 fprintf(ioQQQ, "H2 form grains mean J= %.3f mean v = %.3f ortho/para= %.3f\n",
00852 sumj/sum , sumv/sum , sumo/sump );
00853
00854
00855
00856 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00857 {
00858 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00859 {
00860 H2_X_grain_formation_distribution[ipH2][iVib][iRot] /= sum;
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870 }
00871 }
00872 }
00873
00874 return;
00875 }