27 void iso_level(
const long int ipISO,
const long int nelem)
40 double HighestPColOut = 0.,
45 valarray<int32>
ipiv(numlevels_local);
48 creation(numlevels_local),
49 error(numlevels_local),
50 work(numlevels_local),
51 Save_creation(numlevels_local);
65 z.
alloc(numlevels_local,numlevels_local);
83 for( level=0; level < numlevels_local; level++ )
126 fprintf(
ioQQQ,
" iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n",
164 fprintf(
ioQQQ,
" iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s\n",
193 for(
long n=1; n < numlevels_local; n++ )
202 long SpinUsed[
NISO] = {2,1};
209 fprintf(
ioQQQ,
" iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n",
216 for( level=0; level < numlevels_local; level++ )
230 for( ipLo=0; ipLo < level; ipLo++ )
232 double coll_down = Trans[ipLo].Coll.ColUL *
collider;
235 (Trans[ipLo].Emis->Pesc +
236 Trans[ipLo].Emis->Pelec_esc +
237 Trans[ipLo].Emis->Pdest)*
250 double coll_up = coll_down *
251 (double)StElm[level].g/
252 (
double)StElm[ipLo].g*
255 z[ipLo][ipLo] += coll_up + pump ;
256 z[ipLo][level] = - ( coll_up + pump );
258 double pump_down = pump *
259 (double)StElm[ipLo].g/
260 (
double)StElm[level].g;
262 z[level][level] += RadDecay + pump_down + coll_down;
263 z[level][ipLo] = - (RadDecay + pump_down + coll_down);
265 if( level == indexNmaxP )
267 HighestPColOut += coll_down;
269 if( ipLo == indexNmaxP )
271 HighestPColOut += coll_up;
275 if( (level == 1) && (ipLo==0) )
279 if( (ipLo == 1) && (ipISO==
ipH_LIKE || StElm[level].S==1) )
291 if( HighestPColOut < 1./
StatesElem[ipISO][nelem][indexNmaxP].lifetime )
311 for(
long ion=0; ion<=nelem+1; ++ion )
313 if( ion!=nelem-ipISO )
343 sink +=
mole.
sink[nelem][nelem-ipISO];
348 for( level=0; level < numlevels_local; level++ )
350 z[level][level] += sink;
357 for( level=1; level < numlevels_local; level++ )
359 double RateUp , RateDown;
369 z[level][
ipH1s] -= RateDown;
372 z[
ipH1s][level] -= RateUp;
374 z[level][level] += RateDown;
387 for( ipLo=0; ipLo < numlevels_local; ipLo++ )
388 Save_creation[ipLo] = creation[ipLo];
392 fprintf(
ioQQQ,
" pop level others => (iso_level)\n" );
393 for(
long n=0; n < numlevels_local; n++ )
396 for(
long j=0; j < numlevels_local; j++ )
398 fprintf(
ioQQQ,
"\t%.9e", z[j][n] );
400 fprintf(
ioQQQ,
"\n" );
402 fprintf(
ioQQQ,
" recomb ");
403 for(
long n=0; n < numlevels_local; n++ )
405 fprintf(
ioQQQ,
"\t%.9e", creation[n] );
407 fprintf(
ioQQQ,
"\n" );
408 fprintf(
ioQQQ,
" recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
418 z.
data(),numlevels_local,&ipiv[0],&nerror);
421 &creation[0],numlevels_local,&nerror);
425 fprintf(
ioQQQ,
" iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
431 for( level=
ipH1s; level < numlevels_local; level++ )
435 for(
long n=
ipH1s; n < numlevels_local; n++ )
438 if( fabs(SaveZ[n][level] ) > BigSoln )
439 BigSoln = fabs(SaveZ[n][level]);
441 error[level] += SaveZ[n][level]*creation[n];
446 error[level] = (error[level] - Save_creation[level])/ BigSoln;
458 for( level=
ipH1s; level < numlevels_local; level++ )
461 abserror = fabs( error[level]);
463 if( abserror > BigError )
472 if( BigError > FLT_EPSILON )
475 fprintf(
ioQQQ,
" PROBLEM" );
478 " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
479 "level was %li Search?%c \n",
490 for( level=0; level < numlevels_local; level++ )
492 StatesElem[ipISO][nelem][level].Pop = creation[level];
495 if(
StatesElem[ipISO][nelem][level].Pop < 0. )
498 if(
StatesElem[ipISO][nelem][level].Pop <= 0 )
500 fprintf(
ioQQQ,
" non-positive level pop for iso = %li, nelem = %li = %s, level=%li val=%.3e\n",
508 if(
iso.
PopLTE[ipISO][nelem][level] > 0. )
522 for( level=numlevels_local; level <
iso.
numLevels_max[ipISO][nelem]; level++ )
533 sum_popn_ov_ion = 0.;
538 for( level=0; level < numlevels_local; level++ )
545 sum_popn_ov_ion +=
StatesElem[ipISO][nelem][level].Pop;
551 fprintf(
ioQQQ,
"DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.",
552 nelem+1, nelem-ipISO);
560 if( sum_popn_ov_ion >= 1e-30 )
591 " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\
592 "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
602 fprintf(
ioQQQ,
" level pop are:\n" );
603 for( i=0; i < numlevels_local; i++ )
606 if( (i!=0) && !(i%10) ) fprintf(
ioQQQ,
"\n" );
608 fprintf(
ioQQQ,
"\n" );
621 ratio = rateOutOf2TripS /
646 for( ipHi=1; ipHi<numlevels_local; ++ipHi )
648 for( ipLo=0; ipLo<ipHi; ++ipLo )
654 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc =