42 if( !
exists( SpeciesCurrent ) )
46 fprintf(
ioQQQ,
" PROBLEM dBase_solve did not find molecular species %li\n",ipSpecies);
67 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
73 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
81 for(
long ipHi = 1; ipHi <
dBaseSpecies[ipSpecies].numLevels_max; ipHi++ )
92 (*tr).Emis().xIntensity() = 0.;
93 (*tr).Emis().xObsIntensity() = 0.;
94 (*tr).Coll().col_str() = 0.;
95 (*tr).Coll().cool() = 0.;
96 (*tr).Coll().heat() = 0.;
101 if( !lgMakeInActive )
109 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
118 int ipHi = (*tr).ipHi();
122 (*tr).Coll().rate_coef_ul_set()[k] = 0.f;
131 int ipHi = (*tr).ipHi();
134 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
137 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
140 tr->Coll().is_gbar() = 0;
149 if ((*tr).ipHi() >=
dBaseSpecies[ipSpecies].numLevels_local)
151 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
153 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
156 tr->Coll().is_gbar() = 0;
165 if ((*tr).ipHi() >=
dBaseSpecies[ipSpecies].numLevels_local)
167 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
169 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
172 tr->Coll().is_gbar() = 0;
182 int ipHi = (*tr).ipHi();
234 for(
long intCollNo=0; intCollNo<
ipNCOLLIDER; intCollNo++)
244 if( (*tr).Coll().rate_coef_ul_set()[
ipELECTRON] == 0. )
256 (*tr).Coll().is_gbar() = 1;
263 if( tr->Coll().col_str() == -FLT_MAX && tr->Coll().rate_coef_ul()[
ipELECTRON] > 0.0 )
282 static bool lgFirstPass =
true;
283 static long maxNumLevels = 1;
284 static double *g, *ex, *pops, *
depart, *source, *sink;
285 static double **AulEscp, **AulDest, **AulPump, **CollRate;
289 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
292 g = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
293 ex = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
294 pops = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
295 depart = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
296 source = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
297 sink = (
double *)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double));
299 AulEscp = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
300 AulDest = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
301 AulPump = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
302 CollRate = (
double **)
MALLOC( (
unsigned long)(maxNumLevels)*
sizeof(
double *));
304 AulEscp[0] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels*maxNumLevels)*
sizeof(
double));
305 AulDest[0] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels*maxNumLevels)*
sizeof(
double));
306 AulPump[0] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels*maxNumLevels)*
sizeof(
double));
307 CollRate[0] = (
double *)
MALLOC( (
unsigned long)(maxNumLevels*maxNumLevels)*
sizeof(
double));
308 for(
long j=1; j< maxNumLevels; j++ )
310 AulEscp[j]=AulEscp[j-1]+maxNumLevels;
311 AulDest[j]=AulDest[j-1]+maxNumLevels;
312 AulPump[j]=AulPump[j-1]+maxNumLevels;
313 CollRate[j]=CollRate[j-1]+maxNumLevels;
320 memset( g, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
321 memset( ex, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
322 memset( pops, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
323 memset( depart, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
324 memset( source, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
325 memset( sink, 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
326 for(
long j=0; j< maxNumLevels; j++ )
328 memset( AulEscp[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
329 memset( AulDest[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
330 memset( AulPump[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
331 memset( CollRate[j], 0, (
unsigned long)(maxNumLevels)*
sizeof(
double) );
336 double totalHeating = 0.;
337 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
356 for(
long ipLo = 0; ipLo <
dBaseSpecies[ipSpecies].numLevels_local; ipLo++ )
362 ex[ipLo] =
dBaseStates[ipSpecies][ipLo].energy().WN() -
375 for(
long ipHi= 0; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
377 for(
long ipLo= 0; ipLo<
dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
379 AulEscp[ipHi][ipLo] = 0.;
382 for(
long ipHi= 0; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
384 for(
long ipLo= 0; ipLo<
dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
386 AulDest[ipHi][ipLo] = 0.;
389 for(
long ipLo= 0; ipLo<
dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
391 for(
long ipHi= 0; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
393 AulPump[ipLo][ipHi] = 0.;
397 const bool isO1 = ( strcmp( spName,
"O 1" ) == 0 ),
398 isO3 = (strcmp(spName,
"O 3") == 0),
399 isCa2 = (strcmp(spName,
"Ca 2") == 0),
400 isN1 = (strcmp(spName,
"N 1") == 0),
401 isMg2 = (strcmp(spName,
"Mg 2") == 0);
406 int ipHi = (*tr).ipHi();
407 if (ipHi >=
dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
409 int ipLo = (*tr).ipLo();
410 AulEscp[ipHi][ipLo] = (*tr).Emis().Aul()*(*tr).Emis().Pesc_total();
411 AulDest[ipHi][ipLo] = (*tr).Emis().Aul()*(*tr).Emis().Pdest();
412 AulPump[ipLo][ipHi] = (*tr).Emis().pump();
413 AulPump[ipHi][ipLo] = AulPump[ipLo][ipHi]*g[ipLo]/g[ipHi];
416 double xtraExRate = 0.;
417 double xtraDxRate = 0.;
422 static const double lo1_nrg = 0.;
423 static const double lo2_nrg = 158.265;
424 static const double lo3_nrg = 226.977;
426 static const double hi1_nrg = 97488.378;
427 static const double hi2_nrg = 97488.448;
428 static const double hi3_nrg = 97488.538;
430 if( (
fp_equal( tr->Lo()->energy().WN(),lo1_nrg ) ||
431 fp_equal( tr->Lo()->energy().WN(),lo2_nrg ) ||
432 fp_equal( tr->Lo()->energy().WN(),lo3_nrg ) ) &&
433 (
fp_equal( tr->Hi()->energy().WN(),hi1_nrg ) ||
434 fp_equal( tr->Hi()->energy().WN(),hi2_nrg ) ||
435 fp_equal( tr->Hi()->energy().WN(),hi3_nrg ) ) )
441 if( tr->ipHi() == 3 )
447 if( tr->ipHi() == 3 )
453 if( tr->ipLo() == 0 && tr->ipHi() < 5 )
459 if( tr->ipLo() == 0 && (tr->ipHi() == 1 || tr->ipHi() == 2) )
465 if( tr->ipLo() == 0 && (tr->ipHi() == 4 || tr->ipHi() == 5) )
468 else if( strcmp(spName,
"Fe 2") == 0 )
473 AulPump[ipLo][ipHi] += xtraExRate;
474 AulPump[ipHi][ipLo] += xtraDxRate;
481 int ipHi = (*tr).ipHi();
490 for(
long ipHi= 0; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
492 for(
long ipLo= 0; ipLo<
dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
494 CollRate[ipHi][ipLo] = 0.;
502 int ipHi = (*tr).ipHi();
505 int ipLo = (*tr).ipLo();
506 CollRate[ipHi][ipLo] = (*tr).Coll().ColUL( colld );
509 for(
int ipHi=1; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ++ipHi)
510 arg[ipHi] = -(ex[ipHi]-ex[ipHi-1])*T1CM /
phycon.
te;
512 for (
int ipLo=0; ipLo<
dBaseSpecies[ipSpecies].numLevels_local-1; ++ipLo)
514 double boltz_over_glo=1./ g[ipLo];
515 for (
int ipHi=ipLo+1; ipHi<
dBaseSpecies[ipSpecies].numLevels_local; ++ipHi)
517 boltz_over_glo *= bstep[ipHi];
518 CollRate[ipLo][ipHi] = CollRate[ipHi][ipLo] * g[ipHi] * boltz_over_glo;
536 if( (*tr).ipCont() > 0 )
538 int ipHi = (*tr).ipHi();
541 int ipLo = (*tr).ipLo();
542 (*tr).Coll().rate_lu_nontherm_set() = sfac *
543 ((*tr).Emis().gf()/(*tr).EnergyWN());
545 CollRate[ipLo][ipHi] += (*tr).Coll().rate_lu_nontherm();
546 CollRate[ipHi][ipLo] += (*tr).Coll().rate_lu_nontherm() * g[ipLo] / g[ipHi];
552 enum {DEBUG_LOC=
false};
553 if( DEBUG_LOC && (
nzone>110) )
563 int ipLo = tr->ipLo();
564 int ipHi = tr->ipHi();
565 if( ipLo == 39 && ipHi == 139 )
570 if( ipLo == 36 && ipHi == 133 )
581 double grnd_excit = 0.0;
582 double cooltl, coolder;
584 bool lgZeroPop, lgDeBug =
false;
648 [(*(*
dBaseTrans[ipSpecies].begin()).Lo()).IonStg()-1]
655 fprintf(
ioQQQ,
" PROBLEM in dBase_solve, atom_levelN returned negative population .\n");
660 while( (pops[
dBaseSpecies[ipSpecies].numLevels_local-1]<=0 ) &&
664 for(
long j=0;j<
dBaseSpecies[ipSpecies].numLevels_local; j++ )
667 dBaseStates[ipSpecies][j].DepartCoef() = depart[j];
680 for(
int lvl = 0; lvl <
dBaseSpecies[ipSpecies].numLevels_local; lvl++)
682 for(
int lvl2 = 0; lvl2 <
dBaseSpecies[ipSpecies].numLevels_local; lvl2++)
686 dBaseStates[ipSpecies][lvl].DestCollBB() += CollRate[lvl][lvl2];
687 dBaseStates[ipSpecies][lvl].CreatCollBB() += CollRate[lvl2][lvl]*pops[lvl2];
689 dBaseStates[ipSpecies][lvl].DestPhotoBB() += AulDest[lvl][lvl2];
697 tr != dBaseTrans[ipSpecies].end(); ++tr)
699 int ipHi = (*tr).ipHi();
702 int ipLo = (*tr).ipLo();
703 (*tr).Coll().cool() =
max(Cool[ipHi][ipLo],0.);
704 (*tr).Coll().heat() =
max(-Cool[ipHi][ipLo],0.);
706 if ( (*tr).ipCont() > 0 )
709 (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop()*
710 (*(*tr).Lo()).g()/(*(*tr).Hi()).g();
715 if( CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi] > 0. )
717 (*tr).Emis().ColOvTot() = CollRate[ipLo][ipHi]/
718 (CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi]);
721 (*tr).Emis().ColOvTot() = 0.;
743 tr->Lo()->Pop()*exRate -
744 tr->Hi()->Pop()*dxRate);
763 tr != dBaseTrans[ipSpecies].end(); ++tr)
765 int ipHi = (*tr).ipHi();
782 enum {DEBUG_LOC=
false};
786 fprintf(
ioQQQ,
" Departure coefficients for species %li\n", ipSpecies );
787 for(
long j=0; j<
dBaseSpecies[ipSpecies].numLevels_local; j++ )
789 fprintf(
ioQQQ,
" level %li \t Depar Coef %e\n", j, depart[j] );
821 for(
int j = 0; j < n; j ++)
823 ASSERT( x[j] > 0. && y[j] > 0.);
828 double fupsilon = 0.;
833 else if( ftemp > x[n-1] )
839 fupsilon =
linint(&x[0],&y[0],n,ftemp);
850 tr.
Coll().
col_str() = (rate * (*tr.
Hi()).g()*sqrt(ftemp))/COLL_CONST;
859 rate = (COLL_CONST*fupsilon)/((*tr.
Hi()).g()*sqrt(ftemp));
864 fprintf(
ioQQQ,
"PROBLEM: Stout data format does not support using collision strengths with "
865 "non-electron colliders.\n");
892 rate = (COLL_CONST*fupsilon)/((*tr.
Hi()).g()*sqrt(ftemp));
900 double CHIANTI_Upsilon(
long ipSpecies,
long ipCollider,
long ipHi,
long ipLo,
double ftemp)
902 double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
903 int intxs,inttype,intsplinepts;
907 if(
AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL )
912 intsplinepts =
AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].nSplinePts;
915 fscalingparam =
AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].ScalingParam;
917 fkte = ftemp/fdeltae/1.57888e5;
923 if( inttype ==1 || inttype==4 )
925 fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam)));
927 else if(inttype == 2 || inttype == 3||inttype == 5 || inttype == 6)
929 fxt = fkte/(fkte+fscalingparam);
937 for(intxs=0;intxs<intsplinepts;intxs++)
939 double coeff = (double)1/(intsplinepts-1);
940 xs[intxs] = coeff*intxs;
943 printf(
"The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
948 const bool SPLINE_INTERP=
false;
951 fsups =
linint( xs, spl, intsplinepts, fxt);
962 for(intxs=0;intxs<intsplinepts;intxs++)
964 printf(
"The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]);
969 splint(xs,spl,spl2,intsplinepts,fxt,&fsups);
975 fups = fsups*log(fkte+EE);
977 else if(inttype == 2)
981 else if(inttype == 3)
983 fups = fsups/(fkte+1.0) ;
985 else if(inttype == 4)
987 fups = fsups*log(fkte+fscalingparam) ;
989 else if(inttype == 5)
993 else if(inttype == 6)
995 fups =
exp10(fsups) ;
1004 fprintf(
ioQQQ,
" WARNING: Negative upsilon in species %s, collider %li, indices %4li %4li, Te = %e.\n",
1005 dBaseSpecies[ipSpecies].chLabel, ipCollider, ipHi, ipLo, ftemp );
1028 fixit(
"ticket #78 refers");
1049 " P8446 finds Lbeta, OI widths=%10.2e%10.2e and esc prob=%10.2e%10.2e esAB=%10.2e\n",
1058 xoi = opacoi/(opacoi + opaclb);
1059 xlb = opaclb/(opacoi + opaclb);
1069 " P8446 opac Lb, OI=%10.2e%10.2e X Lb, OI=%10.2e%10.2e FLb, OI=%10.2e%10.2e\n",
1070 opaclb, opacoi, xlb, xoi, flb, foi );
1080 xtraDxRate += (
realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - esab)*foi));
1102 PhotoRate5 = 1.7e-18*hlgam;
1103 PhotoRate4 = 8.4e-19*hlgam;
1104 PhotoRate3 = 7.0e-18*hlgam;
1105 PhotoRate2 = 4.8e-18*hlgam;
1107 if( tr.
ipHi() + 1 == 2)
1109 xtraDxRate += PhotoRate2;
1111 else if( tr.
ipHi() + 1 == 3 )
1113 xtraDxRate += PhotoRate3;
1115 else if( tr.
ipHi() + 1 == 4 )
1117 xtraDxRate += PhotoRate4;
1119 else if( tr.
ipHi() + 1 == 5 )
1121 xtraDxRate += PhotoRate5;
1149 double EnerLyaProf1 = 82259.0 - de*2.0;
1150 double EnerLyaProf2 = 82259.0 - de;
1151 double EnerLyaProf3 = 82259.0 + de;
1152 double EnerLyaProf4 = 82259.0 + de*2.0;
1154 double PhotOccNumLyaCenter = 0.;
1160 PhotOccNumLyaCenter =
1182 if( EnergyWN >= EnerLyaProf1 && EnergyWN <= EnerLyaProf4 && taux > 1 )
1199 double PhotOccNum_at_nu = 0.,
1201 if( EnergyWN < EnerLyaProf2 )
1204 PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnergyWN - EnerLyaProf1)/ de;
1206 else if( EnergyWN < EnerLyaProf3 )
1209 PhotOccNum_at_nu = PhotOccNumLyaCenter;
1214 PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnerLyaProf4 - EnergyWN)/de;
1227 PumpRate = FeIILineWidthHz * PhotOccNum_at_nu * tr.
Emis().
Aul()*
1228 powi(82259.0f/EnergyWN,3);
1231 PumpRate = tr.
Emis().
Aul()* PhotOccNum_at_nu;
1234 xtraDxRate += PumpRate;
1237 xtraExRate += PumpRate * (*tr.
Hi()).g()/(*tr.
Lo()).g();
void CoolAdd(const char *chLabel, realnum lambda, double cool)
void DumpLine(const TransitionProxy &t)
vector< StoutCollArray > StoutCollData
NORETURN void TotalInsanity(void)
double depart(const genericState &gs)
STATIC void setXtraRatesO1(const TransitionProxy &tr, double &xtraExRate, double &xtraDxRate)
void set_xIntensity(const TransitionProxy &t)
double CHIANTI_Upsilon(long, long, long, long, double)
molezone * findspecieslocal(const char buf[])
STATIC double LeidenCollRate(long, long, const TransitionProxy &, double)
multi_arr< CollRateCoeffArray, 2 > AtmolCollRateCoeff
static const bool DEBUGSTATE
void splint(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)
double * rate_coef_ul_set() const
double ** ExcitationGround
t_iso_sp iso_sp[NISO][LIMELM]
void vexp(const double x[], double y[], long nlo, long nhi)
double xIonDense[LIMELM][LIMELM+1]
bool fp_equal(sys_float x, sys_float y, int n=3)
ColliderList colliders(dense)
realnum & EnergyWN() const
double energy(const genericState &gs)
EmissionList::reference Emis() const
STATIC void setXtraRatesFe2(const TransitionProxy &tr, double &xtraExRate, double &xtraDxRate)
STATIC void setXtraRatesCa2(const TransitionProxy &tr, double &xtraDxRate)
double InterpCollRate(const CollRateCoeffArray &rate_table, const long &ipHi, const long &ipLo, const double &ftemp)
qList::iterator Hi() const
double powi(double, long int)
static realnum dBaseAbund(long ipSpecies)
realnum GetDopplerWidth(realnum massAMU)
bool exists(const molecule *m)
vector< multi_arr< CollSplinesArray, 3 > > AtmolCollSplines
TransitionProxy trans(const long ipHi, const long ipLo)
void setHeating(long nelem, long ion, double heating)
realnum AtomicWeight[LIMELM]
qList::iterator Lo() const
void EmLineZero(EmissionList::reference t)
realnum Pesc_total() const
STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy &, double ftemp)
realnum & col_str() const
CollisionProxy Coll() const
vector< vector< long > > ipSpecIon
#define DEBUG_ENTRY(funcname)
double linint(const double x[], const double y[], long n, double xval)
vector< qList > dBaseStates
double elementcool[LIMELM+1]
const double * rate_coef_ul() const
vector< species > dBaseSpecies
int fprintf(const Output &stream, const char *format,...)
STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy &, double ftemp)
t_secondaries secondaries
void CollisionZero(const CollisionProxy &t)
vector< TransitionList > dBaseTrans
void dBaseUpdateCollCoeffs(void)
void MakeCS(const TransitionProxy &t)