/*conorm normalize continuum to proper intensity */ #include #include "cddefines.h" #include "physconst.h" #include "incidspec.h" #include "phycon.h" #include "radius.h" #include "trace.h" #include "sphere.h" #include "fourpi.h" #include "bit32.h" #include "ffun1.h" #include "qintr.h" #include "pintr.h" #include "rfield.h" #include "conorm.h" void conorm(void) { long int i; double ar1, diff, f, flx1, flx2, pentrd, qentrd; # ifdef DEBUG_FUN fputs( "<+>conorm()\n", debug_fp ); # endif ar1 = log10(radius.Radius); /* this is a sanity check, it can't happen */ for( i=0; i < IncidSpec.nspec; i++ ) { if( strcmp(spnorm.chRSpec[i],"UNKN") == 0 ) { fprintf( ioQQQ, " UNKN spectral normalization cannot happen.\n" ); fprintf( ioQQQ, " CONORM punts.\n" ); puts( "[Stop in conorm]" ); exit(1); } else if( strcmp(spnorm.chRSpec[i],"SQCM") != 0 && strcmp(spnorm.chRSpec[i],"4 PI") != 0 ) { fprintf( ioQQQ, " chRSpec must be SQCM or 4 PI, and it was %4.4s. This cannot happen.\n", spnorm.chRSpec[i] ); fprintf( ioQQQ, " CONORM punts.\n" ); puts( "[Stop in conorm]" ); exit(1); } /* this sanity check makes sure that atlas.mod or werner.mod grids * are for the current version of the code */ if( strcmp(spnorm.chSpType[i],"VOLK ") == 0 ) { /* check that wavelength scale is actually defined outside here */ assert( rfield.AnuOrg[NCELL-1]>0. ); diff = fabs(IncidSpec.tnu[NCELL-1][i]-rfield.AnuOrg[NCELL-1])/ rfield.AnuOrg[NCELL-1]; if( diff > 0.001 ) { fprintf( ioQQQ, "%10.2e%10.2e\n", rfield.AnuOrg[NCELL-1], IncidSpec.tnu[NCELL-1][i] ); fprintf( ioQQQ, " CONORM: The energy grid of the stellar atmosphere file does not agree with the grid in this version of the code.\n" ); fprintf( ioQQQ, " A stellar atmosphere grid from an old version of the code is probably in place.\n" ); fprintf( ioQQQ, " A grid for the current version of Cloudy must be generated and used.\n" ); fprintf( ioQQQ, " This is done with the COMPILE STARS command.\n" ); fprintf( ioQQQ, " Sorry.\n" ); puts( "[Stop in conorm]" ); exit(1); } } } /* 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 */ fourpi.pirsq = 0.; fourpi.lgFourPi = FALSE; /* check whether ANY luminosities are present */ for( i=0; i < IncidSpec.nspec; i++ ) { if( strcmp(spnorm.chRSpec[i],"4 PI") == 0 ) { fourpi.pirsq = (float)(1.0992099 + 2.*ar1); fourpi.lgFourPi = TRUE; /* convert down to intensity */ IncidSpec.totpow[i] -= fourpi.pirsq; if( trace.lgTrace ) { fprintf( ioQQQ, " CONORM converts continuum%3ld from luminosity to intensity.\n", i ); } } } /* convert any ionization parameters to q(h) */ for( i=0; i < IncidSpec.nspec; i++ ) { if( strcmp(spnorm.chSpNorm[i],"IONI") == 0 ) { IncidSpec.totpow[i] += (float)(log10(phycon.hden) + log10(SPEEDLIGHT)); strcpy( spnorm.chSpNorm[i], "Q(H)" ); if( trace.lgTrace ) { fprintf( ioQQQ, " CONORM converts continuum%3ld from ionizat par to q(h).\n", i ); } } } /* indicate whether we ended up with luminosity or intensity */ if( trace.lgTrace ) { if( fourpi.lgFourPi ) { 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( !fourpi.lgFourPi ) { sphere.covgeo = 1.; } /* main loop to find continuum normalization starts here */ for( i=0; i < IncidSpec.nspec; i++ ) { IncidSpec.ipspec = i; if( trace.lgTrace ) { fprintf( ioQQQ, " CONORM continuum number%3ld is shape %5.5s range is %11.2e%11.2e\n", i, spnorm.chSpType[i], IncidSpec.range[0][i], IncidSpec.range[1][i] ); } if( strcmp(spnorm.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" ); puts( "[Stop in conorm]" ); exit(1); } /* first find photon flux and Q of prevous continuum */ IncidSpec.ipspec -= 1; flx1 = ffun1(IncidSpec.range[0][i])*IncidSpec.spfac[IncidSpec.ipspec]* IncidSpec.range[0][i]; /* 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" ); puts( "[Stop in conorm]" ); exit(1); } /* return pointer to previous (correct) value, find F, Q */ IncidSpec.ipspec += 1; /* we want a continuum totpow as powerful, flx is now desired flx */ flx1 *= IncidSpec.totpow[i]; /*. find flux of this new continuum at that point */ flx2 = ffun1(IncidSpec.range[1][i])*IncidSpec.range[1][i]; /* this is ratio of desired to actual */ IncidSpec.spfac[i] = (float)(flx1/flx2); if( trace.lgTrace ) { fprintf( ioQQQ, " CONORM ratio will set scale fac to%10.3e at%10.2e Ryd.\n", IncidSpec.totpow[i], IncidSpec.range[0][i] ); } } else if( strcmp(spnorm.chSpNorm[i],"FLUX") == 0 ) { /* specify flux density * option to use arbitrary frequency or range */ f = ffun1(IncidSpec.range[0][i]); f = MAX2(1e-37,f); f = log10(f) + log10(IncidSpec.range[0][i]*EN1RYD/FR1RYD); f = IncidSpec.totpow[i] - f; /* >>chng 96 dec 31, added following test */ if( bit32.lgBit32 ) { if( f > 35. ) { fprintf( ioQQQ, " Continuum source%3ld is too intense for this cpu - is it normalized correctly?\n", i ); puts( "[Stop in conorm]" ); exit(1); } } IncidSpec.spfac[i] = (float)pow(10.,f); if( trace.lgTrace ) { fprintf( ioQQQ, " CONORM will set log fnu to%10.3e at%10.2e Ryd. Factor=%11.4e\n", IncidSpec.totpow[i], IncidSpec.range[0][i], IncidSpec.spfac[i] ); } } else if( strcmp(spnorm.chSpNorm[i],"Q(H)") == 0 || strcmp(spnorm.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", IncidSpec.range[0][i], IncidSpec.range[1][i] , IncidSpec.totpow[i]); } qentrd = qintr(&IncidSpec.range[0][i],&IncidSpec.range[1][i]); diff = IncidSpec.totpow[i] - qentrd; if( bit32.lgBit32 && (diff < -25. || diff > 35.) ) { fprintf( ioQQQ, " Continuum source specified is too extreme for a 32 bit cps.\n" ); puts( "[Stop in conorm]" ); exit(1); } else { IncidSpec.spfac[i] = (float)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", IncidSpec.range[0][i], IncidSpec.range[1][i], qentrd , IncidSpec.spfac[i] ); } } else if( strcmp(spnorm.chSpNorm[i],"LUMI") == 0 ) { /* luminosity entered, special since default is TOTAL lumin */ pentrd = pintr(IncidSpec.range[0][i],IncidSpec.range[1][i]) + log10(EN1RYD); f = IncidSpec.totpow[i] - pentrd; /* >>chng 96 dec 31, added following test */ if( bit32.lgBit32 ) { if( f > 35. ) { fprintf( ioQQQ, " Continuum source%3ld is too intense for this cpu - is it normalized correctly?\n", i ); puts( "[Stop in conorm]" ); exit(1); } } IncidSpec.spfac[i] = (float)pow(10.,f); if( trace.lgTrace ) { fprintf( ioQQQ, " CONORM finds luminosity range is%10.3e to %9.3e Ryd, factor is%11.4e\n", IncidSpec.range[0][i], IncidSpec.range[1][i], IncidSpec.spfac[i] ); } } else { fprintf( ioQQQ, " What chSpNorm label is this? =%s=\n", spnorm.chSpNorm[i]); puts( "[Stop in conorm]" ); exit(1); } /* sec part after .or. added June 93 because sometimes spfac=0 * got past first test */ if( 1./IncidSpec.spfac[i] == 0. || IncidSpec.spfac[i] == 0. ) { fprintf( ioQQQ, " 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 commandsswitched?\n" ); fprintf( ioQQQ, " Cloudy punts. Sorry.\n" ); puts( "[Stop in conorm]" ); exit(1); } } # ifdef DEBUG_FUN fputs( " <->conorm()\n", debug_fp ); # endif return; }