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 /*ffun evaluate total flux for sum of all continuum sources */ 00004 /*ffun1 derive flux at a specific energy, for one continuum */ 00005 /*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */ 00006 #include "cddefines.h" 00007 #include "physconst.h" 00008 #include "rfield.h" 00009 #include "ipoint.h" 00010 #include "opacity.h" 00011 #include "continuum.h" 00012 00013 /*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */ 00014 STATIC void ReadTable(const string& fnam); 00015 00016 /* evaluate sum of all individual continua at one energy, return is 00017 * continuum intensity */ 00018 double ffun( 00019 /* the energy in Rydbergs where the continuum will be evaluated */ 00020 double anu ) 00021 { 00022 double frac_beam_time; 00023 /* fraction of beamed continuum that is constant */ 00024 double frac_beam_const; 00025 /* fraction of continuum that is isotropic */ 00026 double frac_isotropic; 00027 double a; 00028 00029 DEBUG_ENTRY( "ffun()" ); 00030 00031 /* call real function */ 00032 a = ffun( anu , &frac_beam_time , &frac_beam_const , &frac_isotropic ); 00033 return a; 00034 } 00035 00036 /* evaluate sum of all individual continua at one energy, return is 00037 * continuum intensity */ 00038 double ffun( 00039 /* the energy in Rydbergs where the continuum will be evaluated */ 00040 double anu , 00041 /* fraction of beamed continuum that is varies with time */ 00042 double *frac_beam_time, 00043 /* fraction of beamed continuum that is constant */ 00044 double *frac_beam_const, 00045 /* fraction of continuum that is isotropic */ 00046 double *frac_isotropic ) 00047 { 00048 double ffun_v; 00049 static bool lgWarn = false; 00050 double flx_beam_time , flx_beam_const , flx_isotropic; 00051 00052 DEBUG_ENTRY( "ffun()" ); 00053 00054 /* This routine, ffun, returns the sum of photons per unit time, area, energy, 00055 * for all continua in the calculation. ffun1 is called and returns 00056 * a single continuum. We loop over all nspec continuum sources - 00057 * ipspec points to each */ 00058 ffun_v = 0.; 00059 flx_beam_time = 0.; 00060 flx_beam_const = 0.; 00061 flx_isotropic = 0.; 00062 for( rfield.ipspec=0; rfield.ipspec < rfield.nspec; rfield.ipspec++ ) 00063 { 00064 double one = ffun1(anu)*rfield.spfac[rfield.ipspec]; 00065 ffun_v += one; 00066 00067 /* find fraction of total that is constant vs variable and 00068 * isotropic vs beamed */ 00069 if( rfield.lgBeamed[rfield.ipspec] ) 00070 { 00071 if( rfield.lgTimeVary[rfield.ipspec] ) 00072 flx_beam_time += one; 00073 else 00074 flx_beam_const += one; 00075 } 00076 else 00077 flx_isotropic += one; 00078 } 00079 00080 /* at this point rfield.flux is the sum of the following three continua 00081 * now keep track of the three different types */ 00082 if( ffun_v < SMALLFLOAT ) 00083 { 00084 *frac_beam_time = 0.; 00085 *frac_beam_const = 1.; 00086 *frac_isotropic = 0.; 00087 } 00088 else 00089 { 00090 /* fraction of beamed continuum that varies with time */ 00091 *frac_beam_time = flx_beam_time / ffun_v; 00092 /* part of beamed continuum that is constant */ 00093 *frac_beam_const = flx_beam_const / ffun_v; 00094 /* the constant isotropic continuum */ 00095 *frac_isotropic = flx_isotropic / ffun_v; 00096 } 00097 ASSERT( *frac_beam_time >=0. && *frac_beam_time<=1.+3.*DBL_EPSILON ); 00098 ASSERT( *frac_beam_const >=0.&& *frac_beam_const<=1.+3.*DBL_EPSILON ); 00099 ASSERT( *frac_isotropic >=0. && *frac_isotropic<=1.+3.*DBL_EPSILON ); 00100 ASSERT( fabs( 1.-*frac_beam_time-*frac_beam_const-*frac_isotropic)< 00101 10.*DBL_EPSILON); 00102 00103 if( ffun_v > BIGFLOAT && !lgWarn ) 00104 { 00105 lgWarn = true; 00106 fprintf( ioQQQ, " FFUN: The net continuum is very intense.\n" ); 00107 fprintf( ioQQQ, " I will try to press on, but may have problems.\n" ); 00108 } 00109 return( ffun_v ); 00110 } 00111 00112 /*ffun1 derive flux at a specific energy, for one continuum */ 00113 double ffun1(double xnu) 00114 { 00115 char chKey[6]; 00116 long int i; 00117 double fac, 00118 ffun1_v, 00119 y; 00120 static bool lgWarn = false; 00121 00122 DEBUG_ENTRY( "ffun1()" ); 00123 00124 00125 /* confirm that pointer is within range */ 00126 ASSERT( rfield.ipspec >= 0); 00127 ASSERT( rfield.ipspec < rfield.nspec ); 00128 00129 /* FFUN1 returns photons per unit time, area, energy, for one continuum 00130 * ipspec is the pointer to the continuum source, in the order 00131 * entered on the command lines */ 00132 00133 /*begin sanity check */ 00134 ASSERT( xnu >= rfield.emm*0.99 ); 00135 ASSERT( xnu <= rfield.egamry*1.01 ); 00136 /*end sanity check */ 00137 00138 strcpy( chKey, rfield.chSpType[rfield.ipspec] ); 00139 00140 if( strcmp(chKey,"AGN ") == 0 ) 00141 { 00142 /* power law with cutoff at both ends 00143 * nomenclature all screwed up - slope is cutoff energy in Ryd, 00144 * cutoff[1][i] is ratio of two continua from alpha ox 00145 * cutoff[2][i] is slope of bb temp */ 00146 ffun1_v = pow(xnu,-1. + rfield.cutoff[rfield.ipspec][1])* 00147 sexp(xnu/rfield.slope[rfield.ipspec])*sexp(0.01/ 00148 xnu); 00149 /* only add on x-ray for energies above 0.1 Ryd */ 00150 if( xnu > 0.1 ) 00151 { 00152 if( xnu < 7350. ) 00153 { 00154 /* cutoff is actually normalization constant 00155 * below 100keV continuum is nu-1 */ 00156 ffun1_v += pow(xnu,rfield.cutoff[rfield.ipspec][2] - 00157 1.)*rfield.cutoff[rfield.ipspec][0]*sexp(1./ 00158 xnu); 00159 } 00160 else 00161 { 00162 ffun1_v += pow(7350.,rfield.cutoff[rfield.ipspec][2] - 00163 1.)*rfield.cutoff[rfield.ipspec][0]/ 00164 POW3(xnu/7350.); 00165 } 00166 } 00167 00168 } 00169 else if( strcmp(chKey,"POWER") == 0 ) 00170 { 00171 /* power law with cutoff at both ends */ 00172 ffun1_v = pow(xnu,-1. + rfield.slope[rfield.ipspec])* 00173 sexp(xnu/rfield.cutoff[rfield.ipspec][0]+rfield.cutoff[rfield.ipspec][1]/ 00174 xnu); 00175 00176 } 00177 else if( strcmp(chKey,"BLACK") == 0 ) 00178 { 00179 const double db_log = log(DBL_MAX); 00180 00181 /* black body */ 00182 fac = TE1RYD*xnu/rfield.slope[rfield.ipspec]; 00183 /* >>>chng 00 apr 13 from 80 to log(dbl_max) */ 00184 if( fac > db_log ) 00185 { 00186 ffun1_v = 0.; 00187 } 00188 else if( fac > 1.e-5 ) 00189 { 00190 ffun1_v = xnu*xnu/(exp(fac) - 1.); 00191 } 00192 else 00193 { 00194 ffun1_v = xnu*xnu/(fac*(1. + fac/2.)); 00195 } 00196 00197 } 00198 else if( strcmp(chKey,"INTER") == 0 ) 00199 { 00200 /* interpolate on tabulated input spectrum, factor of 1.0001 to 00201 * make sure that requested energy is within bounds of array */ 00202 if( xnu >= rfield.tNuRyd[rfield.ipspec][0]*1.000001 ) 00203 { 00204 /* loop starts at second array element, [1], since want to 00205 * find next continuum energy greater than desired point */ 00206 i = 1; 00207 /* up to NCELL tabulated points may be read in. Very fine 00208 * continuum mesh such as that output by stellar atmospheres 00209 * can have very large number of points */ 00210 while( i< NCELL-1 && rfield.tNuRyd[rfield.ipspec][i]>0. ) 00211 { 00212 if( xnu < rfield.tNuRyd[rfield.ipspec][i] ) 00213 { 00214 /* the energy xnu is between points rfield.tNuRyd[rfield.ipspec][i-1] 00215 * and rfield.tNuRyd[rfield.ipspec][i] - do linear 00216 * interpolation in log log space */ 00217 y = rfield.tFluxLog[rfield.ipspec][i-1] + 00218 rfield.tslop[rfield.ipspec][i-1]* 00219 log10(xnu/rfield.tNuRyd[rfield.ipspec][i-1]); 00220 00221 /* return value is photon density, div by energy */ 00222 ffun1_v = pow(10.,y); 00223 00224 /* this checks that overshoots did not occur - interpolated 00225 * value must be between lowest and highest point */ 00226 # ifndef NDEBUG 00227 double ys1 = MIN2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]); 00228 double ys2 = MAX2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]); 00229 ys1 = pow( 10. , ys1 ); 00230 ys2 = pow( 10. , ys2 ); 00231 ASSERT( ffun1_v >= ys1/(1.+100.*FLT_EPSILON) ); 00232 ASSERT( ffun1_v <= ys2*(1.+100.*FLT_EPSILON) ); 00233 # endif 00234 /* return value is photon density, div by energy */ 00235 return( ffun1_v/xnu ); 00236 } 00237 ++i; 00238 } 00239 /* energy above highest in table */ 00240 ffun1_v = 0.; 00241 } 00242 else 00243 { 00244 /* energy below lowest on table */ 00245 ffun1_v = 0.; 00246 } 00247 } 00248 00249 else if( strcmp(chKey,"BREMS") == 0 ) 00250 { 00251 /* brems continuum, rough gaunt factor */ 00252 fac = TE1RYD*xnu/rfield.slope[rfield.ipspec]; 00253 ffun1_v = sexp(fac)/pow(xnu,1.2); 00254 00255 } 00256 else if( strcmp(chKey,"LASER") == 0 ) 00257 { 00258 const double BIG = 1.e10; 00259 const double SMALL = 1.e-25; 00260 /* a laser, mostly one frequency */ 00261 /* >>chng 01 jul 01, was hard-wired 0.05 rel frac, change to optional 00262 * second parameter, with default of 0.05 */ 00263 /*if( xnu > 0.95*rfield.slope[rfield.ipspec] && xnu < 00264 1.05*rfield.slope[rfield.ipspec] )*/ 00265 if( xnu > (1.-rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] && 00266 xnu < (1.+rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] ) 00267 { 00268 ffun1_v = BIG; 00269 } 00270 else 00271 { 00272 ffun1_v = SMALL; 00273 } 00274 00275 } 00276 else if( strcmp(chKey,"READ ") == 0 ) 00277 { 00278 /* implement the table read command 00279 * if this is first time called then we need to read in file */ 00280 if( rfield.tslop[rfield.ipspec][0] == 0. ) 00281 { 00282 /* >>chng 01 nov 01, array index for below was bogus, only worked for a single continuum */ 00283 ReadTable(rfield.ioTableRead[rfield.ipspec]); 00284 rfield.tslop[rfield.ipspec][0] = 1.; 00285 } 00286 /* use array of values read in on TABLE READ command */ 00287 if( xnu >= rfield.egamry ) 00288 { 00289 ffun1_v = 0.; 00290 } 00291 else 00292 { 00293 i = ipoint(xnu); 00294 if( i > rfield.nupper || i < 1 ) 00295 { 00296 ffun1_v = 0.; 00297 } 00298 else 00299 { 00300 ffun1_v = rfield.ConTabRead[i-1]/rfield.anu[i-1]/rfield.anu[i-1]; 00301 } 00302 } 00303 } 00304 00305 /* >>chng 06 jul 10, retired TABLE STARBURST command, PvH */ 00306 /* >>chng 05 nov 30, retired TABLE TLUSTY command, PvH */ 00307 00308 else if( strcmp(chKey,"VOLK ") == 0 ) 00309 { 00310 /* use array of values read in from Kevin Volk's rebinning of 00311 * large atlas grids */ 00312 if( xnu >= rfield.egamry ) 00313 { 00314 ffun1_v = 0.; 00315 } 00316 else 00317 { 00318 i = ipoint(xnu); 00319 if( i > rfield.nupper ) 00320 { 00321 fprintf( ioQQQ, " ffun1: Too many points - increase ncell\n" ); 00322 fprintf( ioQQQ, " cell needed=%4ld ncell=%4ld\n", 00323 i, rfield.nupper ); 00324 cdEXIT(EXIT_FAILURE); 00325 } 00326 if( i > rfield.nupper || i < 1 ) 00327 { 00328 ffun1_v = 0.; 00329 } 00330 else 00331 { 00332 /* bug fixed Jul 9 93: FFUN1 = TSLOP(IPSPEC,I) / ANU(I) / ANU(I) 00333 * i has value 939 */ 00334 ffun1_v = rfield.tslop[rfield.ipspec][i-1]/ rfield.anu[i-1]; 00335 } 00336 } 00337 } 00338 else 00339 { 00340 fprintf( ioQQQ, " ffun1: I do not understand continuum label \"%s\" for continuum %li.\n", 00341 chKey , rfield.ipspec); 00342 cdEXIT(EXIT_FAILURE); 00343 } 00344 00345 if( ffun1_v > 1e35 && !lgWarn ) 00346 { 00347 lgWarn = true; 00348 fprintf( ioQQQ, " FFUN1: Continuum %ld is very intense.\n", 00349 rfield.ipspec ); 00350 fprintf( ioQQQ, " I will try to press on, but may have problems.\n" ); 00351 } 00352 return( ffun1_v ); 00353 } 00354 00355 /*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */ 00356 STATIC void ReadTable(const string& fnam) 00357 { 00358 char chLine[INPUT_LINE_LENGTH]; 00359 long int i, 00360 nPoints; 00361 double Differ, 00362 EnerLast; 00363 FILE *io; 00364 00365 DEBUG_ENTRY( "ReadTable()" ); 00366 00367 io = open_data( fnam.c_str(), "r", AS_LOCAL_ONLY ); 00368 00369 /* read in first line of header */ 00370 if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL ) 00371 { 00372 fprintf( ioQQQ, " error 1 reading input continuum file.\n" ); 00373 cdEXIT(EXIT_FAILURE); 00374 } 00375 00376 /* now read in the file of numbers */ 00377 i = 0; 00378 /* keep reading until we hit eol or run out of room in the array */ 00379 while( (read_whole_line(chLine, (int)sizeof(chLine),io)!=NULL) && (i<rfield.nupper) ) 00380 { 00381 double help[2]; 00382 sscanf( chLine, "%lf%lf ", &help[0], &help[1] ); 00383 opac.tmn[i] = (realnum)help[0]; 00384 rfield.ConTabRead[i] = (realnum)help[1]; 00385 ++i; 00386 } 00387 /* put pointer at last good value */ 00388 nPoints = i - 1; 00389 00390 /* check that sane number of values entered */ 00391 if( nPoints < 1 ) 00392 { 00393 fprintf( ioQQQ, " ReadTable, file for TABLE READ has too few points, =%5ld\n", 00394 nPoints ); 00395 cdEXIT(EXIT_FAILURE); 00396 } 00397 00398 /* check on units of energy scale, convert to Rydbergs if necessary 00399 * tmn is scale read in from punch file, anu is correct scale 00400 * EnerLast is energy of last point in Rydbergs */ 00401 EnerLast = opac.tmn[nPoints]; 00402 if( fabs(opac.tmn[0]-rfield.anu[0])/rfield.anu[0] > 1e-3 ) 00403 { 00404 /* first energy does not appear to have been in Rydbergs */ 00405 if( opac.tmn[0] <= 0. ) 00406 { 00407 fprintf( ioQQQ, 00408 " DISASTER ReadTable: the first energy in table read file is not positive. Something is wrong.\n" ); 00409 cdEXIT(EXIT_FAILURE); 00410 } 00411 else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0]*1e-4)/opac.tmn[0] < 00412 1e-3 ) 00413 { 00414 /* wavelength in microns */ 00415 EnerLast = RYDLAM/opac.tmn[nPoints]/1e4; 00416 } 00417 else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0])/opac.tmn[0] < 00418 1e-3 ) 00419 { 00420 /* wavelength in Angstroms */ 00421 EnerLast = RYDLAM/opac.tmn[nPoints]; 00422 } 00423 else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD*1e-3)/opac.tmn[0] < 00424 1e-3 ) 00425 { 00426 /* wavelength in keV */ 00427 EnerLast = opac.tmn[nPoints]/EVRYD/1e-3; 00428 } 00429 else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD)/opac.tmn[0] < 00430 1e-3 ) 00431 { 00432 /* wavelength in eV */ 00433 EnerLast = opac.tmn[nPoints]/EVRYD; 00434 } 00435 else 00436 { 00437 fprintf( ioQQQ, " DISASTER ReadTable: the energy scale in the table read file makes no sense to me.\n" ); 00438 cdEXIT(EXIT_FAILURE); 00439 } 00440 } 00441 00442 /* now check now the energies of the highest points agree */ 00443 Differ = fabs(EnerLast-rfield.anu[nPoints])/rfield.anu[nPoints]; 00444 if( Differ > 0.001 ) 00445 { 00446 fprintf( ioQQQ, " DISASTER ReadTable: The energy mesh of the file read in by the TABLE READ command does not agree with this version of Cloudy.\n" ); 00447 fprintf( ioQQQ, " DISASTER ReadTable: Was the file generated by an older version of the code?\n" ); 00448 fprintf( ioQQQ, " DISASTER ReadTable: Did the first model do more than one iteration, but the LAST option was missing on PUNCH LAST TRANSMITTED CONTINUUM?\n"); 00449 fprintf( ioQQQ, " DISASTER ReadTable: Number of points read in=%5ld\n", 00450 nPoints ); 00451 fprintf( ioQQQ, " ReadTable: input, internal energies=%12.4e%12.4e\n", 00452 opac.tmn[nPoints], rfield.anu[nPoints] ); 00453 cdEXIT(EXIT_FAILURE); 00454 } 00455 00456 /* make sure rest of array has valid zeros */ 00457 for( i=nPoints + 1; i < rfield.nupper; i++ ) 00458 { 00459 rfield.ConTabRead[i] = 0.; 00460 } 00461 00462 fclose(io); 00463 return; 00464 } 00465