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