prt_header.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*PrtHeader print program's header, including luminosities and ionization parameters */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "iso.h"
#include "rfield.h"
#include "radius.h"
#include "version.h"
#include "called.h"
#include "thermal.h"
#include "dense.h"
#include "continuum.h"
#include "ipoint.h"
#include "prt.h"
 
void PrtHeader(void)
{
        long int i, 
          ip2500, 
          ip2kev;
        double absbol, 
          absv, 
          alfox, 
          avpow, 
          bolc, 
          gpowl, 
          pballog, 
          pionl, 
          qballog, 
          qgaml, 
          qheiil, 
          qhel, 
          ql, 
          qxl, 
          radpwl, 
          ratio, 
          solar, 
          tcomp, 
          tradio;
 
        DEBUG_ENTRY( "PrtHeader()" );
 
        if( !called.lgTalk )
        { 
                return;
        }
 
        if( prt.lgPrtCitations )
        {
                /* the print citations command */
                CloudyPrintReference();
                fprintf( ioQQQ, "\n" );
                DatabasePrintReference();
                fprintf( ioQQQ, "\n\n" );
        }
 
        fprintf( ioQQQ, "           %4ldCellPeak",rfield.nflux );
        PrintE82(ioQQQ, rfield.anu[prt.ipeak-1] );
        fprintf( ioQQQ, "   Lo");
        fprintf( ioQQQ,PrintEfmt("%9.2e", rfield.anu[0] - rfield.widflx[0]/2. ));
        fprintf( ioQQQ, "=%6.2fcm   Hi-Con:",
          9.117e-6/(rfield.anu[0] - rfield.widflx[0]/2.) );
        PrintE82(ioQQQ,rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.);
        fprintf(ioQQQ," Ryd   E(hi):");
        PrintE82(ioQQQ,rfield.egamry);
        fprintf(ioQQQ,"Ryd     E(hi):  %9.2f MeV\n", rfield.egamry*0.0000136 );
 
        if( prt.xpow > 0. )
        {
                prt.xpow = (realnum)(log10(prt.xpow) + radius.pirsq);
                qxl = log10(prt.qx) + radius.pirsq;
        }
        else
        {
                prt.xpow = 0.;
                qxl = 0.;
        }
 
        if( prt.powion > 0. )
        {
                pionl = log10(prt.powion) + radius.pirsq;
                avpow = prt.powion/rfield.qhtot/EN1RYD;
        }
        else
        {
                pionl = 0.;
                avpow = 0.;
        }
 
        /* >>chng 97 mar 18, break these two out here, so that returns zero
         * when no radiation - had been done in the print statement so pirsq was printed */
        if( prt.pbal > 0. )
        {
                pballog = log10(MAX2(1e-30,prt.pbal)) + radius.pirsq;
                qballog = log10(MAX2(1e-30,rfield.qbal)) + radius.pirsq;
        }
        else
        {
                pballog = 0.;
                qballog = 0.;
        }
 
        if( radius.pirsq > 0. )
        {
                fprintf( ioQQQ, "           L(nu>1ryd):%9.4f   Average nu:",pionl);
                PrintE93(ioQQQ, avpow);
                fprintf( ioQQQ,"   L( X-ray):%9.4f   L(BalC):%9.4f     Q(Balmer C):%9.4f\n", 
                  prt.xpow, pballog, qballog );
                if( pionl > 47. )
                {
                        fprintf(ioQQQ,"\n\n WARNING - the continuum has a luminosity %.2e times greater than the sun.\n",
                                pow( 10. , pionl-log10(SOLAR_LUMINOSITY) ) );
                        fprintf(ioQQQ," WARNING - Is this correct?  Check the luminosity commands.\n\n\n");
                }
        }
        else
        {
                fprintf( ioQQQ, "           I(nu>1ryd):%9.4f   Average nu:",pionl);
                PrintE93(ioQQQ, avpow);
                fprintf( ioQQQ,"   I( X-ray):%9.4f   I(BalC):%9.4f     Phi(BalmrC):%9.4f\n", 
                  prt.xpow, 
                  log10(MAX2(1e-30,prt.pbal)), 
                  log10(MAX2(1e-30,rfield.qbal)) );
        }
 
        if( rfield.qhe > 0. )
        {
                qhel = log10(rfield.qhe) + radius.pirsq;
        }
        else
        {
                qhel = 0.;
        }
 
        if( rfield.qheii > 0. )
        {
                qheiil = log10(rfield.qheii) + radius.pirsq;
        }
        else
        {
                qheiil = 0.;
        }
 
        if( prt.q > 0. )
        {
                ql = log10(prt.q) + radius.pirsq;
        }
        else
        {
                ql = 0.;
        }
 
        if( radius.pirsq != 0. )
        {
                fprintf( ioQQQ, 
                        "           Q(1.0-1.8):%9.4f   Q(1.8-4.0):%9.4f   Q(4.0-20):"
                        "%9.4f   Q(20--):%9.4f     Ion pht flx:", 
                  ql, 
                  qhel, 
                  qheiil, 
                  qxl);
                PrintE93(ioQQQ, rfield.qhtot );
        }
        else
        {
                fprintf( ioQQQ, 
                        "           phi(1.0-1.8):%7.4f   phi(1.8-4.0):%7.3f   phi(4.0-20):"
                        "%7.3f   phi(20--):%7.3f     Ion pht flx:", 
                  ql, 
                  qhel, 
                  qheiil, 
                  qxl );
                PrintE93(ioQQQ, rfield.qhtot );
        }
        fprintf( ioQQQ,"\n");
 
        /* estimate alpha (o-x) */
        if( rfield.anu[rfield.nflux-1] > 150. )
        {
                /* the ratio of fluxes is nearly 403.3^alpha ox */
                ip2kev = ipoint(147.);
                ip2500 = ipoint(0.3645);
                if( rfield.flux[0][ip2500-1] > 1e-28 )
                {
                        ratio = (rfield.flux[0][ip2kev-1]*rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1])/
                          (rfield.flux[0][ip2500-1]*rfield.anu[ip2500-1]/rfield.widflx[ip2500-1]);
                }
                else
                {
                        ratio = 0.;
                }
 
                if( ratio > 0. )
                {
                        alfox = log(ratio)/log(rfield.anu[ip2kev-1]/rfield.anu[ip2500-1]);
                }
                else
                {
                        alfox = 0.;
                }
        }
        else
        {
                alfox = 0.;
        }
 
        if( prt.GammaLumin > 0. )
        {
                gpowl = log10(prt.GammaLumin) + radius.pirsq;
                qgaml = log10(prt.qgam) + radius.pirsq;
        }
        else
        {
                gpowl = 0.;
                qgaml = 0.;
        }
 
        if( prt.pradio > 0. )
        {
                radpwl = log10(prt.pradio) + radius.pirsq;
        }
        else
        {
                radpwl = 0.;
        }
 
        if( radius.pirsq > 0. )
        {
                fprintf( ioQQQ, 
                        "           L(gam ray):%9.4f   Q(gam ray):%9.4f   L(Infred):%9.4f   Alf(ox):%9.4f     Total lumin:%9.4f\n", 
                  gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) + 
                  radius.pirsq );
        }
        else
        {
                fprintf( ioQQQ, 
                        "           I(gam ray):%9.4f   phi(gam r):%9.4f   I(Infred):%9.4f   Alf(ox):%9.4f     Total inten:%9.4f\n", 
                  gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) + 
                  radius.pirsq );
        }
 
        /* magnitudes */
        if( radius.lgPredLumin )
        {
                solar = log10(continuum.TotalLumin) + radius.pirsq - 33.5828;
                absbol = 4.75 - 2.5*solar;
 
                /* absv = 4.79 - 2.5 * (LOG10(MAX(1e-30,continuum.fluxv)) + pirsq - 18.742 -
                 *  1  0.016)
                 * allen 76, page 197 */
                absv = -2.5*(log10(MAX2(1e-30,continuum.fluxv)) + radius.pirsq - 20.654202);
 
                /* >>chng 97 mar 18, following branch so zero returned when no radiation at all */
                if( continuum.fbeta > 0. )
                {
                        continuum.fbeta = (realnum)(log10(MAX2(1e-37,continuum.fbeta)) + radius.pirsq);
                }
                else
                {
                        continuum.fbeta = 0.;
                }
 
                bolc = absbol - absv;
                fprintf( ioQQQ, 
                        "           log L/Lsun:%9.4f   Abs bol mg:%9.4f   Abs V mag:%9.4f   Bol cor:%9.4f     nuFnu(Bbet):%9.4f\n", 
                  solar, 
                  absbol, 
                  absv, 
                  bolc, 
                  continuum.fbeta );
        }
 
        rfield.cmcool = 0.;
        rfield.cmheat = 0.;
 
        for( i=0; i < rfield.nflux; i++ )
        {
                /* CSIGC is Tarter expression times ANU(I)*3.858E-25 */
                rfield.cmcool += (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])*
                  rfield.csigc[i];
 
                /* Compton heating with correction for induced Compton scattering
                 * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25 */
                rfield.cmheat += (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i])*
                  rfield.csigh[i]*(1. + rfield.OccNumbIncidCont[i]);
        }
 
        rfield.cmheat *= dense.eden*1e-15;
 
        /* 6.338E-6 is k in inf mass Rydbergs */
        rfield.cmcool *= dense.eden*4.*6.338e-6*1e-15;
 
        /* all of following need factor of 10^-15 to be true rates */
        if( rfield.cmcool > 0. )
        {
                rfield.lgComUndr = false;
                tcomp = rfield.cmheat/rfield.cmcool;
        }
        else
        {
                /* underflow if Compt cooling rate is zero */
                rfield.lgComUndr = true;
                tcomp = 0.;
        }
 
        /* check whether energy density temp is greater than compton temp */
        if( phycon.TEnerDen > (1.05*tcomp) )
        {
                thermal.lgEdnGTcm = true;
        }
        else
        {
                thermal.lgEdnGTcm = false;
        }
 
        /* say some ionization parameters */
        fprintf( ioQQQ, "           U(1.0----):");
        PrintE93( ioQQQ, rfield.uh);
        fprintf( ioQQQ, "   U(4.0----):");
        PrintE93( ioQQQ,rfield.uheii );
        fprintf( ioQQQ, "   T(En-Den):");
        PrintE93( ioQQQ,phycon.TEnerDen );
        fprintf( ioQQQ, "   T(Comp):");
        PrintE93( ioQQQ,tcomp );
        fprintf( ioQQQ, "     nuJnu(912A):");
        PrintE93( ioQQQ,prt.fx1ryd );
        fprintf( ioQQQ, "\n");
 
        /* some occupation numbers */
        fprintf( ioQQQ, "           Occ(FarIR):");
        PrintE93( ioQQQ, rfield.OccNumbIncidCont[0]);
        fprintf( ioQQQ, "   Occ(H n=6):");
 
        if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_local + iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local >= 6 )
        {
                long ip6p = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[6][1][2];
                PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ip6p].ipIsoLevNIonCon-1]);
        }
        else
        {
                PrintE93( ioQQQ, 0.);
        }
        fprintf( ioQQQ, "   Occ(1Ryd):");
        PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1]);
        fprintf( ioQQQ, "   Occ(4R):");
        PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1]);
        fprintf( ioQQQ, "     Occ (Nu-hi):");
        PrintE93( ioQQQ, rfield.OccNumbIncidCont[rfield.nflux-1] );
        fprintf( ioQQQ, "\n"); 
 
        /* now print brightness temps */
        tradio = rfield.OccNumbIncidCont[0]*TE1RYD*rfield.anu[0];
        fprintf( ioQQQ, "           Tbr(FarIR):");
        PrintE93( ioQQQ, tradio);
        fprintf( ioQQQ, "   Tbr(H n=6):");
        if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_local + iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local >= 6 )
        {
                long ip6p = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[6][1][2];
                PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ip6p].ipIsoLevNIonCon-1]*TE1RYD*rfield.anu[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ip6p].ipIsoLevNIonCon-1]);
        }
        else
        {
                PrintE93( ioQQQ, 0.);
        }
        fprintf( ioQQQ, "   Tbr(1Ryd):");
        PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1]*TE1RYD*rfield.anu[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]);
        fprintf( ioQQQ, "   Tbr(4R):");
        PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1]*TE1RYD*rfield.anu[iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1]);
        fprintf( ioQQQ, "     Tbr (Nu-hi):");
        PrintE93( ioQQQ, rfield.OccNumbIncidCont[rfield.nflux-1]*TE1RYD*rfield.anu[rfield.nflux-1]);
        fprintf( ioQQQ, "\n");
 
        if( tradio > 1e9 )
        {
                fprintf( ioQQQ, 
                        " >>>The radio brightness temperature is very large,%10.2eK at%10.2ecm.  Is this physical???\n", 
                  tradio, 9.115e-6/rfield.anu[0] );
        }
 
        /* skip a line */
        fprintf( ioQQQ, "  \n" );
        return;
}