31void iso_level(
const long int ipISO,
const long int nelem,
double &renorm)
45 double HighestPColOut = 0.;
47 valarray<int32> ipiv(numlevels_local);
63 fabs( (sp->
fb[0].ColIoniz*
dense.EdenHCorr) /
64 SDIV(
ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );
66 if( (
dense.IonHigh[nelem] >= nelem - ipISO) &&
67 (
dense.IonLow[nelem] <= nelem - ipISO) &&
68 ionbal.RateRecomIso[nelem][ipISO] > 0. )
83 if(
trace.lgTrace && (nelem ==
trace.ipIsoTrace[ipISO]) )
85 fprintf(
ioQQQ,
" iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n",
91 sp->
st[0].Pop() =
dense.xIonDense[nelem][nelem-ipISO];
92 for(
long n=1; n < numlevels_local; n++ )
102 creation(numlevels_local),
103 error(numlevels_local),
104 work(numlevels_local),
105 Save_creation(numlevels_local);
107 ASSERT(
dense.xIonDense[nelem][nelem+1-ipISO] >= 0.f );
108 for( level=0; level < numlevels_local; level++ )
111 creation[level] = sp->
fb[level].RateCont2Level *
dense.xIonDense[nelem][nelem+1-ipISO];
114 double ctsource=0.0, ctsink=0.0, ctrec=0.0;
120 ctrec +=
atmdat.CharExcRecTotal[nelem];
121 ctsink +=
atmdat.CharExcIonTotal[nelem];
126 ctrec +=
atmdat.CharExcRecTotal[nelem];
127 ctsink +=
atmdat.CharExcIonTotal[nelem];
134 ctrec +=
atmdat.CharExcRecTo[nelem1][nelem][nelem-ipISO]*
iso_sp[ipISO][nelem1].st[0].Pop();
135 ctsink +=
atmdat.CharExcIonOf[nelem1][nelem][nelem-ipISO]*
dense.xIonDense[nelem1][1];
138 ctsource += ctrec*
dense.xIonDense[nelem][nelem+1-ipISO];
145 ctsink +=
atmdat.CharExcRecTotal[nelem];
146 ction +=
atmdat.CharExcIonTotal[nelem];
153 ctsink +=
atmdat.CharExcRecTo[nelem1][nelem][nelem-(ipISO+1)]*
iso_sp[ipISO1][nelem1].st[0].Pop();
154 ction +=
atmdat.CharExcIonOf[nelem1][nelem][nelem-(ipISO+1)]*
dense.xIonDense[nelem1][1];
157 ctsource += ction *
dense.xIonDense[nelem][nelem-(ipISO+1)];
162 z.
alloc(numlevels_local,numlevels_local);
164 vector<double> Boltzmann(numlevels_local-1);
165 for (
unsigned lev = 0; lev < Boltzmann.size(); ++lev)
166 Boltzmann[lev] = 1.0;
171 long SpinUsed[
NISO] = {2,1};
176 if(
trace.lgTrace && (nelem ==
trace.ipIsoTrace[ipISO]) )
178 fprintf(
ioQQQ,
" iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n",
186 for( level=0; level < numlevels_local; level++ )
190 z[level][level] += sp->
fb[level].RateLevel2Cont;
194 sp->
qTot2S = sp->
fb[level].ColIoniz;
200 double arg = (StElm[level].energy().K()-StElm[level-1].energy().K()) /
203 double bstep =
sexp( arg );
204 for ( ipLo = 0; ipLo < level; ++ipLo )
205 Boltzmann[ipLo] *= bstep;
209 for( ipLo=0; ipLo < level; ipLo++ )
224 coll_down *= sp->
ex[level][ipLo].ErrorFactor[
IPCOLLIS];
225 RadDecay *= sp->
ex[level][ipLo].ErrorFactor[
IPRAD];
226 pump *= sp->
ex[level][ipLo].ErrorFactor[
IPRAD];
229 double coll_up = coll_down *
230 (double)StElm[level].
g()/
231 (double)StElm[ipLo].
g()*
234 z[ipLo][ipLo] += coll_up + pump ;
235 z[ipLo][level] -= coll_up + pump ;
237 double pump_down = pump *
238 (double)StElm[ipLo].
g()/
239 (double)StElm[level].
g();
241 z[level][level] += RadDecay + pump_down + coll_down;
242 z[level][ipLo] -= RadDecay + pump_down + coll_down;
244 if( level == indexNmaxP )
246 HighestPColOut += coll_down;
248 if( ipLo == indexNmaxP )
250 HighestPColOut += coll_up;
254 if( (level == 1) && (ipLo==0) )
258 if( (ipLo == 1) && (ipISO==
ipH_LIKE || StElm[level].
S()==1) )
270 if( HighestPColOut < 1./sp->st[indexNmaxP].lifetime() )
272 iso_ctrl.lgCritDensLMix[ipISO] =
false;
278 for( vector<two_photon>::iterator tnu = sp->
TwoNu.begin(); tnu != sp->
TwoNu.end(); ++tnu )
286 z[tnu->ipHi][tnu->ipLo] -= tnu->AulTotal;
287 z[tnu->ipHi][tnu->ipHi] += tnu->AulTotal;
293 for( vector<two_photon>::iterator tnu = sp->
TwoNu.begin(); tnu != sp->
TwoNu.end(); ++tnu )
299 z[tnu->ipHi][tnu->ipLo] -= tnu->induc_dn;
300 z[tnu->ipHi][tnu->ipHi] += tnu->induc_dn;
303 z[tnu->ipLo][tnu->ipHi] -= tnu->induc_up;
304 z[tnu->ipLo][tnu->ipLo] += tnu->induc_up;
309 for(
long ion=
dense.IonLow[nelem]; ion<=
dense.IonHigh[nelem]; ++ion )
311 if( ion!=nelem-ipISO )
313 source +=
gv.GrainChTrRate[nelem][ion][nelem-ipISO] *
314 dense.xIonDense[nelem][ion];
315 sink +=
gv.GrainChTrRate[nelem][nelem-ipISO][ion];
317 source +=
mole.xMoleChTrRate[nelem][ion][nelem-ipISO] *
319 sink +=
mole.xMoleChTrRate[nelem][nelem-ipISO][ion] *
atmdat.lgCTOn;
326 sink +=
mole.sink[nelem][nelem-ipISO];
346 for( level=0; level < numlevels_local; level++ )
348 creation[level] +=
dynamics.StatesElem[nelem][nelem-ipISO][level];
355 for(
long ion_from=
dense.IonLow[nelem]; ion_from <
MIN2(
dense.IonHigh[nelem], nelem-ipISO ) ; ion_from++ )
357 if(
ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] >= 0. )
360 source +=
ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] *
dense.xIonDense[nelem][ion_from];
363 if( ion_from == nelem-1-ipISO )
371 for( level=0; level < numlevels_local; level++ )
373 z[level][level] +=
sink;
381 for ( level = 0; level < numlevels_local; level++ )
383 partfun += sp->
st[level].Boltzmann()*sp->
st[level].g();
386 for( level=0; level < numlevels_local; level++ )
388 creation[level] +=
source*
389 sp->
st[level].Boltzmann()*sp->
st[level].g();
390 z[level][level] +=
sink;
393 creation[0] += ctsource;
400 for( level=1; level < numlevels_local; level++ )
402 double RateUp , RateDown;
405 RateDown = RateUp * (double)sp->
st[
ipH1s].g() /
406 (double)sp->
st[level].g();
412 z[level][
ipH1s] -= RateDown;
415 z[
ipH1s][level] -= RateUp;
417 z[level][level] += RateDown;
429 for( ipLo=0; ipLo < numlevels_local; ipLo++ )
430 Save_creation[ipLo] = creation[ipLo];
432 if(
trace.lgTrace &&
trace.lgIsoTraceFull[ipISO] && (nelem ==
trace.ipIsoTrace[ipISO]) )
434 const long numlevels_print = numlevels_local;
435 fprintf(
ioQQQ,
" pop level others => (iso_level)\n" );
436 for(
long n=0; n < numlevels_print; n++ )
439 for(
long j=0; j < numlevels_print; j++ )
441 fprintf(
ioQQQ,
"\t%.9e", z[j][n] );
443 fprintf(
ioQQQ,
"\n" );
445 fprintf(
ioQQQ,
" recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
450 fprintf(
ioQQQ,
" recomb ");
451 for(
long n=0; n < numlevels_print; n++ )
453 fprintf(
ioQQQ,
"\t%.9e", creation[n] );
455 fprintf(
ioQQQ,
"\n" );
461 z.
data(),numlevels_local,&ipiv[0],&nerror);
464 &creation[0],numlevels_local,&nerror);
468 fprintf(
ioQQQ,
" iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
474 for( level=
ipH1s; level < numlevels_local; level++ )
476 double qn = 0., qx = 0.;
478 for(
long n=
ipH1s; n < numlevels_local; n++ )
480 double q = SaveZ[n][level]*creation[n];
496 error[level] = (error[level] - Save_creation[level])/qx;
508 for( level=
ipH1s; level < numlevels_local; level++ )
511 abserror = fabs( error[level]);
513 if( abserror > BigError )
522 if( BigError > FLT_EPSILON )
525 fprintf(
ioQQQ,
" PROBLEM" );
528 " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
529 "level was %li Search?%c \n",
544 for ( level = 0; level < numlevels_local; level++ )
546 partfun += sp->
st[level].Boltzmann()*sp->
st[level].g();
548 double scale =
dense.xIonDense[nelem][nelem-ipISO]/partfun;
549 for ( level = 0; level < numlevels_local; level++ )
551 creation[level] = sp->
st[level].Boltzmann()*sp->
st[level].g()*scale;
555 for( level = numlevels_local-1; level > 0; --level )
558 if( creation[level] < 0. )
562 if( lgNegPop &&
dense.lgSetIoniz[nelem] )
567 for( level = 1; level < numlevels_local; ++level )
568 creation[level] = 0.;
569 creation[0] =
dense.xIonDense[nelem][nelem-ipISO];
575 for( level=0; level < numlevels_local; level++ )
577 ASSERT( creation[level] >= 0. );
578 sp->
st[level].Pop() = creation[level];
580 if( sp->
st[level].Pop() <= 0 && !
conv.lgSearch )
583 "PROBLEM non-positive level pop for iso = %li, nelem = "
584 "%li = %s, level=%li val=%.3e nTotalIoniz %li\n",
589 sp->
st[level].Pop() ,
595 for( level=numlevels_local; level < sp->
numLevels_max; level++ )
597 sp->
st[level].Pop() = 0.;
606 double TotalPopExcited = 0.;
608 for( level=1; level < numlevels_local; level++ )
609 TotalPopExcited += sp->
st[level].Pop();
610 ASSERT( TotalPopExcited >= 0. );
611 double TotalPop = TotalPopExcited + sp->
st[0].Pop();
614 if(
dense.lgSetIoniz[nelem] )
616 if( !
fp_equal( TotalPop, (
double)
dense.xIonDense[nelem][nelem-ipISO] ) )
618 if( TotalPopExcited >=
dense.xIonDense[nelem][nelem-ipISO] )
620 for( level=0; level < numlevels_local; level++ )
621 sp->
st[level].Pop() *=
622 dense.xIonDense[nelem][nelem-ipISO] / TotalPop;
627 MAX2( 1e-30 *
dense.xIonDense[nelem][nelem-ipISO],
628 dense.xIonDense[nelem][nelem-ipISO] - TotalPopExcited );
638 if( lgNegPop ||
dense.xIonDense[nelem][nelem-ipISO] < 0. )
641 " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\
642 "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
647 dense.xIonDense[nelem][nelem+1-ipISO]/
SDIV(
dense.xIonDense[nelem][nelem-ipISO]),
652 fprintf(
ioQQQ,
" level pop are:\n" );
653 for( i=0; i < numlevels_local; i++ )
656 if( (i!=0) && !(i%10) ) fprintf(
ioQQQ,
"\n" );
658 fprintf(
ioQQQ,
"\n" );
664 for( ipHi=1; ipHi<numlevels_local; ++ipHi )
666 for( ipLo=0; ipLo<ipHi; ++ipLo )
673 sp->
st[ipLo].Pop() - sp->
st[ipHi].Pop()*
674 sp->
st[ipLo].g()/sp->
st[ipHi].g();