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
00015 static void ReadTable(FILE * io);
00016
00017 double ffun(double anu)
00018 {
00019 double ffun_v;
00020 static bool lgWarn = false;
00021 double flx_time;
00022
00023 DEBUG_ENTRY( "ffun()" );
00024
00025
00026
00027
00028
00029 ffun_v = 0.;
00030 flx_time = 0.;
00031 for( rfield.ipspec=0; rfield.ipspec < rfield.nspec; rfield.ipspec++ )
00032 {
00033 double one = ffun1(anu)*rfield.spfac[rfield.ipspec];
00034 ffun_v += one;
00035
00036
00037 if( rfield.lgTimeVary[rfield.ipspec] )
00038 flx_time += one;
00039 }
00040
00041
00042
00043 rfield.frac_time = (float)(flx_time / SDIV( ffun_v));
00044 ASSERT( rfield.frac_time >=0. && rfield.frac_time<=1. );
00045
00046 if( ffun_v > 1e35 && !lgWarn )
00047 {
00048 lgWarn = true;
00049 fprintf( ioQQQ, " FFUN: The net continuum is very intense.\n" );
00050 fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
00051 }
00052
00053
00054 DEBUG_EXIT( "ffun()" );
00055 return( ffun_v );
00056 }
00057
00058
00059 double ffun1(double xnu)
00060 {
00061 char chKey[6];
00062 long int i;
00063 double fac,
00064 ffun1_v,
00065 y;
00066 static bool lgWarn = false;
00067
00068 DEBUG_ENTRY( "ffun1()" );
00069
00070
00071
00072 ASSERT( rfield.ipspec >= 0);
00073 ASSERT( rfield.ipspec < rfield.nspec );
00074
00075
00076
00077
00078
00079
00080 ASSERT( xnu >= rfield.emm*0.99 );
00081 ASSERT( xnu <= rfield.egamry*1.01 );
00082
00083
00084 strcpy( chKey, rfield.chSpType[rfield.ipspec] );
00085
00086 if( strcmp(chKey,"AGN ") == 0 )
00087 {
00088
00089
00090
00091
00092 ffun1_v = pow(xnu,-1. + rfield.cutoff[rfield.ipspec][1])*
00093 sexp(xnu/rfield.slope[rfield.ipspec])*sexp(0.01/
00094 xnu);
00095
00096 if( xnu > 0.1 )
00097 {
00098 if( xnu < 7350. )
00099 {
00100
00101
00102 ffun1_v += pow(xnu,rfield.cutoff[rfield.ipspec][2] -
00103 1.)*rfield.cutoff[rfield.ipspec][0]*sexp(1./
00104 xnu);
00105 }
00106 else
00107 {
00108 ffun1_v += pow(7350.,rfield.cutoff[rfield.ipspec][2] -
00109 1.)*rfield.cutoff[rfield.ipspec][0]/
00110 POW3(xnu/7350.);
00111 }
00112 }
00113
00114 }
00115 else if( strcmp(chKey,"POWER") == 0 )
00116 {
00117
00118 ffun1_v = pow(xnu,-1. + rfield.slope[rfield.ipspec])*
00119 sexp(xnu/rfield.cutoff[rfield.ipspec][0]+rfield.cutoff[rfield.ipspec][1]/
00120 xnu);
00121
00122 }
00123 else if( strcmp(chKey,"BLACK") == 0 )
00124 {
00125 const double db_log = log(DBL_MAX);
00126
00127
00128 fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
00129
00130 if( fac > db_log )
00131 {
00132 ffun1_v = 0.;
00133 }
00134 else if( fac > 1.e-5 )
00135 {
00136 ffun1_v = xnu*xnu/(exp(fac) - 1.);
00137 }
00138 else
00139 {
00140 ffun1_v = xnu*xnu/(fac*(1. + fac/2.));
00141 }
00142
00143 }
00144 else if( strcmp(chKey,"INTER") == 0 )
00145 {
00146
00147
00148 if( xnu >= rfield.tNuRyd[rfield.ipspec][0]*1.000001 )
00149 {
00150
00151
00152 i = 1;
00153
00154
00155
00156 while( i< NCELL-1 && rfield.tNuRyd[rfield.ipspec][i]>0. )
00157 {
00158 if( xnu < rfield.tNuRyd[rfield.ipspec][i] )
00159 {
00160
00161
00162
00163 y = rfield.tFluxLog[rfield.ipspec][i-1] +
00164 rfield.tslop[rfield.ipspec][i-1]*
00165 log10(xnu/rfield.tNuRyd[rfield.ipspec][i-1]);
00166
00167
00168 ffun1_v = pow(10.,y);
00169
00170
00171
00172 # ifndef NDEBUG
00173 double ys1 = MIN2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]);
00174 double ys2 = MAX2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]);
00175 ys1 = pow( 10. , ys1 );
00176 ys2 = pow( 10. , ys2 );
00177 ASSERT( ffun1_v >= ys1/(1.+100.*FLT_EPSILON) );
00178 ASSERT( ffun1_v <= ys2*(1.+100.*FLT_EPSILON) );
00179 # endif
00180
00181 DEBUG_EXIT( "ffun1()" );
00182
00183 return( ffun1_v/xnu );
00184 }
00185 ++i;
00186 }
00187
00188 ffun1_v = 0.;
00189 }
00190 else
00191 {
00192
00193 ffun1_v = 0.;
00194 }
00195 }
00196
00197 else if( strcmp(chKey,"BREMS") == 0 )
00198 {
00199
00200 fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
00201 ffun1_v = sexp(fac)/pow(xnu,1.2);
00202
00203 }
00204 else if( strcmp(chKey,"LASER") == 0 )
00205 {
00206 const double BIG = 1.e10;
00207 const double SMALL = 1.e-25;
00208
00209
00210
00211
00212
00213 if( xnu > (1.-rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] &&
00214 xnu < (1.+rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] )
00215 {
00216 ffun1_v = BIG;
00217 }
00218 else
00219 {
00220 ffun1_v = SMALL;
00221 }
00222
00223 }
00224 else if( strcmp(chKey,"READ ") == 0 )
00225 {
00226
00227
00228 if( rfield.tslop[rfield.ipspec][0] == 0. )
00229 {
00230
00231 ReadTable(rfield.ioTableRead[rfield.ipspec]);
00232 rfield.tslop[rfield.ipspec][0] = 1.;
00233 }
00234
00235 if( xnu >= rfield.egamry )
00236 {
00237 ffun1_v = 0.;
00238 }
00239 else
00240 {
00241 i = ipoint(xnu);
00242 if( i > rfield.nupper || i < 1 )
00243 {
00244 ffun1_v = 0.;
00245 }
00246 else
00247 {
00248 ffun1_v = rfield.ConTabRead[i-1]/rfield.anu[i-1]/rfield.anu[i-1];
00249 }
00250 }
00251 }
00252
00253
00254
00255
00256 else if( strcmp(chKey,"VOLK ") == 0 )
00257 {
00258
00259
00260 if( xnu >= rfield.egamry )
00261 {
00262 ffun1_v = 0.;
00263 }
00264 else
00265 {
00266 i = ipoint(xnu);
00267 if( i > rfield.nupper )
00268 {
00269 fprintf( ioQQQ, " ffun1: Too many points - increase ncell\n" );
00270 fprintf( ioQQQ, " cell needed=%4ld ncell=%4ld\n",
00271 i, rfield.nupper );
00272 puts( "[Stop in ffun1]" );
00273 cdEXIT(EXIT_FAILURE);
00274 }
00275 if( i > rfield.nupper || i < 1 )
00276 {
00277 ffun1_v = 0.;
00278 }
00279 else
00280 {
00281
00282
00283 ffun1_v = rfield.tslop[rfield.ipspec][i-1]/ rfield.anu[i-1];
00284 }
00285 }
00286 }
00287 else
00288 {
00289 fprintf( ioQQQ, " ffun1: I do not understand continuum label \"%s\" for continuum %li.\n",
00290 chKey , rfield.ipspec);
00291 puts( "[Stop in ffun1]" );
00292 cdEXIT(EXIT_FAILURE);
00293 }
00294
00295 if( ffun1_v > 1e35 && !lgWarn )
00296 {
00297 lgWarn = true;
00298 fprintf( ioQQQ, " FFUN1: Continuum %ld is very intense.\n",
00299 rfield.ipspec );
00300 fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
00301 }
00302
00303 DEBUG_EXIT( "ffun1()" );
00304 return( ffun1_v );
00305 }
00306
00307
00308 static void ReadTable(FILE * io)
00309 {
00310 char chLine[INPUT_LINE_LENGTH];
00311 long int i,
00312 nPoints;
00313 double Differ,
00314 EnerLast;
00315
00316 DEBUG_ENTRY( "ReadTable()" );
00317
00318
00319 ASSERT( io != NULL );
00320
00321
00322 if( NULL==fgets( chLine , (int)sizeof(chLine) , io ) )
00323 {
00324 fprintf( ioQQQ, " error 1 reading input continuum file.\n" );
00325 puts( "[Stop in ReadTable]" );
00326 cdEXIT(EXIT_FAILURE);
00327 }
00328
00329
00330 i = 0;
00331
00332 while( (NULL!=fgets(chLine, (int)sizeof(chLine),io)) && (i<rfield.nupper) )
00333 {
00334 sscanf( chLine, "%f%f ", &opac.tmn[i], &rfield.ConTabRead[i] );
00335 ++i;
00336 }
00337
00338 nPoints = i - 1;
00339
00340
00341 if( nPoints < 1 )
00342 {
00343 fprintf( ioQQQ, " ReadTable, file for TABLE READ has too few points, =%5ld\n",
00344 nPoints );
00345 puts( "[Stop in ReadTable]" );
00346 cdEXIT(EXIT_FAILURE);
00347 }
00348
00349
00350
00351
00352 EnerLast = opac.tmn[nPoints];
00353 if( fabs(opac.tmn[0]-rfield.anu[0])/rfield.anu[0] > 1e-3 )
00354 {
00355
00356 if( opac.tmn[0] <= 0. )
00357 {
00358 fprintf( ioQQQ,
00359 " DISASTER ReadTable: the first energy in table read file is not positive. Something is wrong.\n" );
00360 puts( "[Stop in ReadTable]" );
00361 cdEXIT(EXIT_FAILURE);
00362 }
00363 else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0]*1e-4)/opac.tmn[0] <
00364 1e-3 )
00365 {
00366
00367 EnerLast = RYDLAM/opac.tmn[nPoints]/1e4;
00368 }
00369 else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0])/opac.tmn[0] <
00370 1e-3 )
00371 {
00372
00373 EnerLast = RYDLAM/opac.tmn[nPoints];
00374 }
00375 else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD*1e-3)/opac.tmn[0] <
00376 1e-3 )
00377 {
00378
00379 EnerLast = opac.tmn[nPoints]/EVRYD/1e-3;
00380 }
00381 else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD)/opac.tmn[0] <
00382 1e-3 )
00383 {
00384
00385 EnerLast = opac.tmn[nPoints]/EVRYD;
00386 }
00387 else
00388 {
00389 fprintf( ioQQQ, " DISASTER ReadTable: the energy scale in the table read file makes no sense to me.\n" );
00390 puts( "[Stop in ReadTable]" );
00391 cdEXIT(EXIT_FAILURE);
00392 }
00393 }
00394
00395
00396 Differ = fabs(EnerLast-rfield.anu[nPoints])/rfield.anu[nPoints];
00397 if( Differ > 0.001 )
00398 {
00399 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" );
00400 fprintf( ioQQQ, " DISASTER ReadTable: Was the file generated by an older version of the code?\n" );
00401 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");
00402 fprintf( ioQQQ, " DISASTER ReadTable: Number of points read in=%5ld\n",
00403 nPoints );
00404 fprintf( ioQQQ, " ReadTable: input, internal energies=%12.4e%12.4e\n",
00405 opac.tmn[nPoints], rfield.anu[nPoints] );
00406 puts( "[Stop in ReadTable]" );
00407 cdEXIT(EXIT_FAILURE);
00408 }
00409
00410
00411 for( i=nPoints + 1; i < rfield.nupper; i++ )
00412 {
00413 rfield.ConTabRead[i] = 0.;
00414 }
00415
00416 fclose(io);
00417
00418 DEBUG_EXIT( "ReadTable()" );
00419 return;
00420 }
00421
00422