00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "taulines.h"
00007 #include "iso.h"
00008 #include "conv.h"
00009 #include "trace.h"
00010 #include "secondaries.h"
00011 #include "hmi.h"
00012 #include "h2.h"
00013 #include "atmdat.h"
00014 #include "dense.h"
00015 #include "ionbal.h"
00016 #include "hydrogenic.h"
00017
00018 void Hydrogenic(void)
00019 {
00020 long int ipLo ,
00021 ipHi ,
00022 nelem,
00023 ipISO=ipH_LIKE,
00024 mol;
00025 double coltot,
00026 gamtot;
00027 double sum , error;
00028 static bool lgFinitePop[LIMELM];
00029 static bool lgMustInit=true;
00030 bool lgH_chem_conv;
00031 int loop_H_chem;
00032 double solomon_assumed;
00033 double *PumpSave=NULL;
00034
00035 DEBUG_ENTRY( "Hydrogenic()" );
00036
00037 if( lgMustInit )
00038 {
00039 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
00040 lgFinitePop[nelem] = true;
00041 lgMustInit = false;
00042 }
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 if( hydro.xLymanPumpingScaleFactor!=1.f )
00053 {
00054
00055 PumpSave = (double *)MALLOC( (unsigned)iso.numLevels_local[ipH_LIKE][ipHYDROGEN]*sizeof(double) );
00056 ipLo = 0;
00057
00058
00059
00060 for( ipHi=2; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi )
00061 {
00062 PumpSave[ipHi] = EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].pump;
00063 EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].pump *= hydro.xLymanPumpingScaleFactor;
00064 }
00065 }
00066
00067 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00068 {
00069
00070 if( dense.lgElmtOn[nelem] )
00071 {
00072
00073
00074 if( (dense.IonHigh[nelem] == nelem + 1) )
00075 {
00076
00077
00078
00079 iso_continuum_lower( ipH_LIKE, nelem );
00080
00081
00082 HydroCollid(nelem);
00083
00084
00085 iso_photo( ipH_LIKE , nelem );
00086
00087
00088 HydroRecom(nelem);
00089
00090
00091
00092 iso_ionize_recombine( ipH_LIKE , nelem );
00093
00094
00095 HydroLevel(nelem);
00096
00097
00098 lgFinitePop[nelem] = true;
00099
00100 if( trace.lgTrace )
00101 {
00102 fprintf( ioQQQ, " Hydrogenic Z=%2ld H2OVH1=",nelem);
00103 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.pop_ion_ov_neut[ipH_LIKE][nelem]));
00104 fprintf( ioQQQ, " simple=");
00105 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.pop_ion_ov_neut[ipH_LIKE][nelem]));
00106
00107
00108 fprintf( ioQQQ, "\n");
00109 }
00110
00111 }
00112 else if( lgFinitePop[nelem] )
00113 {
00114
00115
00116 lgFinitePop[nelem] = false;
00117
00118 iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 0.;
00119
00120
00121 for( ipHi=ipH2s; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00122 {
00123 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00124 {
00125
00126 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo = 0.;
00127
00128
00129 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi = 0.;
00130
00131
00132 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc = 0.;
00133 }
00134 }
00135 }
00136
00137 if( dense.lgSetIoniz[nelem] )
00138 {
00139 dense.xIonDense[nelem][nelem+1-ipISO] = dense.SetIoniz[nelem][nelem+1-ipISO]*dense.gas_phase[nelem];
00140 dense.xIonDense[nelem][nelem-ipISO] = dense.SetIoniz[nelem][nelem-ipISO]*dense.gas_phase[nelem];
00141 if( dense.SetIoniz[nelem][nelem+1-ipISO] > SMALLFLOAT )
00142 {
00143 iso.Pop2Ion[ipISO][nelem][ipH1s] = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO];
00144 EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO];
00145 }
00146 else
00147 {
00148 iso.Pop2Ion[ipISO][nelem][ipH1s] = 0.;
00149 EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo = 0.;
00150 }
00151 }
00152 }
00153 }
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166 lgH_chem_conv = false;
00167 loop_H_chem = 0;
00168 while( loop_H_chem < 5 && !lgH_chem_conv )
00169 {
00170
00171
00172 solomon_assumed = hmi.H2_Solomon_dissoc_rate_used_H2g/SDIV(hmi.H2_rate_destroy);
00173
00174
00175 hmole();
00176
00177 H2_LevelPops();
00178
00179
00180 lgH_chem_conv = true;
00181
00182 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
00183 {
00184
00185 if( fabs( solomon_assumed - hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.H2_rate_destroy) ) >
00186 conv.EdenErrorAllowed/5.)
00187 {
00188 lgH_chem_conv = false;
00189 }
00190 }
00191 ++loop_H_chem;
00192 }
00193
00194 {
00195
00196
00197
00198
00199 enum {DEBUG_LOC=false};
00200
00201 if(DEBUG_LOC )
00202 {
00203 fprintf(ioQQQ,"DEBUG\t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
00204 fnzone,
00205 hmi.H2_total ,
00206 dense.xIonDense[ipHYDROGEN][0],
00207 dense.xIonDense[ipHYDROGEN][1],
00208 hmi.H2_Solomon_dissoc_rate_used_H2g,
00209 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
00210 hmi.H2_Solomon_dissoc_rate_TH85_H2g);
00211 }
00212 }
00213
00214 if( dense.lgSetIoniz[ipHYDROGEN] )
00215 {
00216 dense.xIonDense[ipHYDROGEN][1] = dense.SetIoniz[ipHYDROGEN][1]*dense.gas_phase[ipHYDROGEN];
00217 dense.xIonDense[ipHYDROGEN][0] = dense.SetIoniz[ipHYDROGEN][0]*dense.gas_phase[ipHYDROGEN];
00218
00219
00220 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] = dense.SetIoniz[ipHYDROGEN][0] / SDIV( dense.SetIoniz[ipHYDROGEN][1]);
00221 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopLo = dense.SetIoniz[ipHYDROGEN][0] / SDIV( dense.SetIoniz[ipHYDROGEN][1]);
00222 }
00223 else
00224 {
00225
00226
00227
00228
00229
00230
00231 ion_solver( ipHYDROGEN , false );
00232 }
00233
00234
00235
00236
00237
00238
00239 sum = dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1];
00240 for(mol=0;mol<N_H_MOLEC;mol++)
00241 {
00242 sum += hmi.Hmolec[mol]*hmi.nProton[mol];
00243 }
00244
00245 sum -= hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp];
00246
00247
00248
00249
00250
00251 error = ( dense.gas_phase[ipHYDROGEN] - sum )/dense.gas_phase[ipHYDROGEN];
00252
00253 {
00254
00255
00256
00257
00258 enum {DEBUG_LOC=false};
00259
00260 if(DEBUG_LOC && (fabs(error) > 1e-4) )
00261 fprintf(ioQQQ,"PROBLEM hydrogenic zone %li hden %.4e, sum %.4e (h-s)/h %.3e \n", nzone, dense.gas_phase[ipHYDROGEN] , sum ,
00262 error );
00263 }
00264
00265
00266
00267 HydroRenorm();
00268
00269
00270
00271 if( iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 )
00272 {
00273 coltot =
00274 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL*
00275 iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH1s]*
00276 dense.eden*(iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH2p] + iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH2p] +
00277 EmisLines[ipH_LIKE][ipHYDROGEN][3][ipH2p].ColUL*iso.PopLTE[ipH_LIKE][ipHYDROGEN][3]/iso.PopLTE[ipH_LIKE][ipHYDROGEN][ipH2p])/
00278 (EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL*dense.eden +
00279 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul*
00280 (EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pesc +
00281 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pelec_esc +
00282 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest) );
00283
00284
00285
00286 if( coltot > 0.2 )
00287 {
00288 hydro.lgHColionImp = true;
00289 }
00290 }
00291 else
00292 {
00293 hydro.lgHColionImp = false;
00294 }
00295
00296
00297
00298
00299
00300 if( iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]/MAX2(SMALLDOUBLE,iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]) > 0.1 &&
00301 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] > SMALLDOUBLE )
00302 {
00303 hydro.lgHiPop2 = true;
00304 hydro.pop2mx = (float)MAX2(iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH2p]/iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s],
00305 hydro.pop2mx);
00306 }
00307
00308 gamtot = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s] + secondaries.csupra[ipHYDROGEN][0];
00309
00310 coltot = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s] +
00311 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ColUL*4.*
00312 iso.Boltzmann[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s];
00313
00314
00315 if( iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] > SMALLFLOAT )
00316 {
00317 hydro.H_ion_frac_photo =
00318 (float)(iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] );
00319
00320
00321 hydro.H_ion_frac_collis =
00322 (float)(iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.eden/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]);
00323
00324
00325
00326 secondaries.sec2total =
00327 (float)(secondaries.csupra[ipHYDROGEN][0] / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]);
00328
00329
00330 atmdat.HIonFrac = atmdat.HCharExcIonTotal / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s];
00331 }
00332 else
00333 {
00334 hydro.H_ion_frac_collis = 0.;
00335 hydro.H_ion_frac_photo = 0.;
00336 secondaries.sec2total = 0.;
00337 atmdat.HIonFrac = 0.;
00338 }
00339
00340 if( trace.lgTrace )
00341 {
00342 fprintf( ioQQQ, " Hydrogenic retrn %.2f ",fnzone);
00343 fprintf(ioQQQ,"H0:%.3e ", dense.xIonDense[ipHYDROGEN][0]);
00344 fprintf(ioQQQ,"H+:%.3e ", dense.xIonDense[ipHYDROGEN][1]);
00345 fprintf(ioQQQ,"H2:%.3e ", hmi.H2_total);
00346 fprintf(ioQQQ,"H-:%.3e ", hmi.Hmolec[ipMHm]);
00347 fprintf(ioQQQ,"ne:%.3e ", dense.eden);
00348 fprintf( ioQQQ, " REC, COL, GAMT= ");
00349
00350 fprintf(ioQQQ,"%.2e ", iso.RadRec_effec[ipH_LIKE][ipHYDROGEN] );
00351 fprintf(ioQQQ,"%.2e ", coltot);
00352 fprintf(ioQQQ,"%.2e ", gamtot);
00353 fprintf( ioQQQ, " CSUP=");
00354 PrintE82( ioQQQ, secondaries.csupra[ipHYDROGEN][0]);
00355 fprintf( ioQQQ, "\n");
00356 }
00357
00358 if( hydro.xLymanPumpingScaleFactor!=1.f )
00359 {
00360 ipLo = 0;
00361
00362
00363
00364 for( ipHi=2; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi )
00365 {
00366 EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].pump = PumpSave[ipHi];
00367 }
00368 free(PumpSave);
00369 }
00370
00371 DEBUG_EXIT( "Hydrogenic()" );
00372 return;
00373 }
00374
00375
00376 void HydroRenorm( void )
00377 {
00378 double sum_atom_iso , renorm;
00379 long int level,
00380 nelem,
00381 ipHi ,
00382 ipLo;
00383
00384 DEBUG_ENTRY( "HydroRenorm()" );
00385
00386
00387
00388
00389
00390 sum_atom_iso = 0.;
00391
00392 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; level++ )
00393 {
00394 sum_atom_iso += iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][level];
00395 }
00396
00397 sum_atom_iso *= dense.xIonDense[ipHYDROGEN][1];
00398
00399
00400 if( sum_atom_iso > SMALLFLOAT )
00401 {
00402 renorm = dense.xIonDense[ipHYDROGEN][0] / sum_atom_iso;
00403 }
00404 else
00405 {
00406 renorm = 0.;
00407 }
00408 ASSERT( renorm < BIGFLOAT );
00409
00410
00411 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] *= renorm;
00412
00413
00414
00415
00416 nelem = ipHYDROGEN;
00417
00418 for( ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
00419 {
00420 iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipHi] *= renorm;
00421
00422 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00423 {
00424
00425 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo *= renorm;
00426
00427
00428 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi *= renorm;
00429
00430
00431 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc *= renorm;
00432
00433
00434 }
00435 }
00436
00437 DEBUG_EXIT( "HydroRenorm()" );
00438
00439 return;
00440 }