00001
00002
00003
00004
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
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
00051 nelem = ipHYDROGEN;
00052 ipLo = ipH2p;
00053
00054
00055
00056 ASSERT( iso.n_HighestResolved_max[ipH_LIKE][nelem] >=3 );
00057
00058 if( iso.n_HighestResolved_max[ipH_LIKE][nelem] >= 4 )
00059 {
00060 ipHi = ipH4p;
00061 hbetac =
00062 (Transitions[ipH_LIKE][nelem][ipH4p][ipH2s].Emis->Aul *
00063 Transitions[ipH_LIKE][nelem][ipH4p][ipH2s].Emis->Pesc *
00064 StatesElemNEW[nelem][nelem-ipH_LIKE][ipH4p].Pop +
00065 Transitions[ipH_LIKE][nelem][ipH4s][ipH2p].Emis->Aul *
00066 Transitions[ipH_LIKE][nelem][ipH4s][ipH2p].Emis->Pesc *
00067 StatesElemNEW[nelem][nelem-ipH_LIKE][ipH4s].Pop +
00068 Transitions[ipH_LIKE][nelem][ipH4d][ipH2p].Emis->Aul *
00069 Transitions[ipH_LIKE][nelem][ipH4d][ipH2p].Emis->Pesc *
00070 StatesElemNEW[nelem][nelem-ipH_LIKE][ipH4d].Pop) *
00071 Transitions[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg;
00072 }
00073 else
00074 {
00075
00076 ASSERT( iso.n_HighestResolved_max[ipH_LIKE][nelem] == 3 );
00077 ipHi = 6;
00078 hbetac =
00079 (Transitions[ipH_LIKE][nelem][ipHi][ipH2s].Emis->Aul *
00080 Transitions[ipH_LIKE][nelem][ipHi][ipH2s].Emis->Pesc +
00081 Transitions[ipH_LIKE][nelem][ipHi][ipH2p].Emis->Aul *
00082 Transitions[ipH_LIKE][nelem][ipHi][ipH2p].Emis->Pesc ) *
00083 StatesElemNEW[nelem][nelem-ipH_LIKE][ipHi].Pop *
00084 Transitions[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg;
00085 }
00086
00087
00088
00089
00090
00091 rt.fracin = Transitions[ipH_LIKE][nelem][ipHi][ipH2s].Emis->FracInwd;
00092 lindst(hbetac,Transitions[ipH_LIKE][nelem][ipHi][ipH2s].WLAng,"TOTL",
00093 Transitions[ipH_LIKE][nelem][ipHi][ipH2s].ipCont,'i',false,
00094 " H I Balmer beta predicted by model atom " );
00095 rt.fracin = 0.5;
00096
00097 if( iso.n_HighestResolved_max[ipH_LIKE][nelem] < 4 )
00098 {
00099
00100 lindst(hbetac,Transitions[ipH_LIKE][nelem][ipHi][ipH2s].WLAng,"H 1",
00101 Transitions[ipH_LIKE][nelem][ipHi][ipH2s].ipCont,'i',false,
00102 " H I Balmer beta predicted by model atom " );
00103
00104 lindst(hbetac/2.,Transitions[ipH_LIKE][nelem][ipHi][ipH2s].WLAng,"Inwd",
00105 Transitions[ipH_LIKE][nelem][ipHi][ipH2s].ipCont,'i',false,
00106 " H I Balmer beta predicted by model atom " );
00107 }
00108
00109
00110 ipHi = ipH2p;
00111 ipLo = ipH1s;
00112 hlalph =
00113 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul*
00114 StatesElemNEW[nelem][nelem-ipH_LIKE][ipHi].Pop*
00115 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
00116 Transitions[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg;
00117
00118 rt.fracin = Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->FracInwd;
00119 lindst(hlalph,Transitions[ipH_LIKE][nelem][ipHi][ipLo].WLAng,"TOTL",
00120 Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont,'i',false ,
00121 " H I Lya predicted from model atom ");
00122 rt.fracin = 0.5;
00123
00124
00125 if( geometry.iEmissPower == 2 )
00126 {
00127 linadd(continuum.totlsv/radius.dVeffAper,0,"Inci",'i',
00128 "total luminosity in incident continuum");
00129
00130
00131 if( LineSave.ipass > 0 )
00132 {
00133 continuum.totlsv = 0.;
00134 }
00135 }
00136
00137 linadd(thermal.htot,0,"TotH",'i',
00138 " total heating, all forms, information since individuals added later ");
00139
00140 linadd(thermal.ctot,0,"TotC",'i',
00141 " total cooling, all forms, information since individuals added later ");
00142
00143 linadd(thermal.heating[0][0],0,"BFH1",'h',
00144 " hydrogen photoionization heating, ground state only ");
00145
00146 linadd(thermal.heating[0][1],0,"BFHx",'h',
00147 " net hydrogen photoionization heating less rec cooling, all excited states normally zero, positive if excited states are net heating ");
00148
00149 linadd(thermal.heating[0][22],0,"Line",'h',
00150 " heating due to induced lines absorption of continuum ");
00151 if( thermal.htot > 0. )
00152 {
00153 if( thermal.heating[0][22]/thermal.htot > thermal.HeatLineMax )
00154 {
00155 thermal.HeatLineMax = (realnum)(thermal.heating[0][22]/thermal.htot);
00156
00157 GetMaxhLine();
00158 }
00159 }
00160
00161 linadd(thermal.heating[1][0]+thermal.heating[1][1]+thermal.heating[1][2],0,"BFHe",'h',
00162 " total helium photoionization heating, all stages ");
00163
00164 HeatMetal = 0.;
00165
00166 for( nelem=2; nelem<LIMELM; ++nelem)
00167 {
00168
00169 for( i=dense.IonLow[nelem]; i < dense.IonHigh[nelem]; i++ )
00170 {
00171 ASSERT( i < LIMELM );
00172
00173 HeatMetal += thermal.heating[nelem][i];
00174 }
00175 }
00176
00177 linadd(HeatMetal,0,"TotM",'h',
00178 " total heavy element photoionization heating, all stages ");
00179
00180 linadd(thermal.heating[0][21],0,"pair",'h',
00181 " heating due to pair production ");
00182
00183
00184
00185 if( LineSave.ipass > 0 )
00186 {
00187
00188 ionbal.CompHeating_Max = MAX2( ionbal.CompHeating_Max , ionbal.CompRecoilHeatLocal/thermal.htot);
00189 }
00190 else
00191 {
00192 ionbal.CompHeating_Max = 0.;
00193 }
00194
00195 linadd(ionbal.CompRecoilHeatLocal,0,"Cbnd",'h',
00196 " heating due to bound compton scattering ");
00197
00198 linadd(rfield.cmheat,0,"ComH",'h',
00199 " Compton heating ");
00200
00201 linadd(CoolHeavy.tccool,0,"ComC",'c',
00202 " total Compton cooling ");
00203
00204
00205 dynamics.HeatMax = MAX2( dynamics.HeatMax , dynamics.Heat() /thermal.htot );
00206
00207 dynamics.CoolMax = MAX2( dynamics.CoolMax , dynamics.Cool() /thermal.htot );
00208
00209 linadd(dynamics.Cool() , 0 , "advC" , 'i',
00210 " cooling due to advection " );
00211
00212 linadd(dynamics.Heat() , 0 , "advH" , 'i' ,
00213 " heating due to advection ");
00214
00215 linadd( thermal.char_tran_heat ,0,"CT H",'h',
00216 " heating due to charge transfer ");
00217
00218 linadd( thermal.char_tran_cool ,0,"CT C",'c',
00219 " cooling due to charge transfer ");
00220
00221 linadd(thermal.heating[1][6],0,"CR H",'h',
00222 " cosmic ray heating ");
00223
00224 linadd(thermal.heating[0][20],0,"extH",'h',
00225 " extra heat added to this zone, from HEXTRA command ");
00226
00227 linadd(CoolHeavy.cextxx,0,"extC",'c',
00228 " extra cooling added to this zone, from CEXTRA command ");
00229
00230
00231
00232 ee511 = (dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM])*ionbal.PairProducPhotoRate[0]*2.*8.20e-7;
00233 PntForLine(2.427e-2,"e-e+",&ipnt);
00234 lindst(ee511,(realnum)2.427e-2,"e-e+",ipnt,'r',true,
00235 " 511keV annihilation line " );
00236
00237 linadd(CoolHeavy.expans,0,"Expn",'c',
00238 " expansion cooling, only non-zero for wind ");
00239
00240 linadd(iso.RadRecCool[ipH_LIKE][ipHYDROGEN],0,"H FB",'i',
00241 " H radiative recombination cooling ");
00242
00243 linadd(MAX2(0.,iso.FreeBnd_net_Cool_Rate[ipH_LIKE][ipHYDROGEN]),0,"HFBc",'c',
00244 " net free-bound cooling ");
00245
00246 linadd(MAX2(0.,-iso.FreeBnd_net_Cool_Rate[ipH_LIKE][ipHYDROGEN]),0,"HFBh",'h',
00247 " net free-bound heating ");
00248
00249 linadd(iso.RecomInducCool_Rate[ipH_LIKE][ipHYDROGEN],0,"Hind",'c',
00250 " cooling due to induced rec of hydrogen ");
00251
00252 linadd(CoolHeavy.cyntrn,0,"Cycn",'c',
00253 " cyclotron cooling ");
00254
00255
00256 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
00257 {
00258
00259 char chLabel[5];
00260 strncpy( chLabel , Species[ipSpecies].chLabel , 4 );
00261 chLabel[4] = '\0';
00262
00263
00264 linadd(Species[ipSpecies].CoolTotal,0, chLabel,'i',
00265 " net cooling due to database species");
00266 }
00267
00268 return;
00269 }
00270
00271
00272 STATIC void GetMaxhLine(void)
00273 {
00274 long int i;
00275 double strong;
00276
00277 DEBUG_ENTRY( "GetMaxhLine()" );
00278
00279
00280 strong = 0.;
00281
00282
00283 thermal.levlmax = 0;
00284
00285 for( i=1; i <= nLevel1; i++ )
00286 {
00287 if( TauLines[i].Coll.heat > strong )
00288 {
00289 strong = TauLines[i].Coll.heat;
00290 thermal.levlmax = 1;
00291 thermal.ipHeatlmax = i;
00292 }
00293 }
00294
00295 for( i=0; i < nWindLine; i++ )
00296 {
00297 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
00298 {
00299 if( TauLine2[i].Coll.heat > strong )
00300 {
00301 strong = TauLine2[i].Coll.heat;
00302 thermal.levlmax = 2;
00303 thermal.ipHeatlmax = i;
00304 }
00305 }
00306 }
00307
00308 for( i=0; i < nHFLines; i++ )
00309 {
00310 if( HFLines[i].Coll.heat > strong )
00311 {
00312 strong = HFLines[i].Coll.heat;
00313 thermal.levlmax = 3;
00314 thermal.ipHeatlmax = i;
00315 }
00316 }
00317
00318
00319 for( i=0; i < linesAdded2; i++ )
00320 {
00321 if(dBaseLines[i].tran->Coll.heat > strong )
00322 {
00323 strong = dBaseLines[i].tran->Coll.heat;
00324 thermal.levlmax = 4;
00325 thermal.ipHeatlmax = i;
00326 }
00327 }
00328
00329 fixit();
00330
00331
00332 return;
00333 }