00001
00002
00003
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 DEBUG_EXIT( "PrtHeader()" );
00047 return;
00048 }
00049
00050 fprintf( ioQQQ, " %4ldCellPeak",rfield.nflux );
00051 PrintE82(ioQQQ, rfield.anu[prt.ipeak-1] );
00052 fprintf( ioQQQ, " Lo");
00053 fprintf( ioQQQ,PrintEfmt("%9.2e", rfield.anu[0] - rfield.widflx[0]/2. ));
00054 fprintf( ioQQQ, "=%6.2fcm Hi-Con:",
00055 9.117e-6/(rfield.anu[0] - rfield.widflx[0]/2.) );
00056 PrintE82(ioQQQ,rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.);
00057 fprintf(ioQQQ," Ryd E(hi):");
00058 PrintE82(ioQQQ,rfield.egamry);
00059 fprintf(ioQQQ,"Ryd E(hi): %9.2f MeV\n", rfield.egamry*0.0000136 );
00060
00061 if( prt.xpow > 0. )
00062 {
00063 prt.xpow = (float)(log10(prt.xpow) + radius.pirsq);
00064 qxl = log10(prt.qx) + radius.pirsq;
00065 }
00066 else
00067 {
00068 prt.xpow = 0.;
00069 qxl = 0.;
00070 }
00071
00072 if( prt.powion > 0. )
00073 {
00074 pionl = log10(prt.powion) + radius.pirsq;
00075 avpow = prt.powion/rfield.qhtot/EN1RYD;
00076 }
00077 else
00078 {
00079 pionl = 0.;
00080 avpow = 0.;
00081 }
00082
00083
00084
00085 if( prt.pbal > 0. )
00086 {
00087 pballog = log10(MAX2(1e-30,prt.pbal)) + radius.pirsq;
00088 qballog = log10(MAX2(1e-30,rfield.qbal)) + radius.pirsq;
00089 }
00090 else
00091 {
00092 pballog = 0.;
00093 qballog = 0.;
00094 }
00095
00096 if( radius.pirsq > 0. )
00097 {
00098 fprintf( ioQQQ, " L(nu>1ryd):%9.4f Average nu:",pionl);
00099 PrintE93(ioQQQ, avpow);
00100 fprintf( ioQQQ," L( X-ray):%9.4f L(BalC):%9.4f Q(Balmer C):%9.4f\n",
00101 prt.xpow, pballog, qballog );
00102 if( pionl > 47. )
00103 {
00104 fprintf(ioQQQ,"\n\n WARNING - the continuum has a luminosity %.2e times greater than the sun.\n",
00105 pow( 10. , pionl-log10(SOLAR_LUMINOSITY) ) );
00106 fprintf(ioQQQ," WARNING - Is this correct? Check the luminosity commands.\n\n\n");
00107 }
00108 }
00109 else
00110 {
00111 fprintf( ioQQQ, " I(nu>1ryd):%9.4f Average nu:",pionl);
00112 PrintE93(ioQQQ, avpow);
00113 fprintf( ioQQQ," I( X-ray):%9.4f I(BalC):%9.4f Phi(BalmrC):%9.4f\n",
00114 prt.xpow,
00115 log10(MAX2(1e-30,prt.pbal)),
00116 log10(MAX2(1e-30,rfield.qbal)) );
00117 }
00118
00119 if( rfield.qhe > 0. )
00120 {
00121 qhel = log10(rfield.qhe) + radius.pirsq;
00122 }
00123 else
00124 {
00125 qhel = 0.;
00126 }
00127
00128 if( rfield.qheii > 0. )
00129 {
00130 qheiil = log10(rfield.qheii) + radius.pirsq;
00131 }
00132 else
00133 {
00134 qheiil = 0.;
00135 }
00136
00137 if( prt.q > 0. )
00138 {
00139 ql = log10(prt.q) + radius.pirsq;
00140 }
00141 else
00142 {
00143 ql = 0.;
00144 }
00145
00146 if( radius.pirsq != 0. )
00147 {
00148 fprintf( ioQQQ,
00149 " Q(1.0-1.8):%9.4f Q(1.8-4.0):%9.4f Q(4.0-20):"
00150 "%9.4f Q(20--):%9.4f Ion pht flx:",
00151 ql,
00152 qhel,
00153 qheiil,
00154 qxl);
00155 PrintE93(ioQQQ, rfield.qhtot );
00156 }
00157 else
00158 {
00159 fprintf( ioQQQ,
00160 " phi(1.0-1.8):%7.4f phi(1.8-4.0):%7.3f phi(4.0-20):"
00161 "%7.3f phi(20--):%7.3f Ion pht flx:",
00162 ql,
00163 qhel,
00164 qheiil,
00165 qxl );
00166 PrintE93(ioQQQ, rfield.qhtot );
00167 }
00168 fprintf( ioQQQ,"\n");
00169
00170
00171 if( rfield.anu[rfield.nflux-1] > 150. )
00172 {
00173
00174 ip2kev = ipoint(147.);
00175 ip2500 = ipoint(0.3645);
00176 if( rfield.flux[ip2500-1] > 1e-28 )
00177 {
00178 ratio = (rfield.flux[ip2kev-1]*rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1])/
00179 (rfield.flux[ip2500-1]*rfield.anu[ip2500-1]/rfield.widflx[ip2500-1]);
00180 }
00181 else
00182 {
00183 ratio = 0.;
00184 }
00185
00186 if( ratio > 0. )
00187 {
00188 alfox = log(ratio)/log(rfield.anu[ip2kev-1]/rfield.anu[ip2500-1]);
00189 }
00190 else
00191 {
00192 alfox = 0.;
00193 }
00194 }
00195 else
00196 {
00197 alfox = 0.;
00198 }
00199
00200 if( prt.GammaLumin > 0. )
00201 {
00202 gpowl = log10(prt.GammaLumin) + radius.pirsq;
00203 qgaml = log10(prt.qgam) + radius.pirsq;
00204 }
00205 else
00206 {
00207 gpowl = 0.;
00208 qgaml = 0.;
00209 }
00210
00211 if( prt.pradio > 0. )
00212 {
00213 radpwl = log10(prt.pradio) + radius.pirsq;
00214 }
00215 else
00216 {
00217 radpwl = 0.;
00218 }
00219
00220 if( radius.pirsq > 0. )
00221 {
00222 fprintf( ioQQQ,
00223 " L(gam ray):%9.4f Q(gam ray):%9.4f L(Infred):%9.4f Alf(ox):%9.4f Total lumin:%9.4f\n",
00224 gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) +
00225 radius.pirsq );
00226 }
00227 else
00228 {
00229 fprintf( ioQQQ,
00230 " I(gam ray):%9.4f phi(gam r):%9.4f I(Infred):%9.4f Alf(ox):%9.4f Total inten:%9.4f\n",
00231 gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) +
00232 radius.pirsq );
00233 }
00234
00235
00236 if( radius.lgRadiusKnown )
00237 {
00238 solar = log10(continuum.TotalLumin) + radius.pirsq - 33.5828;
00239 absbol = 4.75 - 2.5*solar;
00240
00241
00242
00243
00244 absv = -2.5*(log10(MAX2(1e-30,continuum.fluxv)) + radius.pirsq - 20.654202);
00245
00246
00247 if( continuum.fbeta > 0. )
00248 {
00249 continuum.fbeta = (float)(log10(MAX2(1e-37,continuum.fbeta)) + radius.pirsq);
00250 }
00251 else
00252 {
00253 continuum.fbeta = 0.;
00254 }
00255
00256 bolc = absbol - absv;
00257 fprintf( ioQQQ,
00258 " log L/Lsun:%9.4f Abs bol mg:%9.4f Abs V mag:%9.4f Bol cor:%9.4f nuFnu(Bbet):%9.4f\n",
00259 solar,
00260 absbol,
00261 absv,
00262 bolc,
00263 continuum.fbeta );
00264 }
00265
00266 rfield.cmcool = 0.;
00267 rfield.cmheat = 0.;
00268
00269 for( i=0; i < rfield.nflux; i++ )
00270 {
00271
00272 rfield.cmcool += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])*
00273 rfield.csigc[i];
00274
00275
00276
00277 rfield.cmheat += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])*
00278 rfield.csigh[i]*(1. + rfield.OccNumbIncidCont[i]);
00279 }
00280
00281 rfield.cmheat *= dense.eden*1e-15;
00282
00283
00284 rfield.cmcool *= dense.eden*4.*6.338e-6*1e-15;
00285
00286
00287 if( rfield.cmcool > 0. )
00288 {
00289 rfield.lgComUndr = false;
00290 tcomp = rfield.cmheat/rfield.cmcool;
00291 }
00292 else
00293 {
00294
00295 rfield.lgComUndr = true;
00296 tcomp = 0.;
00297 }
00298
00299
00300 if( phycon.TEnerDen > (1.05*tcomp) )
00301 {
00302 thermal.lgEdnGTcm = true;
00303 }
00304 else
00305 {
00306 thermal.lgEdnGTcm = false;
00307 }
00308
00309
00310 fprintf( ioQQQ, " U(1.0----):");
00311 PrintE93( ioQQQ, rfield.uh);
00312 fprintf( ioQQQ, " U(4.0----):");
00313 PrintE93( ioQQQ,rfield.uheii );
00314 fprintf( ioQQQ, " T(En-Den):");
00315 PrintE93( ioQQQ,phycon.TEnerDen );
00316 fprintf( ioQQQ, " T(Comp):");
00317 PrintE93( ioQQQ,tcomp );
00318 fprintf( ioQQQ, " nuJnu(912A):");
00319 PrintE93( ioQQQ,prt.fx1ryd );
00320 fprintf( ioQQQ, "\n");
00321
00322
00323 fprintf( ioQQQ, " Occ(FarIR):");
00324 PrintE93( ioQQQ, rfield.OccNumbIncidCont[0]);
00325 fprintf( ioQQQ, " Occ(H n=6):");
00326
00327 if( iso.numLevels_local[ipH_LIKE][ipHYDROGEN] > 6 )
00328 {
00329 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][6]-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
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
00349 if( iso.numLevels_local[ipH_LIKE][ipHYDROGEN] > 6 )
00350 {
00351 PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][6]-1]*TE1RYD*rfield.anu[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][6]-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
00373 fprintf( ioQQQ, " \n" );
00374
00375 DEBUG_EXIT( "PrtHeader()" );
00376 return;
00377 }
00378