cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*CO_PopsEmisCool evaluate rotation levels populations, emission, and cooling */ 00004 /*cdCO_colden return column density in H2, negative -1 if cannot find state, 00005 * header is cddrive */ 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 /*lint -e778 constant expression evaluatess to 0 in operation '-' */ 00019 /*=================================================================*/ 00020 00021 /* will be used to save CO column densities */ 00022 static double *col12 , *col13; 00023 00024 /* evaluate CO rotation cooling */ 00025 void CO_PopsEmisCool( 00026 transition ** Rotate , 00027 long int nRotate , 00028 realnum abundan, 00029 const char * chLabel , 00030 realnum * Cooling , 00031 realnum * dCoolingDT ) 00032 { 00033 00034 /* will need to MALLOC space for these but only on first call */ 00035 static double 00036 **AulEscp , 00037 **col_str , 00038 **AulDest, 00039 /* AulPump[low][high] is rate (s^-1) from lower to upper level */ 00040 **AulPump, 00041 **CollRate, 00042 *pops, 00043 *create, 00044 *destroy, 00045 *depart, 00046 /* statistical weight */ 00047 *stat , 00048 /* excitation energies in kelvin */ 00049 *excit; 00050 00051 /*static long int **ipdest;*/ 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 /* will never do this again */ 00070 lgFirst = false; 00071 /* remember how much space we malloced in case ever called with more needed */ 00072 ndimMalloced = nRotate; 00073 /* allocate the 1D arrays*/ 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 /* create space for the 2D arrays */ 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 /* this is test for call with too many rotation levels to handle - logic needs 00096 * for largest rotor to be called first */ 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 cdEXIT(EXIT_FAILURE); 00101 } 00102 00103 /* all elements are used, and must be set to zero if zero */ 00104 for( i=0; i < (nRotate+1); i++ ) 00105 { 00106 create[i] = 0.; 00107 destroy[i] = 0.; 00108 for( j=0; j < (nRotate+1); j++ ) 00109 { 00110 col_str[j][i] = 0.; 00111 AulEscp[j][i] = 0.; 00112 AulDest[j][i] = 0.; 00113 AulPump[j][i] = 0.; 00114 } 00115 } 00116 00117 /* the statistical weights of the levels */ 00118 for( j=0; j < nRotate; j++ ) 00119 { 00120 /* statistical weights for each level */ 00121 stat[j] = (*Rotate)[j].Lo->g; 00122 } 00123 /* this is the highest level, which is one more than the highest line */ 00124 stat[nRotate] = (*Rotate)[nRotate-1].Hi->g; 00125 00126 /* set up the excitation potentials of each level relative to ground - 00127 * the struc saves the energy of the line only */ 00128 excit[0] = 0.; 00129 for( j=1; j < nRotate; j++ ) 00130 { 00131 /* excitation energy of each level relative to ground, in K */ 00132 excit[j] = excit[j-1] + (*Rotate)[j-1].EnergyK; 00133 } 00134 00135 /* this is the highest level, which is one more than the highest line */ 00136 excit[nRotate] = excit[nRotate-1] + (*Rotate)[nRotate-1].EnergyK; 00137 00138 nUsed = nRotate; 00139 00140 /* this determines the largest molecule that can be inverted at this 00141 * temperature. Need Boltzmann factors to be positive for all levels 00142 * nUsed is the index of the highest level, so the number of levels (passed 00143 * to the solver) is nUsed+1 00144 * excit[nUsed]/phycon.te cannot be much larger than 20, or matrix inversion will fail */ 00145 /* >>chng 03 sep 18, allow to go to 1, a two level atom */ 00146 /*while ( (excit[nUsed] > phycon.te*20.) && (nUsed > 1) )*/ 00147 /* >>chng 03 oct 03, chng 20 to 25 */ 00148 /*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/ 00149 /* >>chng 03 oct 03, keep at least 5 levels */ 00150 /*while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 1) )*/ 00151 while ( (excit[nUsed] > phycon.te*25.) && (nUsed > 5) ) 00152 { 00153 --nUsed; 00154 } 00155 00156 for( j=0; j < nRotate; j++ ) 00157 { 00158 /*data[j][j+1] = (*Rotate)[j].Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc);*/ 00159 AulEscp[j+1][j] = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc); 00160 AulDest[j+1][j] = (*Rotate)[j].Emis->Aul*(*Rotate)[j].Emis->Pdest; 00161 /* next 2 not not used by levelN since flag passed saying to use CollRate instead */ 00162 (*Rotate)[j].Coll.cs = 1.; 00163 /*data[j+1][j] = (*Rotate)[j].Coll.cs*hmi.Hmolec[ipMH2g]/dense.eden;*/ 00164 col_str[j+1][j] = (*Rotate)[j].Coll.cs*hmi.H2_total/dense.eden; 00165 /* photon pumping rate */ 00166 AulPump[j][j+1] = (*Rotate)[j].Emis->pump; 00167 /* the continuum indices are on the f, not c, scale, and will be passed to 00168 * atom_levelN, which works on f, not c, scale */ 00169 /*ipdest[j][j+1] = (*Rotate)[j].ipCont;*/ 00170 } 00171 00172 /* next loop for collisional transitions that have delta J >1 */ 00173 for( ilo=0; ilo < nRotate-1; ilo++ ) 00174 { 00175 /* need to have upper limit to to nRotate because 00176 * number ov levels is 1 gt number of lines*/ 00177 for( ihi=ilo+2; ihi <= nRotate; ihi++ ) 00178 { 00179 /* this is not used by levelN since flag passed saying to use CollRate instead */ 00180 /*data[ihi][ilo] = 1. *hmi.Hmolec[ipMH2g]/dense.eden;*/ 00181 /*data[ihi][ilo] = 1. *hmi.H2_total/dense.eden;*/ 00182 col_str[ihi][ilo] = 1. *hmi.H2_total/dense.eden; 00183 /* these are escape and dest rates, which are zero for a rigid rotator */ 00184 /*data[ilo][ihi] = 0;*/ 00185 AulEscp[ihi][ilo] = 0; 00186 AulDest[ihi][ilo] = 0; 00187 } 00188 } 00189 00190 /* now evaluate the H2 collision rates */ 00191 /* recall one more level than lines */ 00192 for( ilo=0; ilo < nRotate; ilo++ ) 00193 { 00194 /* need to have upper limit to to nRotate because 00195 * number of levels is 1 gt number of lines*/ 00196 for( ihi=ilo+1; ihi <= MIN2(nRotate,ilo+5); ihi++ ) 00197 { 00198 /* >>refer CO-H2 collision de Jong, T., Chu, S-I., & Dalgarno, A. 1975, ApJ, 199, 69 */ 00199 double a[5]={1.66, 2.80, 1.19, 1.00, 1.30}; 00200 double b[5]={1.67, 1.47, 1.85, 1.55, 2.24}; 00201 /* >>refer CO-He collision McKee, C.F., Storey, J.W.V., Watson, D.M., & Green, S., 1982, ApJ, 259, 647 */ 00202 /* this reference gives He-CO collisions, 00203 * their table 1 says He collisions are 1.37 slower than H2 collisions*/ 00204 /* >>chng 03 sep 19, from ground H2 to total H2 */ 00205 double collid = hmi.H2_total + dense.xIonDense[ipHELIUM][0]/1.37; 00206 /* de Jong et al. explicitly consider temperatures as low as 20K, 00207 * don't extrapolate significantly lower than this */ 00208 double TeUsed = MAX2(10., phycon.te ); 00209 00210 /* first do deexcitation rate, equation 17 of deJong et al. */ 00211 CollRate[ihi][ilo] = a[ihi-ilo-1]*1.e-10*(*Rotate)[ilo].Lo->g/(*Rotate)[ilo].Hi->g* 00212 (1.+(excit[ihi]-excit[ilo])/TeUsed) * 00213 sexp( b[ihi-ilo-1]*sqrt((excit[ihi]-excit[ilo])/TeUsed) )*collid; 00214 /* this is mainly so that rest of code gets valid cs in case crit dense printed */ 00215 if( ihi == ilo+1 ) 00216 { 00217 /* save downward collision rate */ 00218 (*Rotate)[ilo].Coll.ColUL = (realnum)CollRate[ihi][ilo]; 00219 /* convert into fake electron collision strength */ 00220 LineConvRate2CS( &(*Rotate)[ilo] , (*Rotate)[ilo].Coll.ColUL); 00221 } 00222 00223 /* now get excitation rate */ 00224 CollRate[ilo][ihi] = CollRate[ihi][ilo]* 00225 sexp( (excit[ihi]-excit[ilo])/phycon.te )* 00226 (*Rotate)[ilo].Hi->g / (*Rotate)[ilo].Lo->g; 00227 00228 /* debug print statement */ 00229 /*fprintf(ioQQQ," %li %li %.2e %.2e \n",ilo, ihi, 00230 CollRate[ihi][ilo]/hmi.H2_total , CollRate[ilo][ihi]/hmi.H2_total );*/ 00231 } 00232 /* finish off with zeros */ 00233 for( ihi=ilo+6; ihi <= nRotate; ihi++ ) 00234 { 00235 CollRate[ihi][ilo] = 0.; 00236 CollRate[ilo][ihi] = 0.; 00237 } 00238 } 00239 00240 lgDeBug = false; 00241 00242 atom_levelN( 00243 /* number of levels is number of lines plus one */ 00244 /* set nUsed so that CO is evaluated even at very low temperatures */ 00245 nUsed+1, /*nRotate+1,*/ 00246 abundan, 00247 stat, 00248 excit, 00249 'K', 00250 pops, 00251 depart, 00252 &AulEscp, 00253 &col_str, 00254 &AulDest, 00255 &AulPump, 00256 &CollRate, 00257 create, 00258 destroy, 00259 /* say that we have evaluated the collision rates already */ 00260 true, 00261 /*&ipdest,*/ 00262 &rot_cooling, 00263 &dCoolDT, 00264 chLabel, 00265 /* nNegPop positive if negative pops occured, negative if too cold */ 00266 &nNegPop, 00267 &lgZeroPop, 00268 lgDeBug );/* option to print stuff - set to true for debug printout */ 00269 00270 /* put cooling into place where we can use it later */ 00271 *Cooling = (realnum)rot_cooling; 00272 00273 /* >>chng 03 sep 18, CO rot cooling temp deriv is too small, and temp 00274 * changes too big - incr deriv by 5x 00275 *dCoolingDT = (realnum)(dCoolDT);*/ 00276 *dCoolingDT = (realnum)(rot_cooling/phycon.te); 00277 00278 # if 0 00279 if( lgMainCO ) 00280 fprintf(ioQQQ,"COcool\t%i\t%.2f\t%e\t%e\t%e\t%e\t%e\n", 00281 nUsed, 00282 fnzone, 00283 phycon.te, 00284 *Cooling/abundan, 00285 *dCoolingDT , 00286 *Cooling, 00287 thermal.htot ); 00288 # endif 00289 00290 /* zero out higher populations for case where full CO levels not done */ 00291 for( i=nUsed+1; i<=nRotate; ++i ) 00292 { 00293 pops[i] = 0.; 00294 depart[i] = 0.; 00295 } 00296 00297 /* can only define first LIMLEVELN elements, the vector's length */ 00298 for( i=0; i< MIN2(LIMLEVELN,nRotate); ++i ) 00299 { 00300 atoms.PopLevels[i] = pops[i]; 00301 atoms.DepLTELevels[i] = depart[i]; 00302 } 00303 00304 if( nNegPop > 0 ) 00305 { 00306 fprintf(ioQQQ,"CO_PopsEmisCool called atom_levelN which returned negative populations.\n"); 00307 } 00308 00309 /* now establish information that is passed out to rest of code's infrastructure */ 00310 for( j=0; j<nRotate; ++j ) 00311 { 00312 double EnrLU, EnrUL; 00313 /* lower upper populations, stim emission correct population */ 00314 (*Rotate)[j].Lo->Pop = pops[j]; 00315 (*Rotate)[j].Hi->Pop = pops[j+1]; 00316 (*Rotate)[j].Emis->PopOpc = (pops[j] - pops[j+1]*(*Rotate)[j].Lo->g/(*Rotate)[j].Hi->g); 00317 00318 /* number of photons in the line */ 00319 (*Rotate)[j].Emis->phots = (*Rotate)[j].Emis->Aul*((*Rotate)[j].Emis->Pesc + (*Rotate)[j].Emis->Pelec_esc)*pops[j+1]; 00320 00321 # if 0 00322 /* >>chng 03 oct 04, moved to RT_OTS */ 00323 /* local Emis.ots rates, added to line ots array */ 00324 (*Rotate)[j].Emis->ots = (*Rotate)[j].Aul*(*Rotate)[j].Emis->Pdest*(realnum)(*Rotate)[j].Hi->Pop; 00325 RT_OTS_AddLine( (*Rotate)[j].Emis->ots , (*Rotate)[j].ipCont); 00326 # endif 00327 00328 /* the intensity in the line */ 00329 (*Rotate)[j].Emis->xIntensity = (*Rotate)[j].Emis->phots*(*Rotate)[j].EnergyErg; 00330 00331 /* ratio of collisional to total excitation */ 00332 (*Rotate)[j].Emis->ColOvTot = (realnum)(CollRate[j][j+1]/(CollRate[j][j+1]+(*Rotate)[j].Emis->pump) ); 00333 00334 /* two cases - collisionally excited (usual case) or 00335 * radiatively excited - in which case line can be a heat source 00336 * following are correct heat exchange, if will mix to get correct deriv */ 00337 EnrLU = (*Rotate)[j].Lo->Pop*CollRate[j][j+1]*(*Rotate)[j].EnergyErg; 00338 EnrUL = (*Rotate)[j].Hi->Pop*CollRate[j+1][j]*(*Rotate)[j].EnergyErg; 00339 /* energy exchange due to this level 00340 * net cooling due to excit minus part of de-excit */ 00341 (*Rotate)[j].Coll.cool = EnrLU - EnrUL*(*Rotate)[j].Emis->ColOvTot; 00342 /* net heating is remainder */ 00343 (*Rotate)[j].Coll.heat = EnrUL*(1.f - (*Rotate)[j].Emis->ColOvTot); 00344 /* do not add to cooling, since done with evaluated cooling from atom_levelN */ 00345 /*CoolAdd( chLabel, (long)t10->WLAng , t10->cool);*/ 00346 } 00347 00348 /* generate flag if co cooling important and highest level is fainter 00349 * than second highest level */ 00350 if( rot_cooling / thermal.ctot > 0.1 && 00351 (*Rotate)[nUsed-1].Emis->xIntensity > (*Rotate)[nUsed-2].Emis->xIntensity && 00352 /*>>chng 03 oct 03, add check whether molecule has been trimed down 00353 * due to small Boltzmann factor */ 00354 /* >>chng 03 oct 20, following had; after ), so CoolCaped always true */ 00355 nUsed == nRotate ) 00356 { 00357 co.lgCOCoolCaped = true; 00358 } 00359 00360 return; 00361 } 00362 00363 /*CO_Colden maintain H2 column densities within X */ 00364 void CO_Colden( const char *chLabel ) 00365 { 00366 long int iRot; 00367 00368 DEBUG_ENTRY( "CO_Colden()" ); 00369 00370 if( strcmp(chLabel,"ZERO") == 0 ) 00371 { 00372 static bool lgFIRST=true; 00373 if( lgFIRST ) 00374 { 00375 lgFIRST = false; 00376 col12 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) ); 00377 col13 = (double *)MALLOC( sizeof(double)*(size_t)(nCORotate+1) ); 00378 } 00379 /* the column density (cm-2) of ortho and para H2 */ 00380 /* zero out formation rates and column densites */ 00381 for( iRot=0; iRot<=nCORotate; ++iRot ) 00382 { 00383 /* zero it out */ 00384 col12[iRot] = 0.; 00385 col13[iRot] = 0.; 00386 } 00387 } 00388 else if( strcmp(chLabel,"ADD ") == 0 ) 00389 { 00390 /* add together column densities */ 00391 for( iRot=0; iRot<nCORotate; ++iRot ) 00392 { 00393 /* zero it out */ 00394 col12[iRot] += C12O16Rotate[iRot].Lo->Pop*radius.drad_x_fillfac; 00395 col13[iRot] += C13O16Rotate[iRot].Lo->Pop*radius.drad_x_fillfac; 00396 } 00397 } 00398 00399 /* we will not print column densities so skip that - if not print then we have a problem */ 00400 else if( strcmp(chLabel,"PRIN") != 0 ) 00401 { 00402 fprintf( ioQQQ, " CO_Colden does not understand the label %s\n", 00403 chLabel ); 00404 cdEXIT(EXIT_FAILURE); 00405 } 00406 } 00407 00408 /*cdCO_colden return column density in H2, negative -1 if cannot find state, 00409 * header is cddrive */ 00410 double cdCO_colden( long isotope , long iRot ) 00411 { 00412 00413 /* make sure incoming parameters are ok */ 00414 if( isotope !=12 && isotope != 13 ) 00415 { 00416 fprintf(ioQQQ," cdCO_colden can only do 12CO and 13CO\n"); 00417 return -1.; 00418 } 00419 if( iRot < 0 || iRot > nCORotate ) 00420 { 00421 fprintf(ioQQQ," cdCO_colden - rotation quantum number must be 0 or greater, and less than %li\n", 00422 nCORotate); 00423 return -1.; 00424 } 00425 00426 /* the incoming parameters are fully qualified - return the column density */ 00427 if( isotope == 12 ) 00428 { 00429 return col12[iRot]; 00430 } 00431 else 00432 { 00433 return col13[iRot]; 00434 } 00435 } 00436 00437 /*CO_OTS - add CO lines to ots fields */ 00438 void CO_OTS( void ) 00439 { 00440 00441 long int j; 00442 00443 DEBUG_ENTRY( "CO_OTS()" ); 00444 00445 /* add all CO lines */ 00446 for( j=0; j<nCORotate; ++j ) 00447 { 00448 C12O16Rotate[j].Emis->ots = C12O16Rotate[j].Emis->Aul*C12O16Rotate[j].Emis->Pdest* 00449 C12O16Rotate[j].Hi->Pop; 00450 RT_OTS_AddLine( C12O16Rotate[j].Emis->ots , C12O16Rotate[j].ipCont); 00451 00452 C13O16Rotate[j].Emis->ots = C13O16Rotate[j].Emis->Aul*C13O16Rotate[j].Emis->Pdest* 00453 C13O16Rotate[j].Hi->Pop; 00454 RT_OTS_AddLine( C13O16Rotate[j].Emis->ots , C13O16Rotate[j].ipCont); 00455 } 00456 00457 return; 00458 }