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 /*PrtHeader print program's header, including luminosities and ionization parameters */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "phycon.h" 00007 #include "iso.h" 00008 #include "rfield.h" 00009 #include "radius.h" 00010 #include "called.h" 00011 #include "thermal.h" 00012 #include "dense.h" 00013 #include "continuum.h" 00014 #include "ipoint.h" 00015 #include "prt.h" 00016 00017 void PrtHeader(void) 00018 { 00019 long int i, 00020 ip2500, 00021 ip2kev; 00022 double absbol, 00023 absv, 00024 alfox, 00025 avpow, 00026 bolc, 00027 gpowl, 00028 pballog, 00029 pionl, 00030 qballog, 00031 qgaml, 00032 qheiil, 00033 qhel, 00034 ql, 00035 qxl, 00036 radpwl, 00037 ratio, 00038 solar, 00039 tcomp, 00040 tradio; 00041 00042 DEBUG_ENTRY( "PrtHeader()" ); 00043 00044 if( !called.lgTalk ) 00045 { 00046 return; 00047 } 00048 00049 fprintf( ioQQQ, " %4ldCellPeak",rfield.nflux ); 00050 PrintE82(ioQQQ, rfield.anu[prt.ipeak-1] ); 00051 fprintf( ioQQQ, " Lo"); 00052 fprintf( ioQQQ,PrintEfmt("%9.2e", rfield.anu[0] - rfield.widflx[0]/2. )); 00053 fprintf( ioQQQ, "=%6.2fcm Hi-Con:", 00054 9.117e-6/(rfield.anu[0] - rfield.widflx[0]/2.) ); 00055 PrintE82(ioQQQ,rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.); 00056 fprintf(ioQQQ," Ryd E(hi):"); 00057 PrintE82(ioQQQ,rfield.egamry); 00058 fprintf(ioQQQ,"Ryd E(hi): %9.2f MeV\n", rfield.egamry*0.0000136 ); 00059 00060 if( prt.xpow > 0. ) 00061 { 00062 prt.xpow = (realnum)(log10(prt.xpow) + radius.pirsq); 00063 qxl = log10(prt.qx) + radius.pirsq; 00064 } 00065 else 00066 { 00067 prt.xpow = 0.; 00068 qxl = 0.; 00069 } 00070 00071 if( prt.powion > 0. ) 00072 { 00073 pionl = log10(prt.powion) + radius.pirsq; 00074 avpow = prt.powion/rfield.qhtot/EN1RYD; 00075 } 00076 else 00077 { 00078 pionl = 0.; 00079 avpow = 0.; 00080 } 00081 00082 /* >>chng 97 mar 18, break these two out here, so that returns zero 00083 * when no radiation - had been done in the print statement so pirsq was printed */ 00084 if( prt.pbal > 0. ) 00085 { 00086 pballog = log10(MAX2(1e-30,prt.pbal)) + radius.pirsq; 00087 qballog = log10(MAX2(1e-30,rfield.qbal)) + radius.pirsq; 00088 } 00089 else 00090 { 00091 pballog = 0.; 00092 qballog = 0.; 00093 } 00094 00095 if( radius.pirsq > 0. ) 00096 { 00097 fprintf( ioQQQ, " L(nu>1ryd):%9.4f Average nu:",pionl); 00098 PrintE93(ioQQQ, avpow); 00099 fprintf( ioQQQ," L( X-ray):%9.4f L(BalC):%9.4f Q(Balmer C):%9.4f\n", 00100 prt.xpow, pballog, qballog ); 00101 if( pionl > 47. ) 00102 { 00103 fprintf(ioQQQ,"\n\n WARNING - the continuum has a luminosity %.2e times greater than the sun.\n", 00104 pow( 10. , pionl-log10(SOLAR_LUMINOSITY) ) ); 00105 fprintf(ioQQQ," WARNING - Is this correct? Check the luminosity commands.\n\n\n"); 00106 } 00107 } 00108 else 00109 { 00110 fprintf( ioQQQ, " I(nu>1ryd):%9.4f Average nu:",pionl); 00111 PrintE93(ioQQQ, avpow); 00112 fprintf( ioQQQ," I( X-ray):%9.4f I(BalC):%9.4f Phi(BalmrC):%9.4f\n", 00113 prt.xpow, 00114 log10(MAX2(1e-30,prt.pbal)), 00115 log10(MAX2(1e-30,rfield.qbal)) ); 00116 } 00117 00118 if( rfield.qhe > 0. ) 00119 { 00120 qhel = log10(rfield.qhe) + radius.pirsq; 00121 } 00122 else 00123 { 00124 qhel = 0.; 00125 } 00126 00127 if( rfield.qheii > 0. ) 00128 { 00129 qheiil = log10(rfield.qheii) + radius.pirsq; 00130 } 00131 else 00132 { 00133 qheiil = 0.; 00134 } 00135 00136 if( prt.q > 0. ) 00137 { 00138 ql = log10(prt.q) + radius.pirsq; 00139 } 00140 else 00141 { 00142 ql = 0.; 00143 } 00144 00145 if( radius.pirsq != 0. ) 00146 { 00147 fprintf( ioQQQ, 00148 " Q(1.0-1.8):%9.4f Q(1.8-4.0):%9.4f Q(4.0-20):" 00149 "%9.4f Q(20--):%9.4f Ion pht flx:", 00150 ql, 00151 qhel, 00152 qheiil, 00153 qxl); 00154 PrintE93(ioQQQ, rfield.qhtot ); 00155 } 00156 else 00157 { 00158 fprintf( ioQQQ, 00159 " phi(1.0-1.8):%7.4f phi(1.8-4.0):%7.3f phi(4.0-20):" 00160 "%7.3f phi(20--):%7.3f Ion pht flx:", 00161 ql, 00162 qhel, 00163 qheiil, 00164 qxl ); 00165 PrintE93(ioQQQ, rfield.qhtot ); 00166 } 00167 fprintf( ioQQQ,"\n"); 00168 00169 /* estimate alpha (o-x) */ 00170 if( rfield.anu[rfield.nflux-1] > 150. ) 00171 { 00172 /* the ratio of fluxes is nearly 403.3^alpha ox */ 00173 ip2kev = ipoint(147.); 00174 ip2500 = ipoint(0.3645); 00175 if( rfield.flux[ip2500-1] > 1e-28 ) 00176 { 00177 ratio = (rfield.flux[ip2kev-1]*rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1])/ 00178 (rfield.flux[ip2500-1]*rfield.anu[ip2500-1]/rfield.widflx[ip2500-1]); 00179 } 00180 else 00181 { 00182 ratio = 0.; 00183 } 00184 00185 if( ratio > 0. ) 00186 { 00187 alfox = log(ratio)/log(rfield.anu[ip2kev-1]/rfield.anu[ip2500-1]); 00188 } 00189 else 00190 { 00191 alfox = 0.; 00192 } 00193 } 00194 else 00195 { 00196 alfox = 0.; 00197 } 00198 00199 if( prt.GammaLumin > 0. ) 00200 { 00201 gpowl = log10(prt.GammaLumin) + radius.pirsq; 00202 qgaml = log10(prt.qgam) + radius.pirsq; 00203 } 00204 else 00205 { 00206 gpowl = 0.; 00207 qgaml = 0.; 00208 } 00209 00210 if( prt.pradio > 0. ) 00211 { 00212 radpwl = log10(prt.pradio) + radius.pirsq; 00213 } 00214 else 00215 { 00216 radpwl = 0.; 00217 } 00218 00219 if( radius.pirsq > 0. ) 00220 { 00221 fprintf( ioQQQ, 00222 " L(gam ray):%9.4f Q(gam ray):%9.4f L(Infred):%9.4f Alf(ox):%9.4f Total lumin:%9.4f\n", 00223 gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) + 00224 radius.pirsq ); 00225 } 00226 else 00227 { 00228 fprintf( ioQQQ, 00229 " I(gam ray):%9.4f phi(gam r):%9.4f I(Infred):%9.4f Alf(ox):%9.4f Total inten:%9.4f\n", 00230 gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) + 00231 radius.pirsq ); 00232 } 00233 00234 /* magnitudes */ 00235 if( radius.lgRadiusKnown ) 00236 { 00237 solar = log10(continuum.TotalLumin) + radius.pirsq - 33.5828; 00238 absbol = 4.75 - 2.5*solar; 00239 00240 /* absv = 4.79 - 2.5 * (LOG10(MAX(1e-30,continuum.fluxv)) + pirsq - 18.742 - 00241 * 1 0.016) 00242 * allen 76, page 197 */ 00243 absv = -2.5*(log10(MAX2(1e-30,continuum.fluxv)) + radius.pirsq - 20.654202); 00244 00245 /* >>chng 97 mar 18, following branch so zero returned when no radiation at all */ 00246 if( continuum.fbeta > 0. ) 00247 { 00248 continuum.fbeta = (realnum)(log10(MAX2(1e-37,continuum.fbeta)) + radius.pirsq); 00249 } 00250 else 00251 { 00252 continuum.fbeta = 0.; 00253 } 00254 00255 bolc = absbol - absv; 00256 fprintf( ioQQQ, 00257 " log L/Lsun:%9.4f Abs bol mg:%9.4f Abs V mag:%9.4f Bol cor:%9.4f nuFnu(Bbet):%9.4f\n", 00258 solar, 00259 absbol, 00260 absv, 00261 bolc, 00262 continuum.fbeta ); 00263 } 00264 00265 rfield.cmcool = 0.; 00266 rfield.cmheat = 0.; 00267 00268 for( i=0; i < rfield.nflux; i++ ) 00269 { 00270 /* CSIGC is Tarter expression times ANU(I)*3.858E-25 */ 00271 rfield.cmcool += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])* 00272 rfield.csigc[i]; 00273 00274 /* Compton heating with correction for induced Compton scattering 00275 * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25 */ 00276 rfield.cmheat += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])* 00277 rfield.csigh[i]*(1. + rfield.OccNumbIncidCont[i]); 00278 } 00279 00280 rfield.cmheat *= dense.eden*1e-15; 00281 00282 /* 6.338E-6 is k in inf mass Rydbergs */ 00283 rfield.cmcool *= dense.eden*4.*6.338e-6*1e-15; 00284 00285 /* all of following need factor of 10^-15 to be true rates */ 00286 if( rfield.cmcool > 0. ) 00287 { 00288 rfield.lgComUndr = false; 00289 tcomp = rfield.cmheat/rfield.cmcool; 00290 } 00291 else 00292 { 00293 /* underflow if Compt cooling rate is zero */ 00294 rfield.lgComUndr = true; 00295 tcomp = 0.; 00296 } 00297 00298 /* check whether energy density temp is greater than compton temp */ 00299 if( phycon.TEnerDen > (1.05*tcomp) ) 00300 { 00301 thermal.lgEdnGTcm = true; 00302 } 00303 else 00304 { 00305 thermal.lgEdnGTcm = false; 00306 } 00307 00308 /* say some ionization parameters */ 00309 fprintf( ioQQQ, " U(1.0----):"); 00310 PrintE93( ioQQQ, rfield.uh); 00311 fprintf( ioQQQ, " U(4.0----):"); 00312 PrintE93( ioQQQ,rfield.uheii ); 00313 fprintf( ioQQQ, " T(En-Den):"); 00314 PrintE93( ioQQQ,phycon.TEnerDen ); 00315 fprintf( ioQQQ, " T(Comp):"); 00316 PrintE93( ioQQQ,tcomp ); 00317 fprintf( ioQQQ, " nuJnu(912A):"); 00318 PrintE93( ioQQQ,prt.fx1ryd ); 00319 fprintf( ioQQQ, "\n"); 00320 00321 /* some occupation numbers */ 00322 fprintf( ioQQQ, " Occ(FarIR):"); 00323 PrintE93( ioQQQ, rfield.OccNumbIncidCont[0]); 00324 fprintf( ioQQQ, " Occ(H n=6):"); 00325 00326 if( iso.n_HighestResolved_local[ipH_LIKE][ipHYDROGEN] + iso.nCollapsed_local[ipH_LIKE][ipHYDROGEN] >= 6 ) 00327 { 00328 long ipH6p = iso.QuantumNumbers2Index[ipH_LIKE][ipHYDROGEN][6][1][2]; 00329 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH6p]-1]); 00330 } 00331 else 00332 { 00333 PrintE93( ioQQQ, 0.); 00334 } 00335 fprintf( ioQQQ, " Occ(1Ryd):"); 00336 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]); 00337 fprintf( ioQQQ, " Occ(4R):"); 00338 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1]); 00339 fprintf( ioQQQ, " Occ (Nu-hi):"); 00340 PrintE93( ioQQQ, rfield.OccNumbIncidCont[rfield.nflux-1] ); 00341 fprintf( ioQQQ, "\n"); 00342 00343 /* now print brightness temps */ 00344 tradio = rfield.OccNumbIncidCont[0]*TE1RYD*rfield.anu[0]; 00345 fprintf( ioQQQ, " Tbr(FarIR):"); 00346 PrintE93( ioQQQ, tradio); 00347 fprintf( ioQQQ, " Tbr(H n=6):"); 00348 if( iso.n_HighestResolved_local[ipH_LIKE][ipHYDROGEN] + iso.nCollapsed_local[ipH_LIKE][ipHYDROGEN] >= 6 ) 00349 { 00350 long ipH6p = iso.QuantumNumbers2Index[ipH_LIKE][ipHYDROGEN][6][1][2]; 00351 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH6p]-1]*TE1RYD*rfield.anu[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH6p]-1]); 00352 } 00353 else 00354 { 00355 PrintE93( ioQQQ, 0.); 00356 } 00357 fprintf( ioQQQ, " Tbr(1Ryd):"); 00358 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]*TE1RYD*rfield.anu[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]); 00359 fprintf( ioQQQ, " Tbr(4R):"); 00360 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1]*TE1RYD*rfield.anu[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1]); 00361 fprintf( ioQQQ, " Tbr (Nu-hi):"); 00362 PrintE93( ioQQQ, rfield.OccNumbIncidCont[rfield.nflux-1]*TE1RYD*rfield.anu[rfield.nflux-1]); 00363 fprintf( ioQQQ, "\n"); 00364 00365 if( tradio > 1e9 ) 00366 { 00367 fprintf( ioQQQ, 00368 " >>>The radio brightness temperature is very large,%10.2eK at%10.2ecm. Is this physical???\n", 00369 tradio, 9.115e-6/rfield.anu[0] ); 00370 } 00371 00372 /* skip a line */ 00373 fprintf( ioQQQ, " \n" ); 00374 return; 00375 }