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 /*lines_general put general information and energetics into line intensity stack */ 00004 /*GetMaxhLine find the strongest heating line */ 00005 #include "cddefines.h" 00006 #include "taulines.h" 00007 #include "coolheavy.h" 00008 #include "hydrogenic.h" 00009 #include "dense.h" 00010 #include "thermal.h" 00011 #include "continuum.h" 00012 #include "geometry.h" 00013 #include "dynamics.h" 00014 #include "rt.h" 00015 #include "iso.h" 00016 #include "rfield.h" 00017 #include "trace.h" 00018 #include "ionbal.h" 00019 #include "lines_service.h" 00020 #include "radius.h" 00021 #include "lines.h" 00022 /*GetMaxhLine find the strongest heating line */ 00023 STATIC void GetMaxhLine(void); 00024 00025 void lines_general(void) 00026 { 00027 long int i, 00028 ipHi, 00029 ipLo, 00030 nelem, 00031 ipnt; 00032 00033 double 00034 hbetac, 00035 HeatMetal , 00036 ee511, 00037 hlalph; 00038 00039 DEBUG_ENTRY( "lines_general()" ); 00040 00041 if( trace.lgTrace ) 00042 { 00043 fprintf( ioQQQ, " lines_general called\n" ); 00044 } 00045 00046 i = StuffComment( "general properties" ); 00047 linadd( 0., (realnum)i , "####", 'i', 00048 " start of general properties"); 00049 00050 /* total H-beta from multi-level atom */ 00051 nelem = ipHYDROGEN; 00052 ipHi = ipH4p; 00053 ipLo = ipH2p; 00054 00055 if( iso.n_HighestResolved_max[ipH_LIKE][nelem] >= 4 ) 00056 { 00057 hbetac = 00058 (Transitions[ipH_LIKE][nelem][ipH4p][ipH2s].Emis->Aul * 00059 Transitions[ipH_LIKE][nelem][ipH4p][ipH2s].Emis->Pesc * 00060 StatesElem[ipH_LIKE][nelem][ipH4p].Pop + 00061 Transitions[ipH_LIKE][nelem][ipH4s][ipH2p].Emis->Aul * 00062 Transitions[ipH_LIKE][nelem][ipH4s][ipH2p].Emis->Pesc * 00063 StatesElem[ipH_LIKE][nelem][ipH4s].Pop + 00064 Transitions[ipH_LIKE][nelem][ipH4d][ipH2p].Emis->Aul * 00065 Transitions[ipH_LIKE][nelem][ipH4d][ipH2p].Emis->Pesc * 00066 StatesElem[ipH_LIKE][nelem][ipH4d].Pop) * 00067 Transitions[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg* 00068 dense.xIonDense[ipHYDROGEN][1]; 00069 } 00070 else 00071 { 00072 hbetac = 00073 (Transitions[ipH_LIKE][nelem][ipHi][ipH2s].Emis->Aul * 00074 Transitions[ipH_LIKE][nelem][ipHi][ipH2s].Emis->Pesc + 00075 Transitions[ipH_LIKE][nelem][ipHi][ipH2p].Emis->Aul * 00076 Transitions[ipH_LIKE][nelem][ipHi][ipH2p].Emis->Pesc ) * 00077 StatesElem[ipH_LIKE][nelem][ipHi].Pop * 00078 Transitions[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg* 00079 dense.xIonDense[ipHYDROGEN][1]; 00080 } 00081 00082 /* these lines added to outlin in metdif - following must be false 00083 * this passes array index for line energy in continuum mesh - in rest 00084 * of code this is set by a previous call to PntForLine, this index 00085 * is on the f not c scale */ 00086 rt.fracin = Transitions[ipH_LIKE][nelem][ipHi][ipH2s].Emis->FracInwd; 00087 lindst(hbetac,4861,"TOTL",Transitions[ipH_LIKE][nelem][ipHi][ipH2s].ipCont,'i',false," " ); 00088 rt.fracin = 0.5; 00089 00090 /* total Ly-a from multi-level atom */ 00091 ipHi = ipH2p; 00092 ipLo = ipH1s; 00093 hlalph = 00094 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul* 00095 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Hi->Pop* 00096 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pesc* 00097 Transitions[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg* 00098 dense.xIonDense[ipHYDROGEN][1]; 00099 00100 rt.fracin = Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->FracInwd; 00101 lindst(hlalph,1216,"TOTL",Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont,'i',false , 00102 ""); 00103 rt.fracin = 0.5; 00104 00105 continuum.totlsv = continuum.totlsv/radius.dVeff*geometry.covgeo; 00106 00107 /* H 21 cm, A from 00108 * >>refer H1 A Gould, ApJ 423, 522 00109 * =2.88426e-15 s-1 */ 00110 00111 linadd(continuum.totlsv,0,"Inci",'i', 00112 "total luminosity in incident continuum"); 00113 continuum.totlsv = continuum.totlsv*radius.dVeff/geometry.covgeo; 00114 00115 /* ipass is flag to indicate whether to only set up line array 00116 * (ipass=0) or actually evaluate lines intensities (ipass=1) */ 00117 if( LineSave.ipass > 0 ) 00118 { 00119 continuum.totlsv = 0.; 00120 } 00121 00122 linadd(thermal.htot,0,"TotH",'i', 00123 " total heating, all forms, information since individuals added later "); 00124 00125 linadd(thermal.ctot,0,"TotC",'i', 00126 " total cooling, all forms, information since individuals added later "); 00127 00128 linadd(thermal.heating[0][0],0,"BFH1",'h', 00129 " hydrogen photoionization heating, ground state only "); 00130 00131 linadd(thermal.heating[0][1],0,"BFHx",'h', 00132 " net hydrogen photoionization heating less rec cooling, all excited states normally zero, positive if excited states are net heating "); 00133 00134 linadd(thermal.heating[0][22],0,"Line",'h', 00135 " heating due to induced lines absorption of continuum "); 00136 if( thermal.htot > 0. ) 00137 { 00138 if( thermal.heating[0][22]/thermal.htot > thermal.HeatLineMax ) 00139 { 00140 thermal.HeatLineMax = (realnum)(thermal.heating[0][22]/thermal.htot); 00141 /* finds the strongest heating line */ 00142 GetMaxhLine(); 00143 } 00144 } 00145 00146 linadd(thermal.heating[1][0]+thermal.heating[1][1]+thermal.heating[1][2],0,"BFHe",'h', 00147 " total helium photoionization heating, all stages "); 00148 00149 HeatMetal = 0.; 00150 /* some sums that will be printed in the stack */ 00151 for( nelem=2; nelem<LIMELM; ++nelem) 00152 { 00153 /* we now have final solution for this element */ 00154 for( i=dense.IonLow[nelem]; i < dense.IonHigh[nelem]; i++ ) 00155 { 00156 ASSERT( i < LIMELM ); 00157 /* total metal photo heating for LINES */ 00158 HeatMetal += thermal.heating[nelem][i]; 00159 } 00160 } 00161 00162 linadd(HeatMetal,0,"TotM",'h', 00163 " total heavy element photoionization heating, all stages "); 00164 00165 linadd(thermal.heating[0][21],0,"pair",'h', 00166 " heating due to pair production "); 00167 00168 /* ipass is flag to indicate whether to only set up line array 00169 * (ipass=0) or actually evaluate lines intensities (ipass=1) */ 00170 if( LineSave.ipass > 0 ) 00171 { 00172 /* this will be max local heating due to bound compton */ 00173 ionbal.CompHeating_Max = MAX2( ionbal.CompHeating_Max , ionbal.CompRecoilHeatLocal/thermal.htot); 00174 } 00175 else 00176 { 00177 ionbal.CompHeating_Max = 0.; 00178 } 00179 00180 linadd(ionbal.CompRecoilHeatLocal,0,"Cbnd",'h', 00181 " heating due to bound compton scattering "); 00182 00183 linadd(rfield.cmheat,0,"ComH",'h', 00184 " Compton heating "); 00185 00186 linadd(CoolHeavy.tccool,0,"ComC",'c', 00187 " total Compton cooling "); 00188 00189 /* record max local heating due to advection */ 00190 dynamics.HeatMax = MAX2( dynamics.HeatMax , dynamics.Heat /thermal.htot ); 00191 /* record max local cooling due to advection */ 00192 dynamics.CoolMax = MAX2( dynamics.CoolMax , dynamics.Cool /thermal.htot ); 00193 00194 linadd(dynamics.Cool , 0 , "advC" , 'i', 00195 " cooling due to advection " ); 00196 00197 linadd(dynamics.Heat , 0 , "advH" , 'i' , 00198 " heating due to advection "); 00199 00200 linadd( thermal.char_tran_heat ,0,"CT H",'h', 00201 " heating due to charge transfer "); 00202 00203 linadd( thermal.char_tran_cool ,0,"CT C",'c', 00204 " cooling due to charge transfer "); 00205 00206 linadd(thermal.heating[1][6],0,"CR H",'h', 00207 " cosmic ray heating "); 00208 00209 linadd(thermal.heating[0][20],0,"extH",'h', 00210 " extra heat added to this zone, from HEXTRA command "); 00211 00212 linadd(CoolHeavy.cextxx,0,"extC",'c', 00213 " extra cooling added to this zone, from CEXTRA command "); 00214 00215 00216 ee511 = (dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM])*ionbal.PairProducPhotoRate[0]*2.*8.20e-7; 00217 PntForLine(2.427e-2,"e-e+",&ipnt); 00218 lindst(ee511,0,"e-e+",ipnt,'i',true, 00219 " 511keV annihilation line " ); 00220 00221 linadd(CoolHeavy.expans,0,"Expn",'c', 00222 " expansion cooling, only non-zero for wind "); 00223 00224 linadd(iso.RadRecCool[ipH_LIKE][ipHYDROGEN],0,"H FB",'i', 00225 " H radiative recombination cooling "); 00226 00227 linadd(MAX2(0.,iso.FreeBnd_net_Cool_Rate[ipH_LIKE][ipHYDROGEN]),0,"HFBc",'c', 00228 " net free-bound cooling "); 00229 00230 linadd(MAX2(0.,-iso.FreeBnd_net_Cool_Rate[ipH_LIKE][ipHYDROGEN]),0,"HFBh",'h', 00231 " net free-bound heating "); 00232 00233 linadd(iso.RecomInducCool_Rate[ipH_LIKE][ipHYDROGEN],0,"Hind",'c', 00234 " cooling due to induced rec of hydrogen "); 00235 00236 linadd(CoolHeavy.c3ind,0,"3He2",'c', 00237 " cooling due to induced rec of fully ionized helium "); 00238 00239 linadd(CoolHeavy.cyntrn,0,"Cycn",'c', 00240 " cyclotron cooling "); 00241 return; 00242 } 00243 00244 /*GetMaxhLine find the strongest heating line */ 00245 STATIC void GetMaxhLine(void) 00246 { 00247 long int i; 00248 double strong; 00249 00250 DEBUG_ENTRY( "GetMaxhLine()" ); 00251 00252 /* routine called to find which is the strongest heating line */ 00253 strong = 0.; 00254 00255 /* possible for levlmax to remain 0 if induced processes turned off */ 00256 thermal.levlmax = 0; 00257 00258 for( i=1; i <= nLevel1; i++ ) 00259 { 00260 if( TauLines[i].Coll.heat > strong ) 00261 { 00262 strong = TauLines[i].Coll.heat; 00263 thermal.levlmax = 1; 00264 thermal.ipHeatlmax = i; 00265 } 00266 } 00267 00268 for( i=0; i < nWindLine; i++ ) 00269 { 00270 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00271 { 00272 if( TauLine2[i].Coll.heat > strong ) 00273 { 00274 strong = TauLine2[i].Coll.heat; 00275 thermal.levlmax = 2; 00276 thermal.ipHeatlmax = i; 00277 } 00278 } 00279 } 00280 00281 for( i=0; i < nHFLines; i++ ) 00282 { 00283 if( HFLines[i].Coll.heat > strong ) 00284 { 00285 strong = HFLines[i].Coll.heat; 00286 thermal.levlmax = 3; 00287 thermal.ipHeatlmax = i; 00288 } 00289 } 00290 00291 for( i=0; i < nCORotate; i++ ) 00292 { 00293 if( C12O16Rotate[i].Coll.heat > strong ) 00294 { 00295 strong = C12O16Rotate[i].Coll.heat; 00296 thermal.levlmax = 4; 00297 thermal.ipHeatlmax = i; 00298 } 00299 if( C13O16Rotate[i].Coll.heat > strong ) 00300 { 00301 strong = C13O16Rotate[i].Coll.heat; 00302 thermal.levlmax = 5; 00303 thermal.ipHeatlmax = i; 00304 } 00305 } 00306 00307 /* Atoms & Molecules */ 00308 for( i=0; i < linesAdded2; i++ ) 00309 { 00310 if(atmolEmis[i].tran->Coll.heat > strong ) 00311 { 00312 strong = atmolEmis[i].tran->Coll.heat; 00313 thermal.levlmax = 6; 00314 thermal.ipHeatlmax = i; 00315 } 00316 } 00317 return; 00318 }