00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "taulines.h"
00008 #include "dense.h"
00009 #include "thermal.h"
00010 #include "hmi.h"
00011 #include "radius.h"
00012 #include "atoms.h"
00013 #include "phycon.h"
00014 #include "rt.h"
00015 #include "lines_service.h"
00016 #include "cddrive.h"
00017 #include "mole.h"
00018
00019
00020
00021
00022 static double *col12 , *col13;
00023
00024
00025 void CO_PopsEmisCool(
00026 EmLine ** Rotate ,
00027 long int nRotate ,
00028 float abundan,
00029 const char * chLabel ,
00030 float * Cooling ,
00031 float * dCoolingDT )
00032 {
00033
00034
00035 static double
00036 **AulEscp ,
00037 **col_str ,
00038 **AulDest,
00039
00040 **AulPump,
00041 **CollRate,
00042 *pops,
00043 *create,
00044 *destroy,
00045 *depart,
00046
00047 *stat ,
00048
00049 *excit;
00050
00051
00052
00053 static bool lgFirst=true;
00054 long int i,
00055 j,
00056 ilo ,
00057 ihi;
00058 static long int nUsed;
00059 bool lgDeBug;
00060 int nNegPop;
00061 bool lgZeroPop;
00062 double rot_cooling , dCoolDT;
00063 static long int ndimMalloced = 0;
00064
00065 DEBUG_ENTRY( "CO_PopsEmisCool()" );
00066
00067 if( lgFirst )
00068 {
00069
00070 lgFirst = false;
00071
00072 ndimMalloced = nRotate;
00073
00074 excit = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
00075 stat = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
00076 pops = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
00077 create = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
00078 destroy = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
00079 depart = (double *)MALLOC( sizeof(double)*(size_t)(nRotate+1) );
00080
00081 AulPump = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
00082 CollRate = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
00083 AulDest = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
00084 AulEscp = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
00085 col_str = ((double **)MALLOC((size_t)(nRotate+1)*sizeof(double *)));
00086 for( i=0; i<(nRotate+1); ++i )
00087 {
00088 AulPump[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
00089 CollRate[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
00090 AulDest[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
00091 AulEscp[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
00092 col_str[i] = ((double *)MALLOC((size_t)(nRotate+1)*sizeof(double )));
00093 }
00094 }
00095
00096
00097 if( nRotate > ndimMalloced )
00098 {
00099 fprintf(ioQQQ," CO_PopsEmisCool has been called with the number of rotor levels greater than space allocated.\n");
00100 puts( "[Stop in CO_PopsEmisCool]" );
00101 cdEXIT(EXIT_FAILURE);
00102 }
00103
00104
00105 for( i=0; i < (nRotate+1); i++ )
00106 {
00107 create[i] = 0.;
00108 destroy[i] = 0.;
00109 for( j=0; j < (nRotate+1); j++ )
00110 {
00111 col_str[j][i] = 0.;
00112 AulEscp[j][i] = 0.;
00113 AulDest[j][i] = 0.;
00114 AulPump[j][i] = 0.;
00115 }
00116 }
00117
00118
00119 for( j=0; j < nRotate; j++ )
00120 {
00121
00122 stat[j] = (*Rotate)[j].gLo;
00123 }
00124
00125 stat[nRotate] = (*Rotate)[nRotate-1].gHi;
00126
00127
00128
00129 excit[0] = 0.;
00130 for( j=1; j < nRotate; j++ )
00131 {
00132
00133 excit[j] = excit[j-1] + (*Rotate)[j-1].EnergyK;
00134 }
00135
00136
00137 excit[nRotate] = excit[nRotate-1] + (*Rotate)[nRotate-1].EnergyK;
00138
00139 nUsed = nRotate;
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 5) )
00153 {
00154 --nUsed;
00155 }
00156
00157 for( j=0; j < nRotate; j++ )
00158 {
00159
00160 AulEscp[j+1][j] = (*Rotate)[j].Aul*((*Rotate)[j].Pesc + (*Rotate)[j].Pelec_esc);
00161 AulDest[j+1][j] = (*Rotate)[j].Aul*(*Rotate)[j].Pdest;
00162
00163 (*Rotate)[j].cs = 1.;
00164
00165 col_str[j+1][j] = (*Rotate)[j].cs*hmi.H2_total/dense.eden;
00166
00167 AulPump[j][j+1] = (*Rotate)[j].pump;
00168
00169
00170
00171 }
00172
00173
00174 for( ilo=0; ilo < nRotate-1; ilo++ )
00175 {
00176
00177
00178 for( ihi=ilo+2; ihi <= nRotate; ihi++ )
00179 {
00180
00181
00182
00183 col_str[ihi][ilo] = 1. *hmi.H2_total/dense.eden;
00184
00185
00186 AulEscp[ihi][ilo] = 0;
00187 AulDest[ihi][ilo] = 0;
00188 }
00189 }
00190
00191
00192
00193 for( ilo=0; ilo < nRotate; ilo++ )
00194 {
00195
00196
00197 for( ihi=ilo+1; ihi <= MIN2(nRotate,ilo+5); ihi++ )
00198 {
00199
00200 double a[5]={1.66, 2.80, 1.19, 1.00, 1.30};
00201 double b[5]={1.67, 1.47, 1.85, 1.55, 2.24};
00202
00203
00204
00205
00206 double collid = hmi.H2_total + dense.xIonDense[ipHELIUM][0]/1.37;
00207
00208
00209 double TeUsed = MAX2(10., phycon.te );
00210
00211
00212 CollRate[ihi][ilo] = a[ihi-ilo-1]*1.e-10*(*Rotate)[ilo].gLo/(*Rotate)[ilo].gHi*
00213 (1.+(excit[ihi]-excit[ilo])/TeUsed) *
00214 sexp( b[ihi-ilo-1]*sqrt((excit[ihi]-excit[ilo])/TeUsed) )*collid;
00215
00216 if( ihi == ilo+1 )
00217 {
00218
00219 (*Rotate)[ilo].ColUL = (float)CollRate[ihi][ilo];
00220
00221 LineConvRate2CS( &(*Rotate)[ilo] , (*Rotate)[ilo].ColUL);
00222 }
00223
00224
00225 CollRate[ilo][ihi] = CollRate[ihi][ilo]*
00226 sexp( (excit[ihi]-excit[ilo])/phycon.te )*
00227 (*Rotate)[ilo].gHi / (*Rotate)[ilo].gLo;
00228
00229
00230
00231
00232 }
00233
00234 for( ihi=ilo+6; ihi <= nRotate; ihi++ )
00235 {
00236 CollRate[ihi][ilo] = 0.;
00237 CollRate[ilo][ihi] = 0.;
00238 }
00239 }
00240
00241 lgDeBug = false;
00242
00243 atom_levelN(
00244
00245
00246 nUsed+1,
00247 abundan,
00248 stat,
00249 excit,
00250 'K',
00251 pops,
00252 depart,
00253 &AulEscp,
00254 &col_str,
00255 &AulDest,
00256 &AulPump,
00257 &CollRate,
00258 create,
00259 destroy,
00260
00261 true,
00262
00263 &rot_cooling,
00264 &dCoolDT,
00265 chLabel,
00266
00267 &nNegPop,
00268 &lgZeroPop,
00269 lgDeBug );
00270
00271
00272 *Cooling = (float)rot_cooling;
00273
00274
00275
00276
00277 *dCoolingDT = (float)(rot_cooling/phycon.te);
00278
00279 # if 0
00280 if( lgMainCO )
00281 fprintf(ioQQQ,"COcool\t%i\t%.2f\t%e\t%e\t%e\t%e\t%e\n",
00282 nUsed,
00283 fnzone,
00284 phycon.te,
00285 *Cooling/abundan,
00286 *dCoolingDT ,
00287 *Cooling,
00288 thermal.htot );
00289 # endif
00290
00291
00292 for( i=nUsed+1; i<=nRotate; ++i )
00293 {
00294 pops[i] = 0.;
00295 depart[i] = 0.;
00296 }
00297
00298
00299 for( i=0; i< MIN2(LIMLEVELN,nRotate); ++i )
00300 {
00301 atoms.PopLevels[i] = pops[i];
00302 atoms.DepLTELevels[i] = depart[i];
00303 }
00304
00305 if( nNegPop > 0 )
00306 {
00307 fprintf(ioQQQ,"CO_PopsEmisCool called atom_levelN which returned negative populations.\n");
00308 }
00309
00310
00311 for( j=0; j<nRotate; ++j )
00312 {
00313 double EnrLU, EnrUL;
00314
00315 (*Rotate)[j].PopLo = pops[j];
00316 (*Rotate)[j].PopHi = pops[j+1];
00317 (*Rotate)[j].PopOpc = (pops[j] - pops[j+1]*(*Rotate)[j].gLo/(*Rotate)[j].gHi);
00318
00319
00320 (*Rotate)[j].phots = (*Rotate)[j].Aul*((*Rotate)[j].Pesc + (*Rotate)[j].Pelec_esc)*pops[j+1];
00321
00322 # if 0
00323
00324
00325 (*Rotate)[j].ots = (*Rotate)[j].Aul*(*Rotate)[j].Pdest*(float)(*Rotate)[j].PopHi;
00326 RT_OTS_AddLine( (*Rotate)[j].ots , (*Rotate)[j].ipCont);
00327 # endif
00328
00329
00330 (*Rotate)[j].xIntensity = (*Rotate)[j].phots*(*Rotate)[j].EnergyErg;
00331
00332
00333 (*Rotate)[j].ColOvTot = (float)(CollRate[j][j+1]/(CollRate[j][j+1]+(*Rotate)[j].pump) );
00334
00335
00336
00337
00338 EnrLU = (*Rotate)[j].PopLo*CollRate[j][j+1]*(*Rotate)[j].EnergyErg;
00339 EnrUL = (*Rotate)[j].PopHi*CollRate[j+1][j]*(*Rotate)[j].EnergyErg;
00340
00341
00342 (*Rotate)[j].cool = EnrLU - EnrUL*(*Rotate)[j].ColOvTot;
00343
00344 (*Rotate)[j].heat = EnrUL*(1.f - (*Rotate)[j].ColOvTot);
00345
00346
00347 }
00348
00349
00350
00351 if( rot_cooling / thermal.ctot > 0.1 &&
00352 (*Rotate)[nUsed-1].xIntensity > (*Rotate)[nUsed-2].xIntensity &&
00353
00354
00355
00356 nUsed == nRotate )
00357 {
00358 co.lgCOCoolCaped = true;
00359 }
00360
00361 DEBUG_EXIT( "CO_PopsEmisCool()" );
00362
00363 return;
00364 }
00365
00366
00367 void CO_Colden( const char *chLabel )
00368 {
00369 long int iRot;
00370
00371 if( strcmp(chLabel,"ZERO") == 0 )
00372 {
00373 static bool lgFIRST=true;
00374 if( lgFIRST )
00375 {
00376 lgFIRST = false;
00377 col12 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) );
00378 col13 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) );
00379 }
00380
00381
00382 for( iRot=0; iRot<=nCORotate; ++iRot )
00383 {
00384
00385 col12[iRot] = 0.;
00386 col13[iRot] = 0.;
00387 }
00388 }
00389 else if( strcmp(chLabel,"ADD ") == 0 )
00390 {
00391
00392 for( iRot=0; iRot<nCORotate; ++iRot )
00393 {
00394
00395 col12[iRot] += C12O16Rotate[iRot].PopLo*radius.drad_x_fillfac;
00396 col13[iRot] += C13O16Rotate[iRot].PopLo*radius.drad_x_fillfac;
00397 }
00398 }
00399
00400
00401 else if( strcmp(chLabel,"PRIN") != 0 )
00402 {
00403 fprintf( ioQQQ, " CO_Colden does not understand the label %s\n",
00404 chLabel );
00405 puts( "[Stop in CO_Colden]" );
00406 cdEXIT(EXIT_FAILURE);
00407 }
00408 }
00409
00410
00411
00412 double cdCO_colden( long isotope , long iRot )
00413 {
00414
00415
00416 if( isotope !=12 && isotope != 13 )
00417 {
00418 fprintf(ioQQQ," cdCO_colden can only do 12CO and 13CO\n");
00419 return -1.;
00420 }
00421 if( iRot < 0 || iRot > nCORotate )
00422 {
00423 fprintf(ioQQQ," cdCO_colden - rotation quantum number must be 0 or greater, and less than %li\n",
00424 nCORotate);
00425 return -1.;
00426 }
00427
00428
00429 if( isotope == 12 )
00430 {
00431 return col12[iRot];
00432 }
00433 else
00434 {
00435 return col13[iRot];
00436 }
00437 }
00438
00439
00440 void CO_OTS( void )
00441 {
00442
00443 long int j;
00444
00445 DEBUG_ENTRY( "CO_OTS()" );
00446
00447
00448 for( j=0; j<nCORotate; ++j )
00449 {
00450 C12O16Rotate[j].ots = C12O16Rotate[j].Aul*C12O16Rotate[j].Pdest*
00451 C12O16Rotate[j].PopHi;
00452 RT_OTS_AddLine( C12O16Rotate[j].ots , C12O16Rotate[j].ipCont);
00453
00454 C13O16Rotate[j].ots = C13O16Rotate[j].Aul*C13O16Rotate[j].Pdest*
00455 C13O16Rotate[j].PopHi;
00456 RT_OTS_AddLine( C13O16Rotate[j].ots , C13O16Rotate[j].ipCont);
00457 }
00458
00459 DEBUG_EXIT( "CO_OTS()" );
00460
00461 return;
00462 }
00463