00001
00002
00003
00004 #include "cddefines.h"
00005 #include "hypho.h"
00006
00007 #define NCM 3000
00008 #define NFREQ NCM
00009
00010
00011
00012 static double exp1(double x)
00013 {
00014 double dx,
00015 exp1_v;
00016
00017 DEBUG_ENTRY( "exp1()" );
00018
00019 dx = fabs(x);
00020 if( dx < 1.0e-9 )
00021 {
00022 exp1_v = -x;
00023 }
00024 else if( dx < 1.0e-5 )
00025 {
00026 exp1_v = ((-x*0.5) - 1.0)*x;
00027 }
00028 else if( dx < 1.0e-3 )
00029 {
00030 exp1_v = (((-x*0.1666666667) - 0.5)*x - 1.0)*x;
00031 }
00032 else
00033 {
00034 exp1_v = 1.0 - exp(x);
00035 }
00036
00037
00038 DEBUG_EXIT( "exp1()" );
00039 return( exp1_v );
00040 }
00041
00042
00043 void hypho(
00044
00045 double zed,
00046
00047 long int n,
00048
00049 long int lmin,
00050
00051 long int lmax,
00052
00053 double en,
00054
00055 long int ncell,
00056
00057 float anu[],
00058
00059 float bfnu[])
00060 {
00061 long int i,
00062 ifp,
00063 il,
00064 ilmax,
00065 j,
00066 k,
00067 l,
00068 lg,
00069 ll,
00070 llk,
00071 lll,
00072 llm,
00073 lm,
00074 lmax1,
00075 m,
00076 mulp,
00077 mulr,
00078 muls;
00079
00080
00081
00082 long *mm;
00083 double a,
00084 ai,
00085 al,
00086 alfac,
00087 con2,
00088 e,
00089 ee,
00090 fll,
00091 flm,
00092 fn,
00093 fth,
00094 g11,
00095 g12,
00096 g21,
00097 g22,
00098 g31,
00099 g32,
00100 gl,
00101 gn0,
00102 gn1e,
00103 gne,
00104 gnt,
00105 p,
00106 p1,
00107 p2,
00108 se,
00109 sl,
00110 sl4,
00111 sm,
00112 sm4,
00113 sn,
00114 sn4,
00115 sum,
00116 zed2;
00117
00118 static double ab[NFREQ],
00119 alo[7000],
00120 fal[7000],
00121 freq[NFREQ],
00122 g2[2][NFREQ],
00123 g3[2][NFREQ];
00124
00125 double anl,
00126 con3,
00127 fac,
00128 x;
00129 static int first = true;
00130 static double zero = 0.;
00131 static double one = 1.;
00132 static double two = 2.;
00133 static double four = 4.;
00134 static double ten = 10.;
00135 static double con1 = 0.8559;
00136
00137 DEBUG_ENTRY( "hypho()" );
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 mm = (long*)MALLOC((size_t)(NFREQ*sizeof(long)));
00163
00164 if( ncell > NCM )
00165 {
00166 fprintf( stderr, " increase ncm in hypho to at least%5ld\n",
00167 ncell );
00168 puts( "[Stop]" );
00169 cdEXIT(EXIT_FAILURE);
00170 }
00171
00172 if( first )
00173 {
00174 fal[0] = zero;
00175 for( i=1; i < 7000; i++ )
00176 {
00177 ai = (double)(i);
00178 alo[i-1] = log10(ai);
00179 fal[i] = alo[i-1] + fal[i-1];
00180 }
00181 first = false;
00182 }
00183
00184
00185 ll = -INT_MAX;
00186 lm = -INT_MAX;
00187 anl = -DBL_MAX;
00188
00189
00190
00191
00192
00193
00194
00195
00196 zed2 = zed*zed;
00197 fn = (double)(n);
00198 sn = fn*fn;
00199 sn4 = four*sn;
00200 con2 = con1*POW2(fn/ zed);
00201 fth = one/sn;
00202
00203 gn0 = 2.3052328943 - 2.302585093*fal[n+n-1] - fn*0.61370563888 +
00204 alo[n-1]*(fn + one)*2.30258093;
00205 lmax1 = MIN2(lmax,n-1);
00206 ilmax = n - lmin;
00207
00208
00209
00210 gl = 0.;
00211 if( lmin != 0 )
00212 {
00213 for( i=0; i < lmin; i++ )
00214 {
00215 lg = i;
00216 gl += 2.*lg + 1;
00217 }
00218 }
00219
00220 gnt = 2.*(POW2( (double)n ) - gl);
00221
00222
00223
00224 for( i=0; i < ncell; i++ )
00225 {
00226 mm[i] = 0;
00227 bfnu[i] = (float)zero;
00228 freq[i] = anu[i]*zed2;
00229 for( j=0; j < 2; j++ )
00230 {
00231 g2[j][i] = zero;
00232 g3[j][i] = zero;
00233 }
00234 }
00235
00236
00237
00238 freq[0] = MAX2(freq[0],en);
00239
00240
00241
00242 for( il=1; il <= ilmax; il++ )
00243 {
00244 l = n - il;
00245 m = 0;
00246 al = (double)(l);
00247 k = n - l - 2;
00248 con3 = con2/(two*al + one);
00249
00250
00251 for( ifp=0; ifp < ncell; ifp++ )
00252 {
00253 if( freq[ifp] < fth )
00254 {
00255 if( l <= lmax1 )
00256 anl = 0.;
00257 }
00258 else
00259 {
00260 g11 = zero;
00261 g12 = zero;
00262
00263
00264 se = freq[ifp] - fth;
00265 if( se < 0. )
00266 {
00267 fprintf( stderr, " %4ld%12.4e%12.4e\n", ifp,
00268 freq[ifp], fth );
00269 }
00270
00271
00272 e = sqrt(se);
00273 x = one + sn*se;
00274 if( k < 0 )
00275 {
00276 if( e >= 0.314 )
00277 {
00278 ee = 6.2831853/e;
00279 p = 0.5*log10(exp1(-ee));
00280 }
00281 else
00282 {
00283 p = zero;
00284 }
00285
00286 if( e > 1.0e-6 )
00287 {
00288 a = two*(fn - atan(fn*e)/e);
00289 }
00290 else
00291 {
00292 a = zero;
00293 }
00294
00295 ab[m] = (gn0 + a)/2.302585 - p - (fn + two)*
00296 log10(x);
00297 gne = 0.1;
00298 gn1e = x*gne/(fn + fn);
00299 g3[1][m] = gne;
00300 g3[0][m] = gn1e;
00301 g2[1][m] = gne*fn*x*(fn + fn - one);
00302 g2[0][m] = gn1e*(fn + fn - one)*(four +
00303 (fn - one)*x);
00304 }
00305
00306 g22 = g2[1][m];
00307 g32 = g3[1][m];
00308 g21 = g2[0][m];
00309 g31 = g3[0][m];
00310 muls = mm[m];
00311
00312 if( k < 0 )
00313 {
00314
00315
00316
00317 g11 = g31;
00318 if( l == 0 )
00319 g11 = zero;
00320 g12 = g32;
00321 }
00322
00323 else if( k == 0 )
00324 {
00325
00326
00327
00328 g11 = g21;
00329 if( l == 0 )
00330 g11 = zero;
00331 g12 = g22;
00332 }
00333
00334 else
00335 {
00336
00337
00338
00339 if( k <= 1 )
00340 {
00341 ll = n - 1;
00342 lm = n - 2;
00343 }
00344 sl = POW2( (double)ll );
00345 sl4 = four*sl;
00346 fll = (double)(ll);
00347 g12 = (sn4 - sl4 + (two*sl - fll)*x)*g22 -
00348 sn4*(sn - sl)*(one + POW2(fll + one)*se)*
00349 g32;
00350
00351 if( l != 0 )
00352 {
00353 sm = POW2( (double)lm );
00354 sm4 = four*sm;
00355 flm = (double)(lm);
00356 g11 = (sn4 - sm4 + (two*sm + flm)*x)*g21 -
00357 sn4*(sn - POW2(flm + one))*(one +
00358 sm*se)*g31;
00359 g31 = g21;
00360 g21 = g11;
00361 }
00362
00363 g32 = g22;
00364 g22 = g12;
00365
00366 if( (m+1) == ncell )
00367 {
00368 ll -= 1;
00369 if( l != 0 )
00370 lm -= 1;
00371 }
00372
00373 if( g12 >= 1.e20 )
00374 {
00375 muls += 35;
00376 g22 *= 1.e-35;
00377 g32 *= 1.e-35;
00378 g12 *= 1.e-35;
00379
00380 if( l != 0 )
00381 {
00382 g11 *= 1.e-35;
00383 g21 *= 1.e-35;
00384 g31 *= 1.e-35;
00385 }
00386
00387 }
00388 }
00389
00390 mm[m] = muls;
00391 g2[1][m] = g22;
00392 g3[1][m] = g32;
00393 g2[0][m] = g21;
00394 g3[0][m] = g31;
00395
00396 alfac = fal[n+l] - fal[n-l-1] + two*(al - fn)*
00397 alo[n*2-1];
00398
00399 p1 = one;
00400 lll = l + 1;
00401 llm = l - 1;
00402 mulr = 0;
00403
00404 if( llm >= 1 )
00405 {
00406 for( i=1; i <= llm; i++ )
00407 {
00408 ai = (double)(i);
00409 p1 *= one + ai*ai*se;
00410 if( p1 >= 1.e20 )
00411 {
00412 p1 *= 1.e-10;
00413 mulr += 10;
00414 }
00415 }
00416 }
00417
00418 p2 = p1;
00419 llk = llm + 1;
00420 if( llk < 1 )
00421 llk = 1;
00422
00423 for( i=llk; i <= lll; i++ )
00424 {
00425 ai = (double)(i);
00426 p2 *= one + ai*ai*se;
00427 }
00428
00429 mulp = 0;
00430 while( g12 >= one )
00431 {
00432 mulp -= 10;
00433 g12 *= 1.e-10;
00434 if( l != 0 )
00435 g11 *= 1.e-10;
00436 }
00437
00438 sum = alfac + (float)(mulr) + two*(ab[m] + (float)(muls-
00439 mulp+1));
00440
00441 fac = zero;
00442
00443
00444 if( fabs(sum) < 35. )
00445 fac = pow(ten,sum);
00446 if( l != 0 )
00447 g11 *= g11*p1;
00448 g12 *= g12*p2;
00449
00450
00451
00452
00453 if( l <= lmax1 )
00454 anl = fac*x*con3*(g11*al + g12*(al + 1.));
00455 anl *= 2.*(2.*al + 1.);
00456
00457 bfnu[ifp] += (float)(anl*1e-18);
00458
00459 ++m;
00460 }
00461
00462 }
00463
00464 }
00465
00466
00467 for( i=0; i < ncell; i++ )
00468 {
00469 bfnu[i] /= (float)gnt;
00470 }
00471
00472 free( mm );
00473
00474 DEBUG_EXIT( "hypho()" );
00475 return;
00476 }