00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #define PRT_POPS false
00023
00024 #define LIM_H2_POP_LOOP 10
00025
00026
00027 #define H2_DISS_ALLISON_DALGARNO 6e-19f
00028 #include "cddefines.h"
00029 #include "cddrive.h"
00030 #include "physconst.h"
00031 #include "taulines.h"
00032 #include "atoms.h"
00033 #include "conv.h"
00034 #include "secondaries.h"
00035 #include "pressure.h"
00036 #include "trace.h"
00037 #include "hmi.h"
00038 #include "hextra.h"
00039 #include "mole.h"
00040 #include "rt.h"
00041 #include "radius.h"
00042 #include "ipoint.h"
00043 #include "phycon.h"
00044 #include "thermal.h"
00045 #include "dense.h"
00046 #include "rfield.h"
00047 #include "lines_service.h"
00048 #include "h2.h"
00049 #include "h2_priv.h"
00050
00051 static realnum collider_density[N_X_COLLIDER];
00052 static realnum collider_density_total_not_H2;
00053
00054 void diatomics::H2_X_sink_and_source( void )
00055 {
00056 DEBUG_ENTRY( "diatomics::H2_X_sink_and_source()" );
00057
00058
00059
00060 collider_density_total_not_H2 = collider_density[0] +
00061 collider_density[1] + collider_density[4] +
00062 dense.eden;
00063
00064 for( long ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
00065 {
00066 H2_X_source[ipHi] = 0.;
00067 H2_X_sink[ipHi] = 0.;
00068 }
00069
00070 double source_so_far = 0.;
00071 double sink_so_far = spon_diss_tot * H2_den_s;
00072 double pop_tot = 0.;
00073
00074 for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
00075 {
00076
00077 long iVibHi = ipVib_H2_energy_sort[ipHi];
00078 long iRotHi = ipRot_H2_energy_sort[ipHi];
00079
00080
00081
00082 H2_X_source[ipHi] += H2_X_formation[iVibHi][iRotHi];
00083
00084
00085 H2_X_sink[ipHi] += H2_X_Hmin_back[iVibHi][iRotHi];
00086
00087
00088
00089 H2_X_sink[ipHi] += collider_density_total_not_H2 *
00090 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] * lgColl_deexec_Calc;
00091
00092
00093 H2_X_sink[ipHi] += hmi.H2_total*
00094 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] * lgColl_deexec_Calc;
00095
00096
00097 if( mole_global.lgStancil )
00098 {
00099 H2_X_sink[ipHi] += Cont_Dissoc_Rate[0][iVibHi][iRotHi];
00100 }
00101 else
00102 H2_X_sink[ipHi] += rfield.flux_accum[H2_ipPhoto[iVibHi][iRotHi]-1]*H2_DISS_ALLISON_DALGARNO;
00103
00104 source_so_far += H2_X_source[ipHi];
00105 sink_so_far += H2_X_sink[ipHi] * states[ipHi].Pop();
00106 pop_tot += states[ipHi].Pop();
00107 }
00108
00109
00110 double sink_tot = mole.sink_rate_tot(sp) * pop_tot;
00111
00112 double sink_left = sink_tot - sink_so_far;
00113
00114 ASSERT( pop_tot > 1e-10 * (*dense_total) );
00115 sink_left /= pop_tot;
00116 if( sink_left >= 0. )
00117 {
00118 for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
00119 H2_X_sink[ipHi] += sink_left;
00120 }
00121 else
00122 {
00123 for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
00124 H2_X_sink[ipHi] *= sink_tot / sink_so_far;
00125 }
00126
00127 double sink_so_far_s = 0.;
00128 double pop_tot_s = 0.;
00129 for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
00130 {
00131 if( states[ipHi].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
00132 {
00133 sink_so_far_s += H2_X_sink[ipHi] * states[ipHi].Pop();
00134 pop_tot_s += states[ipHi].Pop();
00135 }
00136 }
00137
00138 fixit();
00139 double sink_tot_s = mole.sink_rate_tot(sp) * pop_tot_s + mole.sink_rate_tot(sp_star) * pop_tot_s;
00140
00141 double sink_left_s = sink_tot_s - sink_so_far_s;
00142
00143 if( pop_tot_s > 1e-30 * (*dense_total) )
00144 sink_left_s /= pop_tot_s;
00145 else
00146 sink_left_s = 0.;
00147 if( sink_left_s >= 0. )
00148 {
00149 for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
00150 {
00151 if( states[ipHi].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
00152 H2_X_sink[ipHi] += sink_left_s;
00153 }
00154 }
00155
00156 fixit();
00157 double source_tot = mole.source_rate_tot(sp) + mole.source_rate_tot(sp_star);
00158 double source_left = source_tot - source_so_far;
00159 if( source_left >= 0. )
00160 {
00161 for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
00162 {
00163 long iElec = states[ipHi].n();
00164 long iVib = states[ipHi].v();
00165 long iRot = states[ipHi].J();
00166 H2_X_source[ipHi] += source_left * H2_populations_LTE[iElec][iVib][iRot];
00167 }
00168 }
00169
00170 return;
00171 }
00172
00173
00174
00175 void diatomics::H2_X_coll_rate_evaluate( void )
00176 {
00177 DEBUG_ENTRY( "diatomics::H2_X_coll_rate_evaluate()" );
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 collider_density[0] = dense.xIonDense[ipHYDROGEN][0];
00188
00189
00190 collider_density[1] = dense.xIonDense[ipHELIUM][0];
00191
00192 collider_density[2] = h2.ortho_density_f;
00193
00194 collider_density[3] = h2.para_density_f;
00195
00196 collider_density[4] = dense.xIonDense[ipHYDROGEN][1];
00197
00198 collider_density[4] += (realnum)findspecieslocal("H3+")->den;
00199
00200 ASSERT( fp_equal( hmi.H2_total_f ,collider_density[2]+collider_density[3]) );
00201
00202 if( nTRACE >= n_trace_full )
00203 {
00204 fprintf(ioQQQ," Collider densities are:");
00205 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
00206 {
00207 fprintf(ioQQQ,"\t%.3e", collider_density[nColl]);
00208 }
00209 fprintf(ioQQQ,"\n");
00210 }
00211
00212 H2_X_coll_rate.zero();
00213
00214 for( long ipHi=0; ipHi < nLevels_per_elec[0]; ++ipHi )
00215 {
00216 if( lgColl_deexec_Calc )
00217 {
00218
00219 for( long ipLo=0; ipLo<ipHi; ++ipLo )
00220 {
00221
00222 double colldown = 0.;
00223 mr3ci CollRate = CollRateCoeff.begin(ipHi, ipLo);
00224 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
00225 {
00226
00227 colldown += CollRate[nColl]*collider_density[nColl];
00228 ASSERT( CollRate[nColl]*collider_density[nColl] >= 0. );
00229 }
00230
00231 H2_X_coll_rate[ipHi][ipLo] += colldown;
00232 }
00233 }
00234 }
00235
00236 return;
00237 }
00238
00239
00240 double diatomics::H2_itrzn( void )
00241 {
00242 if( lgEnabled && nH2_zone>0 )
00243 {
00244 return( (double)nH2_pops / (double)nH2_zone );
00245 }
00246 else
00247 {
00248 return 0.;
00249 }
00250 }
00251
00252
00253 void diatomics::H2_ContPoint( void )
00254 {
00255 if( !lgEnabled )
00256 return;
00257
00258 DEBUG_ENTRY( "H2_ContPoint()" );
00259
00260 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00261 {
00262 ASSERT( (*tr).Emis().Aul() > 0. );
00263 (*tr).ipCont() = ipLineEnergy( (*tr).EnergyRyd(), label.c_str(), 0 );
00264 (*tr).Emis().ipFine() = ipFineCont( (*tr).EnergyRyd());
00265 }
00266 return;
00267 }
00268
00269
00270
00271 double diatomics::H2_Accel(void)
00272 {
00273
00274 if( !lgEnabled )
00275 return 0.;
00276
00277 DEBUG_ENTRY( "H2_Accel()" );
00278
00279
00280
00281 double drive = 0.;
00282
00283 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00284 {
00285 ASSERT( (*tr).ipCont() > 0 );
00286 drive += (*tr).Emis().pump() * (*tr).Emis().PopOpc() * (*tr).EnergyErg();
00287 }
00288
00289 return drive;
00290 }
00291
00292
00293
00294 double diatomics::H2_RadPress(void)
00295 {
00296
00297 realnum smallfloat=SMALLFLOAT*10.f;
00298
00299
00300
00301 if( !lgEnabled || !nCall_this_zone )
00302 return 0.;
00303
00304 DEBUG_ENTRY( "H2_RadPress()" );
00305
00306 realnum doppler_width = GetDopplerWidth( mass_amu );
00307 double press = 0.;
00308
00309 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00310 {
00311 ASSERT( (*tr).ipCont() > 0 );
00312 if( (*(*tr).Hi()).Pop() > smallfloat && (*tr).Emis().PopOpc() > smallfloat )
00313 {
00314 press += PressureRadiationLine( *tr, doppler_width );
00315 }
00316 }
00317
00318 if(nTRACE >= n_trace_full)
00319 fprintf(ioQQQ,
00320 " H2_RadPress returns, radiation pressure is %.2e\n",
00321 press );
00322 return press;
00323 }
00324
00325 #if 0
00326
00327
00328 double diatomics::H2_InterEnergy(void)
00329 {
00330
00331 if( !lgEnabled )
00332 return 0.;
00333
00334 DEBUG_ENTRY( "H2_InterEnergy()" );
00335
00336 double energy = 0.;
00337 for( qList::iterator st = trans.states.begin(); st != trans.states.end(); ++st )
00338 energy += st->Pop() * st->energy();
00339
00340 return energy;
00341 }
00342 #endif
00343
00344
00345 void diatomics::H2_RT_diffuse(void)
00346 {
00347 if( !lgEnabled || !nCall_this_zone )
00348 return;
00349
00350 DEBUG_ENTRY( "H2_RT_diffuse()" );
00351
00352 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00353 {
00354 qList::iterator Hi = (*tr).Hi();
00355 if( (*Hi).n() > 0 )
00356 continue;
00357 (*tr).outline_resonance();
00358 }
00359
00360 return;
00361 }
00362
00363
00364 void diatomics::H2_RTMake( void )
00365 {
00366 if( !lgEnabled )
00367 return;
00368
00369 DEBUG_ENTRY( "H2_RTMake()" );
00370
00371 realnum doppler_width = GetDopplerWidth( mass_amu );
00372 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00373 {
00374
00375
00376
00377
00378 RT_line_one( *tr, false, 0.f, doppler_width );
00379 }
00380
00381 return;
00382 }
00383
00384
00385
00386 void diatomics::H2_RT_tau_inc(void)
00387 {
00388
00389
00390 if( !lgEnabled )
00391 return;
00392
00393 DEBUG_ENTRY( "H2_RT_tau_inc()" );
00394
00395
00396
00397
00398 if( nzone > 0 && nCall_this_iteration>2 )
00399 {
00400 renorm_max = MAX2( H2_renorm_chemistry , renorm_max );
00401 renorm_min = MIN2( H2_renorm_chemistry , renorm_min );
00402 }
00403
00404 realnum doppler_width = GetDopplerWidth( mass_amu );
00405 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00406 {
00407 ASSERT( (*tr).ipCont() > 0 );
00408 RT_line_one_tauinc( *tr,-9, -9, -9, -9, doppler_width );
00409 }
00410
00411 return;
00412 }
00413
00414
00415
00416 void diatomics::H2_LineZero( void )
00417 {
00418 if( !lgEnabled )
00419 return;
00420
00421 DEBUG_ENTRY( "H2_LineZero()" );
00422
00423 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00424 {
00425 (*tr).Zero();
00426 }
00427
00428 return;
00429 }
00430
00431
00432 void diatomics::H2_RT_tau_reset( void )
00433 {
00434 if( !lgEnabled )
00435 return;
00436
00437 DEBUG_ENTRY( "H2_RT_tau_reset()" );
00438
00439 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00440 {
00441 RT_line_one_tau_reset( *tr );
00442 }
00443
00444 return;
00445 }
00446
00447
00448 void diatomics::H2_Level_low_matrix(
00449
00450 realnum abundance )
00451 {
00452
00453 bool lgDoAs;
00454 int nNegPop;
00455 bool lgDeBug,
00456 lgZeroPop;
00457 double rot_cooling , dCoolDT;
00458
00459 DEBUG_ENTRY( "H2_Level_low_matrix()" );
00460
00461
00462 if( nXLevelsMatrix <= 1 )
00463 {
00464 return;
00465 }
00466
00467 if( lgFirst )
00468 {
00469
00470 if( nXLevelsMatrix > nLevels_per_elec[0] )
00471 {
00472
00473 fprintf( ioQQQ,
00474 " The total number of levels used in the matrix solver must be <= %li, the number of levels within X.\n Sorry.\n",
00475 nLevels_per_elec[0]);
00476 cdEXIT(EXIT_FAILURE);
00477 }
00478
00479 lgFirst = false;
00480
00481
00482
00483 ndimMalloced = nLevels_per_elec[0];
00484
00485 excit.resize( ndimMalloced );
00486 stat_levn.resize( ndimMalloced );
00487 pops.resize( ndimMalloced );
00488 create.resize( ndimMalloced );
00489 destroy.resize( ndimMalloced );
00490 depart.resize( ndimMalloced );
00491
00492 AulPump = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
00493 CollRate_levn = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
00494 AulDest = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
00495 AulEscp = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
00496 col_str = ((double **)MALLOC((size_t)(ndimMalloced)*sizeof(double *)));
00497 for( long i=0; i< ndimMalloced; ++i )
00498 {
00499 AulPump[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
00500 CollRate_levn[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
00501 AulDest[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
00502 AulEscp[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
00503 col_str[i] = ((double *)MALLOC((size_t)(ndimMalloced)*sizeof(double )));
00504 }
00505
00506
00507
00508 for( long j=0; j < ndimMalloced; j++ )
00509 {
00510
00511 stat_levn[j] = states[j].g();
00512
00513 excit[j] = states[j].energy().K();
00514 }
00515 }
00516
00517
00518
00519
00520 if( nXLevelsMatrix > ndimMalloced )
00521 {
00522 fprintf(ioQQQ," H2_Level_low_matrix has been called with the number of rotor levels greater than space allocated.\n");
00523 cdEXIT(EXIT_FAILURE);
00524 }
00525
00526
00527 for( long i=0; i < nXLevelsMatrix; i++ )
00528 {
00529 pops[i] = 0.;
00530 depart[i] = 0;
00531 for( long j=0; j < nXLevelsMatrix; j++ )
00532 {
00533 col_str[j][i] = 0.;
00534 }
00535 }
00536
00537
00538 if( nzone!=nzoneAsEval || iteration!=iterationAsEval || nXLevelsMatrix!=levelAsEval)
00539 {
00540 lgDoAs = true;
00541 nzoneAsEval = nzone;
00542 iterationAsEval = iteration;
00543 levelAsEval = nXLevelsMatrix;
00544 ASSERT( levelAsEval <= ndimMalloced );
00545 }
00546 else
00547 {
00548 lgDoAs = false;
00549 }
00550
00551
00552 if( lgDoAs )
00553 {
00554 for( long i=0; i < nXLevelsMatrix; i++ )
00555 {
00556 pops[i] = 0.;
00557 depart[i] = 0;
00558 for( long j=0; j < nXLevelsMatrix; j++ )
00559 {
00560 AulEscp[j][i] = 0.;
00561 AulDest[j][i] = 0.;
00562 AulPump[j][i] = 0.;
00563 CollRate_levn[j][i] = 0.;
00564 }
00565 }
00566 }
00567
00568
00569
00570 for( long ilo=0; ilo < nXLevelsMatrix; ilo++ )
00571 {
00572 long iRot = ipRot_H2_energy_sort[ilo];
00573 long iVib = ipVib_H2_energy_sort[ilo];
00574
00575
00576
00577
00578
00579 destroy[ilo] = H2_X_sink[ilo];
00580
00581
00582 create[ilo] = H2_X_source[ilo];
00583
00584
00585
00586 if( lgDoAs )
00587 {
00588 for( long ihi=ilo+1; ihi<nXLevelsMatrix; ++ihi )
00589 {
00590 long iRotHi = ipRot_H2_energy_sort[ihi];
00591 long iVibHi = ipVib_H2_energy_sort[ihi];
00592 ASSERT( states[ihi].energy().WN() <= states[nXLevelsMatrix-1].energy().WN() );
00593
00594 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
00595 {
00596 if( lgH2_radiative[ihi][ilo] )
00597 {
00598 const TransitionList::iterator&tr = trans.begin()+ ipTransitionSort[ihi][ilo] ;
00599 ASSERT( (*tr).ipCont() > 0 );
00600
00601
00602
00603
00604
00605
00606 AulEscp[ihi][ilo] = (*tr).Emis().Aul()*(
00607 (*tr).Emis().Pesc() +
00608 (*tr).Emis().Pelec_esc());
00609 AulDest[ihi][ilo] = (*tr).Emis().Aul()*(*tr).Emis().Pdest();
00610 AulPump[ilo][ihi] = (*tr).Emis().pump();
00611 }
00612 }
00613 }
00614 }
00615
00616 double rateout = 0.;
00617 double ratein = 0.;
00618
00619
00620 for( long ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
00621 {
00622 long iRotHi = ipRot_H2_energy_sort[ihi];
00623 long iVibHi = ipVib_H2_energy_sort[ihi];
00624 if( (abs(iRotHi-iRot)==2 || (iRotHi-iRot)==0 ) && (iVib<=iVibHi) )
00625 {
00626 if( lgH2_radiative[ihi][ilo] )
00627 {
00628 const TransitionList::iterator&tr = trans.begin()+ ipTransitionSort[ihi][ilo] ;
00629 ASSERT( (*tr).ipCont() > 0 );
00630
00631
00632
00633
00634 ratein +=
00635 (*(*tr).Hi()).Pop() * (
00636 (*tr).Emis().Aul()*( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() ) +
00637 (*tr).Emis().pump() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g() );
00638
00639 rateout +=
00640 (*tr).Emis().pump();
00641 }
00642 }
00643 }
00644
00645
00646 create[ilo] += ratein;
00647
00648
00649 destroy[ilo] += rateout;
00650
00651
00652
00653 create[ilo] += H2_X_rate_from_elec_excited[iVib][iRot];
00654
00655
00656 destroy[ilo] += H2_X_rate_to_elec_excited[iVib][iRot];
00657 }
00658
00659
00660 if( nTRACE >= n_trace_matrix )
00661 lgDeBug = true;
00662 else
00663 lgDeBug = false;
00664
00665
00666 for( long ilo=0; ilo < nXLevelsMatrix; ilo++ )
00667 {
00668 long iRot = ipRot_H2_energy_sort[ilo];
00669 long iVib = ipVib_H2_energy_sort[ilo];
00670 if(lgDeBug)fprintf(ioQQQ,"DEBUG H2_Level_low_matrix, ilo=%li",ilo);
00671 for( long ihi=ilo+1; ihi < nXLevelsMatrix; ihi++ )
00672 {
00673 long iRotHi = ipRot_H2_energy_sort[ihi];
00674 long iVibHi = ipVib_H2_energy_sort[ihi];
00675
00676 CollRate_levn[ihi][ilo] = H2_X_coll_rate[ihi][ilo];
00677
00678 if(lgDeBug)fprintf(ioQQQ,"\t%.1e",CollRate_levn[ihi][ilo]);
00679
00680
00681 CollRate_levn[ilo][ihi] = CollRate_levn[ihi][ilo]*
00682 H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
00683 H2_stat[0][iVibHi][iRotHi] /
00684 H2_stat[0][iVib][iRot];
00685 }
00686 if(lgDeBug)fprintf(ioQQQ,"\n");
00687
00688
00689
00690 for( long ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
00691 {
00692 long iRotHi = ipRot_H2_energy_sort[ihi];
00693 long iVibHi = ipVib_H2_energy_sort[ihi];
00694
00695
00696
00697 double ratein = H2_X_coll_rate[ihi][ilo];
00698 if(lgDeBug)fprintf(ioQQQ,"\t%.1e",ratein);
00699
00700
00701 double rateout = ratein *
00702 H2_Boltzmann[0][iVibHi][iRotHi]/SDIV(H2_Boltzmann[0][iVib][iRot])*
00703 H2_stat[0][iVibHi][iRotHi]/H2_stat[0][iVib][iRot];
00704
00705
00706 create[ilo] += ratein * states[ihi].Pop();
00707 destroy[ilo] += rateout;
00708 }
00709 if(lgDeBug)fprintf(ioQQQ,"\n");
00710 }
00711
00712
00713 {
00714 for( long ihi=2; ihi < nXLevelsMatrix; ihi++ )
00715 {
00716 long iVibHi = ipVib_H2_energy_sort[ihi];
00717 long iRotHi = ipRot_H2_energy_sort[ihi];
00718
00719
00720
00721
00722
00723 CollRate_levn[ihi][H2_lgOrtho[0][iVibHi][iRotHi]] += rate_grain_op_conserve;
00724 }
00725
00726
00727
00728 CollRate_levn[1][0] +=
00729 (realnum)(rate_grain_J1_to_J0);
00730 }
00731
00732
00733 for( long ihi=nXLevelsMatrix; ihi<nLevels_per_elec[0]; ++ihi )
00734 {
00735 long iVibHi = ipVib_H2_energy_sort[ihi];
00736 long iRotHi = ipRot_H2_energy_sort[ihi];
00737
00738
00739
00740 create[H2_lgOrtho[0][iVibHi][iRotHi]] += states[ihi].Pop() * rate_grain_op_conserve;
00741 }
00742
00743
00744 {
00745 enum {DEBUG_LOC=false};
00746 if( DEBUG_LOC || lgDeBug)
00747 {
00748 fprintf(ioQQQ,"DEBUG H2 matexcit");
00749 for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
00750 {
00751 fprintf(ioQQQ,"\t%li",ilo );
00752 }
00753 fprintf(ioQQQ,"\n");
00754 for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
00755 {
00756 fprintf(ioQQQ,"\t%.2e",excit[ihi] );
00757 }
00758 fprintf(ioQQQ,"\n");
00759 for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
00760 {
00761 fprintf(ioQQQ,"\t%.2e",stat_levn[ihi] );
00762 }
00763 fprintf(ioQQQ,"\n");
00764
00765 fprintf(ioQQQ,"AulEscp[n][]\\[][n] = Aul*Pesc\n");
00766 for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
00767 {
00768 fprintf(ioQQQ,"\t%li",ilo );
00769 }
00770 fprintf(ioQQQ,"\n");
00771 for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
00772 {
00773 fprintf(ioQQQ,"%li", ihi);
00774 for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
00775 {
00776 fprintf(ioQQQ,"\t%.2e",AulEscp[ilo][ihi] );
00777 }
00778 fprintf(ioQQQ,"\n");
00779 }
00780
00781 fprintf(ioQQQ,"AulPump [n][]\\[][n]\n");
00782 for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
00783 {
00784 fprintf(ioQQQ,"\t%li",ilo );
00785 }
00786 fprintf(ioQQQ,"\n");
00787 for(long ihi=0; ihi<nXLevelsMatrix;++ihi)
00788 {
00789 fprintf(ioQQQ,"%li", ihi);
00790 for(long ilo=0; ilo<nXLevelsMatrix; ++ilo )
00791 {
00792 fprintf(ioQQQ,"\t%.2e",AulPump[ihi][ilo] );
00793 }
00794 fprintf(ioQQQ,"\n");
00795 }
00796
00797 fprintf(ioQQQ,"CollRate_levn [n][]\\[][n]\n");
00798 for( long ilo=0; ilo<nXLevelsMatrix; ++ilo )
00799 {
00800 fprintf(ioQQQ,"\t%li",ilo );
00801 }
00802 fprintf(ioQQQ,"\n");
00803 for( long ihi=0; ihi<nXLevelsMatrix;++ihi)
00804 {
00805 fprintf(ioQQQ,"%li", ihi);
00806 for( long ilo=0; ilo<nXLevelsMatrix; ++ilo )
00807 {
00808 fprintf(ioQQQ,"\t%.2e",CollRate_levn[ihi][ilo] );
00809 }
00810 fprintf(ioQQQ,"\n");
00811 }
00812 fprintf(ioQQQ,"SOURCE");
00813 for( long ihi=0; ihi<nXLevelsMatrix;++ihi)
00814 {
00815 fprintf(ioQQQ,"\t%.2e",create[ihi]);
00816 }
00817 fprintf(ioQQQ,"\nSINK");
00818 for( long ihi=0; ihi<nXLevelsMatrix;++ihi)
00819 {
00820 fprintf(ioQQQ,"\t%.2e",destroy[ihi]);
00821 }
00822 fprintf(ioQQQ,"\n");
00823 }
00824 }
00825
00826 atom_levelN(
00827
00828 nXLevelsMatrix,
00829 abundance,
00830 &stat_levn[0],
00831 &excit[0],
00832 'K',
00833 &pops[0],
00834 &depart[0],
00835
00836 &AulEscp,
00837
00838 &col_str,
00839 &AulDest,
00840 &AulPump,
00841 &CollRate_levn,
00842 &create[0],
00843 &destroy[0],
00844
00845 true,
00846 &rot_cooling,
00847 &dCoolDT,
00848 " H2 ",
00849
00850 &nNegPop,
00851 &lgZeroPop,
00852 lgDeBug );
00853
00854 for( long i=0; i< nXLevelsMatrix; ++i )
00855 {
00856 states[i].Pop() = pops[i];
00857 }
00858
00859 if( 0 && nTRACE >= n_trace_full)
00860 {
00861
00862 fprintf(ioQQQ,"\n DEBUG H2_Level_lowJ dense_total: %.3e matrix rel pops\n", *dense_total);
00863 fprintf(ioQQQ,"v\tJ\tpop\n");
00864 for( long i=0; i<nXLevelsMatrix; ++i )
00865 {
00866 long iRot = ipRot_H2_energy_sort[i];
00867 long iVib = ipVib_H2_energy_sort[i];
00868 fprintf(ioQQQ,"%3li\t%3li\t%.3e\t%.3e\t%.3e\n",
00869 iVib , iRot , states[i].Pop(), create[i] , destroy[i]);
00870 }
00871 }
00872
00873
00874 if( nNegPop > 0 )
00875 {
00876 fprintf(ioQQQ," H2_Level_low_matrix called atom_levelN which returned negative populations.\n");
00877 ConvFail( "pops" , "H2" );
00878 }
00879 return;
00880 }
00881
00882
00883 void diatomics::H2_LevelPops( bool &lgPopsConverged, double &old_val, double &new_val )
00884 {
00885 DEBUG_ENTRY( "H2_LevelPops()" );
00886
00887
00888
00889 if( !lgEnabled || lgAbort )
00890 {
00891
00892 CalcPhotoionizationRate();
00893 return;
00894 }
00895
00896 double old_solomon_rate=-1.;
00897 long int n_pop_oscil = 0;
00898 int kase=0;
00899 bool lgConv_h2_soln,
00900 lgPopsConv_total,
00901 lgPopsConv_relative,
00902 lgHeatConv,
00903 lgSolomonConv,
00904 lgOrthoParaRatioConv;
00905 double quant_old=-1.,
00906 quant_new=-1.;
00907
00908 bool lgH2_pops_oscil=false,
00909 lgH2_pops_ever_oscil=false;
00910
00911
00912 double PopChgMax_relative=0. , PopChgMaxOld_relative=0., PopChgMax_total=0., PopChgMaxOld_total=0.;
00913 long int iRotMaxChng_relative , iVibMaxChng_relative,
00914 iRotMaxChng_total , iVibMaxChng_total,
00915 nXLevelsMatrix_save;
00916 double popold_relative , popnew_relative , popold_total , popnew_total;
00917
00918 char chReason[100];
00919
00920
00921 double converge_pops_relative=1e-2,
00922 converge_pops_total=1e-3,
00923 converge_ortho_para=1e-2;
00924
00925 double dens_rel_to_lim_react = mole.species[sp->index].xFracLim;
00926
00927 if(nTRACE >= n_trace_full )
00928 {
00929 fprintf(ioQQQ,
00930 "\n***************H2_LevelPops %s call %li this iteration, zone is %.2f, H2/H:%.e Te:%e ne:%e\n",
00931 label.c_str(),
00932 nCall_this_iteration,
00933 fnzone,
00934 dens_rel_to_lim_react,
00935 phycon.te,
00936 dense.eden
00937 );
00938 }
00939 else if( nTRACE >= n_trace_final )
00940 {
00941 static long int nzone_prt=-1;
00942 if( nzone!=nzone_prt )
00943 {
00944 nzone_prt = nzone;
00945 fprintf(ioQQQ,"DEBUG zone %li species %s rel_to_lim:%.3e Te:%.3e *ne:%.3e n(%s):%.3e\n",
00946 nzone,
00947 label.c_str(),
00948 dens_rel_to_lim_react,
00949 phycon.te,
00950 dense.eden,
00951 label.c_str(),
00952 *dense_total );
00953 }
00954 }
00955
00956 CalcPhotoionizationRate();
00957
00958
00959
00960 mole_H2_LTE();
00961
00962
00963
00964
00965 if( (!lgEvaluated && dens_rel_to_lim_react < H2_to_H_limit )
00966 || *dense_total < 1e-20 )
00967 {
00968
00969 if( nTRACE >= n_trace_full )
00970 fprintf(ioQQQ,
00971 " H2_LevelPops %s pops too small, not computing, set to LTE and return, H2/H is %.2e and H2_to_H_limit is %.2e.",
00972 label.c_str(),
00973 dens_rel_to_lim_react,
00974 H2_to_H_limit);
00975 H2_zero_pops_too_low();
00976 fixit();
00977
00978 return;
00979 }
00980
00981
00982
00983
00984
00985
00986
00987 lgEvaluated = true;
00988
00989
00990
00991
00992
00993 nXLevelsMatrix_save = nXLevelsMatrix;
00994 fixit();
00995 if( conv.lgSearch )
00996 {
00997 nXLevelsMatrix = nLevels_per_elec[0];
00998 }
00999
01000
01001
01002
01003
01004
01005
01006
01007 if( !fp_equal(phycon.te,TeUsedColl) )
01008 {
01009 H2_CollidRateEvalAll();
01010 TeUsedColl = phycon.te;
01011 }
01012
01013
01014
01015
01016 if( nCall_this_iteration==0 || lgLTE )
01017 {
01018
01019 if(nTRACE >= n_trace_full )
01020 fprintf(ioQQQ,"%s 1st call - using LTE level pops\n", label.c_str() );
01021
01022 H2_den_s = 0.;
01023 H2_den_g = 0.;
01024 for( qList::iterator st = states.begin(); st != states.end(); ++st )
01025 {
01026 long iElec = (*st).n();
01027 long iVib = (*st).v();
01028 long iRot = (*st).J();
01029
01030
01031 double pop = H2_populations_LTE[iElec][iVib][iRot] * (*dense_total);
01032 H2_old_populations[iElec][iVib][iRot] = pop;
01033 (*st).Pop() = pop;
01034
01035 if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
01036 {
01037 H2_den_s += pop;
01038 }
01039 else
01040 {
01041 H2_den_g += pop;
01042 }
01043 }
01044
01045
01046 ortho_density = 0.75* (*dense_total);
01047 para_density = 0.25* (*dense_total);
01048 {
01049 ortho_density_f = (realnum)ortho_density;
01050 para_density_f = (realnum)para_density;
01051 }
01052 ortho_para_current = ortho_density / SDIV( para_density );
01053 ortho_para_older = ortho_para_current;
01054 ortho_para_old = ortho_para_current;
01055
01056 frac_matrix = 1.;
01057 }
01058
01059
01060 {
01061 pops_per_vib.zero();
01062 fill_n( pops_per_elec, N_ELEC, 0. );
01063 double pop_total = 0.;
01064 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01065 {
01066 long iElec = (*st).n();
01067 long iVib = (*st).v();
01068
01069 pop_total += (*st).Pop();
01070 pops_per_elec[iElec] += (*st).Pop();
01071 pops_per_vib[iElec][iVib] += (*st).Pop();
01072 }
01073 ASSERT( pops_per_elec[0]>SMALLFLOAT );
01074
01075 H2_renorm_chemistry = *dense_total/ SDIV(pop_total);
01076 }
01077
01078 if(nTRACE >= n_trace_full)
01079 fprintf(ioQQQ,
01080 "%s H2_renorm_chemistry is %.4e, *dense_total is %.4e pops_per_elec[0] is %.4e\n",
01081 label.c_str(),
01082 H2_renorm_chemistry ,
01083 *dense_total,
01084 pops_per_elec[0]);
01085
01086
01087 for( qList::iterator st = states.begin(); st != states.end(); ++st )
01088 {
01089 long iElec = (*st).n();
01090 long iVib = (*st).v();
01091 long iRot = (*st).J();
01092
01093 (*st).Pop() *= H2_renorm_chemistry;
01094 H2_old_populations[iElec][iVib][iRot] = (*st).Pop();
01095 }
01096
01097 if(nTRACE >= n_trace_full )
01098 fprintf(ioQQQ,
01099 " H2 entry, old pops sumed to %.3e, renorm to htwo den of %.3e\n",
01100 pops_per_elec[0],
01101 *dense_total);
01102
01103
01104 fixit();
01105 if( conv.lgSearch )
01106 {
01107 converge_pops_relative *= 2.;
01108 converge_pops_total *= 3.;
01109 converge_ortho_para *= 3.;
01110 }
01111
01112 if( !conv.nTotalIoniz )
01113 mole_update_species_cache();
01114
01115
01116 mole_H2_form();
01117
01118
01119 H2_X_coll_rate_evaluate();
01120
01121
01122
01123 lgConv_h2_soln = false;
01124
01125 long loop_h2_pops = 0;
01126 {
01127 if( nzone != nzoneEval )
01128 {
01129 nzoneEval = nzone;
01130
01131 ++nH2_zone;
01132 }
01133 }
01134
01135 if( lgLTE )
01136 lgConv_h2_soln = true;
01137
01138
01139
01140
01141
01142
01143
01144 while( loop_h2_pops < LIM_H2_POP_LOOP-n_pop_oscil && !lgConv_h2_soln && !lgLTE )
01145 {
01146
01147 ++loop_h2_pops;
01148
01149 ++nH2_pops;
01150
01151
01152 H2_X_rate_from_elec_excited.zero();
01153
01154 H2_X_rate_to_elec_excited.zero();
01155 H2_rad_rate_out.zero();
01156 H2_rad_rate_in.zero();
01157 pops_per_vib.zero();
01158 fill_n( pops_per_elec, N_ELEC, 0. );
01159
01160 SolveExcitedElectronicLevels();
01161
01162
01163 H2_X_sink_and_source();
01164
01165
01166
01167 SolveSomeGroundElectronicLevels();
01168
01169
01170
01171 if( nXLevelsMatrix )
01172 {
01173 H2_Level_low_matrix(
01174
01175
01176 (*dense_total) * (realnum)frac_matrix );
01177 }
01178 if(nTRACE >= n_trace_full)
01179 {
01180 long iElecHi = 0;
01181 fprintf(ioQQQ," Rel pop(e=%li)" ,iElecHi);
01182 }
01183
01184
01185
01186 {
01187 pops_per_elec[0] = 0.;
01188 for( md2i it = pops_per_vib.begin(0); it != pops_per_vib.end(0); ++it )
01189 *it = 0.;
01190
01191 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01192 {
01193 long iElec = (*st).n();
01194 if( iElec > 0 ) continue;
01195 long iVib = (*st).v();
01196 pops_per_elec[iElec] += (*st).Pop();
01197 pops_per_vib[iElec][iVib] += (*st).Pop();
01198 }
01199
01200
01201 if(nTRACE >= n_trace_full)
01202 for( md2ci it = pops_per_vib.begin(0); it != pops_per_vib.end(0); ++it )
01203 fprintf(ioQQQ,"\t%.2e", *it/(*dense_total));
01204
01205 ASSERT( pops_per_elec[0] > SMALLFLOAT );
01206 }
01207
01208 if(nTRACE >= n_trace_full)
01209 {
01210 fprintf(ioQQQ,"\n");
01211
01212 fprintf(ioQQQ," Rel pop(0,J)");
01213 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01214 {
01215 long iElec = (*st).n();
01216 if( iElec > 0 ) continue;
01217 long iVib = (*st).v();
01218 if( iVib > 0 ) continue;
01219 fprintf(ioQQQ,"\t%.2e", (*st).Pop()/(*dense_total) );
01220 }
01221 fprintf(ioQQQ,"\n");
01222 }
01223
01224
01225
01226 {
01227 double sum_pops_matrix = 0.;
01228 for( long i=0; i<nXLevelsMatrix; ++i )
01229 {
01230 sum_pops_matrix += states[i].Pop();
01231 }
01232
01233
01234
01235 frac_matrix = sum_pops_matrix / SDIV(pops_per_elec[0]);
01236 }
01237
01238 {
01239 double pop_total = 0.;
01240 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01241 pop_total += (*st).Pop();
01242
01243
01244 double H2_renorm_conserve = *dense_total/ SDIV(pop_total);
01245
01246
01247
01248
01249 pops_per_vib.zero();
01250 fill_n( pops_per_elec, N_ELEC, 0. );
01251 for( qList::iterator st = states.begin(); st != states.end(); ++st )
01252 {
01253 (*st).Pop() *= H2_renorm_conserve;
01254 long iElec = (*st).n();
01255 long iVib = (*st).v();
01256 pops_per_elec[iElec] += (*st).Pop();
01257 pops_per_vib[iElec][iVib] += (*st).Pop();
01258 }
01259 }
01260
01261
01262 PopChgMaxOld_relative = PopChgMax_relative;
01263 PopChgMaxOld_total = PopChgMax_total;
01264 PopChgMax_relative = 0.;
01265 PopChgMax_total = 0.;
01266 iRotMaxChng_relative =-1;
01267 iVibMaxChng_relative = -1;
01268 iRotMaxChng_total =-1;
01269 iVibMaxChng_total = -1;
01270 popold_relative = 0.;
01271 popnew_relative = 0.;
01272 popold_total = 0.;
01273 popnew_total = 0.;
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284 {
01285
01286
01287
01288
01289 double sumold = 0.;
01290
01291 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01292 {
01293 long iElec = (*st).n();
01294 long iVib = (*st).v();
01295 long iRot = (*st).J();
01296 double pop = states[ ipEnergySort[iElec][iVib][iRot] ].Pop();
01297 double rel_change;
01298
01299
01300 if( fabs( pop - H2_old_populations[iElec][iVib][iRot])/
01301
01302
01303
01304 SDIV(pop) > fabs(PopChgMax_relative) &&
01305
01306
01307
01308
01309
01310 pop/SDIV(*dense_total)>1e-6 )
01311 {
01312 PopChgMax_relative =
01313 (pop - H2_old_populations[iElec][iVib][iRot])/SDIV(pop);
01314 iRotMaxChng_relative = iRot;
01315 iVibMaxChng_relative = iVib;
01316 popold_relative = H2_old_populations[iElec][iVib][iRot];
01317 popnew_relative = pop;
01318 }
01319
01320
01321
01322
01323 rel_change = (pop - H2_old_populations[iElec][iVib][iRot])/SDIV(*dense_total);
01324
01325 if( fabs(rel_change) > fabs(PopChgMax_total) )
01326 {
01327 PopChgMax_total = rel_change;
01328 iRotMaxChng_total = iRot;
01329 iVibMaxChng_total = iVib;
01330 popold_total = H2_old_populations[iElec][iVib][iRot];
01331 popnew_total = pop;
01332 }
01333
01334 kase = -1;
01335
01336
01337
01338
01339
01340
01341 rel_change = fabs( H2_old_populations[iElec][iVib][iRot] - pop ) / SDIV( pop );
01342
01343
01344 if( rel_change > 3. && H2_old_populations[iElec][iVib][iRot] * pop > 0 )
01345 {
01346
01347 H2_old_populations[iElec][iVib][iRot] = pow( 10. ,
01348 log10(H2_old_populations[iElec][iVib][iRot])/2. +
01349 log10(pop)/2. );
01350 kase = 2;
01351 }
01352
01353
01354 else if( rel_change> 0.1 )
01355 {
01356 realnum frac_old=0.25f;
01357
01358 H2_old_populations[iElec][iVib][iRot] =
01359 frac_old*H2_old_populations[iElec][iVib][iRot] +
01360 (1.f-frac_old)*pop;
01361 kase = 3;
01362 }
01363 else
01364 {
01365
01366 H2_old_populations[iElec][iVib][iRot] = pop;
01367 kase = 4;
01368 }
01369 sumold += H2_old_populations[iElec][iVib][iRot];
01370 }
01371
01372
01373 double H2_renorm_conserve_init = *dense_total/sumold;
01374
01375
01376
01377
01378 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01379 {
01380 long iElec = (*st).n();
01381 long iVib = (*st).v();
01382 long iRot = (*st).J();
01383 H2_old_populations[iElec][iVib][iRot] *= H2_renorm_conserve_init;
01384 }
01385 }
01386
01387
01388 {
01389 ortho_density = 0.;
01390 para_density = 0.;
01391 H2_den_s = 0.;
01392 H2_den_g = 0.;
01393
01394 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01395 {
01396 long iElec = (*st).n();
01397 long iVib = (*st).v();
01398 long iRot = (*st).J();
01399 const double& pop = (*st).Pop();
01400
01401 if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
01402 {
01403 H2_den_s += pop;
01404 }
01405 else
01406 {
01407 H2_den_g += pop;
01408 }
01409 if( H2_lgOrtho[iElec][iVib][iRot] )
01410 {
01411 ortho_density += pop;
01412 }
01413 else
01414 {
01415 para_density += pop;
01416 }
01417 }
01418 ASSERT( fp_equal_tol( H2_den_s + H2_den_g, *dense_total, 1e-5 * (*dense_total) ) );
01419 }
01420
01421
01422 ortho_para_older = ortho_para_old;
01423 ortho_para_old = ortho_para_current;
01424 ortho_para_current = ortho_density / SDIV( para_density );
01425
01426
01427
01428 old_solomon_rate = Solomon_dissoc_rate_g;
01429
01430
01431
01432 H2_Solomon_rate();
01433
01434
01435
01436
01437 HeatChangeOld = HeatChange;
01438 HeatChange = HeatDexc_old - HeatDexc;
01439 {
01440
01441
01442 if( loop_h2_pops>2 && (
01443 (HeatChangeOld*HeatChange<0. ) ||
01444 (PopChgMax_relative*PopChgMaxOld_relative<0. ) ) )
01445 {
01446 lgH2_pops_oscil = true;
01447 if( loop_h2_pops > 6 )
01448 {
01449 loop_h2_oscil = loop_h2_pops;
01450 lgH2_pops_ever_oscil = true;
01451 ++n_pop_oscil;
01452 }
01453 }
01454 else
01455 {
01456 lgH2_pops_oscil = false;
01457
01458 if( loop_h2_pops - loop_h2_oscil > 4 )
01459 {
01460 lgH2_pops_ever_oscil = false;
01461 }
01462 }
01463 }
01464
01465
01466
01467 HeatDexc_old = HeatDexc;
01468 if(fabs(HeatDexc)/thermal.ctot > conv.HeatCoolRelErrorAllowed/10. ||
01469 HeatDexc==0. )
01470 H2_Cooling("H2lup");
01471
01472
01473 lgConv_h2_soln = true;
01474 lgPopsConv_total = true;
01475 lgPopsConv_relative = true;
01476 lgHeatConv = true;
01477 lgSolomonConv = true;
01478 lgOrthoParaRatioConv = true;
01479
01480
01481
01482 if( fabs(PopChgMax_relative)>converge_pops_relative )
01483 {
01484
01485 lgConv_h2_soln = false;
01486 lgPopsConv_relative = false;
01487
01488
01489 quant_old = PopChgMaxOld_relative;
01490
01491 quant_new = PopChgMax_relative;
01492
01493 strcpy( chReason , "rel pops changed" );
01494 }
01495
01496
01497
01498 else if( fabs(PopChgMax_total)>converge_pops_total)
01499 {
01500 lgConv_h2_soln = false;
01501 lgPopsConv_total = false;
01502
01503
01504 quant_old = PopChgMaxOld_total;
01505
01506 quant_new = PopChgMax_total;
01507
01508 strcpy( chReason , "tot pops changed" );
01509 }
01510
01511
01512
01513
01514
01515 else if( fabs(ortho_para_current-ortho_para_old) / SDIV(ortho_para_current)> converge_ortho_para )
01516
01517
01518 {
01519 lgConv_h2_soln = false;
01520 lgOrthoParaRatioConv = false;
01521 quant_old = ortho_para_old;
01522 quant_new = ortho_para_current;
01523 strcpy( chReason , "ortho/para ratio changed" );
01524 }
01525
01526
01527 else if( !thermal.lgTemperatureConstant &&
01528 fabs(HeatDexc-HeatDexc_old)/MAX2(thermal.ctot,thermal.htot) >
01529 conv.HeatCoolRelErrorAllowed/2.
01530
01531 )
01532 {
01533
01534
01535
01536 lgConv_h2_soln = false;
01537 lgHeatConv = false;
01538 quant_old = HeatDexc_old/MAX2(thermal.ctot,thermal.htot);
01539 quant_new = HeatDexc/MAX2(thermal.ctot,thermal.htot);
01540 strcpy( chReason , "heating changed" );
01541
01542
01543
01544 }
01545
01546
01547
01548
01549
01550 else if( rfield.lgInducProcess &&
01551
01552
01553 hmi.H2_frac_abund_set==0 &&
01554
01555
01556 fabs( Solomon_dissoc_rate_g - old_solomon_rate)/SDIV(hmi.H2_rate_destroy) >
01557 conv.EdenErrorAllowed/5.)
01558 {
01559 lgConv_h2_soln = false;
01560 lgSolomonConv = false;
01561 quant_old = old_solomon_rate;
01562 quant_new = Solomon_dissoc_rate_g;
01563 strcpy( chReason , "Solomon rate changed" );
01564 }
01565
01566
01567 if( !lgConv_h2_soln )
01568 {
01569
01570
01571
01572 if( PRT_POPS || nTRACE >=n_trace_iterations )
01573 {
01574
01575
01576
01577
01578 fprintf(ioQQQ," %s loop %3li no conv oscl?%c why:%s ",
01579 label.c_str(),
01580 loop_h2_pops,
01581 TorF(lgH2_pops_ever_oscil),
01582 chReason );
01583 if( !lgPopsConv_relative )
01584 fprintf(ioQQQ," PopChgMax_relative:%.4e v:%li J:%li old:%.4e new:%.4e",
01585 PopChgMax_relative,
01586 iVibMaxChng_relative,
01587 iRotMaxChng_relative ,
01588 popold_relative ,
01589 popnew_relative );
01590 else if( !lgPopsConv_total )
01591 fprintf(ioQQQ," PopChgMax_total:%.4e v:%li J:%li old:%.4e new:%.4e",
01592 PopChgMax_total,
01593 iVibMaxChng_total,
01594 iRotMaxChng_total ,
01595 popold_total ,
01596 popnew_total );
01597 else if( !lgHeatConv )
01598 fprintf(ioQQQ," heat:%.4e old:%.4e new:%.4e",
01599 (HeatDexc-HeatDexc_old)/MAX2(thermal.ctot,thermal.htot),
01600 quant_old ,
01601 quant_new);
01602
01603 else if( !lgSolomonConv )
01604 fprintf(ioQQQ," d(sol rate)/tot dest\t%2e",(old_solomon_rate - Solomon_dissoc_rate_g)/SDIV(hmi.H2_rate_destroy));
01605 else if( !lgOrthoParaRatioConv )
01606 fprintf(ioQQQ," current, old, older ratios are %.4e %.4e %.4e",
01607 ortho_para_current , ortho_para_old, ortho_para_older );
01608 else
01609 TotalInsanity();
01610 fprintf(ioQQQ,"\n");
01611 }
01612 }
01613
01614
01615 if( trace.nTrConvg >= 5 )
01616 {
01617 fprintf( ioQQQ,
01618 " H2 5lev %li Conv?%c",
01619 loop_h2_pops ,
01620 TorF(lgConv_h2_soln) );
01621
01622 if( fabs(PopChgMax_relative)>0.1 )
01623 fprintf(ioQQQ," pops, rel chng %.3e",PopChgMax_relative);
01624 else
01625 fprintf(ioQQQ," rel heat %.3e rel chng %.3e H2 heat/cool %.2e",
01626 HeatDexc/thermal.ctot ,
01627 fabs(HeatDexc-HeatDexc_old)/thermal.ctot ,
01628 HeatDexc/thermal.ctot);
01629
01630 fprintf( ioQQQ,
01631 " Oscil?%c Ever Oscil?%c",
01632 TorF(lgH2_pops_oscil) ,
01633 TorF(lgH2_pops_ever_oscil) );
01634 fprintf(ioQQQ,"\n");
01635 }
01636
01637 if( nTRACE >= n_trace_full )
01638 {
01639 fprintf(ioQQQ,
01640 "H2 loop\t%li\tkase pop chng\t%i\tchem renorm fac\t%.4e\tortho/para ratio:\t%.3e\tfrac of pop in matrix: %.3f\n",
01641 loop_h2_pops,
01642 kase,
01643 H2_renorm_chemistry,
01644 ortho_density / para_density ,
01645 frac_matrix);
01646
01647
01648 if( iVibMaxChng_relative>=0 && iRotMaxChng_relative>=0 && PopChgMax_relative>1e-10 )
01649 fprintf(ioQQQ,
01650 "end loop %li H2 max rel chng=%.3e from %.3e to %.3e at v=%li J=%li\n\n",
01651 loop_h2_pops,
01652 PopChgMax_relative ,
01653 H2_old_populations[0][iVibMaxChng_relative][iRotMaxChng_relative],
01654 states[ ipEnergySort[0][iVibMaxChng_relative][iRotMaxChng_relative] ].Pop(),
01655 iVibMaxChng_relative , iRotMaxChng_relative
01656 );
01657 }
01658 }
01659
01660
01661
01662 if( !lgConv_h2_soln && !conv.lgSearch )
01663 {
01664 conv.lgConvPops = false;
01665 lgPopsConverged = false;
01666 old_val = quant_old;
01667 new_val = quant_new;
01668 }
01669
01670 for( qList::iterator st = states.begin(); st != states.end(); ++st )
01671 {
01672 ASSERT( (*st).Pop() >= 0. );
01673 }
01674
01675 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
01676 {
01677
01678 (*tr).Coll().cool() = 0.;
01679 (*tr).Coll().heat() = 0.;
01680
01681 (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
01682
01683
01684 (*tr).Emis().phots() = (*tr).Emis().Aul() * ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc()) * (*(*tr).Hi()).Pop();
01685
01686
01687 (*tr).Emis().xIntensity() = (*tr).Emis().phots() * (*tr).EnergyErg();
01688
01689 }
01690
01691 average_energy_g = 0.;
01692 average_energy_s = 0.;
01693
01694 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01695 {
01696 double popTimesE = (*st).Pop() * (*st).energy().WN();
01697 if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
01698 average_energy_s += popTimesE;
01699 else
01700 average_energy_g += popTimesE;
01701 }
01702
01703 average_energy_g /= H2_den_g;
01704 if( H2_den_s > 1e-30 * (*dense_total) )
01705 average_energy_s /= H2_den_s;
01706 else
01707 average_energy_s = 0.;
01708
01709
01710
01711
01712 photodissoc_BigH2_H2s = 0.;
01713 photodissoc_BigH2_H2g = 0.;
01714
01715 Average_collH_dissoc_g = 0.;
01716 Average_collH_dissoc_s = 0.;
01717 Average_collH2_dissoc_g = 0.;
01718 Average_collH2_dissoc_s = 0.;
01719
01720 rel_pop_LTE_g =0.;
01721 rel_pop_LTE_s = 0.;
01722
01723 double exp_disoc = sexp(H2_DissocEnergies[0]/phycon.te_wn);
01724
01725
01726
01727
01728
01729 {
01730 static long ip_cut_off = -1;
01731 if( ip_cut_off < 0 )
01732 {
01733
01734 ip_cut_off = ipoint( 1.14 );
01735 }
01736
01737
01738
01739 double flux_accum_photodissoc_BigH2_H2s = 0;
01740 fixit();
01741 long ip_H2_level = ipoint( 1.07896 - 2.5 / EVRYD);
01742 for( long i= ip_H2_level; i < ip_cut_off; ++i )
01743 {
01744 flux_accum_photodissoc_BigH2_H2s += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
01745 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1] );
01746 }
01747
01748
01749 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
01750 {
01751 long iElec = (*st).n();
01752 if( iElec > 0 ) continue;
01753 long iVib = (*st).v();
01754 long iRot = (*st).J();
01755 const double &pop = (*st).Pop();
01756 fixit();
01757 const double mass_stat_factor = 3.634e-5/(2*2);
01758
01759
01760
01761
01762
01763 if( (*st).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
01764 {
01765 double arg_ratio;
01766 photodissoc_BigH2_H2s += pop * flux_accum_photodissoc_BigH2_H2s;
01767
01768
01769 Average_collH_dissoc_s += pop * H2_coll_dissoc_rate_coef[iVib][iRot];
01770 Average_collH2_dissoc_s += pop * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
01771
01772
01773 arg_ratio = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]);
01774 if( arg_ratio > 0. )
01775 {
01776
01777 rel_pop_LTE_s += SAHA/SDIV(phycon.te32*arg_ratio)*
01778 H2_stat[0][iVib][iRot] * mass_stat_factor;
01779 }
01780 }
01781 else
01782 {
01783 double arg_ratio;
01784
01785 double flux_accum_photodissoc_BigH2_H2g = 0;
01786
01787 ip_H2_level = ipoint( 1.07896 - (*st).energy().Ryd() );
01788
01789 for( long i= ip_H2_level; i < ip_cut_off; ++i )
01790 {
01791 flux_accum_photodissoc_BigH2_H2g += ( rfield.flux[0][i-1] + rfield.ConInterOut[i-1]+
01792 rfield.outlin[0][i-1]+ rfield.outlin_noplot[i-1] );
01793 }
01794
01795 photodissoc_BigH2_H2g += pop * flux_accum_photodissoc_BigH2_H2g;
01796
01797
01798 average_energy_g += (pop * (*st).energy().WN() );
01799
01800
01801 Average_collH_dissoc_g += pop * H2_coll_dissoc_rate_coef[iVib][iRot];
01802 Average_collH2_dissoc_g += pop * H2_coll_dissoc_rate_coef_H2[iVib][iRot];
01803
01804
01805 arg_ratio = exp_disoc/SDIV(H2_Boltzmann[0][iVib][iRot]);
01806 if( arg_ratio > 0. )
01807 {
01808 rel_pop_LTE_g += SAHA/SDIV(phycon.te32*arg_ratio)*
01809 H2_stat[0][iVib][iRot] * mass_stat_factor;
01810 }
01811 }
01812 }
01813 }
01814
01815
01816
01817
01818
01819
01820
01821
01822 if( H2_den_g > SMALLFLOAT )
01823 {
01824 Average_collH_dissoc_g /= SDIV(H2_den_g);
01825 Average_collH2_dissoc_g /= SDIV(H2_den_g);
01826 photodissoc_BigH2_H2g *= H2_DISS_ALLISON_DALGARNO / SDIV(H2_den_g);
01827 }
01828 else
01829 {
01830 Average_collH_dissoc_g = 0.;
01831 Average_collH2_dissoc_g = 0.;
01832 photodissoc_BigH2_H2g = 0.;
01833 }
01834 if( H2_den_s > SMALLFLOAT )
01835 {
01836 Average_collH_dissoc_s /= SDIV(H2_den_s);
01837 Average_collH2_dissoc_s /= SDIV(H2_den_s);
01838 photodissoc_BigH2_H2s *= H2_DISS_ALLISON_DALGARNO / SDIV(H2_den_s);
01839 }
01840 else
01841 {
01842 Average_collH_dissoc_s = 0.;
01843 Average_collH2_dissoc_s = 0.;
01844 photodissoc_BigH2_H2s = 0.;
01845 }
01846
01847
01848
01849 H2_Calc_Average_Rates();
01850
01851 if( nTRACE )
01852 {
01853 fprintf(ioQQQ," H2_LevelPops exit1 %8.2f loops:%3li H2/H:%.3e Sol dis old %.3e new %.3e Sol dis star %.3e g-to-s %.3e photodiss star %.3e\n",
01854 fnzone ,
01855 loop_h2_pops ,
01856 dens_rel_to_lim_react,
01857 old_solomon_rate,
01858 Solomon_dissoc_rate_g,
01859 Solomon_dissoc_rate_s,
01860 gs_rate(),
01861 photodissoc_BigH2_H2s );
01862 }
01863
01864
01865
01866
01867
01868
01869 ++nCall_this_iteration;
01870
01871
01872
01873 ++nCall_this_zone;
01874
01875
01876
01877
01878 nXLevelsMatrix = nXLevelsMatrix_save;
01879
01880
01881
01882
01883 if( nCall_this_iteration && nzone != nzone_nlevel_set )
01884 {
01885
01886 const double FRAC = 0.99999;
01887
01888 double sum_pop = 0.;
01889 long nEner = 0;
01890 long iElec = 0;
01891 const bool PRT = false;
01892 if( PRT ) fprintf(ioQQQ,"DEBUG pops ");
01893 while( nEner < nLevels_per_elec[0] && sum_pop/(*dense_total) < FRAC )
01894 {
01895
01896 ASSERT( iElec == ipElec_H2_energy_sort[nEner] );
01897 long iVib = ipVib_H2_energy_sort[nEner];
01898 long iRot = ipRot_H2_energy_sort[nEner];
01899 sum_pop += H2_old_populations[iElec][iVib][iRot];
01900 if( PRT ) fprintf(ioQQQ,"\t%.3e ", H2_old_populations[iElec][iVib][iRot]);
01901 ++nEner;
01902 }
01903 if( PRT ) fprintf(ioQQQ,"\n");
01904 nzone_nlevel_set = nzone;
01905
01906
01907
01908 }
01909
01910 return;
01911 }
01912
01913
01914 void diatomics::SolveExcitedElectronicLevels( void )
01915 {
01916 DEBUG_ENTRY( "diatomics::SolveExcitedElectronicLevels()" );
01917
01918 multi_arr<double,3> rate_in;
01919 rate_in.alloc( H2_rad_rate_out.clone() );
01920 rate_in.zero();
01921 spon_diss_tot = 0.;
01922 double CosmicRayHILyaExcitationRate = ( hmi.lgLeidenCRHack ) ? secondaries.x12tot : 0.;
01923
01924 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
01925 {
01926 qList::iterator Lo = (*tr).Lo();
01927 long iElecLo = (*Lo).n();
01928 long iVibLo = (*Lo).v();
01929 long iRotLo = (*Lo).J();
01930 qList::iterator Hi = (*tr).Hi();
01931 long iElecHi = (*Hi).n();
01932 if( iElecHi < 1 ) continue;
01933 long iVibHi = (*Hi).v();
01934 long iRotHi = (*Hi).J();
01935
01936
01937
01938
01939
01940
01941 double rate_up = (*tr).Emis().pump() + CosmicRayHILyaExcitationRate * (*tr).Coll().col_str();
01942 double rate_down =
01943 (*tr).Emis().Aul() * ( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() ) +
01944 rate_up * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
01945
01946
01947 ASSERT( H2_lgOrtho[iElecHi][iVibHi][iRotHi] == H2_lgOrtho[iElecLo][iVibLo][iRotLo] );
01948
01949
01950 rate_in[iElecHi][iVibHi][iRotHi] += H2_old_populations[iElecLo][iVibLo][iRotLo]*rate_up;
01951
01952
01953
01954 if( iElecLo==0 )
01955 H2_X_rate_to_elec_excited[iVibLo][iRotLo] += rate_up;
01956 H2_rad_rate_out[iElecLo][iVibLo][iRotLo] += rate_up;
01957
01958
01959
01960 H2_rad_rate_out[iElecHi][iVibHi][iRotHi] += rate_down;
01961 ASSERT( rate_up >= 0. && rate_down >= 0. );
01962 }
01963
01964 for( qList::iterator st = states.begin(); st != states.end(); ++st )
01965 {
01966 if( (*st).n() < 1 )
01967 continue;
01968
01969 long iElec = (*st).n();
01970 long iVib = (*st).v();
01971 long iRot = (*st).J();
01972
01973 H2_rad_rate_out[iElec][iVib][iRot] += H2_dissprob[iElec][iVib][iRot];
01974
01975
01976
01977
01978
01979 double pop = rate_in[iElec][iVib][iRot] / SDIV( H2_rad_rate_out[iElec][iVib][iRot] );
01980 (*st).Pop() = pop;
01981 spon_diss_tot += pop * H2_dissprob[iElec][iVib][iRot];
01982 if( H2_old_populations[iElec][iVib][iRot]==0. )
01983 H2_old_populations[iElec][iVib][iRot] = pop;
01984
01985 pops_per_vib[iElec][iVib] += pop;
01986
01987 pops_per_elec[iElec] += pop;
01988 }
01989
01990 fixit();
01991 if( H2_den_s > 1e-30 * (*dense_total) )
01992 spon_diss_tot /= H2_den_s;
01993 else
01994 spon_diss_tot = 0.;
01995
01996 if(nTRACE >= n_trace_full)
01997 {
01998 for( long iElec=1; iElec<n_elec_states; ++iElec )
01999 {
02000 fprintf(ioQQQ," Pop(e=%li):",iElec);
02001 for( md2i it = pops_per_vib.begin(iElec); it != pops_per_vib.end(iElec); ++it )
02002 fprintf( ioQQQ,"\t%.2e", *it/(*dense_total) );
02003 fprintf(ioQQQ,"\n");
02004 }
02005 }
02006
02007 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
02008 {
02009 qList::iterator Lo = (*tr).Lo();
02010 if( (*Lo).n() != 0 ) continue;
02011 qList::iterator Hi = (*tr).Hi();
02012 if( (*Hi).n() < 1 ) continue;
02013
02014
02015 double rate = (*Hi).Pop() *
02016 ((*tr).Emis().Aul() * ( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() ) +
02017 (*tr).Emis().pump() * (*Lo).g() / (*Hi).g());
02018 H2_X_rate_from_elec_excited[(*Lo).v()][(*Lo).J()] += rate;
02019 H2_rad_rate_in[(*Lo).v()][(*Lo).J()] += rate;
02020 }
02021
02022 return;
02023 }
02024
02025 void diatomics::SolveSomeGroundElectronicLevels( void )
02026 {
02027 DEBUG_ENTRY("diatomics::SolveSomeGroundElectronicLevels()");
02028
02029
02030 H2_col_rate_out.zero();
02031 H2_col_rate_in.zero();
02032
02033
02034 for( long ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi)
02035 {
02036
02037 long iElecHi = ipElec_H2_energy_sort[ipHi];
02038 ASSERT( iElecHi == 0 );
02039 long iVibHi = ipVib_H2_energy_sort[ipHi];
02040 long iRotHi = ipRot_H2_energy_sort[ipHi];
02041
02042 realnum H2stat = H2_stat[iElecHi][iVibHi][iRotHi];
02043 double H2boltz = H2_Boltzmann[iElecHi][iVibHi][iRotHi];
02044
02045 for( long ipLo=0; ipLo<ipHi; ++ipLo )
02046 {
02047 long iVibLo = ipVib_H2_energy_sort[ipLo];
02048 long iRotLo = ipRot_H2_energy_sort[ipLo];
02049
02050
02051 realnum colldn = H2_X_coll_rate[ipHi][ipLo];
02052
02053 realnum collup = colldn *
02054 H2stat / H2_stat[0][iVibLo][iRotLo] *
02055 H2boltz / SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
02056
02057 H2_col_rate_out[iVibHi][iRotHi] += colldn;
02058 H2_col_rate_in[iVibLo][iRotLo] += colldn * H2_old_populations[0][iVibHi][iRotHi];
02059
02060 H2_col_rate_out[iVibLo][iRotLo] += collup;
02061 H2_col_rate_in[iVibHi][iRotHi] += collup * H2_old_populations[0][iVibLo][iRotLo];
02062 }
02063 }
02064
02065
02066
02067
02068
02069
02070
02071 long nEner = nLevels_per_elec[0];
02072 while( (--nEner) >= nXLevelsMatrix )
02073 {
02074
02075
02076 long iElec = ipElec_H2_energy_sort[nEner];
02077 ASSERT( iElec == 0 );
02078 long iVib = ipVib_H2_energy_sort[nEner];
02079 long iRot = ipRot_H2_energy_sort[nEner];
02080
02081 if( nEner+1 < nLevels_per_elec[0] )
02082 ASSERT( states[nEner].energy().WN() < states[nEner+1].energy().WN() ||
02083 fp_equal( states[nEner].energy().WN(), states[nEner+1].energy().WN() ) );
02084
02085
02086
02087
02088 if( nEner >1 )
02089 {
02090 H2_col_rate_out[iVib][iRot] +=
02091
02092
02093 (realnum)(rate_grain_op_conserve);
02094
02095
02096
02097 H2_col_rate_in[0][H2_lgOrtho[0][iVib][iRot]] +=
02098
02099
02100
02101 (realnum)(rate_grain_op_conserve*H2_old_populations[0][iVib][iRot]);
02102 }
02103 else if( nEner == 1 )
02104 {
02105
02106
02107 H2_col_rate_out[0][1] +=
02108
02109
02110
02111 (realnum)(rate_grain_J1_to_J0);
02112
02113 H2_col_rate_in[0][0] +=
02114
02115
02116
02117 (realnum)(rate_grain_J1_to_J0 *H2_old_populations[0][0][1]);
02118 }
02119
02120 double pump_from_below = 0.;
02121 for( long ipLo = 0; ipLo<nEner; ++ipLo )
02122 {
02123 long iElecLo = ipElec_H2_energy_sort[ipLo];
02124 ASSERT( iElecLo == 0 );
02125 long iVibLo = ipVib_H2_energy_sort[ipLo];
02126 long iRotLo = ipRot_H2_energy_sort[ipLo];
02127 const TransitionList::iterator&tr = trans.begin() +ipTransitionSort[nEner][ipLo] ;
02128
02129
02130 if( ( abs(iRotLo-iRot) == 2 || iRotLo == iRot ) && (iVibLo <= iVib) && (*tr).ipCont() > 0 )
02131 {
02132 double rateone = (*tr).Emis().Aul() * ( (*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc() + (*tr).Emis().Pdest() );
02133
02134
02135 double CosmicRayHILyaExcitationRate = ( hmi.lgLeidenCRHack ) ? secondaries.x12tot : 0.;
02136 double pump_up = (*tr).Emis().pump() + CosmicRayHILyaExcitationRate * (*tr).Coll().col_str();
02137 pump_from_below += pump_up * H2_old_populations[iElecLo][iVibLo][iRotLo];
02138 rateone += pump_up * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
02139 ASSERT( rateone >=0 );
02140 H2_rad_rate_out[0][iVib][iRot] += rateone;
02141 H2_rad_rate_in[iVibLo][iRotLo] += rateone * H2_old_populations[iElec][iVib][iRot];
02142 }
02143 }
02144
02145
02146
02147 states[nEner].Pop() =
02148 (H2_col_rate_in[iVib][iRot]+ H2_rad_rate_in[iVib][iRot]+H2_X_source[nEner]+pump_from_below) /
02149 SDIV(H2_col_rate_out[iVib][iRot]+H2_rad_rate_out[0][iVib][iRot]+H2_X_sink[nEner]);
02150
02151 ASSERT( states[nEner].Pop() >= 0. );
02152 }
02153
02154 return;
02155 }
02156
02157
02158
02159
02160
02161 #if defined(__ICC) && defined(__i386)
02162 #pragma optimization_level 1
02163 #endif
02164 void diatomics::H2_Cooling(
02165
02166
02167
02168 const char *chRoutine)
02169 {
02170 DEBUG_ENTRY( "H2_Cooling()" );
02171
02172
02173
02174
02175
02176 if( !lgEnabled || !nCall_this_iteration )
02177 {
02178 HeatDexc = 0.;
02179 HeatDiss = 0.;
02180 HeatDexc_deriv = 0.;
02181 return;
02182 }
02183
02184 if( 0 )
02185 fprintf( ioQQQ, "DEBUG H2_Cooling called by %s.\n", chRoutine );
02186
02187 HeatDiss = 0.;
02188
02189 for( qList::iterator st = states.begin(); st != states.end(); ++st )
02190 {
02191 long iElec = (*st).n();
02192 long iVib = (*st).v();
02193 long iRot = (*st).J();
02194 HeatDiss += (*st).Pop() * H2_dissprob[iElec][iVib][iRot] * H2_disske[iElec][iVib][iRot];
02195 }
02196
02197
02198 HeatDiss *= EN1EV;
02199
02200
02201
02202 HeatDexc = 0.;
02203
02204
02205
02206
02207
02208 HeatDexc_deriv = 0.;
02209 long iElecHi = 0;
02210 long iElecLo = 0;
02211 for( long ipHi=1; ipHi<nLevels_per_elec[iElecHi]; ++ipHi )
02212 {
02213 long iVibHi = ipVib_H2_energy_sort[ipHi];
02214 long iRotHi = ipRot_H2_energy_sort[ipHi];
02215 realnum H2statHi = states[ipHi].g();
02216 double H2boltzHi = H2_Boltzmann[iElecHi][iVibHi][iRotHi];
02217 double H2popHi = states[ipHi].Pop();
02218 double ewnHi = states[ipHi].energy().WN();
02219
02220 for( long ipLo=0; ipLo<ipHi; ++ipLo )
02221 {
02222 long iVibLo = ipVib_H2_energy_sort[ipLo];
02223 long iRotLo = ipRot_H2_energy_sort[ipLo];
02224 double rate_dn_heat = 0.;
02225
02226
02227 mr3ci H2cr = CollRateCoeff.begin(ipHi, ipLo);
02228 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
02229
02230 rate_dn_heat += H2cr[nColl]*collider_density[nColl];
02231
02232
02233 double rate_up_cool = rate_dn_heat * states[ipLo].Pop() *
02234
02235 H2statHi / H2_stat[iElecLo][iVibLo][iRotLo] *
02236 H2boltzHi / SDIV( H2_Boltzmann[iElecLo][iVibLo][iRotLo] );
02237
02238 rate_dn_heat *= H2popHi;
02239
02240
02241
02242
02243
02244 double conversion = (ewnHi - states[ipLo].energy().WN() ) * ERG1CM;
02245 double heatone = rate_dn_heat * conversion;
02246 double coolone = rate_up_cool * conversion;
02247
02248 double oneline = heatone - coolone;
02249 HeatDexc += oneline;
02250
02251
02252
02253
02254 HeatDexc_deriv += (realnum)(oneline * ewnHi);
02255
02256
02257 ASSERT(
02258 (rate_up_cool==0 && rate_dn_heat==0) ||
02259 (states[ipHi].energy().WN() > states[ipLo].energy().WN()) );
02260 }
02261 }
02262
02263
02264 if( PRT_POPS )
02265 fprintf(ioQQQ,
02266 " DEBUG H2 heat fnzone\t%.2f\trenorm\t%.3e\tte\t%.4e\tdexc\t%.3e\theat/tot\t%.3e\n",
02267 fnzone ,
02268 H2_renorm_chemistry ,
02269 phycon.te ,
02270 HeatDexc,
02271 HeatDexc/thermal.ctot);
02272
02273
02274
02275 HeatDexc_deriv /= (realnum)POW2(phycon.te_wn);
02276
02277 if( nTRACE >= n_trace_full )
02278 fprintf(ioQQQ,
02279 " H2_Cooling Ctot\t%.4e\t HeatDiss \t%.4e\t HeatDexc \t%.4e\n" ,
02280 thermal.ctot ,
02281 HeatDiss ,
02282 HeatDexc );
02283
02284
02285
02286
02287
02288
02289
02290
02291 if( conv.lgSearch )
02292 {
02293 HeatDexc = 0.;
02294 HeatDexc_deriv = 0.;
02295 }
02296 return;
02297 }
02298
02299
02300
02301 double cdH2_colden( long iVib , long iRot )
02302 {
02303 diatomics& diatom = h2;
02304
02305
02306
02307
02308
02309
02310 if( iVib < 0 )
02311 {
02312 if( iRot==0 )
02313 {
02314
02315 return( diatom.ortho_colden + diatom.para_colden );
02316 }
02317 else if( iRot==1 )
02318 {
02319
02320 return diatom.ortho_colden;
02321 }
02322 else if( iRot==2 )
02323 {
02324
02325 return diatom.para_colden;
02326 }
02327 else
02328 {
02329 fprintf(ioQQQ," iRot must be 0 (total), 1 (ortho), or 2 (para), returning -1.\n");
02330 return -1.;
02331 }
02332 }
02333 else if( diatom.lgEnabled )
02334 {
02335 return diatom.GetXColden( iVib, iRot );
02336 }
02337
02338 else
02339 return -1;
02340 }
02341
02342 realnum diatomics::GetXColden( long iVib, long iRot )
02343 {
02344 DEBUG_ENTRY( "diatomics::GetXColden()" );
02345
02346
02347
02348 int iElec = 0;
02349 if( iRot <0 || iVib >nVib_hi[iElec] || iRot > nRot_hi[iElec][iVib])
02350 {
02351 fprintf(ioQQQ," iVib and iRot must lie within X, returning -2.\n");
02352 fprintf(ioQQQ," iVib must be <= %li and iRot must be <= %li.\n",
02353 nVib_hi[iElec],nRot_hi[iElec][iVib]);
02354 return -2.;
02355 }
02356 else
02357 return H2_X_colden[iVib][iRot];
02358 }
02359
02360
02361 void diatomics::H2_Colden( const char *chLabel )
02362 {
02363
02364 if( !lgEnabled )
02365 return;
02366
02367 DEBUG_ENTRY( "H2_Colden()" );
02368
02369 if( strcmp(chLabel,"ZERO") == 0 )
02370 {
02371
02372 H2_X_colden.zero();
02373 H2_X_colden_LTE.zero();
02374 }
02375
02376 else if( strcmp(chLabel,"ADD ") == 0 )
02377 {
02378
02379 for( qList::iterator st = states.begin(); st != states.end(); ++st )
02380 {
02381 long iElec = (*st).n();
02382 if( iElec > 0 ) continue;
02383 long iVib = (*st).v();
02384 long iRot = (*st).J();
02385
02386 H2_X_colden[iVib][iRot] += (realnum)( (*st).Pop() * radius.drad_x_fillfac);
02387
02388
02389 H2_X_colden_LTE[iVib][iRot] += (realnum)(H2_populations_LTE[0][iVib][iRot]*
02390 (*dense_total)*radius.drad_x_fillfac);
02391 }
02392 }
02393
02394
02395 else if( strcmp(chLabel,"PRIN") != 0 )
02396 {
02397 fprintf( ioQQQ, " H2_Colden does not understand the label %s\n",
02398 chLabel );
02399 cdEXIT(EXIT_FAILURE);
02400 }
02401
02402 return;
02403 }
02404
02405
02406 double diatomics::H2_DR(void)
02407 {
02408 return BIGFLOAT;
02409 }
02410
02411
02412 void diatomics::H2_RT_OTS( void )
02413 {
02414
02415 if( !lgEnabled || !nCall_this_zone )
02416 return;
02417
02418 DEBUG_ENTRY( "H2_RT_OTS()" );
02419
02420 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
02421 {
02422 qList::iterator Hi = (*tr).Hi();
02423 if( (*Hi).n() > 0 )
02424 continue;
02425
02426
02427 (*tr).Emis().ots() = (*(*tr).Hi()).Pop() * (*tr).Emis().Aul() * (*tr).Emis().Pdest();
02428
02429
02430 RT_OTS_AddLine( (*tr).Emis().ots(), (*tr).ipCont() );
02431 }
02432
02433 return;
02434 }
02435
02436 void diatomics::H2_Calc_Average_Rates( void )
02437 {
02438
02439
02440 double sumpop1 = 0.;
02441 double sumpopA1 = 0.;
02442 double sumpopcollH2O_deexcit = 0.;
02443 double sumpopcollH2p_deexcit = 0.;
02444 double sumpopcollH_deexcit = 0.;
02445 double popH2s = 0.;
02446 double sumpopcollH2O_excit = 0.;
02447 double sumpopcollH2p_excit = 0.;
02448 double sumpopcollH_excit = 0.;
02449 double popH2g = 0.;
02450
02451 for( qList::const_iterator stHi = states.begin(); stHi != states.end(); ++stHi )
02452 {
02453 long iElecHi = (*stHi).n();
02454 if( iElecHi > 0 ) continue;
02455 long iVibHi = (*stHi).v();
02456 long iRotHi = (*stHi).J();
02457 double ewnHi = (*stHi).energy().WN();
02458 for( qList::const_iterator stLo = states.begin(); stLo != stHi; ++stLo )
02459 {
02460 long iVibLo = (*stLo).v();
02461 long iRotLo = (*stLo).J();
02462 double ewnLo2 = (*stLo).energy().WN();
02463 if( ewnHi > ENERGY_H2_STAR && ewnLo2 < ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
02464 {
02465
02466
02467 if( H2_lgOrtho[0][iVibHi][iRotHi] == H2_lgOrtho[0][iVibLo][iRotLo] )
02468 {
02469 long ihi = ipEnergySort[0][iVibHi][iRotHi];
02470 long ilo = ipEnergySort[0][iVibLo][iRotLo];
02471 const TransitionList::iterator&tr =
02472 trans.begin()+ ipTransitionSort[ihi][ilo];
02473 double popHi = (*(*tr).Hi()).Pop();
02474 double popLo = (*(*tr).Lo()).Pop();
02475
02476
02477 popH2s += popHi;
02478 popH2g += popLo;
02479
02480
02481 sumpopcollH_deexcit += popHi * CollRateCoeff[ihi][ilo][0];
02482 sumpopcollH2O_deexcit += popHi * CollRateCoeff[ihi][ilo][2];
02483 sumpopcollH2p_deexcit += popHi * CollRateCoeff[ihi][ilo][3];
02484
02485 double temp = popLo *
02486 H2_stat[0][iVibHi][iRotHi] / H2_stat[0][iVibLo][iRotLo] *
02487 H2_Boltzmann[0][iVibHi][iRotHi] /SDIV( H2_Boltzmann[0][iVibLo][iRotLo] );
02488
02489
02490 sumpopcollH_excit += temp * CollRateCoeff[ihi][ilo][0];
02491 sumpopcollH2O_excit += temp * CollRateCoeff[ihi][ilo][2];
02492 sumpopcollH2p_excit += temp * CollRateCoeff[ihi][ilo][3];
02493
02494 if( lgH2_radiative[ihi][ilo] )
02495 {
02496 sumpop1 += popHi;
02497 sumpopA1 += popHi * (*tr).Emis().Aul();
02498 }
02499 }
02500 }
02501 }
02502 }
02503 Average_A = sumpopA1/SDIV(sumpop1);
02504
02505
02506 Average_collH2_deexcit = (sumpopcollH2O_deexcit+sumpopcollH2p_deexcit)/SDIV(popH2s);
02507 Average_collH2_excit = (sumpopcollH2O_excit+sumpopcollH2p_excit)/SDIV(popH2g);
02508 Average_collH_excit = sumpopcollH_excit/SDIV(popH2g);
02509 Average_collH_deexcit = sumpopcollH_deexcit/SDIV(popH2s);
02510
02511
02512
02513
02514
02515
02516
02517
02518 return;
02519 }
02520
02521 double diatomics::GetExcitedElecDensity( void )
02522 {
02523 double H2_sum_excit_elec_den = 0.;
02524 for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
02525 {
02526 if( iElecHi > 0 )
02527 H2_sum_excit_elec_den += pops_per_elec[iElecHi];
02528 }
02529
02530 return H2_sum_excit_elec_den;
02531 }
02532
02533