cloudy trunk
|
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 /*ParseBackgrd parse options for the BACKGROUND command - this actually enters two continua*/ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "radius.h" 00007 #include "rfield.h" 00008 #include "parse.h" 00009 00010 void ParseBackgrd(long int *nqh, 00011 char *chCard, 00012 realnum *ar1) 00013 { 00014 bool lgEOL; 00015 long int i; 00016 double a, 00017 fac, 00018 rlogl, 00019 z; 00020 char chLocal[INPUT_LINE_LENGTH]; 00021 00022 DEBUG_ENTRY( "ParseBackgrd()" ); 00023 00024 /* check that stack of shape and luminosity specifications 00025 * is parallel, stop if not - this happens is background comes 00026 * BETWEEN another set of shape and luminosity commands */ 00027 if( rfield.nspec != *nqh ) 00028 { 00029 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" ); 00030 fprintf( ioQQQ, " Sorry.\n" ); 00031 cdEXIT(EXIT_FAILURE); 00032 } 00033 00034 /* diffuse x-ray background from 00035 * >>refer cosmic background Ostriker and Ikeuchi ApJL 268, L63. 00036 * parameter on cmnd is redshift, followed by optional J21 00037 * >>refer cosmic background Ikeuchi, S.; Ostriker, J. P., 1986, ApJ 301, 522 */ 00038 00039 /* >>chng 02 aug 12, do not try to enter bare powerlaw, use table power law default instead */ 00040 strcpy( chLocal, "TABLE POWERLAW " ); 00041 ParseTable( nqh, chLocal , ar1); 00042 00043 /* now generate equivalent of luminosity command 00044 * enter phi(h), the number of h-ionizing photons/cm2 */ 00045 00046 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00047 strcpy( rfield.chSpNorm[*nqh], "FLUX" ); 00048 /* this is an isotropic radiation field */ 00049 rfield.lgBeamed[*nqh] = false; 00050 00051 i = 5; 00052 /* this will be the redshift */ 00053 z = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00054 if( lgEOL ) 00055 { 00056 z = 0.; 00057 } 00058 00059 /* optional scale factor */ 00060 fac = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00061 if( lgEOL ) 00062 { 00063 fac = 1.0; 00064 } 00065 00066 /* this equation from Ostriker and Ikeuchi Ap.J.L 268, L63. 00067 * this should be J_nu into 4\pi sr 00068 * >>chng 96 jul 23, from ostriker to vedel, et al mnras 271, 743. */ 00069 rfield.totpow[*nqh] = (log10(PI4*fac*1.e-21/ 00070 (1.+powi(5./(1.+z),4)))); 00071 00072 /* this is at 1 ryd for H 00073 * range(nqh,1) = 1. 00074 *>>chgn 96 dec 18, to H ioniz pot from 1 ryd */ 00075 rfield.range[*nqh][0] = HIONPOT; 00076 00077 /* >>chng 02 aug 12, call table power law routine rather than try to enter spec here, 00078 * so no longer increment nspec */ 00079 /* increment continuum indices 00080 ++rfield.nspec; 00081 if( rfield.nspec >= LIMSPC ) 00082 { 00083 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00084 cdEXIT(EXIT_FAILURE); 00085 } */ 00086 ++*nqh; 00087 if( *nqh >= LIMSPC ) 00088 { 00089 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00090 cdEXIT(EXIT_FAILURE); 00091 } 00092 00093 /* add fireball unless NO FIREBALL or NO CMP in */ 00094 if( !(nMatch("O FI",chCard) || nMatch("O CM",chCard) ) ) 00095 { 00096 00097 /* put in a black body */ 00098 strcpy( rfield.chSpType[rfield.nspec], "BLACK" ); 00099 /* >>chng 03 may 23, CMB temp from 2.756 to 2.725 */ 00100 rfield.slope[rfield.nspec] = (CMB_TEMP*(1. + z)); 00101 rfield.cutoff[rfield.nspec][0] = 0.; 00102 rfield.cutoff[rfield.nspec][1] = 0.; 00103 strcpy( rfield.chSpNorm[*nqh], "LUMI" ); 00104 a = log10(rfield.slope[rfield.nspec]); 00105 rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a; 00106 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00107 rfield.range[*nqh][0] = rfield.emm; 00108 rfield.range[*nqh][1] = rfield.egamry; 00109 rfield.totpow[*nqh] = rlogl; 00110 00111 /* this flag says that CMB has been set */ 00112 rfield.lgCMB_set = true; 00113 00114 ++rfield.nspec; 00115 ++*nqh; 00116 if( *nqh >= LIMSPC ) 00117 { 00118 fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" ); 00119 cdEXIT(EXIT_FAILURE); 00120 } 00121 } 00122 00123 /* set radius to very large value if not already set */ 00124 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00125 if( !radius.lgRadiusKnown ) 00126 { 00127 *ar1 = (realnum)radius.rdfalt; 00128 radius.Radius = pow(10.,radius.rdfalt); 00129 } 00130 return; 00131 }