cloudy  trunk
prt_lines_general.cpp
Go to the documentation of this file.
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 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated for cloudy by doxygen 1.7.6.1