37 static bool lgFirstPass =
true;
38 static long maxNumLevels = 1;
39 double totalHeating = 0.;
46 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
49 g = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
50 ex = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
51 pops = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
52 depart = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
53 source = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
54 sink = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
56 AulEscp = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
57 col_str = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
58 AulDest = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
59 AulPump = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
60 CollRate = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
62 for(
long j=0; j< maxNumLevels; j++ )
64 AulEscp[j] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
65 col_str[j] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
66 AulDest[j] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
67 AulPump[j] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
68 CollRate[j] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
75 memset(
g, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
76 memset(
ex, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
77 memset(
pops, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
78 memset(
depart, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
79 memset(
source, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
80 memset(
sink, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
81 for(
long j=0; j< maxNumLevels; j++ )
83 memset(
AulEscp[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
84 memset(
col_str[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
85 memset(
AulDest[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
86 memset(
AulPump[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
87 memset(
CollRate[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
91 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
94 double cooltl, coolder;
96 bool lgZeroPop, lgDeBug =
false;
114 if( !
conv.nTotalIoniz )
115 fprintf(
ioQQQ,
" PROBLEM dBase_solve did not find molecular species %li\n",ipSpecies);
129 if(
conv.nTotalIoniz == 0)
132 bool lgMakeInActive = (
abund <= 1e-20 *
dense.xNucleiTotal);
133 if( lgMakeInActive &&
dBaseSpecies[ipSpecies].lgActive )
137 for(
long ipHi = 1; ipHi <
dBaseSpecies[ipSpecies].numLevels_max; ipHi++ )
148 (*tr).Emis().phots() = 0.;
149 (*tr).Emis().xIntensity() = 0.;
150 (*tr).Coll().col_str() = 0.;
151 (*tr).Coll().cool() = 0.;
152 (*tr).Coll().heat() = 0.;
157 if( !lgMakeInActive )
166 for(
long ipLo = 0; ipLo <
dBaseSpecies[ipSpecies].numLevels_local; ipLo++ )
180 if(
ex[0] <=
dBaseStates[ipSpecies][0].energy().WN()* 10. *DBL_EPSILON )
185 for(
long ipHi= 0; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
187 for(
long ipLo= 0; ipLo<
dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
207 int ipHi = (*tr).ipHi();
208 if (ipHi >=
dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
210 int ipLo = (*tr).ipLo();
211 AulEscp[ipHi][ipLo] = (*tr).Emis().Aul()*
212 ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc());
213 AulDest[ipHi][ipLo] = (*tr).Emis().Aul()*(*tr).Emis().Pdest();
214 AulPump[ipLo][ipHi] = (*tr).Emis().pump();
218 for(
long ipHi= 0; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
220 for(
long ipLo= 0; ipLo<
dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
231 int ipHi = (*tr).ipHi();
236 (*tr).Coll().rate_coef_ul_set()[k] = 0.f;
246 int ipHi = (*tr).ipHi();
249 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
252 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
258 else if(
dBaseTrans[ipSpecies].chLabel() ==
"Chianti" )
263 if ((*tr).ipHi() >=
dBaseSpecies[ipSpecies].numLevels_local)
265 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
267 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
273 else if(
dBaseTrans[ipSpecies].chLabel() ==
"Stout" )
278 if ((*tr).ipHi() >=
dBaseSpecies[ipSpecies].numLevels_local)
280 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
282 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
294 int ipHi = (*tr).ipHi();
297 int ipLo = (*tr).ipLo();
348 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
358 if( (*tr).Coll().rate_coef_ul_set()[
ipELECTRON] == 0. )
361 if( strcmp(
dBaseSpecies[ipSpecies].chLabel,
"Fe 3") == 0 && ipHi < 14 )
365 else if( strcmp(
dBaseSpecies[ipSpecies].chLabel,
"Fe 4") == 0 && ipHi < 12 )
369 else if( strcmp(
dBaseSpecies[ipSpecies].chLabel,
"Fe 5") == 0 && ipHi < 14 )
373 else if(
atmdat.lgGbarOn )
394 int ipHi = (*tr).ipHi();
397 int ipLo = (*tr).ipLo();
406 if( (*tr).ipCont() > 0 )
411 (*tr).Coll().rate_lu_nontherm_set() =
secondaries.x12tot *
412 ((*tr).Emis().gf()/(*tr).EnergyWN()) /
416 CollRate[ipLo][ipHi] += (*tr).Coll().rate_lu_nontherm();
417 CollRate[ipHi][ipLo] += (*tr).Coll().rate_lu_nontherm() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
482 fprintf(
ioQQQ,
" PROBLEM in dBase_solve, atom_levelN returned negative population .\n");
491 for(
long j=0;j<
dBaseSpecies[ipSpecies].numLevels_local; j++ )
500 for(
int lvl = 0; lvl <
dBaseSpecies[ipSpecies].numLevels_local; lvl++)
502 for(
int lvl2 = 0; lvl2 <
dBaseSpecies[ipSpecies].numLevels_local; lvl2++)
519 int ipHi = (*tr).ipHi();
522 int ipLo = (*tr).ipLo();
523 (*tr).Coll().cool() =
max(Cool[ipHi][ipLo],0.);
524 (*tr).Coll().heat() =
max(-Cool[ipHi][ipLo],0.);
526 if ( (*tr).ipCont() > 0 )
529 (*tr).Emis().phots() = (*tr).Emis().Aul() * (*(*tr).Hi()).Pop() *
530 ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc());
532 (*tr).Emis().xIntensity() = (*tr).Emis().phots() * (*tr).EnergyErg();
535 (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop()*
536 (*(*tr).Lo()).g()/(*(*tr).Hi()).g();
541 (*tr).Emis().ColOvTot() =
CollRate[ipLo][ipHi]/
545 (*tr).Emis().ColOvTot() = 0.;
555 (*tr).Coll().col_str() = 0.;
562 int ipHi = (*tr).ipHi();
579 enum {DEBUG_LOC=
false};
583 fprintf(
ioQQQ,
" Departure coefficients for species %li\n", ipSpecies );
584 for(
long j=0; j<
dBaseSpecies[ipSpecies].numLevels_local; j++ )
586 fprintf(
ioQQQ,
" level %li \t Depar Coef %e\n", j,
depart[j] );
593 thermal.heating[0][27] = totalHeating;
690double CHIANTI_Upsilon(
long ipSpecies,
long ipCollider,
long ipHi,
long ipLo,
double ftemp)
692 double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
693 int intxs,inttype,intsplinepts;
697 if(
AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL )
702 intsplinepts =
AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].nSplinePts;
705 fscalingparam =
AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].ScalingParam;
707 fkte = ftemp/fdeltae/1.57888e5;
713 if( inttype ==1 || inttype==4 )
715 fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam)));
717 else if(inttype == 2 || inttype == 3||inttype == 5 || inttype == 6)
719 fxt = fkte/(fkte+fscalingparam);
724 double xs[9],*spl,*spl2;
726 if(intsplinepts == 5)
729 for(intxs=0;intxs<5;intxs++)
731 xs[intxs] = 0.25*intxs;
734 printf(
"The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
739 else if(intsplinepts == 9)
742 for( intxs=0; intxs<9; intxs++ )
744 xs[intxs] = 0.125*intxs;
747 printf(
"The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
763 for(intxs=0;intxs<intsplinepts;intxs++)
765 printf(
"The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]);
772 fsups =
linint( xs, spl, intsplinepts, fxt);
777 fups = fsups*log(fkte+exp(1.0));
779 else if(inttype == 2)
783 else if(inttype == 3)
785 fups = fsups/(fkte+1.0) ;
787 else if(inttype == 4)
789 fups = fsups*log(fkte+fscalingparam) ;
791 else if(inttype == 5)
795 else if(inttype == 6)
797 fups = pow(10.0,fsups) ;
806 fprintf(
ioQQQ,
" WARNING: Negative upsilon in species %s, collider %li, indices %4li %4li, Te = %e.\n",
807 dBaseSpecies[ipSpecies].chLabel, ipCollider, ipHi, ipLo, ftemp );
void atom_levelN(long int nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double ***AulEscp, double ***col_str, double ***AulDest, double ***AulPump, double ***CollRate, const double source[], const double sink[], bool lgCollRateDone, double *cooltl, double *coolder, const char *chLabel, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE, multi_arr< double, 2 > *Cool, multi_arr< double, 2 > *dCooldT)