00001
00002
00003
00004 #include "cddefines.h"
00005 #include "atmdat.h"
00006 #include "continuum.h"
00007 #include "conv.h"
00008 #include "dense.h"
00009 #include "dynamics.h"
00010 #include "elementnames.h"
00011 #include "grainvar.h"
00012 #include "he.h"
00013 #include "heavy.h"
00014 #include "helike.h"
00015 #include "hmi.h"
00016 #include "hydrogenic.h"
00017 #include "ionbal.h"
00018 #include "iso.h"
00019 #include "mole.h"
00020 #include "opacity.h"
00021 #include "phycon.h"
00022 #include "physconst.h"
00023 #include "rfield.h"
00024 #include "secondaries.h"
00025 #include "taulines.h"
00026 #include "thirdparty.h"
00027 #include "trace.h"
00028 #include "yield.h"
00029
00030
00031 void iso_level( const long int ipISO, const long int nelem, double &renorm)
00032 {
00033 long int ipHi,
00034 ipLo,
00035 i,
00036 level,
00037 level_error;
00038
00039 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00040 const long int numlevels_local = sp->numLevels_local;
00041
00042 double BigError;
00043
00044 int32 nerror;
00045 double HighestPColOut = 0.;
00046 bool lgNegPop=false;
00047 valarray<int32> ipiv(numlevels_local);
00048
00049 double source=0.,
00050 sink=0.;
00051 valarray<double> PopPerN(sp->n_HighestResolved_local+1);
00052
00053 DEBUG_ENTRY( "iso_level()" );
00054 renorm = 1.;
00055
00056
00057 ASSERT( nelem >= ipISO );
00058 ASSERT( nelem < LIMELM );
00059
00060
00061
00062 ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
00063 fabs( (sp->fb[0].ColIoniz* dense.EdenHCorr) /
00064 SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );
00065
00066 if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
00067 (dense.IonLow[nelem] <= nelem - ipISO) &&
00068 ionbal.RateRecomIso[nelem][ipISO] > 0. )
00069 {
00070
00071 sp->xIonSimple = ionbal.RateIonizTot(nelem,nelem-ipISO)/ionbal.RateRecomIso[nelem][ipISO];
00072 }
00073 else
00074 {
00075 sp->xIonSimple = 0.;
00076 }
00077
00078
00079 if( dense.xIonDense[nelem][nelem+1-ipISO] < SMALLFLOAT || sp->xIonSimple < 1e-35 )
00080 {
00081
00082 strcpy( sp->chTypeAtomUsed, "zero " );
00083 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
00084 {
00085 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n",
00086 ipISO, nelem, sp->xIonSimple, sp->chTypeAtomUsed );
00087 }
00088
00089
00090 lgNegPop = false;
00091 sp->st[0].Pop() = dense.xIonDense[nelem][nelem-ipISO];
00092 for( long n=1; n < numlevels_local; n++ )
00093 {
00094 sp->st[n].Pop() = 0.;
00095 }
00096 sp->qTot2S = 0.;
00097 }
00098 else
00099 {
00100
00101 valarray<double>
00102 creation(numlevels_local),
00103 error(numlevels_local),
00104 work(numlevels_local),
00105 Save_creation(numlevels_local);
00106
00107 ASSERT( dense.xIonDense[nelem][nelem+1-ipISO] >= 0.f );
00108 for( level=0; level < numlevels_local; level++ )
00109 {
00110
00111 creation[level] = sp->fb[level].RateCont2Level * dense.xIonDense[nelem][nelem+1-ipISO];
00112 }
00113
00114 {
00115
00116 if( nelem==ipHYDROGEN )
00117 {
00118
00119
00120 source += atmdat.HCharExcRecTotal * dense.xIonDense[nelem][nelem+1-ipISO];
00121 }
00122 else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
00123 {
00124
00125 source += atmdat.HeCharExcRecTotal * dense.xIonDense[nelem][nelem+1-ipISO];
00126 }
00127 else
00128 {
00129 source +=
00130 (atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] +
00131 atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0] )
00132 * dense.xIonDense[nelem][nelem+1-ipISO];
00133 }
00134 if ( nelem > ipISO )
00135 {
00136 if( nelem==ipHELIUM )
00137 sink += atmdat.HeCharExcRecTotal;
00138 else
00139 sink +=
00140 atmdat.HeCharExcRecTo[nelem][nelem-(ipISO+1)]*dense.xIonDense[ipHELIUM][0] +
00141 atmdat.HCharExcRecTo[nelem][nelem-(ipISO+1)]*dense.xIonDense[ipHYDROGEN][0];
00142 }
00143 }
00144
00145 {
00146 if ( nelem==ipHYDROGEN )
00147 sink += atmdat.HCharExcIonTotal;
00148 else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
00149 sink += atmdat.HeCharExcIonTotal;
00150 else
00151 sink +=
00152 atmdat.HeCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][1] +
00153 atmdat.HCharExcIonOf[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][1];
00154
00155 if ( nelem > ipISO )
00156 {
00157 if( nelem==ipHELIUM )
00158 source += atmdat.HeCharExcIonTotal*dense.xIonDense[nelem][nelem-(ipISO+1)];
00159 else
00160 source +=
00161 (atmdat.HeCharExcIonOf[nelem][nelem-(ipISO+1)]*dense.xIonDense[ipHELIUM][1] +
00162 atmdat.HCharExcIonOf[nelem][nelem-(ipISO+1)]*dense.xIonDense[ipHYDROGEN][1]) *
00163 dense.xIonDense[nelem][nelem-(ipISO+1)];
00164 }
00165
00166 }
00167
00168
00169 multi_arr<double,2,C_TYPE> z, SaveZ;
00170 z.alloc(numlevels_local,numlevels_local);
00171 z.zero();
00172 vector<double> Boltzmann(numlevels_local-1);
00173 for (unsigned lev = 0; lev < Boltzmann.size(); ++lev)
00174 Boltzmann[lev] = 1.0;
00175
00176
00177
00178 ASSERT( NISO == 2 );
00179 long SpinUsed[NISO] = {2,1};
00180 long indexNmaxP =
00181 sp->QuantumNumbers2Index[ sp->n_HighestResolved_local ][1][SpinUsed[ipISO]];
00182
00183 strcpy( sp->chTypeAtomUsed, "Popul" );
00184 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
00185 {
00186 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n",
00187 ipISO, nelem, sp->chTypeAtomUsed );
00188 }
00189
00190
00191 qList::const_iterator StElm = sp->st.begin();
00192
00193
00194 for( level=0; level < numlevels_local; level++ )
00195 {
00196
00197
00198 z[level][level] += sp->fb[level].RateLevel2Cont;
00199
00200 if( level == 1 )
00201
00202 sp->qTot2S = sp->fb[level].ColIoniz;
00203
00204
00205
00206 if (level != 0)
00207 {
00208 double arg = (StElm[level].energy().K()-StElm[level-1].energy().K()) /
00209 phycon.te;
00210 arg = MAX2( -38., arg); fixit();
00211 double bstep = sexp( arg );
00212 for ( ipLo = 0; ipLo < level; ++ipLo )
00213 Boltzmann[ipLo] *= bstep;
00214 }
00215
00216
00217 for( ipLo=0; ipLo < level; ipLo++ )
00218 {
00219 double coll_down = sp->trans(level,ipLo).Coll().ColUL( colliders );
00220
00221 double RadDecay = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().Aul()*
00222 (sp->trans(level,ipLo).Emis().Pesc() +
00223 sp->trans(level,ipLo).Emis().Pelec_esc() +
00224 sp->trans(level,ipLo).Emis().Pdest())*
00225 KILL_BELOW_PLASMA(sp->trans(level,ipLo).EnergyRyd()) );
00226
00227 double pump = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().pump() *
00228 KILL_BELOW_PLASMA(sp->trans(level,ipLo).EnergyRyd()) );
00229
00230 if( iso_ctrl.lgRandErrGen[ipISO] )
00231 {
00232 coll_down *= sp->ex[level][ipLo].ErrorFactor[IPCOLLIS];
00233 RadDecay *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
00234 pump *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
00235 }
00236
00237 double coll_up = coll_down *
00238 (double)StElm[level].g()/
00239 (double)StElm[ipLo].g()*
00240 Boltzmann[ipLo];
00241
00242 z[ipLo][ipLo] += coll_up + pump ;
00243 z[ipLo][level] -= coll_up + pump ;
00244
00245 double pump_down = pump *
00246 (double)StElm[ipLo].g()/
00247 (double)StElm[level].g();
00248
00249 z[level][level] += RadDecay + pump_down + coll_down;
00250 z[level][ipLo] -= RadDecay + pump_down + coll_down;
00251
00252 if( level == indexNmaxP )
00253 {
00254 HighestPColOut += coll_down;
00255 }
00256 if( ipLo == indexNmaxP )
00257 {
00258 HighestPColOut += coll_up;
00259 }
00260
00261
00262 if( (level == 1) && (ipLo==0) )
00263 {
00264 sp->qTot2S += coll_down/dense.EdenHCorr;
00265 }
00266 if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S()==1) )
00267 {
00268 sp->qTot2S += coll_up/dense.EdenHCorr;
00269 }
00270 }
00271 }
00272
00273 if( ipISO == nelem )
00274 {
00275
00276
00277
00278 if( HighestPColOut < 1./sp->st[indexNmaxP].lifetime() )
00279 {
00280 iso_ctrl.lgCritDensLMix[ipISO] = false;
00281 }
00282 }
00283
00285 ASSERT( ipISO <= ipHE_LIKE );
00286 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
00287 {
00288
00289
00290
00291
00292 fixit();
00293
00294 z[tnu->ipHi][tnu->ipLo] -= tnu->AulTotal;
00295 z[tnu->ipHi][tnu->ipHi] += tnu->AulTotal;
00296
00297 }
00298
00299 if (iso_ctrl.lgInd2nu_On)
00300 {
00301 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
00302 {
00303
00304
00305
00306
00307 z[tnu->ipHi][tnu->ipLo] -= tnu->induc_dn;
00308 z[tnu->ipHi][tnu->ipHi] += tnu->induc_dn;
00309
00310
00311 z[tnu->ipLo][tnu->ipHi] -= tnu->induc_up;
00312 z[tnu->ipLo][tnu->ipLo] += tnu->induc_up;
00313 }
00314 }
00315
00316
00317 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00318 {
00319 if( ion!=nelem-ipISO )
00320 {
00321 source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] *
00322 dense.xIonDense[nelem][ion];
00323 sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
00324
00325 source += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] *
00326 dense.xIonDense[nelem][ion] * atmdat.lgCTOn;
00327 sink += mole.xMoleChTrRate[nelem][nelem-ipISO][ion] * atmdat.lgCTOn;
00328
00329 }
00330 }
00331
00332
00333 source += mole.source[nelem][nelem-ipISO];
00334 sink += mole.sink[nelem][nelem-ipISO];
00335
00336 #if 1
00337
00338
00339 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
00340 dynamics.Rate != 0.0 &&
00341 !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
00342 {
00343
00344 source += dynamics.Source[nelem][nelem-ipISO];
00345
00346 sink += dynamics.Rate;
00347 }
00348 #else
00349
00350 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
00351 dynamics.Rate != 0.0 &&
00352 !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
00353 {
00354 for( level=0; level < numlevels_local; level++ )
00355 {
00356 creation[level] += dynamics.StatesElem[nelem][nelem-ipISO][level];
00357 z[level][level] += dynamics.Rate;
00358 }
00359 }
00360 #endif
00361
00362
00363 for(long ion_from=dense.IonLow[nelem]; ion_from < MIN2( dense.IonHigh[nelem], nelem-ipISO ) ; ion_from++ )
00364 {
00365 if( ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] >= 0. )
00366 {
00367
00368 source += ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] * dense.xIonDense[nelem][ion_from];
00369 }
00370
00371 if( ion_from == nelem-1-ipISO )
00372 sink += ionbal.RateRecomTot[nelem][ion_from];
00373 }
00374
00375 ASSERT( source >= 0.f );
00376 if (0)
00377 {
00378 creation[0] += source;
00379 for( level=0; level < numlevels_local; level++ )
00380 {
00381 z[level][level] += sink;
00382 }
00383 }
00384 else
00385 {
00386
00387 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00388 double partfun=0.0;
00389 for ( level = 0; level < numlevels_local; level++ )
00390 {
00391 partfun += sp->st[level].Boltzmann()*sp->st[level].g();
00392 }
00393 source /= partfun;
00394 for( level=0; level < numlevels_local; level++ )
00395 {
00396 creation[level] += source*
00397 sp->st[level].Boltzmann()*sp->st[level].g();
00398 z[level][level] += sink;
00399 }
00400 }
00401
00402
00403 if( sp->trans(iso_ctrl.nLyaLevel[ipISO],0).Coll().rate_lu_nontherm() * iso_ctrl.lgColl_excite[ipISO] > 0. )
00404 {
00405
00406 for( level=1; level < numlevels_local; level++ )
00407 {
00408 double RateUp , RateDown;
00409
00410 RateUp = sp->trans(level,0).Coll().rate_lu_nontherm();
00411 RateDown = RateUp * (double)sp->st[ipH1s].g() /
00412 (double)sp->st[level].g();
00413
00414
00415 z[ipH1s][ipH1s] += RateUp;
00416
00417
00418 z[level][ipH1s] -= RateDown;
00419
00420
00421 z[ipH1s][level] -= RateUp;
00422
00423 z[level][level] += RateDown;
00424 }
00425 }
00426
00427
00428
00429
00430
00431
00432
00433 SaveZ = z;
00434
00435 for( ipLo=0; ipLo < numlevels_local; ipLo++ )
00436 Save_creation[ipLo] = creation[ipLo];
00437
00438 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
00439 {
00440 const long numlevels_print = numlevels_local;
00441 fprintf( ioQQQ, " pop level others => (iso_level)\n" );
00442 for( long n=0; n < numlevels_print; n++ )
00443 {
00444 fprintf( ioQQQ, " %s %s %2ld", iso_ctrl.chISO[ipISO], elementnames.chElementNameShort[nelem], n );
00445 for( long j=0; j < numlevels_print; j++ )
00446 {
00447 fprintf( ioQQQ,"\t%.9e", z[j][n] );
00448 }
00449 fprintf( ioQQQ, "\n" );
00450 }
00451 fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
00452 atmdat.HeCharExcRecTotal,
00453 findspecieslocal("CO")->den ,
00454 atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] ,
00455 atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0] );
00456 fprintf(ioQQQ," recomb ");
00457 for( long n=0; n < numlevels_print; n++ )
00458 {
00459 fprintf( ioQQQ,"\t%.9e", creation[n] );
00460 }
00461 fprintf( ioQQQ, "\n" );
00462 }
00463
00464 nerror = 0;
00465
00466 getrf_wrapper(numlevels_local,numlevels_local,
00467 z.data(),numlevels_local,&ipiv[0],&nerror);
00468
00469 getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0],
00470 &creation[0],numlevels_local,&nerror);
00471
00472 if( nerror != 0 )
00473 {
00474 fprintf( ioQQQ, " iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
00475 cdEXIT(EXIT_FAILURE);
00476 }
00477
00478
00479
00480 for( level=ipH1s; level < numlevels_local; level++ )
00481 {
00482 double qn = 0., qx = 0.;
00483 error[level] = 0.;
00484 for( long n=ipH1s; n < numlevels_local; n++ )
00485 {
00486 double q = SaveZ[n][level]*creation[n];
00487
00488
00489 if ( q > qx )
00490 qx = q;
00491 else if (q < qn)
00492 qn = q;
00493
00494 error[level] += q;
00495 }
00496
00497 if (-qn > qx)
00498 qx = -qn;
00499
00500 if( qx > 0. )
00501 {
00502 error[level] = (error[level] - Save_creation[level])/qx;
00503 }
00504 else
00505 {
00506 error[level] = 0.;
00507 }
00508 }
00509
00510
00511 BigError = -1.;
00512 level_error = -1;
00513
00514 for( level=ipH1s; level < numlevels_local; level++ )
00515 {
00516 double abserror;
00517 abserror = fabs( error[level]);
00518
00519 if( abserror > BigError )
00520 {
00521 BigError = abserror;
00522 level_error = level;
00523 }
00524 }
00525
00526
00527
00528 if( BigError > FLT_EPSILON )
00529 {
00530 if( !conv.lgSearch )
00531 fprintf(ioQQQ," PROBLEM" );
00532
00533 fprintf(ioQQQ,
00534 " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
00535 "level was %li Search?%c \n",
00536 fnzone,
00537 ipISO,
00538 elementnames.chElementName[nelem],
00539 nelem ,
00540 BigError ,
00541 level_error,
00542 TorF(conv.lgSearch) );
00543 }
00544
00545
00546 if ( iso_ctrl.lgLTE_levels[ipISO] )
00547 {
00548 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00549 double partfun=0.0;
00550 for ( level = 0; level < numlevels_local; level++ )
00551 {
00552 partfun += sp->st[level].Boltzmann()*sp->st[level].g();
00553 }
00554 double scale = dense.xIonDense[nelem][nelem-ipISO]/partfun;
00555 for ( level = 0; level < numlevels_local; level++ )
00556 {
00557 creation[level] = sp->st[level].Boltzmann()*sp->st[level].g()*scale;
00558 }
00559 }
00560
00561 for( level = numlevels_local-1; level > 0; --level )
00562 {
00563
00564 if( creation[level] < 0. )
00565 lgNegPop = true;
00566 }
00567
00568 if( lgNegPop && dense.lgSetIoniz[nelem] )
00569 {
00570
00571
00572
00573 for( level = 1; level < numlevels_local; ++level )
00574 creation[level] = 0.;
00575 creation[0] = dense.xIonDense[nelem][nelem-ipISO];
00576
00577 lgNegPop = false;
00578 }
00579
00580
00581 for( level=0; level < numlevels_local; level++ )
00582 {
00583 ASSERT( creation[level] >= 0. );
00584 sp->st[level].Pop() = creation[level];
00585
00586 if( sp->st[level].Pop() <= 0 && !conv.lgSearch )
00587 {
00588 fprintf(ioQQQ,
00589 "PROBLEM non-positive level pop for iso = %li, nelem = "
00590 "%li = %s, level=%li val=%.3e nTotalIoniz %li\n",
00591 ipISO,
00592 nelem ,
00593 elementnames.chElementSym[nelem],
00594 level,
00595 sp->st[level].Pop() ,
00596 conv.nTotalIoniz);
00597 }
00598 }
00599
00600
00601 for( level=numlevels_local; level < sp->numLevels_max; level++ )
00602 {
00603 sp->st[level].Pop() = 0.;
00604
00605
00606 }
00607
00608
00609
00610 iso_renorm( nelem, ipISO, renorm );
00611
00612 double TotalPopExcited = 0.;
00613
00614 for( level=1; level < numlevels_local; level++ )
00615 TotalPopExcited += sp->st[level].Pop();
00616 ASSERT( TotalPopExcited >= 0. );
00617 double TotalPop = TotalPopExcited + sp->st[0].Pop();
00618
00619
00620 if( dense.lgSetIoniz[nelem] )
00621 {
00622 if( !fp_equal( TotalPop, (double)dense.xIonDense[nelem][nelem-ipISO] ) )
00623 {
00624 if( TotalPopExcited >= dense.xIonDense[nelem][nelem-ipISO] )
00625 {
00626 for( level=0; level < numlevels_local; level++ )
00627 sp->st[level].Pop() *=
00628 dense.xIonDense[nelem][nelem-ipISO] / TotalPop;
00629 }
00630 else
00631 {
00632 sp->st[0].Pop() =
00633 MAX2( 1e-30 * dense.xIonDense[nelem][nelem-ipISO],
00634 dense.xIonDense[nelem][nelem-ipISO] - TotalPopExcited );
00635 }
00636 sp->lgPopsRescaled = true;
00637 }
00638 ASSERT( sp->st[0].Pop() >= 0. );
00639 }
00640 }
00641
00642
00643
00644 if( lgNegPop || dense.xIonDense[nelem][nelem-ipISO] < 0. )
00645 {
00646 fprintf( ioQQQ,
00647 " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\
00648 "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
00649 ipISO,
00650 nelem,
00651 elementnames.chElementSym[nelem],
00652 sp->chTypeAtomUsed,
00653 dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO]),
00654 sp->xIonSimple,
00655 phycon.te,
00656 nzone );
00657
00658 fprintf( ioQQQ, " level pop are:\n" );
00659 for( i=0; i < numlevels_local; i++ )
00660 {
00661 fprintf( ioQQQ,PrintEfmt("%8.1e", sp->st[i].Pop() ));
00662 if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
00663 }
00664 fprintf( ioQQQ, "\n" );
00665 ContNegative();
00666 ShowMe();
00667 cdEXIT(EXIT_FAILURE);
00668 }
00669
00670 for( ipHi=1; ipHi<numlevels_local; ++ipHi )
00671 {
00672 for( ipLo=0; ipLo<ipHi; ++ipLo )
00673 {
00674 if( sp->trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
00675 continue;
00676
00677
00678 sp->trans(ipHi,ipLo).Emis().PopOpc() =
00679 sp->st[ipLo].Pop() - sp->st[ipHi].Pop()*
00680 sp->st[ipLo].g()/sp->st[ipHi].g();
00681
00682
00683 if( N_(ipHi) > sp->n_HighestResolved_local || opac.lgCaseB )
00684 sp->trans(ipHi,ipLo).Emis().PopOpc() = MAX2( 0., sp->trans(ipHi,ipLo).Emis().PopOpc() );
00685 }
00686 }
00687
00688 return;
00689 }
00690
00691 void iso_set_ion_rates( long ipISO, long nelem)
00692 {
00693 DEBUG_ENTRY("iso_set_ion_rates()");
00694 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00695 const long int numlevels_local = sp->numLevels_local;
00696
00697
00698 double TotalPop = 0.;
00699 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = 0.;
00700 for( long level=0; level < numlevels_local; level++ )
00701 {
00702
00703
00704 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] +=
00705 sp->st[level].Pop() * sp->fb[level].RateLevel2Cont;
00706 TotalPop += sp->st[level].Pop();
00707 }
00708
00709 if( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] > BIGDOUBLE )
00710 {
00711 fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.",
00712 nelem+1, nelem-ipISO);
00713 cdEXIT(EXIT_FAILURE);
00714 }
00715
00716 if (TotalPop <= SMALLFLOAT)
00717 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = sp->fb[0].RateLevel2Cont;
00718 else
00719 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] /= TotalPop;
00720
00721 if( ionbal.RateRecomIso[nelem][ipISO] > 0. )
00722 sp->xIonSimple = ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1]/ionbal.RateRecomIso[nelem][ipISO];
00723 else
00724 sp->xIonSimple = 0.;
00725
00726 ASSERT( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] >= 0. );
00727
00728 if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 )
00729 {
00730
00731 double ratio;
00732 double rateOutOf2TripS = sp->st[ipHe2s3S].Pop() * sp->fb[ipHe2s3S].RateLevel2Cont;
00733 if( rateOutOf2TripS > SMALLFLOAT )
00734 {
00735 ratio = rateOutOf2TripS /
00736 ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
00737 }
00738 else
00739 {
00740 ratio = 0.;
00741 }
00742 if( ratio > he.frac_he0dest_23S )
00743 {
00744
00745 he.nzone = nzone;
00746 he.frac_he0dest_23S = ratio;
00747 rateOutOf2TripS = sp->st[ipHe2s3S].Pop() *sp->fb[ipHe2s3S].gamnc;
00748 if( rateOutOf2TripS > SMALLFLOAT )
00749 {
00750 he.frac_he0dest_23S_photo = rateOutOf2TripS /
00751 ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
00752 }
00753 else
00754 {
00755 he.frac_he0dest_23S_photo = 0.;
00756 }
00757 }
00758 }
00759 }