19 #define DEBUGSTATE false
32 double *g, *ex, *pops, *depart;
33 double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate,fcolldensity[9];
34 double *source, *sink;
35 double cooltl, coolder;
37 long ipLo,ipHi,intCollNo;
39 bool lgZeroPop, lgDeBug =
false;
51 fprintf(
ioQQQ,
" PROBLEM atmol_popsolve did not find molecular species %li\n",i);
53 abund = SpeciesCurrent->
hevmol;
67 for(
long ipHi = 1; ipHi < nlev; ipHi++ )
70 for(
long ipLo = 0; ipLo < ipHi; ipLo++ )
82 g = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
83 ex = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
84 pops = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
85 depart = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
86 source = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
87 sink = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
89 AulEscp = (
double **)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double *));
90 col_str = (
double **)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double *));
91 AulDest = (
double **)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double *));
92 AulPump = (
double **)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double *));
93 CollRate = (
double **)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double *));
95 for(
long j=0; j< nlev; j++ )
97 AulEscp[j] = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
98 col_str[j] = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
99 AulDest[j] = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
100 AulPump[j] = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
101 CollRate[j] = (
double *)
MALLOC( (
unsigned long)(nlev)*
sizeof(
double));
104 for( ipLo = 0; ipLo < nlev; ipLo++ )
112 AulEscp[ipLo][ipLo] = 0.;
113 AulDest[ipLo][ipLo] = 0.;
114 AulPump[ipLo][ipLo] = 0.;
117 for( ipHi=1; ipHi<nlev; ipHi++ )
119 for( ipLo=0; ipLo<ipHi; ipLo++ )
137 AulEscp[ipLo][ipHi] = 0.;
138 AulDest[ipLo][ipHi] = 0.;
139 AulPump[ipHi][ipLo] = 0.;
144 for( ipHi= 0; ipHi<nlev; ipHi++)
146 for( ipLo= 0; ipLo<nlev; ipLo++)
148 col_str[ipHi][ipLo] = 0.;
149 CollRate[ipHi][ipLo] = 0.;
157 for( ipHi=1; ipHi<nlev; ipHi++)
159 for( ipLo=0; ipLo<ipHi; ipLo++)
217 fcolldensity[6] = (fcolldensity[8])*(0.75);
218 fcolldensity[7] = (fcolldensity[8])*(0.25);
222 for( ipHi=0; ipHi<nlev; ipHi++ )
224 for( ipLo=0; ipLo<nlev; ipLo++ )
228 CollRate[ipHi][ipLo] = 0.;
236 CollRate[ipHi][ipLo] +=
CollRatesArray[i][intCollNo][ipHi][ipLo]*fcolldensity[intCollNo];
247 CollRate[ipHi][ipLo] += 0.7 *
CollRatesArray[i][8][ipHi][ipLo]*fcolldensity[3];
256 while ( ((ex[nlev-1]*
T1CM) >
phycon.
te*25.) && (nlev > 1) )
317 ASSERT( lgZeroPop ==
false );
319 for(
long j=0; j< nlev; j++ )
324 else if( nNegPop > 0 )
329 fprintf(
ioQQQ,
" PROBLEM, atom_levelN returned negative population .\n");
337 for(
long j=
Species[i].numLevels_max-1; j>0; j--)
362 atom_levelN(nlev,abund,g,ex,
'w',pops,depart,&AulEscp,&col_str,
363 &AulDest, &AulPump, &CollRate, source, sink,
true,&cooltl,
364 &coolder, spName, &nNegPop, &lgZeroPop, lgDeBug );
366 for(
long j=0; j< nlev; j++ )
379 for(
long ipHi = 1; ipHi < nlev; ipHi++ )
381 for(
long ipLo = 0; ipLo < ipHi; ipLo++ )
410 enum {DEBUG_LOC=
false};
414 fprintf(
ioQQQ,
" Departure coefficients for species %li\n", i );
415 for(
long j=0; j< nlev; j++ )
417 fprintf(
ioQQQ,
" level %li \t Depar Coef %e\n", j, depart[j] );
431 for(
long j=0; j< nlev; j++ )
464 double LeidenCollRate(
long ipSpecies,
long ipCollider,
long ipHi,
long ipLo,
double ftemp)
493 printf(
"\nThe collision rate at temperature %f is %e",ftemp,ret_collrate);
496 return(ret_collrate);
501 double Chianti_Upsilon(
long ipSpecies,
long ipCollider,
long ipHi,
long ipLo,
double ftemp)
503 double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
504 double *xs,*spl,*spl2;
505 int i,intxs,inttype,intsplinepts;
509 if(
AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL )
519 fkte = ftemp/fdeltae/1.57888e5;
525 if(inttype ==1 || inttype == 4)
527 fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam)));
529 else if(inttype == 2 || inttype == 3||inttype == 5 || inttype == 6)
531 fxt = fkte/(fkte+fscalingparam);
537 if(intsplinepts == 5)
539 xs = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
540 spl = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
541 spl2 = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
542 for(intxs=0;intxs<5;intxs++)
544 xs[intxs] = 0.25*intxs;
548 printf(
"The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
553 else if(intsplinepts == 9)
555 xs = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
556 spl = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
557 spl2 = (
double *)
MALLOC((
unsigned long)intsplinepts*
sizeof(double));
558 for( intxs=0; intxs<9; intxs++ )
560 xs[intxs] = 0.125*intxs;
564 printf(
"The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
575 for( i=0; i<intsplinepts; i++)
583 for(intxs=0;intxs<intsplinepts;intxs++)
585 printf(
"The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]);
590 splint(xs,spl,spl2,intsplinepts,fxt,&fsups);
595 fups = fsups*log(fkte+exp(1.0));
597 else if(inttype == 2)
601 else if(inttype == 3)
603 fups = fsups/(fkte+1.0) ;
605 else if(inttype == 4)
607 fups = fsups*log(fkte+fscalingparam) ;
609 else if(inttype == 5)
613 else if(inttype == 6)
615 fups = pow(10.0,fsups) ;