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