/* 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 */ /*ContSetIntensity derive intensity of incident continuum */ /*extin do extinction of incident continuum as set by extinguish command */ /*sumcon sums L and Q for net incident continuum */ /*ptrcer show continuum pointers in real time following drive pointers command */ /*conorm normalize continuum to proper intensity */ /*qintr integrates Q for any continuum between two limits, used for normalization */ /*pintr integrates L for any continuum between two limits, used for normalization */ #include "cddefines.h" #include "physconst.h" #include "iso.h" #include "noexec.h" #include "ionbal.h" #include "hextra.h" #include "trace.h" #include "dense.h" #include "oxy.h" #include "conv.h" #include "prt.h" #include "heavy.h" #include "rfield.h" #include "phycon.h" #include "called.h" #include "hydrogenic.h" #include "timesc.h" #include "secondaries.h" #include "opacity.h" #include "thermal.h" #include "ipoint.h" #include "atmdat.h" #include "rt.h" #include "radius.h" #include "geometry.h" #include "grainvar.h" #include "continuum.h" #include "taulines.h" /* these are weights used for continuum integration */ static const double aweigh[4]={-0.4305682,-0.1699905, 0.1699905, 0.4305682}; static const double fweigh[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274}; /*conorm normalize continuum to proper intensity */ STATIC void conorm(); /*pintr integrates L for any continuum between two limits, used for normalization */ STATIC double pintr(double penlo, double penhi); /*qintr integrates Q for any continuum between two limits, used for normalization */ STATIC double qintr(double *qenlo, double *qenhi); /*sumcon sums L and Q for net incident continuum */ STATIC void sumcon(long int il, long int ih, realnum *q, realnum *p, realnum *panu); /*extin do extinction of incident continuum as set by extinguish command */ STATIC void extin(realnum *ex1ryd); /*ptrcer show continuum pointers in real time following drive pointers command */ STATIC void ptrcer(); void IncidentContinuumHere() { DEBUG_ENTRY( "IncidentContinuumHere()" ); double frac_beam_time; /* fraction of beamed continuum that is constant */ double frac_beam_const; /* fraction of continuum that is isotropic */ double frac_isotropic; double BigLog = 0.; for( long i = 0; iSMALLFLOAT && flux_now>SMALLFLOAT ) { ratio = fabs( log10( flux_now / flux_org ) ); BigLog = max( ratio , BigLog ); } } if( BigLog > 0.01 ) fprintf(ioQQQ , "DEBUG diff continua %.2e\n", BigLog ); return; } /* called by Cloudy to set up continuum */ void ContSetIntensity() { bool lgCheckOK; long int i, ip, j, n; realnum EdenHeav, ex1ryd, factor, occ1, p, p1, p2, p3, p4, p5, p6, p7, p8, pgn, phe, pheii, qgn; realnum xIoniz; double HCaseBRecCoeff, wanu[4], alf, bet, fntest, fsum, ecrit, tcompr, tcomp, RatioIonizToRecomb, r3ov2; double amean, amean2, amean3, peak, wfun[4]; /* fraction of beamed continuum that is varies with time */ double frac_beam_time , frac_beam_time1; /* fraction of beamed continuum that is constant */ double frac_beam_const , frac_beam_const1; /* fraction of continuum that is isotropic */ double frac_isotropic , frac_isotropic1; long int nelem , ion; DEBUG_ENTRY( "ContSetIntensity()" ); /* set continuum */ if( trace.lgTrace ) { fprintf( ioQQQ, " ContSetIntensity called.\n" ); } /* find normalization factors for the continua - this decides whether continuum is * per unit area of luminosity, and which is desired final product */ conorm(); /* define factors to convert rfeld.flux array into photon occupation array OCCNUM * by multiplication */ factor = (realnum)(EN1RYD/PI4/FR1RYD/HNU3C2); /*------------------------------------------------------------- */ lgCheckOK = true; for( i=0; i < rfield.nupper; i++ ) { /* this was original anu array with no continuum averaging */ rfield.anu[i] = rfield.AnuOrg[i]; rfield.ContBoltz[i] = 0.; fsum = 0.; amean = 0.; amean2 = 0.; amean3 = 0.; frac_beam_time = 0.; frac_beam_const = 0.; frac_isotropic = 0.; for( j=0; j < 4; j++ ) { /* aweigh is symmetric about 0.5 widflx */ wanu[j] = rfield.anu[i] + rfield.widflx[i]*aweigh[j]; /* >>chng 02 jul 16, add test on continuum limits - * this was exceeded when resolution set very large */ wanu[j] = MAX2( wanu[j] , rfield.emm ); wanu[j] = MIN2( wanu[j] , rfield.egamry ); /* >>chng 06 feb 03, the continuum binning can change dramatically * at some energies - make sure that this cell does not overextend the * boundaries of its neighbors */ if( i > 0 && i < rfield.nupper-1 ) { wanu[j] = MAX2( wanu[j] , rfield.anu[i-1] + 0.5*rfield.widflx[i-1] ); wanu[j] = MIN2( wanu[j] , rfield.anu[i+1] - 0.5*rfield.widflx[i+1] ); } wfun[j] = fweigh[j]*ffun(wanu[j] , &frac_beam_time1 , &frac_beam_const1 , &frac_isotropic1 ); /*if( i==76 ) fprintf(ioQQQ,"DEBUG ffun %li %.4e at %.4e R\n", j , ffun(wanu[j]) , wanu[j] );*/ fsum += wfun[j]; amean += wanu[j]*wfun[j]; amean2 += wanu[j]*wanu[j]*wfun[j]; amean3 += wanu[j]*wanu[j]*wanu[j]*wfun[j]; frac_beam_time += fweigh[j]*frac_beam_time1; frac_beam_const += fweigh[j]*frac_beam_const1; frac_isotropic += fweigh[j]*frac_isotropic1; } ASSERT( fabs( 1.-frac_beam_time-frac_beam_const-frac_isotropic)< 10.*FLT_EPSILON); /* This is a fix for ticket #1 */ if( fsum*rfield.widflx[i] > BIGFLOAT ) { fprintf( ioQQQ, "\n Cannot continue. The continuum is far too intense.\n" ); for( j=0; j < rfield.nShape; j++ ) { if( (wfun[j]*rfield.widflx[i] > BIGFLOAT) && ( rfield.nShape > 1 ) ) { fprintf( ioQQQ, " Problem is with source number %li\n", j ); break; } } cdEXIT(EXIT_FAILURE); } rfield.flux[0][i] = (realnum)(fsum*rfield.widflx[i]); /* save separation into isotropic constant and beamed, and possibly * time-variable beamed continua */ rfield.flux_beam_const[i] = rfield.flux[0][i] * (realnum)frac_beam_const; rfield.flux_beam_time[i] = rfield.flux[0][i] * (realnum)frac_beam_time; rfield.flux_isotropic[i] = rfield.flux[0][i] * (realnum)frac_isotropic; if( rfield.flux[0][i] > 0. ) { /*fprintf(ioQQQ,"DEBUG %4li %e %e %.3e %.3e\n", i , rfield.anu[i] , (amean2/amean) , amean2 , amean );*/ rfield.anu[i] = (realnum)(amean2/amean); rfield.anu2[i] = (realnum)(amean3/amean); /* mesh must be strictly monotonically increasing - make it so */ if( i > 0 && rfield.anu[i] <= rfield.anu[i-1] ) { /* prevent roundoff from allowing i cell to lie below i-1 * cell when continuum mesh is very fine. */ /* use 2*epsilon to protect against unusual rounding modes */ rfield.anu[i] = rfield.anu[i-1]*(1.f+2.f*FLT_EPSILON); rfield.anu2[i] = pow2(rfield.anu[i]); } ASSERT( i==0 || rfield.anu[i] > rfield.anu[i-1] ); /* define array of LOG10( nu(ryd) ) these are not guaranteed to * be monotonically increasing due to loss of precision in log * of float */ rfield.anulog[i] = (realnum)log10(rfield.anu[i]); } else if( rfield.flux[0][i] == 0. ) { rfield.anu2[i] = rfield.anu[i]*rfield.anu[i]; rfield.anulog[i] = (realnum)log10(rfield.anu[i]); } else { rfield.anu2[i] = rfield.anu[i]*rfield.anu[i]; fprintf( ioQQQ, " negative continuum returned at%6ld%10.2e%10.2e\n", i, rfield.anu[i], rfield.flux[0][i] ); lgCheckOK = false; } rfield.anu3[i] = rfield.anu2[i]*rfield.anu[i]; rfield.ConEmitReflec[0][i] = 0.; rfield.ConEmitOut[0][i] = 0.; rfield.convoc[i] = factor/rfield.widflx[i]/rfield.anu2[i]; /* following are Compton exchange factors from Tarter */ alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i])); bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/ 4.; rfield.csigh[i] = (realnum)(alf*rfield.anu2[i]*3.858e-25); rfield.csigc[i] = (realnum)(alf*bet*rfield.anu[i]*3.858e-25); rfield.lgMeshSetUp = true; } /*i = 76; fprintf(ioQQQ,"\nDEBUG %4li %e \n", i , rfield.anu[i] );*/ /* >>chng 05 jul 01 add this * finished with stored continua - return these vectors */ #if 0 /* commented out since we must conserve energy, and continuum was set with old widflx */ /* now fix widflx array so that it is correct */ for( i=1; i (realnum)peak ) { /* prt.ipeak points to largest f_nu at H-ionizing energies * and is passed to other parts of code */ /* i+1 to keep ipeak on fortran version of energy array */ prt.ipeak = i+1; peak = rfield.flux[0][i]*rfield.anu[i]/rfield.widflx[i]; } } /* find highest energy to consider in continuum flux array * peak is the peak product nu*flux */ peak = rfield.flux[0][prt.ipeak-1]/rfield.widflx[prt.ipeak-1]* rfield.anu2[prt.ipeak-1]; /* say what type of cpu this is, if desired */ if( trace.lgTrace ) { fprintf( ioQQQ, " ContSetIntensity: The peak of the H-ion continuum is at %.3e Ryd - its value is %.2e\n", rfield.anu[prt.ipeak-1] , peak); } if( peak > 1e38 ) { fprintf( ioQQQ, " PROBLEM DISASTER The continuum is too intense to compute. Use a fainter continuum. (This is the nu*f_nu test)\n" ); fprintf( ioQQQ, " Sorry.\n" ); cdEXIT(EXIT_FAILURE); } /* FluxFaint set in zero.c; normally 1e-10 */ /* this will be faintest level of continuum we want to consider. * peak was set above, is peak of hydrogen ionizing radiation field, * and is zero if no H-ionizing radiation */ fntest = peak*rfield.FluxFaint; { enum {DEBUG_LOC=false}; /* print flux array then quit */ if( DEBUG_LOC ) { for( i=0; i 0. ) { /* this test is not done in pdr conditions where NO H-ionizing radiation, * since fntest is zero*/ i = rfield.nupper; /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs * where continuum barely goes into H-ionizing radiation */ while( i > prt.ipeak+3 && rfield.flux[0][i-1]*rfield.anu2[i-1]/rfield.widflx[i-1] < (realnum)fntest ) { --i; } } else { /* when no H-ionizing radiation set to Lyman edge */ /* >>chng 05 aug 16, change ipeak to ipeak+3, to avoid possible one-off bugs * where continuum barely goes into H-ionizing radiation */ i = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+3; } /* * this line of code dates from 1979 and IOA Cambridge. removed july 19 95 * I think it was the last line of the original Cambridge Fortran source nflux = MAX( ineon(1)+4 , i ) */ /* >>chng 99 apr 28, reinstate the rfield.FluxFaint limit with nflux */ rfield.nflux = i; /* trim down nflux, was set to rfield.nupper, the dimension of all vectors, in zero.c, * in ContCreatePointers was set to nupper, the number of cells needed to get up to the * high-energy limit of the code */ while( rfield.flux[0][rfield.nflux-1] < SMALLFLOAT && rfield.nflux > 1 ) { --rfield.nflux; } /* make sure we go high enough, avoid 1-off bugs as in draine field */ ++rfield.nflux; if( rfield.nflux == 1 ) { fprintf( ioQQQ, " PROBLEM DISASTER This incident continuum appears to have no radiation.\n" ); fprintf( ioQQQ, " Sorry.\n" ); cdEXIT(EXIT_FAILURE); } /* >>chng 04 oct 10, add this limit - arrays will malloc to nupper, but will add unit * continuum to [nflux] - this must be within array */ rfield.nflux = MIN2( rfield.nupper-1 , rfield.nflux ); /* >>chng 05 mar 01, nfine should not extend beyond continuum * make sure fine opacity scale does not extend beyond continuum we will use */ rfield.nfine = rfield.nfine_malloc; while( rfield.nfine > 0 && rfield.fine_anu[rfield.nfine-1] > rfield.anu[rfield.nflux-1] ) { --rfield.nfine; } /* check that continuum defined everywhere - look for zero's and comment if present */ continuum.lgCon0 = false; ip = 0; for( i=0; i < rfield.nflux; i++ ) { if( rfield.flux[0][i] == 0. ) { if( called.lgTalk && !continuum.lgCon0 ) { fprintf( ioQQQ, " NOTE Setcon: continuum has zero intensity starting at %11.4e Ryd.\n", rfield.anu[i] ); continuum.lgCon0 = true; } ++ip; } } if( continuum.lgCon0 && called.lgTalk ) { fprintf( ioQQQ, "%6ld cells in the incident continuum have zero intensity. Problems???\n\n", ip ); } /* check for devastating error in the continuum mesh or intensity */ lgCheckOK = true; for( i=1; i < rfield.nflux; i++ ) { if( rfield.flux[0][i] < 0. ) { fprintf( ioQQQ, " PROBLEM DISASTER Continuum has negative intensity at %.4e Ryd=%.2e %4.4s %4.4s\n", rfield.anu[i], rfield.flux[0][i], rfield.chLineLabel[i], rfield.chContLabel[i] ); lgCheckOK = false; } else if( rfield.anu[i] <= rfield.anu[i-1] ) { fprintf( ioQQQ, " PROBLEM DISASTER cont_setintensity - internal error - continuum energies not in increasing order: energies follow\n" ); fprintf( ioQQQ, "%ld %e %ld %e %ld %e\n", i -1 , rfield.anu[i-1], i, rfield.anu[i], i +1, rfield.anu[i+1] ); lgCheckOK = false; } } /* either of the ways lgCheckOK would be set true would be a major internal error */ if( !lgCheckOK ) { ShowMe(); cdEXIT(EXIT_FAILURE); } /* turn off recoil ionization if high energy < 190R */ if( rfield.anu[rfield.nflux-1] <= 190 ) { ionbal.lgCompRecoil = false; } /* sum photons and energy, save mean */ /* sum from low energy to Balmer edge */ sumcon(1,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon-1,&rfield.qrad,&prt.pradio,&p1); /* sum over Balmer continuum */ sumcon(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2p].ipIsoLevNIonCon,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1,&rfield.qbal,&prt.pbal,&p2); /* sum from Lyman edge to HeI edge */ sumcon(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon, iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1,&prt.q,&p,&p3); /* sum from HeI to HeII edges */ sumcon(iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon, iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1,&rfield.qhe,&phe,&p4); /* sum from Lyman edge to carbon k-shell */ sumcon(iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon,opac.ipCKshell-1,&rfield.qheii,&pheii,&p5); /* sum from c k-shell to gamma ray - where pairs start */ sumcon( opac.ipCKshell , rfield.ipEnerGammaRay-1 , &prt.qx, &prt.xpow , &p6); /* complete sum up to high-energy limit */ sumcon(rfield.ipEnerGammaRay,rfield.nflux,&prt.qgam,&prt.GammaLumin, &p7); /* find to estimate photoerosion timescale */ n = ipoint(7.35e5); sumcon(n,rfield.nflux,&qgn,&pgn,&p8); timesc.TimeErode = qgn; /* find Compton temp */ tcompr = (p1 + p2 + p3 + p4 + p5 + p6 + p7)/(prt.pradio + prt.pbal + p + phe + pheii + prt.xpow + prt.GammaLumin); tcomp = tcompr/(4.*6.272e-6); if( trace.lgTrace ) { fprintf( ioQQQ, " mean photon energy=%10.3eR =%10.3eK low, high nu=%12.4e%12.4e\n", tcompr, tcomp, rfield.anu[0] - rfield.widflx[0]/2., rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2. ); } /* this is total power in ionizing radiation, per unit area */ prt.powion = p + phe + pheii + prt.xpow + prt.GammaLumin; /* this is the total radiation intensity, erg cm-2 s-1 */ continuum.TotalLumin = prt.pradio + prt.powion + prt.pbal; /* this is placed into the line stack on the first zone, then * reset to zero, to end up with luminosity in the emission lines array. * at end of iteration it is reset to TotalLumin */ continuum.totlsv = continuum.TotalLumin; /* total H-ionizing photon number, */ rfield.qhtot = prt.q + rfield.qhe + rfield.qheii + prt.qx + prt.qgam; /* ftotal photon number, all energies */ rfield.qtot = rfield.qhtot + rfield.qbal + rfield.qrad; if( prt.powion <= 0. && called.lgTalk ) { rfield.lgHionRad = true; fprintf( ioQQQ, " NOTE There is no hydrogen-ionizing radiation.\n" ); fprintf( ioQQQ, " Was this intended?\n\n" ); /* check if any Balmer ionizing radiation since even metals will be * totally neutral */ if( prt.pbal <=0. && called.lgTalk ) { fprintf( ioQQQ, " NOTE There is no Balmer continuum radiation.<<<<\n" ); fprintf( ioQQQ, " Was this intended?\n\n" ); } } else { rfield.lgHionRad = false; } /* option to add energy deposition due to fast neutrons, * entered as fraction of total photon luminosity */ if( hextra.lgNeutrnHeatOn ) { /* hextra.totneu is erg cm-2 s-1 in neutrons * hextra.effneu - efficiency default is unity */ hextra.totneu = (realnum)(hextra.effneu*continuum.TotalLumin* pow((realnum)10.f,hextra.frcneu)); } else { hextra.totneu = (realnum)0.; } /* temp correspond to energy density, printed in STARTR */ phycon.TEnerDen = pow(continuum.TotalLumin/SPEEDLIGHT/7.56464e-15,0.25); /* sanity check for single blackbody, that energy density temperature * is smaller than black body temperature */ if( rfield.nShape==1 && strcmp( rfield.chSpType[rfield.nShape-1], "BLACK" )==0 ) { /* single black body, now confirm that TEnerDen is less than this temperature, * >>>chng 99 may 02, * in lte these are very very close, factor of 1.00001 allows for numerical * errors, and apparently slightly different atomic coef in different parts * of code. eventaully all mustuse physonst.h and agree exactly */ if( phycon.TEnerDen > 1.0001f*rfield.slope[rfield.nShape-1] ) { fprintf( ioQQQ, "\n WARNING: The energy density temperature (%g) is greater than the" " black body temperature (%g). This is unphysical.\n\n", phycon.TEnerDen , rfield.slope[rfield.nShape-1]); } } /* incident continuum nu*f_nu at Hbeta and Ly-alpha */ continuum.cn4861 = (realnum)(ffun(0.1875)*HPLANCK*FR1RYD*0.1875*0.1875); continuum.cn1216 = (realnum)(ffun(0.75)*HPLANCK*FR1RYD*0.75*0.75); continuum.sv4861 = continuum.cn4861; continuum.sv1216 = continuum.cn1216; /* flux density in V, erg / s / cm2 / hz */ continuum.fluxv = (realnum)(ffun(0.1643)*HPLANCK*0.1643); continuum.fbeta = (realnum)(ffun(0.1875)*HPLANCK*0.1875*6.167e14); /* flux density nu*Fnu = erg / s / cm2 * EX1RYD is optional extinction factor at 1 ryd */ prt.fx1ryd = (realnum)(ffun(1.000)*HPLANCK*ex1ryd*FR1RYD); realnum plsFrqConstant = (realnum)(ELEM_CHARGE_ESU/sqrt(PI*ELECTRON_MASS)/FR1RYD); ASSERT( plsFrqConstant > 2.7e-12 && plsFrqConstant < 2.8e-12 ); /* check for plasma frequency - then zero out incident continuum * for energies below this * this is critical electron density, shielding of incident continuum * if electron density is greater than this */ ecrit = POW2(rfield.anu[0]/plsFrqConstant); if( dense.gas_phase[ipHYDROGEN]*1.2 > ecrit ) { rfield.lgPlasNu = true; // This should be electron density, but we don't know that yet. // We use n_H (with 1.2 factor for He) as an ansatz. rfield.plsfrq = plsFrqConstant*sqrt(dense.gas_phase[ipHYDROGEN]*1.2f); rfield.plsfrqmax = rfield.plsfrq; rfield.ipPlasma = ipoint(rfield.plsfrq); /* save max pointer too */ rfield.ipPlasmax = rfield.ipPlasma; /* now loop over all affected energies, setting incident continuum * to zero there, and counting all as reflected */ /* >>chng 01 jul 14, from i < ipPlasma to ipPlasma-1 - * when ipPlasma is 1 plasma freq is not on energy scale */ for( i=0; i < rfield.ipPlasma-1; i++ ) { /* count as reflected incident continuum */ rfield.ConRefIncid[0][i] = rfield.flux[0][i]; /* set continuum to zero there */ rfield.flux_beam_const[i] = 0.; rfield.flux_beam_time[i] = 0.; rfield.flux_isotropic[i] = 0.; rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] + rfield.flux_isotropic[i]; } } else { rfield.lgPlasNu = false; /* >>chng 01 jul 14, from 0 to 1 - 1 is the first array element on the F scale, * ipoint would return this, so rest of code assumes ipPlasma is 1 plus correct index */ rfield.ipPlasma = 1; rfield.plsfrqmax = 0.; rfield.plsfrq = 0.; } if( rfield.ipPlasma > 1 && called.lgTalk ) { fprintf( ioQQQ, " !The plasma frequency is %.2e Ryd. The incident continuum is set to 0 below this.\n", rfield.plsfrq ); } rfield.occmax = 0.; rfield.tbrmax = 0.; for( i=0; i < rfield.nupper; i++ ) { /* set up occupation number array */ rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i]; if( rfield.OccNumbIncidCont[i] > rfield.occmax ) { rfield.occmax = rfield.OccNumbIncidCont[i]; rfield.occmnu = rfield.anu[i]; } /* following product is continuum brightness temperature */ if( rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i] > rfield.tbrmax ) { rfield.tbrmax = (realnum)(rfield.OccNumbIncidCont[i]*TE1RYD*rfield.anu[i]); rfield.tbrmnu = rfield.anu[i]; } /* save continuum for next iteration */ rfield.flux_total_incident[0][i] = rfield.flux[0][i]; rfield.flux_beam_const_save[i] = rfield.flux_beam_const[i]; rfield.flux_time_beam_save[i] = rfield.flux_beam_time[i]; rfield.flux_isotropic_save[i] = rfield.flux_isotropic[i]; /*fprintf(ioQQQ,"DEBUG type cont %li\t%.3e\t%.2e\t%.2e\t%.2e\t%.2e\n", i, rfield.anu[i], rfield.flux[0][i],rfield.flux_beam_const[i],rfield.flux_beam_time[i], rfield.flux_isotropic[i]); fflush(ioQQQ);*/ } /* if continuum brightness temp is large, where does it fall below * 1e4k??? */ if( rfield.tbrmax > 1e4 ) { i = ipoint(rfield.tbrmnu)-1; while( i < rfield.nupper-1 && (rfield.OccNumbIncidCont[i]*TE1RYD* rfield.anu[i] > 1e4) ) { ++i; } rfield.tbr4nu = rfield.anu[i]; } else { rfield.tbr4nu = 0.; } /* if continuum occ num is large, where does it fall below 1? */ if( rfield.occmax > 1. ) { i = ipoint(rfield.occmnu)-1; while( i < rfield.nupper && (rfield.OccNumbIncidCont[i] > 1.) ) { ++i; } rfield.occ1nu = rfield.anu[i]; } else { rfield.occ1nu = 0.; } /* remember if incident radiation field is less than 10*Habing ISM */ /* >>chng 06 aug 01, change this test from continuum.TotalLumin to * energy in balmer and ionizing continuum, since this is the true habing field * and is the continuum that interacts with gas. When CMB set this * tests on total did not trigger due to cold blackbody, which has little * effect on gas, other than compton */ if( (prt.powion + prt.pbal) < 1.8e-2 ) { /* thermal.ConstTemp def is zero, set pos when constant temperature is set */ rfield.lgHabing = true; /* >>chng 06 aug 01 also print warning if substantially below Habing, this may be a mistake */ if( ((prt.powion + prt.pbal) < 1.8e-12) && /* this is test for not constant temperature */ (!thermal.lgTemperatureConstant) ) { fprintf( ioQQQ, "\n >>>\n" " >>> NOTE The incident continuum is surprisingly faint.\n" ); fprintf( ioQQQ, " >>> The total energy in the Balmer and Lyman continua is %.2e erg cm-2 s-1.\n" ,(prt.powion + prt.pbal)); fprintf( ioQQQ, " >>> This is many orders of magnitude fainter than the ISM galactic background.\n" ); fprintf( ioQQQ, " >>> This seems unphysical - please check that the continuum intensity has been properly set.\n" ); fprintf( ioQQQ, " >>> YOU MAY BE MAKING A BIG MISTAKE!!\n >>>\n\n\n\n" ); } } /* fix ionization parameter (per hydrogen) at inner edge */ rfield.uh = (realnum)(rfield.qhtot/ (dense.gas_phase[ipHYDROGEN]*SPEEDLIGHT)); rfield.uheii = (realnum)((rfield.qheii + prt.qx)/ (dense.gas_phase[ipHYDROGEN]*SPEEDLIGHT)); if( rfield.uh > 1e10 ) { fprintf( ioQQQ, "\n\n" " CAUTION The incident radiation field is surprisingly intense.\n" ); fprintf( ioQQQ, " The dimensionless hydrogen ionization parameter is %.2e.\n" , rfield.uh ); fprintf( ioQQQ, " This is many orders of magnitude brighter than commonly seen.\n" ); fprintf( ioQQQ, " This seems unphysical - please check that the radiation field intensity has been properly set.\n" ); fprintf( ioQQQ, " YOU MAY BE MAKING A BIG MISTAKE!!\n\n\n\n\n" ); } /* guess first temperature and neutral h density */ double TeNew; if( thermal.ConstTemp > 0. ) { TeNew = thermal.ConstTemp; } else { if( rfield.uh > 0. ) { TeNew = (20000.+log10(rfield.uh)*5000.); TeNew = MAX2(8000. , TeNew ); } else { TeNew = 1000.; } } TempChange( TeNew ); /* this is an option to stop after printing header only */ if( noexec.lgNoExec ) return; /* estimate secondary ionization rate - probably 0, but possible extra * SetCsupra set with "set csupra" */ /* >>>chng 99 apr 29, added cosmic ray ionization since this is used to get * helium ionization fraction, and was zero in pdr, so He turned off at start, * and never turned back on */ /* coef on cryden is from highen.c */ for( nelem=ipHYDROGEN; nelemSMALLFLOAT ) PhotonEnergy = prt.powion/rfield.qhtot/EN1RYD; double OtherIonization = rfield.qhtot*6.3e-18/POW3(PhotonEnergy) + secondaries.csupra[ipHYDROGEN][0]; double newEden = dense.eden; long loopCount = 0; do { /* update electron density */ EdenChange( newEden ); double RatioIoniz = (CollIoniz*dense.eden+OtherIonization)/(HCaseBRecCoeff*dense.eden); if( RatioIoniz<1e-3 ) { /* very low ionization solution */ dense.xIonDense[ipHYDROGEN][1] = (realnum)( dense.gas_phase[ipHYDROGEN]*RatioIoniz); dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1]; ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. && dense.xIonDense[ipHYDROGEN][0]<=dense.gas_phase[ipHYDROGEN]); ASSERT( dense.xIonDense[ipHYDROGEN][1]>=0. && dense.xIonDense[ipHYDROGEN][1]1e3 ) { /* very high ionization solution */ dense.xIonDense[ipHYDROGEN][0] = (realnum)( dense.gas_phase[ipHYDROGEN]/RatioIoniz); dense.xIonDense[ipHYDROGEN][1] = dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][0]; ASSERT( dense.xIonDense[ipHYDROGEN][0]>=0. && dense.xIonDense[ipHYDROGEN][0]=0. && dense.xIonDense[ipHYDROGEN][1]<=dense.gas_phase[ipHYDROGEN]); //fprintf(ioQQQ,"DEBUG br 2 H0 %.2e rat %.2e\n", dense.xIonDense[ipHYDROGEN][0], // dense.gas_phase[ipHYDROGEN]/RatioIoniz); } else { /* intermediate ionization - solve quadratic */ double alpha = HCaseBRecCoeff + CollIoniz ; double beta = HCaseBRecCoeff*EdenExtraLocal + OtherIonization + (EdenExtraLocal - dense.gas_phase[ipHYDROGEN])*CollIoniz; double gamma = -dense.gas_phase[ipHYDROGEN]*(OtherIonization+EdenExtraLocal*CollIoniz); double discriminant = POW2(beta) - 4.*alpha*gamma; if( discriminant <0 ) { /* oops */ fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found negative discriminant.\n"); TotalInsanity(); } dense.xIonDense[ipHYDROGEN][1] = (realnum)((-beta + sqrt(discriminant))/(2.*alpha)); if( dense.xIonDense[ipHYDROGEN][1]> dense.gas_phase[ipHYDROGEN] ) { /* oops */ fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H+)>n(H).\n"); TotalInsanity(); } dense.xIonDense[ipHYDROGEN][0] = dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1]; if( dense.xIonDense[ipHYDROGEN][0]<= 0 ) { /* oops */ fprintf(ioQQQ," DISASTER PROBLEM cont_initensity found n(H0)<0.\n"); TotalInsanity(); } //fprintf(ioQQQ,"DEBUG br 3 H0 %.2e\n", dense.xIonDense[ipHYDROGEN][0]); } if( dense.xIonDense[ipHYDROGEN][1] > 1e-30 ) { iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = dense.xIonDense[ipHYDROGEN][0]; } else { iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = 0.; } /* now save estimates of whether induced recombination is going * to be important -this is a code pacesetter since GammaBn is slower * than GammaK */ hydro.lgHInducImp = false; for( i=ipH1s; i < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max; i++ ) { if( rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[i].ipIsoLevNIonCon-1] > 0.01 ) hydro.lgHInducImp = true; } /******************************************************************* * * * estimate helium's level of ionization * * * *******************************************************************/ /* only if helium is turned on */ if( dense.lgElmtOn[ipHELIUM] ) { /* next estimate level of helium singly ionized */ xIoniz = (realnum)t_ADfA::Inst().coll_ion_wrapper(1,0,phycon.te); /* >>chng 05 jan 05, add cosmic rays */ xIoniz = (realnum)(xIoniz*dense.eden + rfield.qhe*1e-18 + secondaries.csupra[ipHELIUM][1] ); double RecTot = HCaseBRecCoeff*dense.eden; RatioIonizToRecomb = xIoniz/RecTot; /* now estimate level of helium doubly ionized */ xIoniz = (realnum)t_ADfA::Inst().coll_ion_wrapper(1,1,phycon.te); /* >>chng 05 jan 05, add cosmic rays */ xIoniz = (realnum)(xIoniz*dense.eden + rfield.qheii*1e-18 + ionbal.CosRayIonRate ); /* rough charge dependence */ RecTot *= 4.; r3ov2 = xIoniz/RecTot; /* now set level of helium ionization */ if( RatioIonizToRecomb > 0. ) { double r1 = dense.gas_phase[ipHELIUM]/(1./RatioIonizToRecomb + 1. + r3ov2); dense.xIonDense[ipHELIUM][1] = (realnum)(r1); dense.xIonDense[ipHELIUM][0] = (realnum)(r1/RatioIonizToRecomb); dense.xIonDense[ipHELIUM][2] = (realnum)(r1*r3ov2); } else { /* no He ionizing radiation */ dense.xIonDense[ipHELIUM][1] = dense.xIonDense[ipHELIUM][2] = 0.; dense.xIonDense[ipHELIUM][0] = dense.gas_phase[ipHELIUM]; } if( dense.xIonDense[ipHELIUM][2] > 1e-30 ) { iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = dense.xIonDense[ipHELIUM][1]; } else { iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = 0.; } } else { /* case where helium is turned off */ dense.xIonDense[ipHELIUM][1] = 0.; dense.xIonDense[ipHELIUM][0] = 0.; dense.xIonDense[ipHELIUM][2] = 0.; iso_sp[ipH_LIKE][ipHELIUM].st[ipH1s].Pop() = 0.; } /* update electron density */ newEden = dense.xIonDense[ipHYDROGEN][1] + EdenExtraLocal + dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2]; loopCount++; } /* repeat above until guessed and calculated eden agree to at least 0.1%. */ while( loopCount < 10 && fabs(newEden/dense.eden - 1.) > 0.001 ); if( dense.xIonDense[ipHYDROGEN][0]<0.) TotalInsanity(); else if( dense.xIonDense[ipHYDROGEN][0] == 0. ) { fprintf(ioQQQ,"PROBLEM the derived atomic hydrogen density is zero.\n"); if( dense.gas_phase[ipHYDROGEN]<1e-5 && rfield.uh > 1e10) { fprintf(ioQQQ,"This is almost certainly due to floating point " "limits on this computer.\nThe ionization parameter is very large," " the density is very small,\nand the H^0 density cannot be" " stored in a float.\n"); //cdEXIT( EXIT_FAILURE ); } } //ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.); /* update electron density */ EdenChange( newEden ); if( dense.eden <= SMALLFLOAT ) { /* no electrons! */ fprintf(ioQQQ,"\n PROBLEM DISASTER - this simulation has no source" " of ionization. The electron density is zero. Consider " "adding a source of ionization such as cosmic rays.\n\n"); cdEXIT(EXIT_FAILURE); } /* fix range of stages of ionization */ for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) { if( dense.lgElmtOn[nelem] ) { // IonHigh[n] is the highest stage of ionization present // the IonHigh array index is on the C scale, so [0] is hydrogen // the value is also on the C scale, so element [nelem] can range // from 0 to nelem+1 dense.IonHigh[nelem] = nelem + 1; dense.IonLow[nelem] = 0; // for very intense radiation fields very heavy elements, N>Fe, will fail // in ion_solver due to ill conditioned matrix. all populations are in // fully stripped state. Start will fully stripped ion distribution in this case. if( rfield.uh > 1e15 ) { //trim down highest stage to be within incident radiation field while ( rfield.anu[Heavy.ipHeavy[nelem][dense.IonHigh[nelem]-1]] > rfield.anu[rfield.nflux] && dense.IonHigh[nelem]>1 ) --dense.IonHigh[nelem]; dense.IonLow[nelem] = max( 0 , dense.IonHigh[nelem]-1 ); } /* >>chng 04 jan 13, add this test, caught by Orly Gnat */ /* check on actual zero abundances of lower stages - this will only * happen when ionization is set with element ionization command */ if( dense.lgSetIoniz[nelem] ) { while( dense.SetIoniz[nelem][dense.IonLow[nelem]] < dense.density_low_limit ) ++dense.IonLow[nelem]; while( dense.SetIoniz[nelem][dense.IonHigh[nelem]] < dense.density_low_limit ) --dense.IonHigh[nelem]; } // make low-stage populations zero for( ion=0; iondense.IonHigh[nelem]; --ion ) { dense.xIonDense[nelem][dense.IonHigh[nelem]] += dense.xIonDense[nelem][ion]; dense.xIonDense[nelem][ion] = 0.; } } else { /* this element is turned off, make stages impossible */ dense.IonLow[nelem] = -1; dense.IonHigh[nelem] = -1; } } // make first estimate of iso continuum lowering for( long ipISO=ipH_LIKE; ipISO=1) { realnum low2dens = dense.xIonDense[nelem][0]+dense.xIonDense[nelem][1]; dense.xIonDense[nelem][0] = low2dens * atomFrac; dense.xIonDense[nelem][1] = low2dens * firstIonFrac; } for (long ion=1;ion>chng 05 jan 05, insure positive eden */ EdenChange( MAX2( SMALLFLOAT , dense.eden ) ); if( dense.EdenSet > 0. ) { EdenChange( dense.EdenSet ); } dense.EdenHCorr = dense.eden; dense.EdenHCorr_f = (realnum)dense.EdenHCorr; if( dense.eden < 0. ) { fprintf( ioQQQ, " PROBLEM DISASTER Negative electron density results in ContSetIntensity.\n" ); fprintf( ioQQQ, "%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", dense.eden, dense.xIonDense[ipHYDROGEN][1], dense.xIonDense[ipHELIUM][1], dense.xIonDense[ipHELIUM][2], dense.gas_phase[ipCARBON], dense.EdenExtra ); TotalInsanity(); } if( dense.EdenSet > 0. ) { dense.EdenTrue = dense.EdenSet; } else dense.EdenTrue = dense.eden; if( trace.lgTrace ) { fprintf( ioQQQ, " ContSetIntensity sets initial EDEN to %.4e, contributors H+=%.2e He+, ++= %.2e %.2e Heav %.2e extra %.2e\n", dense.eden , dense.xIonDense[ipHYDROGEN][1], dense.xIonDense[ipHELIUM][1], 2.*dense.xIonDense[ipHELIUM][2], EdenHeav, dense.EdenExtra); } // photon occupation number at 1 Ryd - used for printout occ1 = (realnum)(prt.fx1ryd/HNU3C2/PI4/FR1RYD); if( occ1 > 1. ) rfield.lgOcc1Hi = true; else rfield.lgOcc1Hi = false; if( trace.lgTrace && trace.lgConBug ) { fprintf(ioQQQ,"\ntrace continuum print of %li incident spectral " "components\n", rfield.nShape); fprintf(ioQQQ," # type Illum Beam? 1/cos TimeVary?\n"); for( rfield.ipSpec=0; rfield.ipSpec0 ); for( i=rfield.ipG0_DB96_lo; i < rfield.ipG0_DB96_hi; i++ ) sum += rfield.flux[0][i]; fprintf(ioQQQ,"\n Sum DB96 photons %.2e", sum); sum = 0.; for( i=(Heavy.ipHeavy[ipHYDROGEN][0] - 1); i continuum.filbnd[continuum.nrange] ) { fprintf( ioQQQ," ipoint would crash since highest desired energy of %e ryd is above limit of %e\n", rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2.*0.9 , continuum.filbnd[continuum.nrange-1] ); fprintf( ioQQQ," width of cell is %e\n",rfield.widflx[rfield.nflux]); fprintf( ioQQQ," this, previous cells are %e %e\n", rfield.anu[rfield.nflux-1],rfield.anu[rfield.nflux-2]); cdEXIT(EXIT_FAILURE); } /* this branch check everything, write errors to error file */ fprintf( ioQQQ, " errors output on errors.txt\n"); fprintf( ioQQQ, " IP(cor),IP(fount),nu lower, upper of found, desired cell.\n" ); /* error file not open, set to null so we can check later */ ioERRORS = NULL; for( i=0; i < rfield.nflux-1; i++ ) { t1 = rfield.anu[i] - rfield.widflx[i]/2.*0.9; t2 = rfield.anu[i] + rfield.widflx[i]/2.*0.9; j = ipoint(t1); if( j != i+1 ) { /* open file for errors if not already open */ if( ioERRORS == NULL ) { ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY ); fprintf( ioQQQ," created errors.txt file with error summary\n"); } fprintf( ioQQQ, " Pointers do not agree for lower bound of cell%4ld, %e\n", i, rfield.anu[i]); fprintf( ioERRORS, " Pointers do not agree for lower bound of cell%4ld, %e\n", i, rfield.anu[i] ); } j = ipoint(t2); if( j != i+1 ) { /* open file for errors if not already open */ if( ioERRORS == NULL ) { ioERRORS = open_data( "errors.txt", "w", AS_LOCAL_ONLY ); fprintf( ioQQQ," created errors.txt file with error summary\n"); } fprintf( ioQQQ, " Pointers do not agree for upper bound of cell%4ld, %e\n", i , rfield.anu[i]); fprintf( ioERRORS, " Pointers do not agree for upper bound of cell%4ld, %e\n", i , rfield.anu[i]); } } } else { fprintf( ioQQQ, "I do not understand this key, sorry. %c\n", chKey ); cdEXIT(EXIT_FAILURE); } if( ioERRORS!=NULL ) fclose( ioERRORS ); cdEXIT(EXIT_FAILURE); } /*extin do extinction of incident continuum as set by extinguish command */ STATIC void extin(realnum *ex1ryd) { DEBUG_ENTRY( "extin()" ); /* modify input continuum by leaky absorber * power law fit to * >>refer XUV extinction Cruddace, R., Paresce, F., Bowyer, S., & Lampton, M. 1974, ApJ, 187, 497. */ if( rfield.ExtinguishColumnDensity == 0. ) { *ex1ryd = 1.; for( long i=0; iRSFCheck, continuum.ResolutionScaleFactor ) ) { fprintf( ioQQQ,"\n\n PROBLEM DISASTER At least one of the grain opacity files " "has been compiled with a different energy grid resolution factor.\n" ); fprintf( ioQQQ, " Please recompile this file using the COMPILE GRAINS command " "and make sure that you use the correct SET CONTINUUM RESOLUTION factor.\n" ); cdEXIT(EXIT_FAILURE); } } /* default is is to predict line intensities, * but if any continuum specified as luminosity then would override this - * following two values are correct for intensities */ radius.pirsq = 0.; radius.lgPredLumin = false; /* check whether ANY luminosities are present */ for( i=0; i < rfield.nShape; i++ ) { if( strcmp(rfield.chRSpec[i],"4 PI") == 0 ) { radius.pirsq = (realnum)(1.0992099 + 2.*xLog_radius_inner); radius.lgPredLumin = true; /* convert down to intensity */ rfield.totpow[i] -= radius.pirsq; if( trace.lgTrace ) { fprintf( ioQQQ, " conorm converts continuum %ld from luminosity to intensity.\n", i ); } } } /* if total luminosities are present, must have specified a starting radius */ if( radius.lgPredLumin && !radius.lgRadiusKnown ) { fprintf(ioQQQ,"PROBLEM DISASTER conorm: - A continuum source was specified as a luminosity, but the inner radius of the cloud was not set.\n"); fprintf(ioQQQ,"Please set an inner radius.\nSorry.\n"); cdEXIT(EXIT_FAILURE); } /* convert ionization parameters to number of photons, called "q(h)" * at this stage q(h) and "PHI " are the same */ for( i=0; i < rfield.nShape; i++ ) { if( strcmp(rfield.chSpNorm[i],"IONI") == 0 ) { /* the log of the ionization parameter was stored here, this converts * it to the log of the number of photons per sq cm */ rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) + log10(SPEEDLIGHT); strcpy( rfield.chSpNorm[i], "Q(H)" ); if( trace.lgTrace ) { fprintf( ioQQQ, " conorm converts continuum %ld from ionizat par to q(h).\n", i ); } } } /* convert x-ray ionization parameter xi to intensity */ for( i=0; i < rfield.nShape; i++ ) { if( strcmp(rfield.chSpNorm[i],"IONX") == 0 ) { /* this converts it to an intensity */ rfield.totpow[i] += log10(dense.gas_phase[ipHYDROGEN]) - log10(PI4); strcpy( rfield.chSpNorm[i], "LUMI" ); if( trace.lgTrace ) { fprintf( ioQQQ, " conorm converts continuum%3ld from x-ray ionizat par to I.\n", i ); } } } /* indicate whether we ended up with luminosity or intensity */ if( trace.lgTrace ) { if( radius.lgPredLumin ) { fprintf( ioQQQ, " Cloudy will predict lumin into 4pi\n" ); } else { fprintf( ioQQQ, " Cloudy will do surface flux for lumin\n" ); } } /* if intensity per unit area is predicted then geometric * covering factor must be unity * variable can also be set elsewhere */ if( !radius.lgPredLumin ) { geometry.covgeo = 1.; } /* main loop over all continuum shapes to find continuum normalization * for each one */ for( i=0; i < rfield.nShape; i++ ) { rfield.ipSpec = i; /* check that, if laser, bounds include laser energy */ if( strcmp(rfield.chSpType[rfield.ipSpec],"LASER") == 0 ) { if( !( rfield.range[rfield.ipSpec][0] < rfield.slope[rfield.ipSpec] && rfield.range[rfield.ipSpec][1] > rfield.slope[rfield.ipSpec]) ) { fprintf(ioQQQ,"PROBLEM DISASTER, a continuum source is a laser at %f Ryd, but the intensity was specified over a range from %f to %f Ryd.\n", rfield.slope[rfield.ipSpec], rfield.range[rfield.ipSpec][0], rfield.range[rfield.ipSpec][1]); fprintf(ioQQQ,"Please specify the continuum flux where the laser is active.\n"); cdEXIT(EXIT_FAILURE); } } if( trace.lgTrace ) { long int jj; fprintf( ioQQQ, " conorm continuum number %ld is shape %s range is %.2e %.2e\n", i, rfield.chSpType[i], rfield.range[i][0], rfield.range[i][1] ); fprintf( ioQQQ, "the continuum points follow\n"); jj = 0; if( rfield.lgContMalloc[rfield.ipSpec] ) { while( rfield.tNu[rfield.ipSpec][jj].Ryd() != 0. && jj < 100 ) { fprintf( ioQQQ, "%li %e %e\n", jj , rfield.tNu[rfield.ipSpec][jj].Ryd(), rfield.tslop[rfield.ipSpec][jj] ); ++jj; } } } if( strcmp(rfield.chSpNorm[i],"RATI") == 0 ) { /* option to scale relative to previous continua * this must come first since otherwise may be too late * BUT ratio cannot be the first continuum source */ if( trace.lgTrace ) { fprintf( ioQQQ, " conorm this is ratio to 1st con\n" ); } /* check that this is not the first continuum source, we must ratio */ if( i == 0 ) { fprintf( ioQQQ, " I cant form a ratio if continuum is first source.\n" ); cdEXIT(EXIT_FAILURE); } /* first find photon flux and Q of previous continuum */ rfield.ipSpec -= 1; flx1 = ffun1(rfield.range[i][0])*rfield.spfac[rfield.ipSpec]* rfield.range[i][0]; /* check that previous continua were not zero where ratio is formed */ if( flx1 <= 0. ) { fprintf( ioQQQ, " Previous continua were zero where ratio is desired.\n" ); cdEXIT(EXIT_FAILURE); } /* return pointer to previous (correct) value, find F, Q */ rfield.ipSpec += 1; /* we want a continuum totpow as powerful, flx is now desired flx */ flx1 *= rfield.totpow[i]; /*. find flux of this new continuum at that point */ flx2 = ffun1(rfield.range[i][1])*rfield.range[i][1]; /* this is ratio of desired to actual */ rfield.spfac[i] = flx1/flx2; if( trace.lgTrace ) { fprintf( ioQQQ, " conorm ratio will set scale fac to%10.3e at%10.2e Ryd.\n", rfield.totpow[i], rfield.range[i][0] ); } } else if( strcmp(rfield.chSpNorm[i],"FLUX") == 0 ) { /* specify flux density * option to use arbitrary frequency or range */ f = ffun1(rfield.range[i][0]); /* make sure this is positive, could be zero if out of range of table, * or for various forms of insanity */ if( f<=SMALLDOUBLE ) { fprintf( ioQQQ, "\n\n PROBLEM DISASTER\n The intensity of continuum source %ld is non-positive at the energy used to normalize it (%.3e Ryd). Something is seriously wrong.\n", i , rfield.range[i][0]); /* is this a table? if so, is ffun within its bounds? */ if( strcmp(rfield.chSpType[i],"INTER") == 0 ) fprintf( ioQQQ, " This continuum shape given by a table of points - check that the intensity is specified at an energy within the range of that table.\n"); fprintf( ioQQQ, " Also check that the numbers used to specify the shape and intensity do not under or overflow on this cpu.\n\n"); cdEXIT(EXIT_FAILURE); } /* now convert to log in continuum units we shall use */ f = log10(f) + log10(rfield.range[i][0]*EN1RYD/FR1RYD); f = rfield.totpow[i] - f; rfield.spfac[i] = pow(10.,f); if( trace.lgTrace ) { fprintf( ioQQQ, " conorm will set log fnu to%10.3e at%10.2e Ryd. Factor=%11.4e\n", rfield.totpow[i], rfield.range[i][0], rfield.spfac[i] ); } } else if( strcmp(rfield.chSpNorm[i],"Q(H)") == 0 || strcmp(rfield.chSpNorm[i],"PHI ") == 0 ) { /* some type of photon density entered */ if( trace.lgTrace ) { fprintf( ioQQQ, " conorm calling qintr range=%11.3e %11.3e desired val is %11.3e\n", rfield.range[i][0], rfield.range[i][1] , rfield.totpow[i]); } /* the total number of photons over the specified range in * the arbitrary system of units that the code save the continuum shape */ ASSERT( rfield.range[i][0] < rfield.range[i][1] ); qentrd = qintr(&rfield.range[i][0],&rfield.range[i][1]); /* this is the log of the scale factor that must multiply the * continuum shape to get the final set of numbers */ diff = rfield.totpow[i] - qentrd; /* >>chng 03 mar 13, from diff < -25 to <-35, * tripped for very low U models used for H2 simulations */ /*if( diff < -25. || diff > 35. )*/ if( diff < -35. || diff > 35. ) { fprintf( ioQQQ, " PROBLEM DISASTER Continuum source specified is too extreme.\n" ); fprintf( ioQQQ, " The integral over the continuum shape gave (log) %.3e photons, and the command requested (log) %.3e.\n" , qentrd , rfield.totpow[i]); fprintf( ioQQQ, " The difference in the log is %.3e.\n" , diff ); if( diff>0. ) { fprintf( ioQQQ, " The continuum source is too bright.\n" ); } else { fprintf( ioQQQ, " The continuum source is too faint.\n" ); } /* explain what happened */ fprintf( ioQQQ, " The usual cause for this problem is an incorrect continuum intensity/luminosity or radius command.\n" ); fprintf( ioQQQ, " There were a total of %li continuum shape commands entered - the problem is with number %li.\n", rfield.nShape , i+1 ); cdEXIT(EXIT_FAILURE); } else { rfield.spfac[i] = pow(10.,diff); } if( trace.lgTrace ) { fprintf( ioQQQ, " conorm finds Q over range from%11.4e-%11.4e Ryd, integral= %10.4e Factor=%11.4e\n", rfield.range[i][0], rfield.range[i][1], qentrd , rfield.spfac[i] ); } } else if( strcmp(rfield.chSpNorm[i],"LUMI") == 0 ) { /* luminosity entered, special since default is TOTAL lumin */ /*pintr integrates L for any continuum between two limits, used for normalization, * return units are log of ryd cm-2 s-1, last log conv to ergs */ pentrd = pintr(rfield.range[i][0],rfield.range[i][1]) + log10(EN1RYD); f = rfield.totpow[i] - pentrd; rfield.spfac[i] = pow(10.,f); if( trace.lgTrace ) { fprintf( ioQQQ, " conorm finds luminosity range is %10.3e to %9.3e Ryd, factor is %11.4e\n", rfield.range[i][0], rfield.range[i][1], rfield.spfac[i] ); } } else { fprintf( ioQQQ, "PROBLEM DISASTER What chSpNorm label is this? =%s=\n", rfield.chSpNorm[i]); TotalInsanity(); } /* spfac used to renormalize SED into flux */ if( fabs(rfield.spfac[i]) <=SMALLDOUBLE ) { fprintf( ioQQQ, "PROBLEM DISASTER conorm finds infinite continuum scale factor.\n" ); fprintf( ioQQQ, "The continuum is too intense to compute with this cpu.\n" ); fprintf( ioQQQ, "Were the intensity and luminosity commands switched?\n" ); fprintf( ioQQQ, "Sorry, but I cannot go on.\n" ); cdEXIT(EXIT_FAILURE); } } /* this is conversion factor for final units of line intensities or luminosities in printout, * will be intensities (==0) unless luminosity is to be printed, or flux at Earth * pirsq is the log of 4 pi r_in^2 */ radius.Conv2PrtInten = radius.pirsq; /* >>chng 02 apr 25, add option for slit on aperture command */ if( geometry.iEmissPower == 1 ) { if( radius.lgPredLumin ) { /* factor should be divided by 2 r_in (so that 2pi*r_in remains) */ radius.Conv2PrtInten -= (log10(2.) + xLog_radius_inner); } else if( !radius.lgPredLumin ) { /* this is an error - slit requested but radius is not known */ fprintf( ioQQQ, "PROBLEM DISASTER conorm: Aperture slit specified, but not predicting luminosity.\n" ); fprintf( ioQQQ, "conorm: Please specify an inner radius to determine L.\nSorry\n" ); cdEXIT(EXIT_FAILURE); } } if( geometry.iEmissPower == 0 && radius.lgPredLumin ) { /* leave Conv2PrtInten at zero if not predicting luminosity */ radius.Conv2PrtInten = log10(2.); } /* this is option to give final absolute results as flux observed at Earth */ if( radius.distance > 0. && radius.lgRadiusKnown && prt.lgPrintFluxEarth ) { /* this implements the conversion from Q_alpha -> F_alpha * described in the section on the APERTURE command in Hazy 1 */ if( geometry.iEmissPower == 0 ) radius.Conv2PrtInten += log10(double(geometry.size)/SQAS_SKY); else if( geometry.iEmissPower == 1 ) radius.Conv2PrtInten += log10(double(geometry.size)/(PI4*AS1RAD*radius.distance)); else if( geometry.iEmissPower == 2 ) radius.Conv2PrtInten -= log10( PI4*pow2(radius.distance) ); else TotalInsanity(); } /* normally lines are into 4pi, this is option to do per sr or arcsec^2 */ if( prt.lgSurfaceBrightness ) { if( radius.pirsq != 0. ) { /* make sure we are predicting line intensities, not luminosity */ fprintf( ioQQQ, " PROBLEM DISASTER Sorry, but both luminosity and surface brightness have been requested for lines.\n" ); fprintf( ioQQQ, " the PRINT LINE SURFACE BRIGHTNESS command can only be used when lines are predicted per unit cloud area.\n" ); cdEXIT(EXIT_FAILURE); } if( prt.lgSurfaceBrightness_SR ) { /* we want final units to be per sr */ radius.Conv2PrtInten -= log10( PI4 ); } else { /* we want final units to be per square arcsec */ radius.Conv2PrtInten -= log10(SQAS_SKY); } } return; } /*qintr integrates Q for any continuum between two limits, used for normalization */ STATIC double qintr(double *qenlo, double *qenhi) { long int i, ipHi, ipLo, j; double qintr_v, sum, wanu; DEBUG_ENTRY( "qintr()" ); /* returns LOG of number of photons over energy interval */ /* this is copy of logic that occurs three times across code */ ASSERT(*qenhi > *qenlo); ipLo = ipoint(*qenlo); ipHi = ipoint(*qenhi); /* this is actual sum of photons within band */ sum = 0.; for( i=ipLo-1; i < (ipHi - 1); i++ ) { /*sum += ffun1(rfield.anu[i])*rfield.widflx[i];*/ for( j=0; j < 4; j++ ) { wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j]; /* >>chng 02 jul 16, add test on continuum limits - * this was exceeded when resolution set very large */ wanu = MAX2( wanu , rfield.emm ); wanu = MIN2( wanu , rfield.egamry ); sum += fweigh[j]*ffun1(wanu)*rfield.widflx[i]; } } if( sum <= 0. ) { fprintf( ioQQQ, " PROBLEM DISASTER Photon number sum in QINTR is %.3e\n", sum ); fprintf( ioQQQ, " This source has no ionizing radiation, and the number of ionizing photons was specified.\n" ); fprintf( ioQQQ, " This was continuum source number%3ld\n", rfield.ipSpec ); fprintf( ioQQQ, " Sorry, but I cannot go on. ANU and FLUX arrays follow. Enjoy.\n" ); fprintf( ioQQQ, "\n\n This error is also caused by an old table read file whose energy mesh does not agree with the code.\n" ); for( i=0; i < rfield.nupper; i++ ) { fprintf( ioQQQ, "%.2e\t%.2e\n", rfield.anu[i], rfield.flux[0][i] ); } cdEXIT(EXIT_FAILURE); } else { qintr_v = log10(sum); } return qintr_v; } /*pintr integrates L for any continuum between two limits, used for normalization, * return units are log of ryd cm-2 s-1 */ STATIC double pintr(double penlo, double penhi) { long int i, j; double fsum, pintr_v, sum, wanu, wfun; long int ip1 , ip2; DEBUG_ENTRY( "pintr()" ); /* computes log of luminosity in radiation over some integral * answer is in Ryd per sec */ sum = 0.; /* >>chng 02 oct 27, do not call qg32, do same type sum as * final integration */ /* laser is special since delta function, this is center of laser */ /* >>chng 01 jul 01, was +-21 cells, change to call to ipoint */ ip1 = ipoint( penlo ); ip2 = ipoint( penhi ); for( i=ip1-1; i < ip2-1; i++ ) { fsum = 0.; for( j=0; j < 4; j++ ) { wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j]; /*++iiii;fprintf(ioQQQ,"DEBUG iii %li %e \n",iiii, wanu );*/ wfun = fweigh[j]*ffun1(wanu)*wanu; fsum += wfun; } sum += fsum*rfield.widflx[i]; } if( sum > 0. ) { pintr_v = log10(sum); } else { pintr_v = -38.; } return pintr_v; }