21 double frac_beam_time;
23 double frac_beam_const;
25 double frac_isotropic;
31 a =
ffun( anu , &frac_beam_time , &frac_beam_const , &frac_isotropic );
41 double *frac_beam_time,
43 double *frac_beam_const,
45 double *frac_isotropic )
48 static bool lgWarn =
false;
49 double flx_beam_time , flx_beam_const , flx_isotropic;
73 flx_beam_const += one;
84 *frac_beam_const = 1.;
90 *frac_beam_time = flx_beam_time / ffun_v;
92 *frac_beam_const = flx_beam_const / ffun_v;
94 *frac_isotropic = flx_isotropic / ffun_v;
96 ASSERT( *frac_beam_time >=0. && *frac_beam_time<=1.+3.*DBL_EPSILON );
97 ASSERT( *frac_beam_const >=0.&& *frac_beam_const<=1.+3.*DBL_EPSILON );
98 ASSERT( *frac_isotropic >=0. && *frac_isotropic<=1.+3.*DBL_EPSILON );
99 ASSERT( fabs( 1.-*frac_beam_time-*frac_beam_const-*frac_isotropic)<
105 fprintf(
ioQQQ,
" FFUN: The net continuum is very intense.\n" );
106 fprintf(
ioQQQ,
" I will try to press on, but may have problems.\n" );
119 static bool lgWarn =
false;
139 if( strcmp(chKey,
"AGN ") == 0 )
145 ffun1_v = pow(xnu,-1. +
rfield.cutoff[
rfield.ipSpec][1])*
168 else if( strcmp(chKey,
"POWER") == 0 )
176 else if( strcmp(chKey,
"DISKB") == 0 )
179 double TempHi, TempLo;
197 ASSERT( TempLo < TempHi );
198 double LogDeltaT = (log10(TempHi) - log10(TempLo))/(numSteps-1.);
200 for(
long i=0; i<numSteps; i++ )
202 double Temp = pow( 10., log10(TempHi) - i * LogDeltaT );
203 double relativeWeight = pow( TempHi/Temp, 2.6666 ) * pow( 10., LogDeltaT );
208 else if( strcmp(chKey,
"BLACK") == 0 )
212 else if( strcmp(chKey,
"INTER") == 0 )
216 if( xnu >=
rfield.tNu[
rfield.ipSpec][0].Ryd()*1.000001 )
236 ffun1_v = pow(10.,y);
243 ys1 = pow( 10. , ys1 );
244 ys2 = pow( 10. , ys2 );
245 ASSERT( ffun1_v >= ys1/(1.+100.*FLT_EPSILON) );
246 ASSERT( ffun1_v <= ys2*(1.+100.*FLT_EPSILON) );
249 return( ffun1_v/xnu );
263 else if( strcmp(chKey,
"BREMS") == 0 )
267 ffun1_v =
sexp(fac)/pow(xnu,1.2);
270 else if( strcmp(chKey,
"LASER") == 0 )
272 const double BIG = 1.e10;
273 const double SMALL = 1.e-25;
290 else if( strcmp(chKey,
"READ ") == 0 )
293 if( xnu >=
rfield.egamry )
300 if( i >
rfield.nupper || i < 1 )
321 else if( strcmp(chKey,
"VOLK ") == 0 )
325 if( xnu >=
rfield.egamry )
334 fprintf(
ioQQQ,
" ffun1: Too many points - increase ncell\n" );
335 fprintf(
ioQQQ,
" cell needed=%4ld ncell=%4ld\n",
339 if( i >
rfield.nupper || i < 1 )
353 fprintf(
ioQQQ,
" ffun1: I do not understand continuum label \"%s\" for continuum %li.\n",
358 if( ffun1_v > 1e35 && !lgWarn )
361 fprintf(
ioQQQ,
" FFUN1: Continuum %ld is very intense.\n",
363 fprintf(
ioQQQ,
" I will try to press on, but may have problems.\n" );
370 const double db_log = log(DBL_MAX);
381 else if( fac > 1.e-5 )
383 ffun1_v = E_Ryd*E_Ryd/(exp(fac) - 1.);
387 ffun1_v = E_Ryd*E_Ryd/(fac*(1. + fac/2.));
393void outsum(
double *outtot,
double *outin,
double *outout)
401 for( i=0; i <
rfield.nflux; i++ )
409 *outtot = *outin + *outout;
sys_float sexp(sys_float x)
bool fp_equal(sys_float x, sys_float y, int n=3)
#define DEBUG_ENTRY(funcname)
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
void outsum(double *outtot, double *outin, double *outout)
double PlanckFunction(double Temp, double E_Ryd)
long ipoint(double energy_ryd)
UNUSED const realnum BIGFLOAT
UNUSED const double EN1RYD
UNUSED const double TE1RYD