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