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