/*ParseBlackbody parse parameters off black body command */ #include "cddefines.h" #include #include "physconst.h" #include "incidspec.h" #include "bounds.h" #include "radius.h" #include "ffmtread.h" #include "lgmatch.h" #include "parse.h" void ParseBlackbody(char *chCard, /* input command line, already changed to caps */ long int *nqh, /* counter for which continuum source this is */ float *ar1) /* optional area that might be set here */ { int lgEOL; long int i; double a, dil, rlogl; # ifdef DEBUG_FUN fputs( "<+>ParseBlackbody()\n", debug_fp ); # endif /* increment number of continus */ IncidSpec.nspec = MIN2(IncidSpec.nspec+1,LIMSPC); /* type is blackbody */ strcpy( spnorm.chSpType[IncidSpec.nspec-1], "BLACK" ); /* these two are not used for this continuum shape */ IncidSpec.cutoff[0][IncidSpec.nspec-1] = 0.; IncidSpec.cutoff[1][IncidSpec.nspec-1] = 0.; /* get the temperature */ i = 5; IncidSpec.slope[IncidSpec.nspec-1] = (float)FFmtRead(chCard,&i,76,&lgEOL); /* there must be at least one number, the temperature */ if( lgEOL ) { fprintf( ioQQQ, " There must be at least 1 number on a BLACKBODY command line. Sorry.\n" ); puts( "[Stop in ParseBlackbody]" ); exit(1); } /* this is the temperature - make sure its linear in the end * there are two keys, LINEAR and LOG, that could be here, * else choose which is here by which side of 10 */ if( (IncidSpec.slope[IncidSpec.nspec-1] <= 10. && !lgMatch("LINE",chCard)) || lgMatch(" LOG",chCard) ) { /* log option */ IncidSpec.slope[IncidSpec.nspec-1] = (float)pow(10.,IncidSpec.slope[IncidSpec.nspec-1]); } /* check that temp is not too low - could happen if log misused */ if( IncidSpec.slope[IncidSpec.nspec-1] < 1e4 ) { fprintf( ioQQQ, " Is T(star)=%10.2e correct???\n", IncidSpec.slope[IncidSpec.nspec-1] ); } /* now check that temp not too low - would peak below low * energy limit of the code * factor is temperature of 1 Ryd, egamry is high energy limit of code */ if( IncidSpec.slope[IncidSpec.nspec-1]/TE1RYD < bounds.emm ) { fprintf( ioQQQ, " This temperature is very low - the blackbody will have significant flux low the low energy limit of the code, presently %10.2e Ryd.\n", bounds.emm ); fprintf( ioQQQ, " Was this intended?\n" ); } /* now check that temp not too high - would extend beyond high * energy limit of the code * factor is temperature of 1 Ryd, egamry is high energy limit of code */ if( IncidSpec.slope[IncidSpec.nspec-1]/TE1RYD*2. > bounds.egamry ) { fprintf( ioQQQ, " This temperature is very high - the blackbody will have significant flux above the high energy limit of the code,%10.2e Ryd.\n", bounds.egamry ); fprintf( ioQQQ, " Was this intended?\n" ); } /* also possible to input log(total luminosity)=real log(l) */ a = FFmtRead(chCard,&i,76,&lgEOL); if( lgEOL ) { /* there was not a second number on the line; check if LTE */ if( lgMatch(" LTE",chCard) ) { /* set energy density to LTE value */ *nqh = MIN2(*nqh+1,LIMSPC); /* check that stack of shape and luminosity specifications * is parallel, stop if not - this happens is background comes * BETWEEN another set of shape and luminosity commands */ if( IncidSpec.nspec != *nqh ) { fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete eachcontinuum specification before starting another.\n" ); fprintf( ioQQQ, " Sorry.\n" ); puts( "[Stop in ParseBlackbody]" ); exit(1); } strcpy( spnorm.chSpNorm[*nqh-1], "LUMI" ); a = log10(IncidSpec.slope[IncidSpec.nspec-1]); rlogl = log10(2.99792e10*7.56464e-15) + 4.*a; if( radius.Radius == 0. ) { *ar1 = (float)radius.rdfalt; radius.Radius = pow(10.,*ar1); } strcpy( spnorm.chRSpec[*nqh-1], "SQCM" ); IncidSpec.range[0][*nqh-1] = bounds.emm; IncidSpec.range[1][*nqh-1] = bounds.egamry; IncidSpec.totpow[*nqh-1] = (float)rlogl; } # ifdef DEBUG_FUN fputs( " <->ParseBlackbody()\n", debug_fp ); # endif return; } /* the second number will be some sort of luminosity */ *nqh = MIN2(*nqh+1,LIMSPC); /* check that stack of shape and luminosity specifications * is parallel, stop if not - this happens is background comes * BETWEEN another set of shape and luminosity commands */ if( IncidSpec.nspec != *nqh ) { fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete eachcontinuum specification before starting another.\n" ); fprintf( ioQQQ, " Sorry.\n" ); puts( "[Stop in ParseBlackbody]" ); exit(1); } strcpy( spnorm.chSpNorm[*nqh-1], "LUMI" ); /* a second number was entered, what was it? */ if( lgMatch("LUMI",chCard) ) { rlogl = a; strcpy( spnorm.chRSpec[*nqh-1], "4 PI" ); } else if( lgMatch("RADI",chCard) ) { /* radius was entered, convert to total lumin */ rlogl = -3.147238 + 2.*a + 4.*log10(IncidSpec.slope[IncidSpec.nspec-1]); strcpy( spnorm.chRSpec[*nqh-1], "4 PI" ); } else if( lgMatch("DENS",chCard) ) { /* number was temperature to deduce energy density * number is linear if greater than 10, or if LINEAR appears on line * want number to be log of temperature at end of this */ if( lgMatch("LINE",chCard) || a > 10. ) { a = log10(a); } rlogl = log10(2.99792e10*7.56464e-15) + 4.*a; if( radius.Radius == 0. ) { *ar1 = (float)radius.rdfalt; radius.Radius = pow(10.,*ar1); } strcpy( spnorm.chRSpec[*nqh-1], "SQCM" ); } else if( lgMatch("DILU",chCard) ) { /* number is dilution factor, if negative then its log */ if( a < 0. ) { dil = a; } else { dil = log10(a); } if( dil > 0. ) { fprintf( ioQQQ, " Is a dilution factor greater than one physical?\n" ); } /* this is LTE energy density */ a = log10(IncidSpec.slope[IncidSpec.nspec-1]); rlogl = log10(2.99792e10*7.56464e-15) + 4.*a; /* add on dilution factor */ rlogl += dil; if( radius.Radius == 0. ) { *ar1 = (float)radius.rdfalt; radius.Radius = pow(10.,*ar1); } strcpy( spnorm.chRSpec[*nqh-1], "SQCM" ); } else { fprintf( ioQQQ, " I dont understand what the second number was- keywords are LUMINOSITY, RADIUS, DILUTION, and DENSITY.\n" ); puts( "[Stop in ParseBlackbody]" ); exit(1); } IncidSpec.range[0][*nqh-1] = bounds.emm; IncidSpec.range[1][*nqh-1] = bounds.egamry; IncidSpec.totpow[*nqh-1] = (float)rlogl; # ifdef DEBUG_FUN fputs( " <->ParseBlackbody()\n", debug_fp ); # endif return; }