25 void solveions(
double *ion,
double *rec,
double *snk,
double *src,
26 long int nlev,
long int nmax);
43 double dennew, rateone;
45 static bool lgMustMalloc=
true;
46 static double *sink, *source , *sourcesave;
49 long ion_low, ion_range, i, j, ion_to , ion_from;
50 static double sum_dense;
54 double abund_total, renormnew;
55 bool lgHomogeneous =
true;
56 static double *xmat , *xmatsave;
133 abund_total = renorm * sum;
137 if( abund_total < 0. )
149 "PROBLEM ion_solver - neg net atomic abundance zero for nelem= %li, rel val= %.2e conv.nTotalIoniz=%li, fixed\n",
155 abund_total = -abund_total/2.;
176 for( ion=0; ion<nelem+2; ++ion )
203 ASSERT( ion_range <= nelem+2 );
213 lgMustMalloc =
false;
218 xmat=(
double*)
MALLOC( (
sizeof(
double)*(
unsigned)(ncell*ncell) ));
219 xmatsave=(
double*)
MALLOC( (
sizeof(
double)*(unsigned)(ncell*ncell) ));
220 sink=(
double*)
MALLOC( (
sizeof(
double)*(unsigned)ncell ));
221 source=(
double*)
MALLOC( (
sizeof(
double)*(unsigned)ncell ));
222 sourcesave=(
double*)
MALLOC( (
sizeof(
double)*(unsigned)ncell ));
223 ipiv=(int32*)
MALLOC( (
sizeof(int32)*(unsigned)ncell ));
227 auger=(
double*)
MALLOC( (
sizeof(
double)*(unsigned)(
LIMELM+1) ));
230 for( i=0; i<ion_range;i++ )
237 # define MAT(M_,I_,J_) (*((M_)+(I_)*(ion_range)+(J_)))
243 for( ion=0; ion <= limit; ion++ )
251 for( i=0; i<
LIMELM+1; ++i )
257 for( i=0; i< ion_range; ++i )
259 for( j=0; j< ion_range; ++j )
261 MAT( xmat, i, j ) = 0.;
268 enum {DEBUG_LOC=
false};
277 for( ion=
dense.
IonLow[nelem]; ion <= limit; ion++ )
296 for( ion_to=
dense.
IonLow[nelem]; ion_to <= limit; ion_to++ )
301 if( ion_to != ion_from )
307 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
308 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
333 for( ion_to=low; ion_to <=
dense.
IonHigh[nelem]; ion_to++ )
338 if( ion_to != ion_from )
343 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
344 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
350 for( ion=
dense.
IonLow[nelem]; ion <= limit; ion++ )
360 if( ion+1-ion_low < ion_range )
363 MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
364 MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
368 if( ion-1-ion_low >= 0 )
388 if( ion+1-ion_low < ion_range )
413 MAT( xmat, ion-ion_low, IonProduced-ion_low ) += rateone;
418 auger[IonProduced-1] += rateone;
428 if( ion+1-ion_low < ion_range )
430 MAT( xmat, ion-ion_low, ion-ion_low ) -= rateone;
431 MAT( xmat, ion-ion_low, ion+1-ion_low ) += rateone;
440 j =
MAX2( ion_low , j );
443 ASSERT( ion>=0 && ion<nelem+2 );
445 if( ion+1-ion_low < ion_range )
452 if( ion-1-ion_low >= 0 )
469 for( i=0; i<ion_range;i++ )
487 lgHomogeneous =
false;
494 for( i=0; i<ion_range;i++ )
504 lgHomogeneous =
false;
507 for( i=0; i< ion_range; ++i )
509 for( j=0; j< ion_range; ++j )
511 MAT( xmatsave, i, j ) =
MAT( xmat, i, j );
513 sourcesave[i] = source[i];
523 if( !lgHomogeneous && ion_range==2 )
526 double a =
MAT( xmatsave, 0, 0 ),
527 b =
MAT( xmatsave, 1, 0 ) ,
528 c =
MAT( xmatsave, 0, 1 ) ,
529 d =
MAT( xmatsave, 1, 1 );
532 double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
533 if( fabs(ratio1-ratio2)/
MAX2(fratio1,fratio2) <DBL_EPSILON*1e4 )
536 lgHomogeneous =
true;
544 lgHomogeneous =
true;
553 double rate_ioniz=1., rate_recomb ;
556 for( i=0; i<ion_range-1;i++)
565 if( rate_recomb <= 0.)
568 rate_ioniz /= rate_recomb;
578 for( i=0; i<ion_range;i++ )
580 MAT(xmat,i,jmax) = 1.;
582 source[jmax] = abund_total;
585 for( i=0; i< ion_range; ++i )
587 for( j=0; j< ion_range; ++j )
589 MAT( xmatsave, i, j ) =
MAT( xmat, i, j );
591 sourcesave[i] = source[i];
605 enum {DEBUG_LOC=
false};
607 if( DEBUG_LOC &&
nzone==1 && lgPrintIt )
610 " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
613 fprintf(
ioQQQ ,
"Homogeneous \n");
614 for( i=0; i<ion_range; ++i )
616 for( j=0;j<ion_range;j++ )
618 fprintf(
ioQQQ,
"%e\t",
MAT(xmatsave,i,j));
622 fprintf(
ioQQQ,
"source follows\n");
623 for( i=0; i<ion_range;i++ )
625 fprintf(
ioQQQ,
"%e\t",source[i]);
634 bool lgLapack=false , lgTriOptimized=
true;
646 d = (
double *)
MALLOC((
unsigned)ion_range*
sizeof(double));
647 du = (
double *)
MALLOC((
unsigned)(ion_range-1)*
sizeof(
double));
648 dl = (
double *)
MALLOC((
unsigned)(ion_range-1)*
sizeof(
double));
650 for(i=0;i<ion_range-1;i++)
652 du[i] =
MAT(xmat,i+1,i);
653 d[i] =
MAT(xmat,i,i);
654 dl[i] =
MAT(xmat,i,i+1);
656 d[i] =
MAT(xmat,i,i);
662 lgBad = dgtsv_wrapper(&i1, &i2, dl, d, du, source, &i2);
665 fprintf(
ioQQQ,
" dgtsz error.\n");
667 free(dl);free(du);free(d);
669 else if(lgTriOptimized)
676 for(i=0;i<ion_range;i++)
678 source[i] = sink[i] = 0.;
682 for(i=0;i<ion_range;i++)
691 for(i=0;i<ion_range;i++)
700 sink,source,ion_range,jmax);
708 getrf_wrapper(ion_range, ion_range, xmat, ion_range, ipiv, &nerror);
712 " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, limit=%li, nConv %li xmat follows\n",
718 for( i=0; i<ion_range; ++i )
720 for( j=0;j<ion_range;j++ )
722 fprintf(
ioQQQ,
"%e\t",
MAT(xmatsave,j,i));
726 fprintf(
ioQQQ,
"source follows\n");
727 for( i=0; i<ion_range;i++ )
729 fprintf(
ioQQQ,
"%e\t",source[i]);
734 getrs_wrapper(
'N', ion_range, 1, xmat, ion_range, ipiv, source, ion_range, &nerror);
737 fprintf(
ioQQQ,
" DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
746 enum {DEBUG_LOC=
false};
750 fprintf(
ioQQQ,
"debuggg ion_solver1 \t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n",
758 fprintf(
ioQQQ,
" Poprat %.3e nomol %.3e\n",source[1]/source[0],
763 for( i=0; i<ion_range; i++ )
774 for(i=0; i<ion_range; i++) {
776 for(j=0; j<ion_range; j++) {
777 test = test+source[j]*
MAT(xmatsave,j,i);
779 fprintf(
ioQQQ,
"%e\t",test);
784 fprintf(
ioQQQ,
" ion %li abundance %.3e\n",nelem,abund_total);
788 fprintf(
ioQQQ,
" %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low],
789 source[ion-ion_low+1]/source[ion-ion_low],
792 fprintf(
ioQQQ,
" %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]);
793 test += source[ion-ion_low];
795 test += source[ion-ion_low];
821 for( i=1;i < ion_range; i++ )
844 for( i=0;i < ion_range; i++ )
853 renormnew = sum_dense / dennew;
868 fprintf(
ioQQQ,
"PROBLEM impossible value of renorm \n");
876 for( i=0; i < ion_range; i++ )
888 if( source[i]<-2e-9 )
890 " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
900 if( ion == nelem+1-
NISO )
902 long int ipISO = nelem - ion;
932 fprintf(
ioQQQ,
" PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n",
935 fprintf(
ioQQQ,
" Populations were" );
940 fprintf(
ioQQQ,
"\n" );
942 fprintf(
ioQQQ,
" destroy vector =" );
947 fprintf(
ioQQQ,
"\n" );
952 fprintf(
ioQQQ,
" CTHeavy vector =" );
957 fprintf(
ioQQQ,
"\n" );
959 fprintf(
ioQQQ,
" HCharExcIonOf vtr=" );
964 fprintf(
ioQQQ,
"\n" );
966 fprintf(
ioQQQ,
" CollidRate vtr=" );
971 fprintf(
ioQQQ,
"\n" );
974 fprintf(
ioQQQ,
" photo rates per subshell, ion\n" );
977 fprintf(
ioQQQ,
"%3ld", ion );
982 fprintf(
ioQQQ,
"\n" );
987 fprintf(
ioQQQ,
" create vector =" );
992 fprintf(
ioQQQ,
"\n" );
1005 "\n %s ion_solver DEBUG ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n",
1020 fprintf(
ioQQQ,
"\n" );
1023 fprintf(
ioQQQ,
" %s mole sink ",
1029 fprintf(
ioQQQ,
"\n" );
1034 fprintf(
ioQQQ,
" %s dynm sink ",
1040 fprintf(
ioQQQ,
"\n" );
1049 fprintf(
ioQQQ,
"\n" );
1052 fprintf(
ioQQQ,
" %s mole source ",
1058 fprintf(
ioQQQ,
"\n" );
1063 fprintf(
ioQQQ,
" %s dynm sour ",
1069 fprintf(
ioQQQ,
"\n" );
1078 fprintf(
ioQQQ,
"\n" );
1086 fprintf(
ioQQQ,
"\n" );
1092 fprintf(
ioQQQ,
" %9.2e",
1095 fprintf(
ioQQQ,
"\n" );
1101 fprintf(
ioQQQ,
" %9.2e",
1104 fprintf(
ioQQQ,
"\n" );
1110 fprintf(
ioQQQ,
" %9.2e",
1113 fprintf(
ioQQQ,
"\n" );
1121 fprintf(
ioQQQ,
"\n" );
1130 fprintf(
ioQQQ,
"\n" );
1148 fprintf(
ioQQQ,
" %9.2e", sum );
1150 fprintf(
ioQQQ,
"\n" );
1158 fprintf(
ioQQQ,
"\n" );
1166 fprintf(
ioQQQ,
"\n" );
1174 fprintf(
ioQQQ,
"\n" );
1183 fprintf(
ioQQQ,
"\n" );
1192 fprintf(
ioQQQ,
"\n" );
1215 fprintf(
ioQQQ,
" %9.2e", sum );
1217 fprintf(
ioQQQ,
"\n" );
1225 fprintf(
ioQQQ,
"\n" );
1267 void solveions(
double *ion,
double *rec,
double *snk,
double *src,
1268 long int nlev,
long int nmax)
1279 for(i=nmax;i<nlev-1;i++)
1280 src[i+1] = src[i]*ion[i]/rec[i];
1281 for(i=nmax-1;i>=0;i--)
1282 src[i] = src[i+1]*rec[i]/ion[i];
1288 for(i=0;i<nlev-1;i++)
1293 fprintf(
ioQQQ,
"Ionization solver error\n");
1298 src[i+1] += ion[i]*src[i];
1299 snk[i] = bet*rec[i];
1300 kap = kap*snk[i]+snk[i+1];
1305 fprintf(
ioQQQ,
"Ionization solver error\n");
1310 for(i=nlev-2;i>=0;i--)
1312 src[i] += snk[i]*src[i+1];