53 static bool lgFirst=
true;
58 static long int nUsed;
62 double rot_cooling , dCoolDT;
63 static long int ndimMalloced = 0;
72 ndimMalloced = nRotate;
74 excit = (
double *)
MALLOC(
sizeof(
double)*(size_t)(nRotate+1) );
75 stat = (
double *)
MALLOC(
sizeof(
double)*(size_t)(nRotate+1) );
76 pops = (
double *)
MALLOC(
sizeof(
double)*(size_t)(nRotate+1) );
77 create = (
double *)
MALLOC(
sizeof(
double)*(size_t)(nRotate+1) );
78 destroy = (
double *)
MALLOC(
sizeof(
double)*(size_t)(nRotate+1) );
79 depart = (
double *)
MALLOC(
sizeof(
double)*(size_t)(nRotate+1) );
81 AulPump = ((
double **)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double *)));
82 CollRate = ((
double **)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double *)));
83 AulDest = ((
double **)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double *)));
84 AulEscp = ((
double **)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double *)));
85 col_str = ((
double **)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double *)));
86 for( i=0; i<(nRotate+1); ++i )
88 AulPump[i] = ((
double *)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double )));
89 CollRate[i] = ((
double *)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double )));
90 AulDest[i] = ((
double *)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double )));
91 AulEscp[i] = ((
double *)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double )));
92 col_str[i] = ((
double *)
MALLOC((
size_t)(nRotate+1)*
sizeof(
double )));
97 if( nRotate > ndimMalloced )
99 fprintf(
ioQQQ,
" CO_PopsEmisCool has been called with the number of rotor levels greater than space allocated.\n");
104 for( i=0; i < (nRotate+1); i++ )
108 for( j=0; j < (nRotate+1); j++ )
118 for( j=0; j < nRotate; j++ )
121 stat[j] = (*Rotate)[j].Lo->g;
124 stat[nRotate] = (*Rotate)[nRotate-1].Hi->g;
129 for( j=1; j < nRotate; j++ )
132 excit[j] = excit[j-1] + (*Rotate)[j-1].EnergyK;
136 excit[nRotate] = excit[nRotate-1] + (*Rotate)[nRotate-1].EnergyK;
151 while ( (excit[nUsed] >
phycon.
te*25.) && (nUsed > 5) )
156 for( j=0; j < nRotate; j++ )
159 AulEscp[j+1][j] = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc);
160 AulDest[j+1][j] = (*Rotate)[j].Emis->Aul*(*Rotate)[j].Emis->Pdest;
162 (*Rotate)[j].Coll.cs = 1.;
166 AulPump[j][j+1] = (*Rotate)[j].Emis->pump;
173 for( ilo=0; ilo < nRotate-1; ilo++ )
177 for( ihi=ilo+2; ihi <= nRotate; ihi++ )
185 AulEscp[ihi][ilo] = 0;
186 AulDest[ihi][ilo] = 0;
192 for( ilo=0; ilo < nRotate; ilo++ )
196 for( ihi=ilo+1; ihi <=
MIN2(nRotate,ilo+5); ihi++ )
199 double a[5]={1.66, 2.80, 1.19, 1.00, 1.30};
200 double b[5]={1.67, 1.47, 1.85, 1.55, 2.24};
211 CollRate[ihi][ilo] = a[ihi-ilo-1]*1.e-10*(*Rotate)[ilo].Lo->g/(*Rotate)[ilo].Hi->g*
212 (1.+(excit[ihi]-excit[ilo])/TeUsed) *
213 sexp( b[ihi-ilo-1]*sqrt((excit[ihi]-excit[ilo])/TeUsed) )*collid;
218 (*Rotate)[ilo].Coll.ColUL = (
realnum)CollRate[ihi][ilo];
224 CollRate[ilo][ihi] = CollRate[ihi][ilo]*
226 (*Rotate)[ilo].Hi->g / (*Rotate)[ilo].Lo->g;
233 for( ihi=ilo+6; ihi <= nRotate; ihi++ )
235 CollRate[ihi][ilo] = 0.;
236 CollRate[ilo][ihi] = 0.;
271 *Cooling = (
realnum)rot_cooling;
280 fprintf(
ioQQQ,
"COcool\t%i\t%.2f\t%e\t%e\t%e\t%e\t%e\n",
291 for( i=nUsed+1; i<=nRotate; ++i )
306 fprintf(
ioQQQ,
"CO_PopsEmisCool called atom_levelN which returned negative populations.\n");
310 for( j=0; j<nRotate; ++j )
314 (*Rotate)[j].Lo->Pop = pops[j];
315 (*Rotate)[j].Hi->Pop = pops[j+1];
316 (*Rotate)[j].Emis->PopOpc = (pops[j] - pops[j+1]*(*Rotate)[j].Lo->g/(*Rotate)[j].Hi->g);
319 (*Rotate)[j].Emis->phots = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc)*pops[j+1];
324 (*Rotate)[j].Emis->ots = (*Rotate)[j].Aul*(*Rotate)[j].Emis->Pdest*(
realnum)(*Rotate)[j].Hi->Pop;
329 (*Rotate)[j].Emis->xIntensity = (*Rotate)[j].Emis->phots*(*Rotate)[j].EnergyErg;
332 (*Rotate)[j].Emis->ColOvTot = (
realnum)(CollRate[j][j+1]/(CollRate[j][j+1]+(*Rotate)[j].Emis->pump) );
337 EnrLU = (*Rotate)[j].Lo->Pop*CollRate[j][j+1]*(*Rotate)[j].EnergyErg;
338 EnrUL = (*Rotate)[j].Hi->Pop*CollRate[j+1][j]*(*Rotate)[j].EnergyErg;
341 (*Rotate)[j].Coll.cool = EnrLU - EnrUL*(*Rotate)[j].Emis->ColOvTot;
343 (*Rotate)[j].Coll.heat = EnrUL*(1.f - (*Rotate)[j].Emis->ColOvTot);
351 (*Rotate)[nUsed-1].Emis->xIntensity > (*Rotate)[nUsed-2].Emis->xIntensity &&
370 if( strcmp(chLabel,
"ZERO") == 0 )
372 static bool lgFIRST=
true;
388 else if( strcmp(chLabel,
"ADD ") == 0 )
400 else if( strcmp(chLabel,
"PRIN") != 0 )
402 fprintf(
ioQQQ,
" CO_Colden does not understand the label %s\n",
414 if( isotope !=12 && isotope != 13 )
416 fprintf(
ioQQQ,
" cdCO_colden can only do 12CO and 13CO\n");
421 fprintf(
ioQQQ,
" cdCO_colden - rotation quantum number must be 0 or greater, and less than %li\n",