00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "taulines.h"
00008 #include "iso.h"
00009 #include "dense.h"
00010 #include "secondaries.h"
00011 #include "trace.h"
00012 #include "elementnames.h"
00013 #include "conv.h"
00014 #include "phycon.h"
00015 #include "dynamics.h"
00016 #include "opacity.h"
00017 #include "ionbal.h"
00018 #include "continuum.h"
00019 #include "hydrogenic.h"
00020
00021
00022
00023
00024 static void PrtHydroTrace1a(long nelem )
00025 {
00026 double colfrc,
00027 phtfrc,
00028 secfrc,
00029 collider;
00030
00031
00032
00033
00034
00035
00036 if( nelem==ipHYDROGEN )
00037 {
00038
00039 collider = dense.EdenHontoHCorr;
00040 }
00041 else
00042 {
00043 collider = dense.EdenHCorr;
00044 }
00045
00046
00047 if( iso.xIonSimple[ipH_LIKE][nelem] > 0. )
00048 {
00049
00050 phtfrc = iso.gamnc[ipH_LIKE][nelem][ipH1s]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] +
00051 ionbal.CotaRate[nelem]) )*
00052 iso.xIonSimple[ipH_LIKE][nelem]);
00053
00054
00055
00056 colfrc = (iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*collider )/
00057 ((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] +
00058 ionbal.CotaRate[0]) )*
00059 iso.xIonSimple[ipH_LIKE][nelem]);
00060
00061
00062 secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso.RadRec_effec[ipH_LIKE][nelem] +
00063 ionbal.CotaRate[0]) )*
00064 iso.xIonSimple[ipH_LIKE][nelem]);
00065 }
00066 else
00067 {
00068 phtfrc = 0.;
00069 colfrc = 0.;
00070 secfrc = 0.;
00071 }
00072
00073 fprintf( ioQQQ, " HydroLevel Z=%2ld called, simple II/I=",nelem);
00074 PrintE93( ioQQQ, iso.xIonSimple[ipH_LIKE][nelem]);
00075 fprintf( ioQQQ," PhotFrc:");
00076 PrintE93( ioQQQ,phtfrc);
00077 fprintf(ioQQQ," ColFrc:");
00078 PrintE93( ioQQQ,colfrc);
00079 fprintf( ioQQQ," SecFrc");
00080 PrintE93( ioQQQ, secfrc);
00081 fprintf( ioQQQ," Te:");
00082 PrintE93( ioQQQ,phycon.te);
00083 fprintf( ioQQQ," eden:");
00084 PrintE93( ioQQQ,dense.eden);
00085 fprintf( ioQQQ,"\n");
00086 return;
00087 }
00088
00089
00090 static void PrtHydroTrace1(long nelem )
00091 {
00092 long int ipHi , ipLo , i;
00093
00094 fprintf( ioQQQ,
00095 " HydroLevel%3ld finds arrays, with optical depths defined? %1c induced 2ph=%12.3e\n",
00096 nelem, TorF(opac.lgTauOutOn), EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].pump );
00097
00098 for( ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
00099 {
00100 fprintf( ioQQQ, "up:%2ld", ipHi );
00101 fprintf( ioQQQ, "lo" );
00102 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00103 {
00104 fprintf( ioQQQ, "%9ld", ipLo );
00105 }
00106 fprintf( ioQQQ, "\n" );
00107
00108 fprintf( ioQQQ, "%3ld", ipHi );
00109 fprintf( ioQQQ, " A*esc" );
00110 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00111 {
00112 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
00113 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc ));
00114 }
00115 fprintf( ioQQQ, "\n" );
00116
00117 fprintf( ioQQQ, "%3ld", ipHi );
00118 fprintf( ioQQQ, " A*ees" );
00119 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00120 {
00121 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
00122 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pelec_esc ));
00123 }
00124 fprintf( ioQQQ, "\n" );
00125
00126 fprintf( ioQQQ, "%3ld", ipHi );
00127 fprintf( ioQQQ, " tauin" );
00128 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00129 {
00130 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn ));
00131 }
00132 fprintf( ioQQQ, "\n" );
00133
00134 fprintf( ioQQQ, "%3ld", ipHi );
00135 fprintf( ioQQQ, " t tot" );
00136 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00137 {
00138 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauTot ));
00139 }
00140 fprintf( ioQQQ, "\n" );
00141
00142 fprintf( ioQQQ, "%3ld", ipHi );
00143 fprintf( ioQQQ, " Esc " );
00144 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00145 {
00146 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pesc ));
00147 }
00148 fprintf( ioQQQ, "\n" );
00149
00150 fprintf( ioQQQ, "%3ld", ipHi );
00151 fprintf( ioQQQ, " Eesc " );
00152 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00153 {
00154 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pelec_esc ));
00155 }
00156 fprintf( ioQQQ, "\n" );
00157
00158 fprintf( ioQQQ, "%3ld", ipHi );
00159 fprintf( ioQQQ, " Dest " );
00160 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00161 {
00162 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pdest) );
00163 }
00164 fprintf( ioQQQ, "\n" );
00165
00166 fprintf( ioQQQ, "%3ld", ipHi );
00167 fprintf( ioQQQ, " A*dst" );
00168 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00169 {
00170 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
00171 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Pdest ));
00172 }
00173 fprintf( ioQQQ, "\n" );
00174
00175 fprintf( ioQQQ, "%3ld", ipHi );
00176 fprintf( ioQQQ, " StrkE" );
00177 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00178 {
00179 fprintf( ioQQQ,PrintEfmt("%9.2e", hydro.pestrk[ipLo][ipHi] ));
00180 }
00181 fprintf( ioQQQ, "\n" );
00182
00183 fprintf( ioQQQ, "%3ld", ipHi );
00184 fprintf( ioQQQ, " B(ul)" );
00185 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00186 {
00187 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].pump*
00188 iso.stat[ipH_LIKE][nelem][ipLo]/iso.stat[ipH_LIKE][nelem][ipHi] ));
00189 }
00190 fprintf( ioQQQ, "\n" );
00191
00192 fprintf( ioQQQ, "%3ld", ipHi );
00193 fprintf( ioQQQ, " tcont" );
00194 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00195 {
00196 fprintf( ioQQQ,PrintEfmt("%9.2e", EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauCon ));
00197 }
00198 fprintf( ioQQQ, "\n" );
00199
00200 fprintf( ioQQQ, "%3ld", ipHi );
00201 fprintf( ioQQQ, " C(ul)" );
00202 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00203 {
00204 fprintf( ioQQQ,PrintEfmt("%9.2e",
00205 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ColUL*dense.eden ));
00206 }
00207 fprintf( ioQQQ, "\n" );
00208
00209 if( ipHi == 2 )
00210 {
00211 fprintf( ioQQQ, " FeIIo");
00212 fprintf( ioQQQ,PrintEfmt("%9.2e",
00213 hydro.dstfe2lya* EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul ));
00214 fprintf( ioQQQ, "\n");
00215 }
00216 }
00217
00218 fprintf( ioQQQ, " " );
00219
00220 for( i=1; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00221 {
00222 fprintf( ioQQQ, "%9ld", i );
00223 }
00224 fprintf( ioQQQ, "\n" );
00225 return;
00226 }
00227
00228
00229
00230
00231 void HydroLevel(long int nelem)
00232 {
00233 long int i,
00234 ipHi,
00235 ipLo,
00236 level;
00237 double
00238 sum_popn_ov_ion,
00239 collider;
00240
00241 double TooSmall;
00242
00243 DEBUG_ENTRY( "HydroLevel()" );
00244
00245
00246 ASSERT( nelem >= 0);
00247 ASSERT( nelem < LIMELM );
00248
00249
00250
00251
00252
00253 if( nelem==ipHYDROGEN )
00254 {
00255
00256 collider = dense.EdenHontoHCorr;
00257 }
00258 else
00259 {
00260 collider = dense.EdenHCorr;
00261 }
00262
00263
00264
00265 ASSERT( ionbal.CollIonRate_Ground[nelem][nelem][0]< SMALLFLOAT ||
00266 fabs( iso.ColIoniz[ipH_LIKE][nelem][0]* collider /
00267 SDIV(ionbal.CollIonRate_Ground[nelem][nelem][0] ) - 1.) < 0.001 );
00268
00269
00270 if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
00271 PrtHydroTrace1(nelem);
00272
00273 if( trace.lgHBug && trace.lgTrace )
00274 PrtHydroTrace1a(nelem);
00275
00276
00277
00278 if( nelem==ipHYDROGEN )
00279 TooSmall = 1e-28;
00280 else if( nelem==ipHELIUM )
00281 TooSmall = 1e-18;
00282 else
00283 TooSmall = 2e-18;
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293 if( iso.xIonSimple[ipH_LIKE][nelem] < 1e-35 )
00294 {
00295
00296 strcpy( iso.chTypeAtomUsed[ipH_LIKE][nelem], "zero " );
00297 if( trace.lgHBug && trace.lgTrace )
00298 {
00299 fprintf( ioQQQ, " HydroLevel Z=%2ld simple II/I=%10.2e so not doing equilibrium.\n",
00300 nelem, iso.xIonSimple[ipH_LIKE][nelem] );
00301 }
00302
00303
00304 for( i=ipH1s; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
00305 {
00306 iso.DepartCoef[ipH_LIKE][nelem][i] = 1.;
00307 iso.Pop2Ion[ipH_LIKE][nelem][i] = 0.;
00308 }
00309
00310 ionbal.RateRecomTot[nelem][nelem] = 0.;
00311 iso.xIonSimple[ipH_LIKE][nelem] = 0.;
00312 iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 0.;
00313
00314 }
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 else if( (strcmp( iso.chTypeAtomSet[ipH_LIKE] , "LOWT" )==0) ||
00336 (iso.xIonSimple[ipH_LIKE][nelem] < TooSmall) )
00337 # if 0
00338 ( (nelem==ipHYDROGEN) && (iso.xIonSimple[ipH_LIKE][nelem] < 1e-28) ) ||
00339 ( (nelem==ipHELIUM) && (iso.xIonSimple[ipH_LIKE][nelem] < 1e-18) ) ||
00340
00341
00342 ( (nelem>ipHELIUM) && (iso.xIonSimple[ipH_LIKE][nelem] < 2e-18) ) )
00343 # endif
00344 {
00345 strcpy( iso.chTypeAtomUsed[ipH_LIKE][nelem], "Low T" );
00346
00347
00348
00349
00350
00351
00352 HydroT2Low( nelem );
00353 }
00354
00355
00356
00357
00358 else
00359 {
00360
00361 strcpy( iso.chTypeAtomUsed[ipH_LIKE][nelem], "Popul" );
00362 HydroLevelPop(nelem );
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372 for( ipLo=ipH1s; ipLo < (iso.numLevels_local[ipH_LIKE][nelem] - 1); ipLo++ )
00373 {
00374 for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
00375 {
00376
00377 double PopOpcSave = EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc;
00378
00379
00380 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo =
00381 iso.Pop2Ion[ipH_LIKE][nelem][ipLo];
00382
00383
00384 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi =
00385 iso.Pop2Ion[ipH_LIKE][nelem][ipHi];
00386
00387
00388 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc =
00389 iso.Pop2Ion[ipH_LIKE][nelem][ipLo] - EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi*
00390 iso.stat[ipH_LIKE][nelem][ipLo]/iso.stat[ipH_LIKE][nelem][ipHi];
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400 if( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].TauIn > 0.01 &&
00401 ipLo>2 && EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc /
00402 SDIV( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo ) <0.9 )
00403 {
00404 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc = (0.75*PopOpcSave+
00405 0.25*EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc);
00406 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc = MIN2(
00407 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc ,
00408 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo );
00409
00410 }
00411 }
00412 }
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 for( ipHi=iso.nTopOff[ipH_LIKE][nelem]-1; ipHi<iso.numLevels_local[ipH_LIKE][nelem]; ++ipHi )
00432 {
00433 for( ipLo=ipH1s; ipLo < ipHi; ++ipLo )
00434 {
00435
00436 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc =
00437
00438 EmisLines[ipH_LIKE][nelem][iso.numLevels_local[ipH_LIKE][nelem]-1][ipLo].PopLo;
00439 }
00440 }
00441
00442
00443
00444
00445 for( ipHi=ipH2p+1; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
00446 {
00447
00448 double sum = EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].PopOpc +
00449 EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].PopOpc;
00450
00451
00452
00453
00454 EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].PopOpc = sum;
00455 EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].PopOpc = sum;
00456 }
00457
00458
00459
00460
00461 sum_popn_ov_ion = 0.;
00462
00463
00464 ionbal.RateIonizTot[nelem][nelem] = 0.;
00465
00466
00467 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][nelem]; level++ )
00468 {
00469
00470 ionbal.RateIonizTot[nelem][nelem] +=
00471 iso.RateLevel2Cont[ipH_LIKE][nelem][level]*iso.Pop2Ion[ipH_LIKE][nelem][level];
00472
00473 sum_popn_ov_ion += iso.Pop2Ion[ipH_LIKE][nelem][level];
00474 }
00475
00476
00477 if( ( ionbal.RateIonizTot[nelem][nelem]/MAX2(SMALLDOUBLE , sum_popn_ov_ion) ) > BIGDOUBLE )
00478 {
00479 fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.",
00480 nelem+1, nelem);
00481 cdEXIT(EXIT_FAILURE);
00482 }
00483 else
00484 ionbal.RateIonizTot[nelem][nelem] /= MAX2(SMALLDOUBLE , sum_popn_ov_ion);
00485
00486
00487 if( sum_popn_ov_ion < 0. )
00488 {
00489 fprintf( ioQQQ,
00490 " HydroLevel finds negative H-like ion fraction for nelem=%2ld %s using routine %s, val= %.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
00491 nelem,
00492 elementnames.chElementSym[nelem],
00493 iso.chTypeAtomUsed[ipH_LIKE][nelem] ,
00494 sum_popn_ov_ion,
00495 iso.xIonSimple[ipH_LIKE][nelem],
00496 phycon.te,
00497 nzone );
00498 fprintf( ioQQQ, " level pop are:" );
00499
00500 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00501 {
00502 fprintf( ioQQQ,PrintEfmt("%8.1e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
00503 }
00504 fprintf( ioQQQ, "\n" );
00505 ContNegative();
00506 ShowMe();
00507 puts( "[Stop in HydroLevel]" );
00508 cdEXIT(EXIT_FAILURE);
00509 }
00510
00511
00512
00513 else if( sum_popn_ov_ion >= 0. && sum_popn_ov_ion < 1e-30 )
00514 {
00515 iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 0.;
00516
00517
00518 dense.xIonDense[nelem][nelem+1] = 0.;
00519
00520
00521
00522 dense.IonHigh[nelem] = nelem;
00523 ionbal.RateIonizTot[nelem][nelem] = 0.;
00524
00525
00526
00527 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
00528 {
00529 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00530 {
00531
00532 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo = 0.;
00533
00534
00535 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopHi = 0.;
00536
00537
00538 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc = 0.;
00539
00540
00541 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ots = 0.;
00542 }
00543 }
00544 }
00545
00546 else
00547 {
00548
00549
00550
00551
00552 iso.pop_ion_ov_neut[ipH_LIKE][nelem] = 1./sum_popn_ov_ion;
00553
00554
00555
00556
00557 # ifndef NDEBUG
00558 {
00559 double pop_ion_ov_neut = ionbal.RateIonizTot[nelem][nelem] /
00560 SDIV(ionbal.RateRecomTot[nelem][nelem]);
00561
00562 double ratio_error = fabs(iso.pop_ion_ov_neut[ipH_LIKE][nelem] - pop_ion_ov_neut ) /
00563
00564
00565
00566
00567
00568
00569
00570
00571 MAX2(1e-7,iso.pop_ion_ov_neut[ipH_LIKE][nelem]);
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583 ASSERT( conv.lgSearch || dynamics.lgAdvection ||
00584
00585
00586 ionbal.RateIonizTot[nelem][nelem] < SMALLFLOAT ||
00587
00588 iso.pop_ion_ov_neut[ipH_LIKE][nelem] < 1e-5 ||
00589 ratio_error < 0.001 );
00590 # if 0
00591 ||
00592
00593
00594 hydro.xLymanPumpingScaleFactor > 1. );
00595 # endif
00596 }
00597 # endif
00598 }
00599
00600
00601 if( (trace.lgIsoTraceFull[ipH_LIKE] && trace.lgTrace) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
00602 {
00603 fprintf( ioQQQ, " HLEV HGAMNC" );
00604 PrintE93( ioQQQ, iso.gamnc[ipH_LIKE][nelem][ipH1s] );
00605
00606 for( i=ipH2s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00607 {
00608 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.gamnc[ipH_LIKE][nelem][i] ));
00609 }
00610 fprintf( ioQQQ, "\n" );
00611
00612 fprintf( ioQQQ, " HLEV TOTCAP" );
00613
00614 for( i=1; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00615 {
00616 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RateCont2Level[ipH_LIKE][nelem][i]/dense.eden ));
00617 }
00618 fprintf( ioQQQ," tot");
00619 fprintf( ioQQQ,PrintEfmt("%10.2e", ionbal.RateRecomTot[ipH_LIKE][nelem]/dense.eden ) );
00620 fprintf( ioQQQ, "\n" );
00621
00622 fprintf( ioQQQ, " HLEV IND Rc" );
00623
00624 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00625 {
00626 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.RecomInducRate[ipH_LIKE][nelem][i]*iso.PopLTE[ipH_LIKE][nelem][i] ));
00627 }
00628 fprintf( ioQQQ, "\n" );
00629
00630
00631 fprintf( ioQQQ, " IND Rc LTE " );
00632
00633 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00634 {
00635 fprintf(ioQQQ,PrintEfmt("%9.2e",
00636 iso.gamnc[ipH_LIKE][nelem][i]*iso.PopLTE[ipH_LIKE][nelem][i] ));
00637 }
00638 fprintf( ioQQQ, "\n" );
00639
00640
00641 fprintf( ioQQQ, " HLEV HLTE" );
00642
00643 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00644 {
00645 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.PopLTE[ipH_LIKE][nelem][i] ));
00646 }
00647 fprintf( ioQQQ, "\n" );
00648
00649
00650 fprintf( ioQQQ, " HLEVfr cion" );
00651
00652 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00653 {
00654 fprintf(ioQQQ,PrintEfmt("%9.2e",
00655 iso.ColIoniz[ipH_LIKE][nelem][i]*
00656 iso.Pop2Ion[ipH_LIKE][nelem][i]*collider/MAX2(1e-30,iso.RateLevel2Cont[ipH_LIKE][nelem][i]) ) );
00657 }
00658 fprintf( ioQQQ, "\n" );
00659
00660
00661 if( ionbal.RateRecomTot[nelem][nelem]> 0. )
00662 {
00663 fprintf( ioQQQ, " HLEVfrPhIon" );
00664
00665 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00666 {
00667 fprintf(ioQQQ,PrintEfmt("%9.2e",
00668 iso.gamnc[ipH_LIKE][nelem][i]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][i]/MAX2(1e-30,iso.RateLevel2Cont[ipH_LIKE][nelem][i]) ) );
00669 }
00670 fprintf( ioQQQ, "\n" );
00671 }
00672
00673 fprintf( ioQQQ, " HLEV HN" );
00674
00675 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00676 {
00677 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
00678 }
00679 fprintf( ioQQQ, "\n" );
00680
00681 fprintf( ioQQQ, " HLEV b(n)" );
00682
00683 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00684 {
00685 fprintf(ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipH_LIKE][nelem][i] ));
00686 }
00687 fprintf( ioQQQ, "\n" );
00688
00689 fprintf( ioQQQ, " HLEV X12tot");
00690 fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.x12tot ));
00691 fprintf( ioQQQ," Grn dest:");
00692 fprintf(ioQQQ,PrintEfmt("%9.2e",
00693 ionbal.RateIonizTot[nelem][nelem] ));
00694 fprintf(ioQQQ, "\n");
00695 }
00696
00697
00698
00699
00700 if( iso.pop_ion_ov_neut[ipH_LIKE][nelem] > 0. )
00701 {
00702
00703 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][nelem]; level++ )
00704 {
00705
00706
00707
00708
00709 if( iso.Pop2Ion[ipH_LIKE][nelem][level] < 0. )
00710 {
00711 fprintf( ioQQQ,
00712 " HydroLevel finds negative hydrogen level population for %s nelem=%ld level %ld value=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
00713 elementnames.chElementName[nelem],
00714 nelem, level,
00715 iso.Pop2Ion[ipH_LIKE][nelem][level],
00716 iso.xIonSimple[ipH_LIKE][nelem],
00717 phycon.te, nzone );
00718 fprintf( ioQQQ, " level pop are:" );
00719
00720 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00721 {
00722 fprintf( ioQQQ,PrintEfmt("%8.1e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
00723 }
00724 fprintf( ioQQQ, "\n" );
00725 ContNegative();
00726 ShowMe();
00727 puts( "[Stop in HydroLevel]" );
00728 cdEXIT(EXIT_FAILURE);
00729 }
00730 }
00731 }
00732 else if( iso.pop_ion_ov_neut[ipH_LIKE][nelem] < 0.)
00733 {
00734 ionbal.RateIonizTot[nelem][nelem] = 0.;
00735
00736 fprintf( ioQQQ,
00737 " HydroLevel finds negative hydrogen ion frac for nelem=%ld (%s) value=%.3e simple= %.3e TE=%.3e ZONE=%ld\n",
00738 nelem,
00739 elementnames.chElementSym[nelem],
00740 iso.Pop2Ion[ipH_LIKE][nelem][level],
00741 iso.xIonSimple[ipH_LIKE][nelem],
00742 phycon.te, nzone );
00743 fprintf( ioQQQ, " level pop are:" );
00744
00745 for( i=ipH1s; i < iso.numLevels_local[ipH_LIKE][nelem]; i++ )
00746 {
00747 fprintf( ioQQQ,PrintEfmt("%8.1e", iso.Pop2Ion[ipH_LIKE][nelem][i] ));
00748 }
00749 fprintf( ioQQQ, "\n" );
00750 ContNegative();
00751 ShowMe();
00752 puts( "[Stop in HydroLevel]" );
00753 cdEXIT(EXIT_FAILURE);
00754 }
00755
00756
00757
00758
00759
00760
00761
00762 HydroRenorm();
00763
00764 if( trace.lgTrace )
00765 {
00766
00767
00768
00769 fprintf( ioQQQ, " HydroLevel%3ld retrn %s te=",
00770 nelem,
00771 iso.chTypeAtomUsed[ipH_LIKE][nelem] );
00772 PrintE93( ioQQQ,phycon.te);
00773 fprintf( ioQQQ," ion/atom=%.4e",iso.pop_ion_ov_neut[ipH_LIKE][nelem]);
00774
00775 fprintf( ioQQQ," simple=%.4e",iso.xIonSimple[ipH_LIKE][nelem]);
00776
00777 fprintf( ioQQQ," b1=%.2e",iso.DepartCoef[ipH_LIKE][nelem][ipH1s]);
00778
00779 fprintf( ioQQQ," ion rate=%.4e",ionbal.RateIonizTot[nelem][nelem]);
00780
00781 fprintf( ioQQQ," TotRec=%.4e",ionbal.RateRecomTot[nelem][nelem]);
00782
00783 fprintf( ioQQQ," RadRec=%.4e",iso.RadRec_effec[ipH_LIKE][nelem]);
00784 fprintf( ioQQQ, "\n");
00785 }
00786 #if 0
00787 {
00788 static int ncall=0;
00789 if( nelem==0)
00790 ++ncall;
00791 if( ncall > 11 )
00792 {
00793 fprintf(ioQQQ,"debugggg exit\n");
00794 exit(1);
00795 }
00796 }
00797 #endif
00798
00799 ASSERT( fabs( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].PopLo - iso.Pop2Ion[ipH_LIKE][nelem][0] ) /
00800 SDIV( iso.Pop2Ion[ipH_LIKE][nelem][0] ) < 1e-4 || EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].PopLo < 1e-30 );
00801 # if !defined(NDEBUG)
00802
00803 for( ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
00804 {
00805 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00806 {
00807
00808
00809 ASSERT( (ipLo==ipH2s || ipLo==ipH2p) ||
00810 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc <=
00811 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopLo );
00812 ASSERT( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].PopOpc <BIGFLOAT );
00813 }
00814 }
00815 # endif
00816
00817 DEBUG_EXIT( "HydroLevel()" );
00818 return;
00819 }
00820
00821
00822