35STATIC void fill_array(
long int nelem,
long ion_range, valarray<double> &xmat, valarray<double> &
source, valarray<double> &auger,
double *abund_total );
38STATIC void combine_arrays( valarray<double> &xmat,
const valarray<double> &xmat1,
const valarray<double> &xmat2,
long ion_range1,
long ion_range2 );
43STATIC void HomogeneousSource(
long nelem,
long ion_low,
long ion_range, valarray<double> &xmat, valarray<double> &
source,
double abund_total );
47STATIC void PrintRates(
long nelem,
bool lgNegPop,
double abund_total, valarray<double> &auger,
bool lgPrintIt );
49void solveions(
double *ion,
double *rec,
double *snk,
double *src,
50 long int nlev,
long int nmax);
54# define MAT(M_,I_,J_) ((M_)[(I_)*(ion_range)+(J_)])
57# define MAT1(M_,I_,J_) ((M_)[(I_)*(ion_range1)+(J_)])
60# define MAT2(M_,I_,J_) ((M_)[(I_)*(ion_range2)+(J_)])
64 double abund_total = 0.0;
65 long ion_low =
dense.IonLow[nelem];
66 bool lgNegPop =
false;
73 long ion_range =
dense.IonHigh[nelem]-
dense.IonLow[nelem]+1;
74 valarray<double> xmat(ion_range*ion_range);
75 valarray<double>
source(ion_range);
76 valarray<double> auger(
LIMELM+1);
83 for (
long it=0; it<4; ++it)
88 if( (
dense.IonHigh[nelem] >= nelem - ipISO) &&
89 (
dense.IonLow[nelem] <= nelem - ipISO))
116 if (ion_range !=
dense.IonHigh[nelem]-
dense.IonLow[nelem]+1)
118 ion_range =
dense.IonHigh[nelem]-
dense.IonLow[nelem]+1;
119 xmat.resize(ion_range*ion_range);
131 if (thiserror > error)
146 if(
prt.lgPrtArry[nelem] || lgPrintIt )
147 PrintRates( nelem, lgNegPop, abund_total, auger, lgPrintIt );
154void ion_solver(
long int nelem1,
long int nelem2,
bool lgPrintIt)
156 bool lgHomogeneous =
true;
157 bool lgNegPop1 =
false;
158 bool lgNegPop2 =
false;
162 long ion_low1 =
dense.IonLow[nelem1];
163 long ion_low2 =
dense.IonLow[nelem2];
164 long ion_range1 =
dense.IonHigh[nelem1]-
dense.IonLow[nelem1]+1;
165 long ion_range2 =
dense.IonHigh[nelem2]-
dense.IonLow[nelem2]+1;
166 long ion_range = ion_range1 + ion_range2;
169 valarray<double> xmat(ion_range*ion_range);
170 valarray<double> xmat1(ion_range1*ion_range1);
171 valarray<double> xmat2(ion_range2*ion_range2);
172 valarray<double>
source(ion_range);
173 valarray<double> auger1(
LIMELM+1);
174 valarray<double> auger2(
LIMELM+1);
176#define ENABLE_SIMULTANEOUS_SOLUTION 0
183 else if(
lgTrivialSolution( nelem2, abund_total2 ) || !ENABLE_SIMULTANEOUS_SOLUTION )
193 fprintf(
ioQQQ,
"This routine is currently intended to do only O and H when both have significant neutral fractions.\n" );
194 fprintf(
ioQQQ,
"It should be generalized further for other cases. Exiting. Sorry.\n" );
198 ASSERT( nelem1 != nelem2 );
205 combine_arrays( xmat, xmat1, xmat2, ion_range1, ion_range2 );
208 for(
long i=0; i< ion_range; ++i )
212 for(
long i=0; i< ion_range1; ++i )
213 MAT( xmat, i, 0 ) = 1.;
215 for(
long i=ion_range1; i< ion_range; ++i )
216 MAT( xmat, i, ion_range1 ) = 1.;
219 source[ion_range1] =
dense.xIonDense[nelem2][0] +
dense.xIonDense[nelem2][1];
222 lgHomogeneous =
true;
230 ASSERT( lgNegPop2 ==
false );
232 if(
prt.lgPrtArry[nelem1] || lgPrintIt )
233 PrintRates( nelem1, lgNegPop1, abund_total1, auger1, lgPrintIt );
243 if(
dense.IonHigh[nelem] ==
dense.IonLow[nelem] )
246 dense.xIonDense[nelem][
dense.IonHigh[nelem]] = abund_total;
251 else if(
dense.lgSetIoniz[nelem] )
254 for(
long ion=0; ion<nelem+2; ++ion )
255 dense.xIonDense[nelem][ion] =
dense.SetIoniz[nelem][ion]*abund_total;
267 valarray<double> xmatsave(ion_range*ion_range);
268 valarray<double> sourcesave(ion_range);
269 valarray<int32> ipiv(ion_range);
274 for(
long i=0; i< ion_range; ++i )
276 sourcesave[i] =
source[i];
277 for(
long j=0; j< ion_range; ++j )
279 MAT( xmatsave, i, j ) =
MAT( xmat, i, j );
291 fprintf(
ioQQQ,
"Error for nelem %ld, active ion range %ld--%ld\n",
292 nelem,
dense.IonLow[nelem],
dense.IonHigh[nelem]);
293 fprintf(
ioQQQ,
"Initial ion abundances:");
294 for(
long j=0; j<nelem+2; ++j )
295 fprintf(
ioQQQ,
" %g",
dense.xIonDense[nelem][j]);
297 for(
long j=0; j<ion_range; ++j )
306 getrf_wrapper(ion_range, ion_range, &xmat[0], ion_range, &ipiv[0], &nerror);
310 " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, nConv %li IonLow %li IonHi %li\n",
316 fprintf(
ioQQQ,
" xmat follows\n");
317 for(
long i=0; i<ion_range; ++i )
319 for(
long j=0;j<ion_range;j++ )
321 fprintf(
ioQQQ,
"%e\t",
MAT(xmatsave,j,i));
325 fprintf(
ioQQQ,
"source follows\n");
326 for(
long i=0; i<ion_range;i++ )
328 fprintf(
ioQQQ,
"%e\t",sourcesave[i]);
333 getrs_wrapper(
'N', ion_range, 1, &xmat[0], ion_range, &ipiv[0], &
source[0], ion_range, &nerror);
336 fprintf(
ioQQQ,
" DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
344 enum {DEBUG_LOC=
false};
347 fprintf(
ioQQQ,
"debuggg ion_solver1 %ld\t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n",
352 ionbal.RateIonizTot(nelem,0) ,
353 ionbal.RateRecomTot[nelem][0]);
354 fprintf(
ioQQQ,
" Msrc %.3e %.3e\n",
mole.source[nelem][0],
mole.source[nelem][1]);
355 fprintf(
ioQQQ,
" Msnk %.3e %.3e\n",
mole.sink[nelem][0],
mole.sink[nelem][1]);
357 ionbal.RateIonizTot(nelem,0)/
ionbal.RateRecomTot[nelem][0]);
361 for(
long i=0; i<ion_range; i++ )
372#define THRESHOLD 0.25
374 bool lgDominant =
false;
379 lgDominant = lgDominant ||
385 lgDominant = lgDominant ||
393STATIC void combine_arrays( valarray<double> &xmat,
const valarray<double> &xmat1,
const valarray<double> &xmat2,
long ion_range1,
long ion_range2 )
397 long ion_range = ion_range1 + ion_range2;
399 for(
long i=0; i<ion_range1; i++ )
400 for(
long j=0; j<ion_range1; j++ )
401 MAT( xmat, i, j ) =
MAT1( xmat1, i, j );
403 for(
long i=0; i<ion_range2; i++ )
404 for(
long j=0; j<ion_range2; j++ )
405 MAT( xmat, i+ion_range1, j+ion_range1 ) =
MAT2( xmat2, i, j );
408 bool lgNoDice =
false;
409 for(
long i=0; i<ion_range1; i++ )
420 for(
long i=ion_range1; i<ion_range; i++ )
421 MAT( xmat, i, 0 ) = 1.0;
435 ASSERT( ion_range <= nelem + 2 );
437 ASSERT( ion_low <= nelem + 1 );
445 for(
long i=0; i<ion_range; i++) {
447 for(
long j=0; j<ion_range; j++) {
450 fprintf(
ioQQQ,
"%e\t",test);
455 fprintf(
ioQQQ,
" ion %li abundance %.3e\n",nelem,abund_total);
456 for(
long ion=
dense.IonLow[nelem]; ion <
dense.IonHigh[nelem]; ion++ )
458 if(
ionbal.RateRecomTot[nelem][ion] != 0 &&
source[ion-ion_low] != 0 )
459 fprintf(
ioQQQ,
" %li %.3e %.3e : %.3e\n",ion,
source[ion-ion_low],
461 ionbal.RateIonizTot(nelem,ion)/
ionbal.RateRecomTot[nelem][ion]);
463 fprintf(
ioQQQ,
" %li %.3e [One ratio infinity]\n",ion,
source[ion-ion_low]);
464 test +=
source[ion-ion_low];
492 for(
long i=0; i < ion_range; i++ )
494 long ion = i+ion_low;
505 " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
523 if( ion > nelem-
NISO && ion < nelem + 1 )
525 long int ipISO = nelem - ion;
527 for(
long level = 0; level <
iso_sp[ipISO][nelem].numLevels_max; level++ )
528 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
534 double renormnew = 1.;
539 for(
long i=0;i < ion_range; i++ )
546 renormnew = abund_total / dennew;
558 for(
long i=0;i < ion_range; i++ )
568 fprintf(
ioQQQ,
"PROBLEM impossible value of renorm \n");
572 for(
long i=0; i < ion_range; i++ )
574 long ion = i+ion_low;
578 fprintf(
ioQQQ,
"PROBLEM DISASTER Huge density in ion_solver, nelem %ld ion %ld source %e renormnew %e\n",
579 nelem, ion,
source[i], renormnew );
587 while(
dense.IonHigh[nelem] >
dense.IonLow[nelem] &&
588 dense.xIonDense[nelem][
dense.IonHigh[nelem]] < 1e-25*abund_total &&
589 dense.IonHigh[nelem] > 1)
593 dense.xIonDense[nelem][
dense.IonHigh[nelem]] = 0.;
596 --
dense.IonHigh[nelem];
608 (
dense.IonLow[nelem]==0 &&
dense.IonHigh[nelem]==0 ) ||
610 (
dense.IonLow[nelem]==nelem+1 &&
dense.IonHigh[nelem]==nelem+1 ) );
624 ionbal.elecsnk[nelem] = 0.;
625 ionbal.elecsrc[nelem] = 0.;
627 double abund_total = 0.;
628 for(
long ion=
dense.IonLow[nelem]; ion<=
dense.IonHigh[nelem]; ++ion )
630 abund_total +=
dense.xIonDense[nelem][ion];
637 if (fabs(tot2-tot1) >
conv.GasPhaseAbundErrorAllowed*tot1)
638 fprintf(
ioQQQ,
"%ld %13.8g %13.8g %13.8g %13.8g\n",nelem,tot1,tot2,abund_total,tot2/tot1 - 1.0);
648STATIC void fill_array(
long int nelem,
long ion_range, valarray<double> &xmat, valarray<double> &
source, valarray<double> &auger,
double *abund_total )
655 valarray<double>
sink(ion_range);
656 valarray<int32> ipiv(ion_range);
682 ASSERT( ion_range <= nelem+2 );
684 ion_low =
dense.IonLow[nelem];
685 for(
long i=0; i<ion_range;i++ )
694 for(
long ion_from=0; ion_from <= limit; ion_from++ )
696 for(
long ion_to=0; ion_to < nelem+2; ion_to++ )
698 ionbal.RateIoniz[nelem][ion_from][ion_to] = 0.;
705 for(
long i=0; i<
LIMELM+1; ++i )
711 for(
long i=0; i< ion_range; ++i )
713 for(
long j=0; j< ion_range; ++j )
715 MAT( xmat, i, j ) = 0.;
722 enum {DEBUG_LOC=
false};
726 dense.IonLow[nelem] = 0;
727 dense.IonHigh[nelem] = 3;
730 for(
long ion=
dense.IonLow[nelem]; ion <= limit; ion++ )
735 ionbal.RateRecomTot[nelem][ion] = 100.;
736 for(
long ns=0; ns <
Heavy.nsShells[nelem][ion]; ns++ )
739 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
745 bool lgGrainsOn =
gv.lgDustOn() &&
ionbal.lgGrainIonRecom &&
gv.lgGrainPhysicsOn;
751 for(
long ion_to=
dense.IonLow[nelem]; ion_to <=
dense.IonHigh[nelem]; ion_to++ )
753 for(
long ion_from=
dense.IonLow[nelem]; ion_from <=
dense.IonHigh[nelem]; ++ion_from )
756 if( ion_to != ion_from )
759 rateone =
mole.xMoleChTrRate[nelem][ion_from][ion_to] *
atmdat.lgCTOn;
760 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
761 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
763 rateone =
gv.GrainChTrRate[nelem][ion_from][ion_to]*lgGrainsOn;
764 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
765 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
770 for(
long ion=
dense.IonLow[nelem]; ion <= limit; ion++ )
773 ionbal.RateIoniz[nelem][ion][ion+1] +=
774 ionbal.CollIonRate_Ground[nelem][ion][0] +
777 ionbal.UTA_ionize_rate[nelem][ion];
780 for(
long ns=0; ns <
Heavy.nsShells[nelem][ion]; ns++ )
789 if( ion+1-ion_low < ion_range )
803 IonProduced =
MIN2(ion+nej,
dense.IonHigh[nelem]);
804 rateone =
ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
808 ionbal.RateIoniz[nelem][ion][IonProduced] += rateone;
813 auger[IonProduced-1] += rateone;
819 for(
long ion=
dense.IonLow[nelem]; ion <
dense.IonHigh[nelem]; ion++ )
823 long ipISO=nelem-ion;
827 ction =
atmdat.CharExcIonTotal[nelem] *
iso_sp[ipISO][nelem].st[0].Pop() /
SDIV(
dense.xIonDense[nelem][nelem-ipISO]);
833 ction +=
atmdat.CharExcIonOf[nelem1][nelem][ion]*
dense.xIonDense[nelem1][1];
837 MAT( xmat, ion-ion_low, ion-ion_low ) -= ction;
838 MAT( xmat, ion-ion_low, ion+1-ion_low ) += ction;
841 for(
long ion=
dense.IonLow[nelem]; ion <
dense.IonHigh[nelem]; ion++ )
845 long ipISO=nelem-ion;
848 ctrec =
atmdat.CharExcRecTotal[nelem];
850 ctrec =
atmdat.CharExcRecTotal[nelem];
858 atmdat.CharExcRecTo[nelem1][nelem][ion]*
iso_sp[ipISO][nelem1].st[0].Pop();
862 MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -= ctrec;
863 MAT( xmat, ion+1-ion_low, ion-ion_low ) += ctrec;
868 for(
long ion_from=
dense.IonLow[nelem]; ion_from <
dense.IonHigh[nelem]; ion_from++ )
870 for(
long ion_to=ion_from+1; ion_to <=
dense.IonHigh[nelem]; ion_to++ )
872 ionbal.elecsrc[nelem] +=
ionbal.RateIoniz[nelem][ion_from][ion_to]*
dense.xIonDense[nelem][ion_from]*
875 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -=
ionbal.RateIoniz[nelem][ion_from][ion_to];
876 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) +=
ionbal.RateIoniz[nelem][ion_from][ion_to];
880 for(
long ion=
dense.IonLow[nelem]; ion<
dense.IonHigh[nelem]; ion++ )
883 MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -=
ionbal.RateRecomTot[nelem][ion];
884 MAT( xmat, ion+1-ion_low, ion-ion_low ) +=
ionbal.RateRecomTot[nelem][ion];
885 ionbal.elecsnk[nelem] +=
ionbal.RateRecomTot[nelem][ion]*
dense.xIonDense[nelem][ion+1];
888 for(
long i=0; i<ion_range;i++ )
890 long ion = i+ion_low;
895 MAT( xmat, i, i ) -=
mole.sink[nelem][ion];
903 bool lgHomogeneous =
true;
912 totsrc = totscl = maxdiag = 0.;
913 for(
long i=0; i<ion_range;i++ )
915 long ion = i + ion_low;
922 double diag = -(
MAT( xmat, i, i )+
mole.sink[nelem][ion]);
925 totscl += diag*
dense.xIonDense[nelem][ion];
930 const double CONDITION_SCALE = 1e8;
931 double conserve_scale = maxdiag/CONDITION_SCALE;
932 ASSERT( conserve_scale > 0.0 );
935 if( totsrc > 1e-10*totscl )
936 lgHomogeneous =
false;
944 for(
long i=0; i<ion_range;i++ )
946 long ion = i+ion_low;
952 lgHomogeneous =
false;
962 if( !lgHomogeneous && ion_range==2 )
965 double a =
MAT( xmat, 0, 0 ),
966 b =
MAT( xmat, 1, 0 ) ,
967 c =
MAT( xmat, 0, 1 ) ,
968 d =
MAT( xmat, 1, 1 );
971 double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
972 if( fabs(ratio1-ratio2) <= DBL_EPSILON*1e4*
MAX2(fratio1,fratio2) )
974 lgHomogeneous =
true;
980 fabs(
dense.xMolecules(nelem)) <DBL_EPSILON*100.*
dense.gas_phase[nelem] )
982 lgHomogeneous =
true;
989 if( 1 || lgHomogeneous )
991 double rate_ioniz=1., rate_recomb ;
994 for(
long i=0; i<ion_range-1;i++)
996 long ion = i+ion_low;
997 rate_ioniz *=
ionbal.RateIonizTot(nelem,ion);
998 rate_recomb =
ionbal.RateRecomTot[nelem][ion];
1000 if( rate_ioniz == 0)
1003 if( rate_recomb <= 0.)
1006 rate_ioniz /= rate_recomb;
1007 if( rate_ioniz > 1.)
1016 for(
long i=0; i<ion_range;i++ )
1018 MAT(xmat,i,jmax) -= conserve_scale;
1020 source[jmax] -= abund_total*conserve_scale;
1027 fprintf(
ioQQQ,
" %.3e %.3e\n",
ionbal.RateIonizTot(nelem,0),
ionbal.RateIonizTot(nelem,1) );
1028 fprintf(
ioQQQ,
" %.3e %.3e\n",
ionbal.RateRecomTot[nelem][0],
ionbal.RateRecomTot[nelem][1]);
1034 enum {DEBUG_LOC=
false};
1035 if( DEBUG_LOC &&
nzone==1 && lgPrintIt )
1038 " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
1039 nelem , ion_range,limit ,
conv.nTotalIoniz );
1041 fprintf(
ioQQQ ,
"Homogeneous \n");
1042 for(
long i=0; i<ion_range; ++i )
1044 for(
long j=0;j<ion_range;j++ )
1046 fprintf(
ioQQQ,
"%e\t",
MAT(xmat,i,j));
1048 fprintf(
ioQQQ,
"\n");
1050 fprintf(
ioQQQ,
"source follows\n");
1051 for(
long i=0; i<ion_range;i++ )
1055 fprintf(
ioQQQ,
"\n");
1056 fprintf(
ioQQQ,
"totsrc/totscl %g %g\n",totsrc,totscl);
1063STATIC void PrintRates(
long nelem,
bool lgNegPop,
double abund_total, valarray<double> &auger,
bool lgPrintIt )
1072 fprintf(
ioQQQ,
" PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n",
1075 fprintf(
ioQQQ,
" Populations were" );
1076 for( ion=1; ion <=
dense.IonHigh[nelem]+1; ion++ )
1078 fprintf(
ioQQQ,
"%9.1e",
dense.xIonDense[nelem][ion-1] );
1080 fprintf(
ioQQQ,
"\n" );
1082 fprintf(
ioQQQ,
" destroy vector =" );
1083 for( ion=1; ion <=
dense.IonHigh[nelem]; ion++ )
1085 fprintf(
ioQQQ,
"%9.1e",
ionbal.RateIonizTot(nelem,ion-1) );
1087 fprintf(
ioQQQ,
"\n" );
1089 fprintf(
ioQQQ,
" CTHeavy vector =" );
1090 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1094 fprintf(
ioQQQ,
"\n" );
1096 fprintf(
ioQQQ,
" CharExcIonOf[ipHYDROGEN] vtr=" );
1097 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1101 fprintf(
ioQQQ,
"\n" );
1103 fprintf(
ioQQQ,
" CollidRate vtr=" );
1104 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1106 fprintf(
ioQQQ,
"%9.1e",
ionbal.CollIonRate_Ground[nelem][ion][0] );
1108 fprintf(
ioQQQ,
"\n" );
1111 fprintf(
ioQQQ,
" photo rates per subshell, ion\n" );
1112 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1114 fprintf(
ioQQQ,
"%3ld", ion );
1115 for(
long ns=0; ns <
Heavy.nsShells[nelem][ion]; ns++ )
1117 fprintf(
ioQQQ,
"%9.1e",
ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
1119 fprintf(
ioQQQ,
"\n" );
1123 fprintf(
ioQQQ,
" create vector =" );
1124 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1126 fprintf(
ioQQQ,
"%9.1e",
ionbal.RateRecomTot[nelem][ion] );
1128 fprintf(
ioQQQ,
"\n" );
1133 if( lgPrintIt ||
prt.lgPrtArry[nelem] || lgNegPop )
1135 const char* format =
" %9.2e";
1138 "\n %s ion_solver ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n",
1144 dense.gas_phase[nelem],
1146 dense.xMolecules(nelem) );
1149 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1151 fprintf(
ioQQQ, format,
ionbal.RateIonizTot(nelem,ion) );
1153 fprintf(
ioQQQ,
"\n" );
1157 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1159 fprintf(
ioQQQ, format,
ionbal.RateRecomTot[nelem][ion] );
1161 fprintf(
ioQQQ,
"\n" );
1165 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1167 double ColIoniz =
ionbal.CollIonRate_Ground[nelem][ion][0];
1168 if( ion > nelem -
NISO )
1170 long ipISO = nelem-ion;
1172 if(
dense.xIonDense[nelem][nelem-ipISO] > 0. )
1174 ColIoniz *=
iso_sp[ipISO][nelem].st[0].Pop();
1175 for(
long ipLevel=1; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
1177 ColIoniz +=
iso_sp[ipISO][nelem].fb[ipLevel].ColIoniz *
dense.EdenHCorr *
iso_sp[ipISO][nelem].st[ipLevel].Pop();
1179 ColIoniz /=
dense.xIonDense[nelem][nelem-ipISO];
1182 fprintf(
ioQQQ, format, ColIoniz );
1184 fprintf(
ioQQQ,
"\n" );
1188 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1190 fprintf(
ioQQQ, format,
ionbal.UTA_ionize_rate[nelem][ion] );
1192 fprintf(
ioQQQ,
"\n" );
1196 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1198 double PhotIoniz = 0.;
1199 for(
long ipShell = 0; ipShell <
Heavy.nsShells[nelem][ion]; ipShell++ )
1200 PhotIoniz +=
ionbal.PhotoRate_Shell[nelem][ion][ipShell][0];
1203 if( ion > nelem -
NISO )
1205 long ipISO = nelem-ion;
1207 if(
dense.xIonDense[nelem][nelem-ipISO]>0 )
1209 PhotIoniz *=
iso_sp[ipISO][nelem].st[0].Pop();
1210 for(
long ipLevel=1; ipLevel <
iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
1212 PhotIoniz +=
iso_sp[ipISO][nelem].fb[ipLevel].gamnc *
iso_sp[ipISO][nelem].st[ipLevel].Pop();
1214 PhotIoniz /=
dense.xIonDense[nelem][nelem-ipISO];
1217 fprintf(
ioQQQ, format, PhotIoniz );
1219 fprintf(
ioQQQ,
"\n" );
1223 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1231 fprintf(
ioQQQ,
"\n" );
1235 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1243 fprintf(
ioQQQ,
"\n" );
1247 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1249 fprintf(
ioQQQ, format,
1252 fprintf(
ioQQQ,
"\n" );
1256 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1258 fprintf(
ioQQQ, format,
gv.GrainChTrRate[nelem][ion][ion+1] );
1260 fprintf(
ioQQQ,
"\n" );
1264 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1266 fprintf(
ioQQQ, format,
mole.xMoleChTrRate[nelem][ion][ion+1] );
1268 fprintf(
ioQQQ,
"\n" );
1272 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1276 long ipISO=nelem-ion;
1280 sum =
atmdat.CharExcIonTotal[nelem] *
iso_sp[ipISO][nelem].st[0].Pop() /
SDIV(
dense.xIonDense[nelem][nelem-ipISO]);
1286 sum +=
atmdat.CharExcIonOf[nelem1][nelem][ion] *
dense.xIonDense[nelem1][1];
1289 fprintf(
ioQQQ, format, sum );
1291 fprintf(
ioQQQ,
"\n" );
1295 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1299 fprintf(
ioQQQ,
"\n" );
1303 for( ion=0; ion <
min(nelem-1,
dense.IonHigh[nelem]); ion++ )
1305 fprintf(
ioQQQ, format,
dense.eden*
ionbal.DR_Badnell_rate_coef[nelem][ion] );
1307 fprintf(
ioQQQ,
"\n" );
1311 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1315 fprintf(
ioQQQ,
"\n" );
1320 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1322 fprintf(
ioQQQ, format,
gv.GrainChTrRate[nelem][ion+1][ion] );
1324 fprintf(
ioQQQ,
"\n" );
1328 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1330 fprintf(
ioQQQ, format,
mole.xMoleChTrRate[nelem][ion+1][ion] );
1332 fprintf(
ioQQQ,
"\n" );
1336 for( ion=0; ion <
dense.IonHigh[nelem]; ion++ )
1342 sum =
atmdat.CharExcRecTotal[nelem];
1346 sum =
atmdat.CharExcRecTotal[nelem];
1354 sum +=
atmdat.CharExcRecTo[nelem1][nelem][ion] *
iso_sp[ipISO][nelem1].st[0].Pop();
1357 fprintf(
ioQQQ, format, sum );
1359 fprintf(
ioQQQ,
"\n" );
1363 fprintf(
ioQQQ,
" %s molecule src",
1365 for( ion=0; ion <=
dense.IonHigh[nelem]; ion++ )
1369 fprintf(
ioQQQ,
"\n" );
1372 fprintf(
ioQQQ,
" %s molecule snk",
1374 for( ion=0; ion <=
dense.IonHigh[nelem]; ion++ )
1376 fprintf(
ioQQQ,format,
mole.sink[nelem][ion] );
1378 fprintf(
ioQQQ,
"\n" );
1383 fprintf(
ioQQQ,
" %s dynamics src",
1385 for( ion=0; ion <=
dense.IonHigh[nelem]; ion++ )
1389 fprintf(
ioQQQ,
"\n" );
1392 fprintf(
ioQQQ,
" %s dynamics snk",
1394 for( ion=0; ion <=
dense.IonHigh[nelem]; ion++ )
1398 fprintf(
ioQQQ,
"\n" );
1404 for( ion=0; ion <=
dense.IonHigh[nelem]; ion++ )
1406 fprintf(
ioQQQ, format,
dense.xIonDense[nelem][ion] );
1408 fprintf(
ioQQQ,
"\n" );
1448void solveions(
double *ion,
double *rec,
double *snk,
double *src,
1449 long int nlev,
long int nmax)
1460 for(i=nmax;i<nlev-1;i++)
1461 src[i+1] = src[i]*ion[i]/rec[i];
1462 for(i=nmax-1;i>=0;i--)
1463 src[i] = src[i+1]*rec[i]/ion[i];
1469 for(i=0;i<nlev-1;i++)
1474 fprintf(
ioQQQ,
"Ionization solver error\n");
1479 src[i+1] += ion[i]*src[i];
1480 snk[i] = bet*rec[i];
1481 kap = kap*snk[i]+snk[i+1];
1486 fprintf(
ioQQQ,
"Ionization solver error\n");
1491 for(i=nlev-2;i>=0;i--)
1493 src[i] += snk[i]*src[i+1];
1521 fprintf(
ioQQQ,
" ion_wrapper returns; %s fracs = ",
elementnames.chElementSym[nelem] );
1522 for(
long ion = 0; ion<=nelem+1; ion++ )
1523 fprintf(
ioQQQ,
"%10.3e ",
dense.xIonDense[nelem][ion]/
dense.gas_phase[nelem] );
1524 fprintf(
ioQQQ,
"\n" );
bool fp_equal(sys_float x, sys_float y, int n=3)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
long nelec_eject(long n, long i, long ns) const
realnum elec_eject_frac(long n, long i, long ns, long ne) const
bool lgElemsConserved(void)
t_elementnames elementnames
void IonNelem(bool lgPrintIt, long int nelem)
STATIC void store_new_densities(long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop)
STATIC double get_total_abundance_ions(long int nelem)
STATIC bool lgTrivialSolution(long nelem, double abund_total)
void ion_solver(long int nelem, bool lgPrintIt)
STATIC void HomogeneousSource(long nelem, long ion_low, long ion_range, valarray< double > &xmat, valarray< double > &source, double abund_total)
STATIC void PrintRates(long nelem, bool lgNegPop, double abund_total, valarray< double > &auger, bool lgPrintIt)
STATIC void fill_array(long int nelem, long ion_range, valarray< double > &xmat, valarray< double > &source, valarray< double > &auger, double *abund_total)
void ion_wrapper(long nelem)
STATIC void find_solution(long nelem, long ion_range, valarray< double > &xmat, valarray< double > &source)
bool lgOH_ChargeTransferDominant(void)
void solveions(double *ion, double *rec, double *snk, double *src, long int nlev, long int nmax)
t_iso_sp iso_sp[NISO][LIMELM]
void iso_departure_coefficients(long ipISO, long nelem)
void iso_renorm(long nelem, long ipISO, double &renorm)
void iso_satellite_update(long nelem)
void iso_set_ion_rates(long ipISO, long nelem)
void iso_charge_transfer_update(long nelem)
void iso_solve(long ipISO, long nelem, double &maxerr)
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
t_secondaries secondaries
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)