cloudy: /home66/gary/public_html/cloudy/c08_branch/source/prt

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 }
/home66/gary/public_html/cloudy/c08_branch/source/prt_header.cpp
