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 /*atom_level2 do level population and cooling for two level atom, 00004 * side effects: 00005 * set elements of transition struc 00006 * cooling via CoolAdd( chLab, (long)t->WLAng , t->cool); 00007 * cooling derivative */ 00008 #include "cddefines.h" 00009 #include "phycon.h" 00010 #include "lines_service.h" 00011 #include "dense.h" 00012 #include "rfield.h" 00013 #include "thermal.h" 00014 #include "cooling.h" 00015 #include "atoms.h" 00016 00017 void atom_level2(transition * t) 00018 { 00019 char chLab[5]; 00020 long int ion, 00021 ip, 00022 nelem; 00023 double AbunxIon, 00024 a21, 00025 boltz, 00026 col12, 00027 col21, 00028 coolng, 00029 g1, 00030 g2, 00031 omega, 00032 pfs1, 00033 pfs2, 00034 plower, 00035 r, 00036 rate12, 00037 ri21; 00038 00039 DEBUG_ENTRY( "atom_level2()" ); 00040 00041 /* result is density (cm-3) of excited state times A21 00042 * result normalized to N1+N2=ABUND 00043 * routine increments dCooldT, call CoolAdd 00044 * CDSQTE is EDEN / SQRTE * 8.629E-6 00045 */ 00046 00047 ion = t->Hi->IonStg; 00048 nelem = t->Hi->nelem; 00049 00050 /* dense.xIonDense[nelem][i] is density of ith ionization stage (cm^-3) */ 00051 AbunxIon = dense.xIonDense[nelem-1][ion-1]; 00052 00053 /* continuum pointer */ 00054 ip = t->ipCont; 00055 00056 /* approximate Boltzmann factor to see if results zero */ 00057 boltz = rfield.ContBoltz[ip-1]; 00058 00059 /* collision strength for this transition, omega is zero for hydrogenic 00060 * species which are done in special hydro routines */ 00061 omega = t->Coll.cs; 00062 00063 /* ROUGH check whether upper level has significant population,*/ 00064 r = (boltz*dense.cdsqte + t->Emis->pump)/(dense.cdsqte + t->Emis->Aul); 00065 00066 /* following first test needed for 32 bit cpu on search phase 00067 * >>chng 96 jul 02, was 1e-30 crashed on dec, change to 1e-25 00068 * if( AbunxIon.lt.1e-25 .or. boltz.gt.30. ) then 00069 * >>chng 96 jul 11, to below, since can be strong pumping when 00070 * Boltzmann factor essentially zero */ 00071 /* omega in following is zero for hydrogenic species, since done 00072 * in hydro routines, so this should cause us to quit on this test */ 00073 /* >>chng 99 nov 29, from 1e-25 to 1e-30 to keep same result for 00074 * very low density models, where AbunxIon is very small but still significant*/ 00075 /*if( omega*AbunxIon < 1e-25 || r < 1e-25 )*/ 00076 if( omega*AbunxIon < 1e-30 || r < 1e-25 ) 00077 { 00078 /* put in pop since possible just too cool */ 00079 t->Lo->Pop = AbunxIon; 00080 t->Emis->PopOpc = AbunxIon; 00081 t->Hi->Pop = 0.; 00082 t->Emis->xIntensity = 0.; 00083 t->Coll.cool = 0.; 00084 t->Emis->ots = 0.; 00085 t->Emis->phots = 0.; 00086 t->Emis->ColOvTot = 0.; 00087 t->Coll.heat = 0.; 00088 /* level populations */ 00089 atoms.PopLevels[0] = AbunxIon; 00090 atoms.PopLevels[1] = 0.; 00091 atoms.DepLTELevels[0] = 1.; 00092 atoms.DepLTELevels[1] = 0.; 00093 return; 00094 } 00095 00096 /* net rate down A21*(escape + destruction) */ 00097 a21 = t->Emis->Aul*(t->Emis->Pesc+ t->Emis->Pdest + t->Emis->Pelec_esc); 00098 00099 /* put null terminated line label into chLab using line array*/ 00100 chIonLbl(chLab,t); 00101 00102 /* statistical weights of upper and lower levels */ 00103 g1 = t->Lo->g; 00104 g2 = t->Hi->g; 00105 00106 /* now get real Boltzmann factor */ 00107 boltz = t->EnergyK/phycon.te; 00108 00109 ASSERT( boltz > 0. ); 00110 boltz = sexp(boltz); 00111 00112 ASSERT( g1 > 0. && g2 > 0. ); 00113 00114 /* this lacks the upper statistical weight */ 00115 col21 = dense.cdsqte*omega; 00116 /* upward coll rate s-1 */ 00117 col12 = col21/g1*boltz; 00118 /* downward coll rate s-1 */ 00119 col21 /= g2; 00120 00121 /* rate 1 to 2 is both collisions and pumping */ 00122 /* the total excitation rate from lower to upper, collisional and pumping */ 00123 rate12 = col12 + t->Emis->pump; 00124 00125 /* induced emissions down */ 00126 ri21 = t->Emis->pump*g1/g2; 00127 00128 /* this is the ratio of lower to upper level populations */ 00129 r = (a21 + col21 + ri21)/rate12; 00130 00131 /* upper level pop */ 00132 pfs2 = AbunxIon/(r + 1.); 00133 atoms.PopLevels[1] = pfs2; 00134 t->Hi->Pop = pfs2; 00135 00136 /* pop of ground */ 00137 pfs1 = pfs2*r; 00138 atoms.PopLevels[0] = pfs1; 00139 00140 /* compute ratio Aul/(Aul+Cul) */ 00141 /* level population with correction for stimulated emission */ 00142 t->Lo->Pop = atoms.PopLevels[0]; 00143 00144 t->Emis->PopOpc = (atoms.PopLevels[0] - atoms.PopLevels[1]*g1/g2 ); 00145 00146 /* departure coef of excited state rel to ground */ 00147 atoms.DepLTELevels[0] = 1.; 00148 if( boltz > 1e-20 && atoms.PopLevels[1] > 1e-20 ) 00149 { 00150 /* this line could have zero boltz factor but radiatively excited 00151 * dec alpha does not obey () in fast mode?? */ 00152 atoms.DepLTELevels[1] = (atoms.PopLevels[1]/atoms.PopLevels[0])/ 00153 (boltz*g2/g1); 00154 } 00155 else 00156 { 00157 atoms.DepLTELevels[1] = 0.; 00158 } 00159 00160 /* number of escaping line photons, used elsewhere for outward beam */ 00161 t->Emis->phots = t->Emis->Aul*(t->Emis->Pesc + t->Emis->Pelec_esc)*pfs2; 00162 00163 /* intensity of line */ 00164 t->Emis->xIntensity = t->Emis->phots*t->EnergyErg; 00165 plower = AbunxIon - pfs2; 00166 00167 /* ratio of collisional to total (collisional + pumped) excitation */ 00168 t->Emis->ColOvTot = col12/rate12; 00169 00170 /* two cases - collisionally excited (usual case) or 00171 * radiatively excited - in which case line can be a heat source 00172 * following are correct heat exchange, they will mix to get correct deriv 00173 * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to 00174 * keep stable solution by effectively dividing up heating and cooling, 00175 * so that negative cooling does not occur */ 00176 00177 //double Enr12 = plower*col12*t->EnergyErg; 00178 //double Enr21 = pfs2*col21*t->EnergyErg; 00179 00180 /* energy exchange due to this level 00181 * net cooling due to excit minus part of de-excit - 00182 * note that ColOvTot cancels out in the sum heat - cool */ 00183 //coolng = Enr12 - Enr21*t->Emis->ColOvTot; 00184 00185 /* this form of coolng is guaranteed to be non-negative */ 00186 coolng = t->EnergyErg*AbunxIon*col12*(a21 + ri21)/(a21 + col21 + ri21 + rate12); 00187 ASSERT( coolng >= 0. ); 00188 00189 t->Coll.cool = coolng; 00190 00191 /* net heating is remainder */ 00192 //t->Coll.heat = Enr21*(1. - t->Emis->ColOvTot); 00193 t->Coll.heat = t->EnergyErg*AbunxIon*col21*(t->Emis->pump)/(a21 + col21 + ri21 + rate12); 00194 00195 00196 /* expression pre jul 3 95, changed for case where line heating dominates 00197 * coolng = (plower*col12 - pfs2*col21)*t->t(ipLnEnrErg) 00198 * t->t(ipLnCool) = cooling */ 00199 00200 /* add to cooling - heating part was taken out above, 00201 * and is not added in here - it will be added to thermal.heating[0][22] 00202 * in CoolSum */ 00203 CoolAdd( chLab, t->WLAng , t->Coll.cool); 00204 00205 /* derivative of cooling function */ 00206 thermal.dCooldT += coolng * (t->EnergyK * thermal.tsq1 - thermal.halfte ); 00207 return; 00208 }