00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "physconst.h"
00007 #include "mole.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
00017
00018
00019 enum {DEBUG_ENER=false};
00020
00021
00022
00023
00024
00025
00026 static double XVIB[H2_TOP] = { 0.70 , 0.60 , 0.20 };
00027 static double Xdust[H2_TOP] = { 0.04 , 0.10 , 0.40 };
00028
00029
00030
00031
00032 static const double energy_off = 0.273*FREQ_1EV/SPEEDLIGHT;
00033
00034 STATIC double EH2_eval( long int iVib , int ipH2 )
00035 {
00036 double EH2_here;
00037 double Evm = H2_DissocEnergies[0]* XVIB[ipH2] + energy_off;
00038
00039 double Ev = (energy_wn[0][iVib][0]+energy_off);
00040
00041
00042
00043
00044
00045 double Edust = H2_DissocEnergies[0] * Xdust[ipH2] *
00046 ( 1. - ( (Ev - Evm) / (H2_DissocEnergies[0]+energy_off-Evm)) *
00047 ( (1.-Xdust[ipH2])/2.) );
00048 ASSERT( Edust >= 0. );
00049
00050
00051
00052 EH2_here = H2_DissocEnergies[0]+energy_off - Edust;
00053 ASSERT( EH2_here >= 0.);
00054
00055 return EH2_here;
00056 }
00057
00058
00059 STATIC double H2_vib_dist( long int iVib , int ipH2 , double EH2)
00060 {
00061 double G1[H2_TOP] = { 0.3 , 0.4 , 0.9 };
00062 double G2[H2_TOP] = { 0.6 , 0.6 , 0.4 };
00063 double Evm = H2_DissocEnergies[0]* XVIB[ipH2] + energy_off;
00064 double Fv;
00065 if( (energy_wn[0][iVib][0]+energy_off) <= Evm )
00066 {
00067
00068 Fv = sexp( POW2( (energy_wn[0][iVib][0]+energy_off - Evm)/(G1[ipH2]* Evm ) ) );
00069 }
00070 else
00071 {
00072
00073 Fv = sexp( POW2( (energy_wn[0][iVib][0]+energy_off - Evm)/(G2[ipH2]*(EH2 - Evm ) ) ) );
00074 }
00075 return Fv;
00076 }
00077
00078
00079
00080
00081 void H2_Create(void)
00082 {
00083 long int i , iElecHi , iElecLo;
00084 long int iVibHi , iVibLo;
00085 long int iRotHi , iRotLo;
00086 long int iElec, iVib , iRot;
00087 long int nColl,
00088 nlines;
00089 int ier;
00090 int nEner;
00091
00092 int ipH2;
00093 realnum sum , sumj , sumv , sumo , sump;
00094
00095
00096
00097
00098 if( lgH2_READ_DATA || !h2.lgH2ON )
00099 return;
00100
00101 DEBUG_ENTRY( "H2_Create()" );
00102
00103
00104 if( mole.nH2_TRACE )
00105 fprintf(ioQQQ," H2_Create called in DEBUG mode.\n");
00106
00107
00108
00109
00110 if( h2.nElecLevelOutput < 1 )
00111 h2.nElecLevelOutput = mole.n_h2_elec_states;
00112
00113
00114 lgH2_READ_DATA = true;
00115
00116
00117
00118
00119
00120 iElecHi = 0;
00121
00122
00123 CollRateFit.reserve(h2.nVib_hi[iElecHi]+1);
00124 H2_CollRate.reserve(h2.nVib_hi[iElecHi]+1);
00125 for( iVibHi = 0; iVibHi <= h2.nVib_hi[iElecHi]; ++iVibHi )
00126 {
00127 CollRateFit.reserve(iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
00128 H2_CollRate.reserve(iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
00129 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00130 {
00131 CollRateFit.reserve(iVibHi,iRotHi,h2.nVib_hi[iElecHi]+1);
00132 H2_CollRate.reserve(iVibHi,iRotHi,h2.nVib_hi[iElecHi]+1);
00133 for( iVibLo=0; iVibLo<(h2.nVib_hi[iElecHi]+1); ++iVibLo )
00134 {
00135 CollRateFit.reserve(iVibHi,iRotHi,iVibLo,h2.nRot_hi[iElecHi][iVibLo]+1);
00136 H2_CollRate.reserve(iVibHi,iRotHi,iVibLo,h2.nRot_hi[iElecHi][iVibLo]+1);
00137 for( iRotLo=0; iRotLo<=h2.nRot_hi[iElecHi][iVibLo]; ++iRotLo )
00138 {
00139 H2_CollRate.reserve(iVibHi,iRotHi,iVibLo,iRotLo,N_X_COLLIDER);
00140 CollRateFit.reserve(iVibHi,iRotHi,iVibLo,iRotLo,N_X_COLLIDER);
00141 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
00142 {
00143
00144 CollRateFit.reserve(iVibHi,iRotHi,iVibLo,iRotLo,nColl,3);
00145 }
00146 }
00147 }
00148 }
00149 }
00150
00151 CollRateFit.alloc();
00152 H2_CollRate.alloc();
00153
00154
00155 CollRateFit.zero();
00156 H2_CollRate.zero();
00157
00158
00159 H2_populations.reserve(mole.n_h2_elec_states);
00160 pops_per_vib.reserve(mole.n_h2_elec_states);
00161 H2_dissprob.reserve(mole.n_h2_elec_states);
00162
00163 for( iElec = 0; iElec<mole.n_h2_elec_states; ++iElec )
00164 {
00165
00166 if( mole.nH2_TRACE >= mole.nH2_trace_full)
00167 fprintf(ioQQQ,"elec %li highest vib= %li\n", iElec , h2.nVib_hi[iElec] );
00168
00169 ASSERT( h2.nVib_hi[iElec] > 0 );
00170
00171
00172
00173 H2_populations.reserve(iElec,h2.nVib_hi[iElec]+1);
00174 pops_per_vib.reserve(iElec,h2.nVib_hi[iElec]+1);
00175 if( iElec > 0 )
00176 H2_dissprob.reserve(iElec,h2.nVib_hi[iElec]+1);
00177
00178
00179
00180
00181 for( iVib = 0; iVib <= h2.nVib_hi[iElec]; ++iVib )
00182 {
00183
00184 H2_populations.reserve(iElec,iVib,h2.nRot_hi[iElec][iVib]+1);
00185 if( iElec > 0 )
00186 H2_dissprob.reserve(iElec,iVib,h2.nRot_hi[iElec][iVib]+1);
00187 }
00188 }
00189
00190 H2_populations.alloc();
00191 H2_populations_LTE.alloc( H2_populations.clone() );
00192 H2_old_populations.alloc( H2_populations.clone() );
00193 H2_Boltzmann.alloc( H2_populations.clone() );
00194 H2_stat.alloc( H2_populations.clone() );
00195 energy_wn.alloc( H2_populations.clone() );
00196 H2_rad_rate_out.alloc( H2_populations.clone() );
00197 H2_lgOrtho.alloc( H2_populations.clone() );
00198
00199 pops_per_vib.alloc();
00200
00201 H2_dissprob.alloc();
00202 H2_disske.alloc( H2_dissprob.clone() );
00203
00204
00205 H2_rad_rate_out.zero();
00206
00207
00208 H2_ipPhoto.reserve(h2.nVib_hi[0]+1);
00209
00210
00211 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00212 {
00213
00214 H2_ipPhoto.reserve(iVib,h2.nRot_hi[0][iVib]+1);
00215 }
00216
00217 H2_ipPhoto.alloc();
00218 H2_col_rate_in.alloc( H2_ipPhoto.clone() );
00219 H2_col_rate_out.alloc( H2_ipPhoto.clone() );
00220 H2_rad_rate_in.alloc( H2_ipPhoto.clone() );
00221 H2_coll_dissoc_rate_coef.alloc( H2_ipPhoto.clone() );
00222 H2_coll_dissoc_rate_coef_H2.alloc( H2_ipPhoto.clone() );
00223 H2_X_colden.alloc( H2_ipPhoto.clone() );
00224 H2_X_rate_from_elec_excited.alloc( H2_ipPhoto.clone() );
00225 H2_X_rate_to_elec_excited.alloc( H2_ipPhoto.clone() );
00226 H2_X_colden_LTE.alloc( H2_ipPhoto.clone() );
00227 H2_X_formation.alloc( H2_ipPhoto.clone() );
00228 H2_X_Hmin_back.alloc( H2_ipPhoto.clone() );
00229
00230 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00231 {
00232 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
00233 {
00234
00235 H2_rad_rate_in[iVib][iRot] = -BIGFLOAT;
00236 H2_coll_dissoc_rate_coef[iVib][iRot] = -BIGFLOAT;
00237 H2_coll_dissoc_rate_coef_H2[iVib][iRot] = -BIGFLOAT;
00238 }
00239 }
00240
00241 H2_X_colden.zero();
00242 H2_X_colden_LTE.zero();
00243 H2_X_formation.zero();
00244 H2_X_Hmin_back.zero();
00245
00246 H2_X_rate_from_elec_excited.zero();
00247
00248 H2_X_rate_to_elec_excited.zero();
00249
00250
00251 H2_X_hminus_formation_distribution.reserve(nTE_HMINUS);
00252 for( i=0; i<nTE_HMINUS; ++i )
00253 {
00254 H2_X_hminus_formation_distribution.reserve(i,h2.nVib_hi[0]+1);
00255
00256 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00257 {
00258 H2_X_hminus_formation_distribution.reserve(i,iVib,h2.nRot_hi[0][iVib]+1);
00259 }
00260 }
00261 H2_X_hminus_formation_distribution.alloc();
00262 H2_X_hminus_formation_distribution.zero();
00263 H2_Read_hminus_distribution();
00264
00265
00266
00267
00268
00269
00270
00271 H2_X_grain_formation_distribution.reserve(H2_TOP);
00272 for( ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
00273 {
00274 H2_X_grain_formation_distribution.reserve(ipH2,h2.nVib_hi[0]+1);
00275
00276
00277 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00278 {
00279 H2_X_grain_formation_distribution.reserve(ipH2,iVib,h2.nRot_hi[0][iVib]+1);
00280 }
00281 }
00282 H2_X_grain_formation_distribution.alloc();
00283 H2_X_grain_formation_distribution.zero();
00284
00285
00286
00287 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00288 {
00289
00290 H2_ReadEnergies(iElec);
00291
00292
00293 if( iElec > 0 )
00294 H2_ReadDissprob(iElec);
00295 }
00296
00297
00298
00299
00300 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00301 {
00302 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
00303 {
00304
00305
00306
00307
00308
00309
00310
00311
00313
00314 double thresh = (H2_DissocEnergies[1] - energy_wn[0][iVib][iRot])*WAVNRYD;
00315
00316
00317
00318 ASSERT( thresh > 0.74 );
00319 H2_ipPhoto[iVib][iRot] = ipoint(thresh);
00320 }
00321 }
00322
00323 nH2_energies = 0;
00324 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec)
00325 {
00326
00327 nH2_energies += nLevels_per_elec[iElec];
00328 }
00329
00330 if( mole.nH2_TRACE >= mole.nH2_trace_full )
00331 {
00332 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec)
00333 {
00334
00335 fprintf(ioQQQ,"\t(%li %li)", iElec , nLevels_per_elec[iElec] );
00336 }
00337 fprintf(ioQQQ,
00338 " H2_Create: there are %li electronic levels, in each level there are",
00339 mole.n_h2_elec_states);
00340 fprintf(ioQQQ,
00341 " for a total of %li levels.\n", nH2_energies );
00342 }
00343
00344
00345 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
00346 {
00347
00348 H2_CollidRateRead(nColl);
00349 }
00350
00351 if( mole.lgH2_NOISE )
00352 {
00353 iElecHi = 0;
00354
00355 for( iVibHi = 0; iVibHi <= VIB_COLLID; ++iVibHi )
00356 {
00357 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00358 {
00359 for( iVibLo=0; iVibLo<(VIB_COLLID+1); ++iVibLo )
00360 {
00361 for( iRotLo=0; iRotLo<=h2.nRot_hi[iElecHi][iVibLo]; ++iRotLo )
00362 {
00363
00364
00365
00366
00367
00368 for( nColl=0; nColl<N_X_COLLIDER-2; ++nColl )
00369 {
00370
00371
00372
00373 if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] != 0. )
00374 {
00375
00376 realnum r = (realnum)RandGauss( mole.xMeanNoise , mole.xSTDNoise );
00377
00378
00379 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] += r;
00380 }
00381 }
00382
00383
00384 if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][N_X_COLLIDER-2][0] != 0. )
00385 {
00386
00387 realnum r = (realnum)RandGauss( mole.xMeanNoise , mole.xSTDNoise );
00388
00389
00390 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][N_X_COLLIDER-2][0] *= pow((realnum)10.f,r);
00391 }
00392 }
00393 }
00394 }
00395 }
00396 }
00397
00398
00399
00400 H2_energies = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nH2_energies );
00401 H2_ipX_ener_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
00402 ipElec_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
00403 ipVib_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
00404 ipRot_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies );
00405
00406
00407 H2_X_source = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nLevels_per_elec[0] );
00408 H2_X_sink = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nLevels_per_elec[0] );
00409
00410 H2_X_coll_rate.reserve(nLevels_per_elec[0]);
00411
00412 for( i=1; i<nLevels_per_elec[0]; ++i )
00413 {
00414 H2_X_coll_rate.reserve(i,i);
00415 }
00416 H2_X_coll_rate.alloc();
00417
00418
00419 ipEnergySort.reserve(mole.n_h2_elec_states);
00420 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00421 {
00422 ipEnergySort.reserve(iElecHi,h2.nVib_hi[iElecHi]+1);
00423 for( iVibHi = 0; iVibHi <= h2.nVib_hi[iElecHi]; ++iVibHi )
00424 {
00425 ipEnergySort.reserve(iElecHi,iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
00426 }
00427 }
00428 ipEnergySort.alloc();
00429
00430 nEner = 0;
00431 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00432 {
00433
00434 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00435 {
00436 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00437 {
00438 H2_energies[nEner] = (realnum)energy_wn[iElecHi][iVibHi][iRotHi];
00439 ipElec_H2_energy_sort[nEner] = iElecHi;
00440 ipVib_H2_energy_sort[nEner] = iVibHi;
00441 ipRot_H2_energy_sort[nEner] = iRotHi;
00442 ipEnergySort[iElecHi][iVibHi][iRotHi] = -1;
00443 ++nEner;
00444 }
00445 }
00446 }
00447
00448 ASSERT( nH2_energies == nEner );
00449
00450
00451
00452 spsort(
00453
00454 H2_energies,
00455
00456 nH2_energies,
00457
00458 H2_ipX_ener_sort,
00459
00460
00461 1,
00462
00463 &ier);
00464
00465
00466 for( nEner=0; nEner<nH2_energies; ++nEner )
00467 {
00468 if( nEner+1 < nLevels_per_elec[0] )
00469 ASSERT( H2_energies[H2_ipX_ener_sort[nEner]] <
00470 H2_energies[H2_ipX_ener_sort[nEner+1]] );
00471
00472
00473
00474
00475
00476 i = H2_ipX_ener_sort[nEner];
00477 iElec = ipElec_H2_energy_sort[i];
00478 iRot = ipRot_H2_energy_sort[i];
00479 iVib = ipVib_H2_energy_sort[i];
00480
00481 ipEnergySort[iElec][iVib][iRot] = nEner;
00482 }
00483
00484
00485 for( nEner=0; nEner<nLevels_per_elec[0]; ++nEner )
00486 {
00487 i = H2_ipX_ener_sort[nEner];
00488 iRot = ipRot_H2_energy_sort[i];
00489 iVib = ipVib_H2_energy_sort[i];
00490 if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR )
00491 break;
00492 nEner_H2_ground = nEner;
00493 }
00494
00495
00496 ++nEner_H2_ground;
00497
00498
00499
00500
00501
00502 if( nXLevelsMatrix<0 )
00503 {
00504 nXLevelsMatrix = nLevels_per_elec[0];
00505 }
00506 else if( nXLevelsMatrix > nLevels_per_elec[0] )
00507 {
00508 fprintf( ioQQQ,
00509 " 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",
00510 nXLevelsMatrix ,
00511 nLevels_per_elec[0]);
00512 cdEXIT(EXIT_FAILURE);
00513 }
00514
00515
00516 H2Lines.reserve(mole.n_h2_elec_states);
00517
00518 nlines = 0;
00519 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00520 {
00521 H2Lines.reserve(iElecHi,h2.nVib_hi[iElecHi]+1);
00522 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00523 {
00524 H2Lines.reserve(iElecHi,iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1);
00525 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00526 {
00527
00528
00529
00530
00531 long int lim_elec_lo = 0;
00532 H2Lines.reserve(iElecHi,iVibHi,iRotHi,1);
00533 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00534 {
00535
00536
00537
00538 long int nv = h2.nVib_hi[iElecLo];
00539
00540 if( iElecLo==iElecHi )
00541 nv = iVibHi;
00542 H2Lines.reserve(iElecHi,iVibHi,iRotHi,iElecLo,nv+1);
00543 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00544 {
00545 long nr = h2.nRot_hi[iElecLo][iVibLo];
00546 if( iElecLo==iElecHi && iVibHi==iVibLo )
00547
00548 nr = MAX2(1,iRotHi-1);
00549 H2Lines.reserve(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,nr+1);
00550 nlines += nr+1;
00551 }
00552 }
00553 }
00554 }
00555 }
00556
00557 H2Lines.alloc();
00558 H2_SaveLine.alloc( H2Lines.clone() );
00559 lgH2_line_exists.alloc( H2Lines.clone() );
00560
00561
00562 H2_SaveLine.zero();
00563
00564 lgH2_line_exists.zero();
00565
00566 if( mole.nH2_TRACE >= mole.nH2_trace_full )
00567 fprintf(ioQQQ," There are a total of %li lines in the entire H2 molecule.\n", nlines );
00568
00569
00570 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00571 {
00572 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00573 {
00574 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00575 {
00576
00577
00578
00579 long int lim_elec_lo = 0;
00580 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00581 {
00582
00583
00584 long int nv = h2.nVib_hi[iElecLo];
00585 if( iElecLo==iElecHi )
00586 nv = iVibHi;
00587 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00588 {
00589 long nr = h2.nRot_hi[iElecLo][iVibLo];
00590 if( iElecLo==iElecHi && iVibHi==iVibLo )
00591 nr = iRotHi-1;
00592 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00593 {
00594 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Junk();
00595 }
00596 }
00597 }
00598 }
00599 }
00600 }
00601
00602
00603 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00604 {
00605 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00606 {
00607 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00608 {
00609 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi = AddState2Stack();
00610
00611
00612
00613
00614 long int lim_elec_lo = 0;
00615 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00616 {
00617
00618
00619 long int nv = h2.nVib_hi[iElecLo];
00620 if( iElecLo==iElecHi )
00621 nv = iVibHi;
00622 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00623 {
00624 long nr = h2.nRot_hi[iElecLo][iVibLo];
00625 if( iElecLo==iElecHi && iVibHi==iVibLo )
00626 nr = iRotHi-1;
00627 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00628 {
00629 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi =
00630 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi;
00631
00632 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo =
00633 H2Lines[iElecLo][iVibLo][iRotLo][0][0][0].Hi;
00634
00635
00636 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Zero();
00637 }
00638 }
00639 }
00640 }
00641 }
00642 }
00643
00644
00645 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
00646 {
00647
00648 H2_ReadTransprob(iElec);
00649 }
00650
00651
00652
00653 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00654 {
00655 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00656 {
00657 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00658 {
00659
00660
00661
00662 if( is_odd(iRotHi+H2_nRot_add_ortho_para[iElecHi]) )
00663 {
00664
00665 H2_lgOrtho[iElecHi][iVibHi][iRotHi] = true;
00666 H2_stat[iElecHi][iVibHi][iRotHi] = 3.f*(2.f*iRotHi+1.f);
00667 }
00668 else
00669 {
00670
00671 H2_lgOrtho[iElecHi][iVibHi][iRotHi] = false;
00672 H2_stat[iElecHi][iVibHi][iRotHi] = (2.f*iRotHi+1.f);
00673 }
00674 }
00675 }
00676 }
00677
00678
00679 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi )
00680 {
00681 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi )
00682 {
00683 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00684 {
00685
00686
00687
00688
00689 long int lim_elec_lo = 0;
00690 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
00691 {
00692
00693
00694 long int nv = h2.nVib_hi[iElecLo];
00695 if( iElecLo==iElecHi )
00696 nv = iVibHi;
00697 for( iVibLo=0; iVibLo<=nv; ++iVibLo )
00698 {
00699 long nr = h2.nRot_hi[iElecLo][iVibLo];
00700
00701 if( iElecLo==iElecHi && iVibHi==iVibLo )
00702 nr = iRotHi-1;
00703 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo )
00704 {
00705
00706
00707 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->nelem = LIMELM+3;
00708
00709 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->IonStg = 1;
00710
00711
00712 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g = H2_stat[iElecLo][iVibLo][iRotLo];
00713 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g = H2_stat[iElecHi][iVibHi][iRotHi];
00714
00715
00716 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN =
00717 (realnum)(energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]);
00718
00719
00720 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN > SMALLFLOAT)
00721 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng =
00722 (realnum)(1.e8f/H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN /
00723 RefIndex( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN));
00724
00725 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK =
00726 (realnum)(T1CM)*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN;
00727
00728
00729 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg =
00730 (realnum)(ERG1CM)*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN;
00731
00732
00733 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] )
00734 {
00735
00736
00737
00738 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->iRedisFun = ipCRDW;
00739
00740
00741 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauIn = opac.taumin;
00742 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauCon = opac.taumin;
00743
00744 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauTot = 1e20f;
00745
00746
00747 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->dampXvel =
00748 (realnum)(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul/
00749 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN/PI4);
00750
00751 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->gf =
00752 (realnum)(GetGF(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul,
00753 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN,
00754 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g ) );
00755
00756
00757 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->opacity = (realnum)(
00758 abscf(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->gf,
00759 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN,
00760 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g));
00761
00762 if( iElecHi > 0 )
00763 {
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.col_str = (realnum)(
00775 pow3(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng*1e-8)*
00776 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g/
00777 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g)*
00778 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul*
00779 log(3.)*HPLANCK/(160.f*pow3(PI)*0.5*1e-8*EN1EV)/6.e-17);
00780 ASSERT(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.col_str>0.);
00781 }
00782 else
00783 {
00784
00785 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.col_str = 0.;
00786 }
00787 }
00788 else
00789 {
00790
00791
00792 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.col_str = 0.;
00793 }
00794 }
00795 }
00796 }
00797 }
00798 }
00799 }
00800
00801
00802
00803
00804
00805
00806
00807
00808 for( ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
00809 {
00810 sum = 0.;
00811 sumj = 0.;
00812 sumv = 0.;
00813 sumo = 0.;
00814 sump = 0.;
00815 iElec = 0;
00816
00817 if( hmi.chGrainFormPump == 'D' )
00818 {
00819
00820
00821
00822 double T_H2_FORM = 50000.;
00823 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00824 {
00825 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
00826 {
00827
00828 H2_X_grain_formation_distribution[ipH2][iVib][iRot] =
00829
00830 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * (1.f+iVib) *
00831 (realnum)sexp( energy_wn[iElec][iVib][iRot]*T1CM/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 if( hmi.chGrainFormPump == 'T' )
00848 {
00849
00850 double Xrot[H2_TOP] = { 0.14 , 0.15 , 0.15 };
00851 double Xtrans[H2_TOP] = { 0.12 , 0.15 , 0.25 };
00852
00853 double sumvib = 0.;
00854 double EH2;
00855
00856 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00857 {
00858 double vibdist;
00859 EH2 = EH2_eval( iVib , ipH2 );
00860 vibdist = H2_vib_dist( iVib , ipH2 , EH2);
00861 sumvib += vibdist;
00862 }
00863
00864
00865 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00866 {
00867 double Ev = (energy_wn[iElec][iVib][0]+energy_off);
00868 double Fv;
00869
00870
00871 double Erot;
00872
00873
00874 EH2 = EH2_eval( iVib , ipH2 );
00875
00876
00877 Erot = (EH2 - Ev) * Xrot[ipH2] / (Xrot[ipH2] + Xtrans[ipH2]);
00878
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902 if( Erot > 0. )
00903 {
00904
00905 Fv = H2_vib_dist( iVib , ipH2 , EH2) / sumvib;
00906
00907
00908 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
00909 {
00910
00911 double gaussian =
00912 sexp( POW2( (energy_wn[iElec][iVib][iRot] - energy_wn[iElec][iVib][0] - Erot) /
00913 (0.5 * Erot) ) );
00914
00915 double thermal_dist =
00916 sexp( (energy_wn[iElec][iVib][iRot] - energy_wn[iElec][iVib][0]) /
00917 Erot );
00918
00919
00920 double aver = ( gaussian + thermal_dist ) / 2.;
00921
00922
00923
00924
00925
00926
00927 ASSERT( gaussian <= 1. );
00928
00929 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = (realnum)(
00930
00931 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * Fv * (2.*iRot+1.) * aver );
00932
00933 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00934 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00935 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00936 if( H2_lgOrtho[iElec][iVib][iRot] )
00937 {
00938 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00939 }
00940 else
00941 {
00942 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00943 }
00944
00945 }
00946 }
00947 else
00948 {
00949
00950 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
00951 {
00952 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 0.;
00953 }
00954 }
00955 }
00956 }
00957 else if( hmi.chGrainFormPump == 't' )
00958 {
00959
00960
00961
00962
00963
00964 double T_H2_FORM = 17329.;
00965 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00966 {
00967 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
00968 {
00969
00970 H2_X_grain_formation_distribution[ipH2][iVib][iRot] =
00971
00972 H2_stat[0][iVib][iRot] *
00973 (realnum)sexp( energy_wn[0][iVib][iRot]*T1CM/T_H2_FORM );
00974 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00975 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00976 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00977 if( H2_lgOrtho[iElec][iVib][iRot] )
00978 {
00979 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00980 }
00981 else
00982 {
00983
00984 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
00985 }
00986 }
00987 }
00988 }
00989 else
00990 TotalInsanity();
00991
00992 if( mole.nH2_TRACE >= mole.nH2_trace_full )
00993 fprintf(ioQQQ, "H2 form grains mean J= %.3f mean v = %.3f ortho/para= %.3f\n",
00994 sumj/sum , sumv/sum , sumo/sump );
00995
00996 iElec = 0;
00997
00998 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib )
00999 {
01000 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot )
01001 {
01002 H2_X_grain_formation_distribution[ipH2][iVib][iRot] /= sum;
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012 }
01013 }
01014 }
01015
01016
01017
01018 {
01019
01020 if( DEBUG_ENER )
01021 {
01022
01023
01024 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec )
01025 {
01026
01027 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib )
01028 {
01029 for( iRot=0; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot )
01030 {
01031 fprintf(ioQQQ,"%li\t%li\t%li\t%.5e\n",
01032 iElec, iVib, iRot ,
01033 energy_wn[iElec][iVib][iRot]);
01034 }
01035 }
01036 }
01037
01038 cdEXIT(EXIT_SUCCESS);
01039 }
01040 }
01041
01042 return;
01043 }