00001
00002
00003
00004
00005
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
00014 STATIC void ReadTable(const string& fnam);
00015
00016
00017
00018 double ffun(
00019
00020 double anu )
00021 {
00022 double frac_beam_time;
00023
00024 double frac_beam_const;
00025
00026 double frac_isotropic;
00027 double a;
00028
00029 DEBUG_ENTRY( "ffun()" );
00030
00031
00032 a = ffun( anu , &frac_beam_time , &frac_beam_const , &frac_isotropic );
00033 return a;
00034 }
00035
00036
00037
00038 double ffun(
00039
00040 double anu ,
00041
00042 double *frac_beam_time,
00043
00044 double *frac_beam_const,
00045
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
00055
00056
00057
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
00068
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
00081
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
00091 *frac_beam_time = flx_beam_time / ffun_v;
00092
00093 *frac_beam_const = flx_beam_const / ffun_v;
00094
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
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
00126 ASSERT( rfield.ipspec >= 0);
00127 ASSERT( rfield.ipspec < rfield.nspec );
00128
00129
00130
00131
00132
00133
00134 ASSERT( xnu >= rfield.emm*0.99 );
00135 ASSERT( xnu <= rfield.egamry*1.01 );
00136
00137
00138 strcpy( chKey, rfield.chSpType[rfield.ipspec] );
00139
00140 if( strcmp(chKey,"AGN ") == 0 )
00141 {
00142
00143
00144
00145
00146 ffun1_v = pow(xnu,-1. + rfield.cutoff[rfield.ipspec][1])*
00147 sexp(xnu/rfield.slope[rfield.ipspec])*sexp(0.01/
00148 xnu);
00149
00150 if( xnu > 0.1 )
00151 {
00152 if( xnu < 7350. )
00153 {
00154
00155
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
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
00182 fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
00183
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
00201
00202 if( xnu >= rfield.tNuRyd[rfield.ipspec][0]*1.000001 )
00203 {
00204
00205
00206 i = 1;
00207
00208
00209
00210 while( i< NCELL-1 && rfield.tNuRyd[rfield.ipspec][i]>0. )
00211 {
00212 if( xnu < rfield.tNuRyd[rfield.ipspec][i] )
00213 {
00214
00215
00216
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
00222 ffun1_v = pow(10.,y);
00223
00224
00225
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
00235 return( ffun1_v/xnu );
00236 }
00237 ++i;
00238 }
00239
00240 ffun1_v = 0.;
00241 }
00242 else
00243 {
00244
00245 ffun1_v = 0.;
00246 }
00247 }
00248
00249 else if( strcmp(chKey,"BREMS") == 0 )
00250 {
00251
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
00261
00262
00263
00264
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
00279
00280 if( rfield.tslop[rfield.ipspec][0] == 0. )
00281 {
00282
00283 ReadTable(rfield.ioTableRead[rfield.ipspec]);
00284 rfield.tslop[rfield.ipspec][0] = 1.;
00285 }
00286
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
00306
00307
00308 else if( strcmp(chKey,"VOLK ") == 0 )
00309 {
00310
00311
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
00333
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
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
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
00377 i = 0;
00378
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
00388 nPoints = i - 1;
00389
00390
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
00399
00400
00401 EnerLast = opac.tmn[nPoints];
00402 if( fabs(opac.tmn[0]-rfield.anu[0])/rfield.anu[0] > 1e-3 )
00403 {
00404
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
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
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
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
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
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
00457 for( i=nPoints + 1; i < rfield.nupper; i++ )
00458 {
00459 rfield.ConTabRead[i] = 0.;
00460 }
00461
00462 fclose(io);
00463 return;
00464 }
00465