00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /************************************************************************************************/ 00004 /*H_photo_cs_lin returns hydrogenic photoionization cross section in cm-2 */ 00005 /*H_Einstein_A_lin calculates Einstein A for any nlz */ 00006 /*hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */ 00007 /************************************************************************************************/ 00008 /*************************** LOG version of h_bauman.c ****************************************/ 00009 /* In this version, quantities that would normally cause a 64-bit floating point processor */ 00010 /* to either underflow or overflow are evaluated using logs instead of floating point math. */ 00011 /* This allows us to use an upper principal quantum number `n' greater than which the */ 00012 /* other version begins to fail. The trade-off is, of course, lower accuracy */ 00013 /* ( or is it precision ). We use LOG_10 for convenience. */ 00014 /************************************************************************************************/ 00015 #include "cddefines.h" 00016 #include "physconst.h" 00017 #include "thirdparty.h" 00018 #include "hydro_bauman.h" 00019 00020 struct t_mx 00021 { 00022 double m; 00023 long int x; 00024 }; 00025 00026 typedef struct t_mx mx; 00027 00028 struct t_mxq 00029 { 00030 struct t_mx mx; 00031 long int q; 00032 }; 00033 00034 typedef struct t_mxq mxq; 00035 00036 /************************************************************************************************/ 00037 /* these routines were written by Robert Bauman */ 00038 /* The main reference for this section of code is */ 00039 /* M. Brocklehurst */ 00040 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 00041 /* */ 00042 /* The recombination coefficient is obtained from the */ 00043 /* photoionization cross-section (see Burgess 1965). */ 00044 /* We have: */ 00045 /* */ 00046 /* - - l'=l+1 */ 00047 /* | 2 pi^(1/2) alpha^4 a_o^2 c | 2 y^(1/2) --- */ 00048 /* alpha(n,l,Z,Te)=|-----------------------------| --------- Z > I(n, l, l', t) */ 00049 /* | 3 | n^2 --- */ 00050 /* - - l'=l-1 */ 00051 /* */ 00052 /* where OO */ 00053 /* - */ 00054 /* | */ 00055 /* I(n, l, l', t) = max(l,l') y | (1 + n^2 K^2)^2 Theta(n,l; K, l') exp( -K^2 y ) d(K^2) */ 00056 /* | */ 00057 /* - */ 00058 /* 0 */ 00059 /* */ 00060 /* Here K = k / Z */ 00061 /* */ 00062 /* */ 00063 /* and */ 00064 /* */ 00065 /* */ 00066 /* y = Z^2 Rhc/(k Te)= 15.778/t */ 00067 /* */ 00068 /* where "t" is the scaled temperature, and "Te" is the electron Temperature */ 00069 /* */ 00070 /* t = Te/(10^4 Z^2) */ 00071 /* Te in kelvin */ 00072 /* */ 00073 /* | |^2 */ 00074 /* Theta(n,l; K, l') = (1 + n^2 K^2) | g(n,l; K,l') | */ 00075 /* | | */ 00076 /* */ 00077 /* */ 00078 /* ---- Not sure if this is K or k */ 00079 /* OO / but think it is K */ 00080 /* where - v */ 00081 /* K^2 | */ 00082 /* g(n,l; K,l') = ----- | R_nl(r) F_(K,l) r dr */ 00083 /* n^2 | */ 00084 /* - */ 00085 /* 0 */ 00086 /* */ 00087 /* */ 00088 /* - - */ 00089 /* | 2 pi^(1/2) alpha^4 a_o^2 c | */ 00090 /* |-----------------------------| */ 00091 /* | 3 | */ 00092 /* - - */ 00093 /* */ 00094 /* = 2 * (3.141592654)^1/2 * (7.29735308e-3)^4 */ 00095 /* * (0.529177249e-10)^2 * (2.99792458e8) / 3 */ 00096 /* */ 00097 /* = 2.8129897e-21 */ 00098 /* Mathematica gives 2.4764282710571237e-21 */ 00099 /* */ 00100 /* The photoionization cross-section (see also Burgess 1965) */ 00101 /* is given by; */ 00102 /* _ _ l'=l+1 */ 00103 /* |4 PI alpha a_o^2 | n^2 --- max(l,l') */ 00104 /* a(Z';n,l,k)=|---------------- | --- > --------- Theta(n,l; K, l') */ 00105 /* | 3 | Z^2 --- (2l + 1) */ 00106 /* _ _ l'=l-1 */ 00107 /* */ 00108 /* */ 00109 /* where Theta(n,l; K, l') is defined above */ 00110 /************************************************************************************************/ 00111 /************************************************************************************************/ 00112 /* For the transformation: */ 00113 /* Z -> rZ = Z' */ 00114 /* */ 00115 /* k -> rk = k' */ 00116 /* then */ 00117 /* */ 00118 /* K -> K = K' */ 00119 /* */ 00120 /* and the cross-sections satisfy; */ 00121 /* 1 */ 00122 /* a(Z'; n,l,k') = --- a(Z; n,l,k) */ 00123 /* r^2 */ 00124 /* */ 00125 /* Similiarly, the recombination coefficient satisfies */ 00126 /************************************************************************************************/ 00127 00128 /************************************************************************/ 00129 /* IN THE FOLLOWING WE HAVE n > n' */ 00130 /************************************************************************/ 00131 /* returns hydrogenic photoionization cross section in cm-2 */ 00132 /* this routine is called by H_photo_cs when n is small */ 00133 STATIC double H_photo_cs_lin( 00134 /* photon energy relative to threshold energy */ 00135 double rel_photon_energy, 00136 /* principal quantum number, 1 for ground */ 00137 long int n, 00138 /* angular momentum, 0 for s */ 00139 long int l, 00140 /* charge, 1 for H+, 2 for He++, etc */ 00141 long int iz ); 00142 00167 double H_photo_cs_log10( 00168 double photon_energy, 00169 long int n, 00170 long int l, 00171 long int iz 00172 ); 00173 00174 /****************************************************************************/ 00175 /* Calculates the Einstein A's for hydrogen */ 00176 STATIC double H_Einstein_A_lin( /* IN THE FOLLOWING WE HAVE n > n' */ 00177 /* principal quantum number, 1 for ground, upper level */ 00178 long int n, 00179 /* angular momentum, 0 for s */ 00180 long int l, 00181 /* principal quantum number, 1 for ground, lower level */ 00182 long int np, 00183 /* angular momentum, 0 for s */ 00184 long int lp, 00185 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 00186 long int iz 00187 ); 00188 00202 double H_Einstein_A_log10( 00203 long int n, 00204 long int l, 00205 long int np, 00206 long int lp, 00207 long int iz 00208 ); 00209 00230 inline double OscStr_f( 00231 long int n, 00232 long int l, 00233 long int np, 00234 long int lp, 00235 long int iz 00236 ); 00237 00258 inline double OscStr_f_log10( 00259 long int n, 00260 long int l, 00261 long int np, 00262 long int lp, 00263 long int iz 00264 ); 00265 00266 /****************************************************************************** 00267 * F21() 00268 * Calculates the Hyper_Spherical_Function 2_F_1(a,b,c;y) 00269 * Here a,b, and c are (long int) 00270 * y is of type (double) 00271 * A is of type (char) and specifies whether the recursion is over 00272 * a or b. It has values A='a' or A='b'. 00273 ******************************************************************************/ 00274 00275 STATIC double F21( 00276 long int a, 00277 long int b, 00278 long int c, 00279 double y, 00280 char A 00281 ); 00282 00283 STATIC double F21i( 00284 long int a, 00285 long int b, 00286 long int c, 00287 double y, 00288 double *yV 00289 ); 00290 00291 /****************************************************************************************/ 00292 /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */ 00293 /* In the following, we have n > n' */ 00294 /****************************************************************************************/ 00295 00296 inline double hv( 00297 /* returns energy in ergs */ 00298 /* principal quantum number, 1 for ground, upper level */ 00299 long int n, 00300 /* principal quantum number, 1 for ground, lower level */ 00301 long int nprime, 00302 long int iz 00303 ); 00304 00305 /********************************************************************************/ 00306 /* In the following, we have n > n' */ 00307 /********************************************************************************/ 00308 00309 STATIC double fsff( 00310 /* principal quantum number, 1 for ground, upper level */ 00311 long int n, 00312 /* angular momentum, 0 for s */ 00313 long int l, 00314 /* principal quantum number, 1 for ground, lower level */ 00315 long int np 00316 ); 00317 00318 STATIC double log10_fsff( 00319 /* principal quantum number, 1 for ground, upper level */ 00320 long int n, 00321 /* angular momentum, 0 for s */ 00322 long int l, 00323 /* principal quantum number, 1 for ground, lower level */ 00324 long int np 00325 ); 00326 00327 /********************************************************************************/ 00328 /* F21_mx() */ 00329 /* Calculates the Hyper_Spherical_Function 2_F_1(a,b,c;y) */ 00330 /* Here a,b, and c are (long int) */ 00331 /* y is of type (double) */ 00332 /* A is of type (char) and specifies whether the recursion is over */ 00333 /* a or b. It has values A='a' or A='b'. */ 00334 /********************************************************************************/ 00335 00336 STATIC mx F21_mx( 00337 long int a, 00338 long int b, 00339 long int c, 00340 double y, 00341 char A 00342 ); 00343 00344 STATIC mx F21i_log( 00345 long int a, 00346 long int b, 00347 long int c, 00348 double y, 00349 mxq *yV 00350 ); 00351 00365 inline double hri( 00366 long int n, 00367 long int l, 00368 long int np, 00369 long int lp, 00370 long int iz 00371 ); 00372 00385 inline double hri_log10( 00386 long int n, 00387 long int l, 00388 long int np, 00389 long int lp, 00390 long int iz 00391 ); 00392 00393 /******************************************************************************/ 00394 /* In the following, we have n > n' */ 00395 /******************************************************************************/ 00396 STATIC double hrii( 00397 /* principal quantum number, 1 for ground, upper level */ 00398 long int n, 00399 /* angular momentum, 0 for s */ 00400 long int l, 00401 /* principal quantum number, 1 for ground, lower level */ 00402 long int np, 00403 /* angular momentum, 0 for s */ 00404 long int lp 00405 ); 00406 00407 STATIC double hrii_log( 00408 /* principal quantum number, 1 for ground, upper level */ 00409 long int n, 00410 /* angular momentum, 0 for s */ 00411 long int l, 00412 /* principal quantum number, 1 for ground, lower level */ 00413 long int np, 00414 /* angular momentum, 0 for s */ 00415 long int lp 00416 ); 00417 00418 STATIC double bh( 00419 double k, 00420 long int n, 00421 long int l, 00422 double *rcsvV 00423 ); 00424 00425 STATIC double bh_log( 00426 double k, 00427 long int n, 00428 long int l, 00429 mxq *rcsvV_mxq 00430 ); 00431 00432 STATIC double bhintegrand( 00433 double k, 00434 long int n, 00435 long int l, 00436 long int lp, 00437 double *rcsvV 00438 ); 00439 00440 STATIC double bhintegrand_log( 00441 double k, 00442 long int n, 00443 long int l, 00444 long int lp, 00445 mxq *rcsvV_mxq 00446 ); 00447 00448 STATIC double bhG( 00449 double K, 00450 long int n, 00451 long int l, 00452 long int lp, 00453 double *rcsvV 00454 ); 00455 00456 STATIC mx bhG_mx( 00457 double K, 00458 long int n, 00459 long int l, 00460 long int lp, 00461 mxq *rcsvV_mxq 00462 ); 00463 00464 STATIC double bhGp( 00465 long int q, 00466 double K, 00467 long int n, 00468 long int l, 00469 long int lp, 00470 double *rcsvV, 00471 double GK 00472 ); 00473 00474 STATIC mx bhGp_mx( 00475 long int q, 00476 double K, 00477 long int n, 00478 long int l, 00479 long int lp, 00480 mxq *rcsvV_mxq, 00481 const mx& GK_mx 00482 ); 00483 00484 STATIC double bhGm( 00485 long int q, 00486 double K, 00487 long int n, 00488 long int l, 00489 long int lp, 00490 double *rcsvV, 00491 double GK 00492 ); 00493 00494 STATIC mx bhGm_mx( 00495 long int q, 00496 double K, 00497 long int n, 00498 long int l, 00499 long int lp, 00500 mxq *rcsvV_mxq, 00501 const mx& GK_mx 00502 ); 00503 00504 STATIC double bhg( 00505 double K, 00506 long int n, 00507 long int l, 00508 long int lp, 00509 double *rcsvV 00510 ); 00511 00512 STATIC double bhg_log( 00513 double K, 00514 long int n, 00515 long int l, 00516 long int lp, 00517 mxq *rcsvV_mxq 00518 ); 00519 00520 inline void normalize_mx( mx& target ); 00521 inline mx add_mx( const mx& a, const mx& b ); 00522 inline mx sub_mx( const mx& a, const mx& b ); 00523 inline mx mxify( double a ); 00524 inline double unmxify( const mx& a_mx ); 00525 inline mx mxify_log10( double log10_a ); 00526 inline mx mult_mx( const mx& a, const mx& b ); 00527 00528 inline double local_product( double K , long int lp ); 00529 inline double log10_prodxx( long int lp, double Ksqrd ); 00530 00531 /****************************************************************************************/ 00532 /* 64 pi^4 (e a_o)^2 64 pi^4 (a_o)^2 e^2 1 1 */ 00533 /* ----------------- = ----------------- -------- ---- = 7.5197711e-38 ----- */ 00534 /* 3 h c^3 3 c^2 hbar c 2 pi sec */ 00535 /****************************************************************************************/ 00536 00537 static const double CONST_ONE = 32.*pow3(PI)*pow2(BOHR_RADIUS_CM)*FINE_STRUCTURE/(3.*pow2(SPEEDLIGHT)); 00538 00539 /************************************************************************/ 00540 /* (4.0/3.0) * PI * FINE_STRUCTURE_CONSTANT * BOHRRADIUS * BOHRRADIUS */ 00541 /* */ 00542 /* */ 00543 /* 4 PI alpha a_o^2 */ 00544 /* ---------------- */ 00545 /* 3 */ 00546 /* */ 00547 /* where alpha = Fine Structure Constant */ 00548 /* a_o = Bohr Radius */ 00549 /* */ 00550 /* = 3.056708^-02 (au Length)^2 */ 00551 /* = 8.56x10^-23 (meters)^2 */ 00552 /* = 8.56x10^-19 (cm)^2 */ 00553 /* = 8.56x10^+05 (barns) */ 00554 /* = 0.856 (MB or megabarns) */ 00555 /* */ 00556 /* */ 00557 /* 1 barn = 10^-28 (meter)^2 */ 00558 /************************************************************************/ 00559 00560 static const double PHYSICAL_CONSTANT_TWO = 4./3.*PI*FINE_STRUCTURE*pow2(BOHR_RADIUS_CM); 00561 00562 /************************Start of program***************************/ 00563 00564 double H_photo_cs( 00565 /* incident photon energy relative to threshold */ 00566 double rel_photon_energy, 00567 /* principal quantum number, 1 for ground */ 00568 long int n, 00569 /* angular momentum, 0 for s */ 00570 long int l, 00571 /* charge, 1 for H+, 2 for He++, etc */ 00572 long int iz ) 00573 { 00574 DEBUG_ENTRY( "H_photo_cs()" ); 00575 00576 double result; 00577 if( n<= 25 ) 00578 { 00579 result = H_photo_cs_lin( rel_photon_energy , n , l , iz ); 00580 } 00581 else 00582 { 00583 result = H_photo_cs_log10( rel_photon_energy , n , l , iz ); 00584 } 00585 return result; 00586 } 00587 00588 /************************************************************************/ 00589 /* IN THE FOLLOWING WE HAVE n > n' */ 00590 /************************************************************************/ 00591 00592 /* returns hydrogenic photoionization cross section in cm-2 */ 00593 /* this routine is called by H_photo_cs when n is small */ 00594 STATIC double H_photo_cs_lin( 00595 /* photon energy relative to threshold energy */ 00596 double rel_photon_energy, 00597 /* principal quantum number, 1 for ground */ 00598 long int n, 00599 /* angular momentum, 0 for s */ 00600 long int l, 00601 /* charge, 1 for H+, 2 for He++, etc */ 00602 long int iz ) 00603 { 00604 DEBUG_ENTRY( "H_photo_cs_lin()" ); 00605 00606 long int dim_rcsvV; 00607 00608 /* >>chng 02 sep 15, make rcsvV always NPRE_FACGTORIAL+3 long */ 00609 double rcsvV[NPRE_FACTORIAL+3]; 00610 int i; 00611 00612 double electron_energy; 00613 double result = 0.; 00614 double xn_sqrd = (double)(n*n); 00615 double z_sqrd = (double)(iz*iz); 00616 double Z = (double)iz; 00617 double K = 0.; /* K = k / Z */ 00618 double k = 0.; /* k^2 = ejected-electron-energy (Ryd) */ 00619 00620 /* expressions blow up at precisely threshold */ 00621 if( rel_photon_energy < 1.+FLT_EPSILON ) 00622 { 00623 /* below or very close to threshold, return zero */ 00624 return 0.; 00625 } 00626 00627 if( n < 1 || l >= n ) 00628 { 00629 fprintf(ioQQQ," The quantum numbers are impossible.\n"); 00630 cdEXIT(EXIT_FAILURE); 00631 } 00632 00633 if( ((2*n) - 1) >= NPRE_FACTORIAL ) 00634 { 00635 fprintf(ioQQQ," This value of n is too large.\n"); 00636 cdEXIT(EXIT_FAILURE); 00637 } 00638 00639 /* k^2 is the ejected photoelectron energy in ryd */ 00640 /*electron_energy = SDIV( (photon_energy/ryd) - (z_sqrd/xn_sqrd) );*/ 00641 00642 electron_energy = (rel_photon_energy-1.) * (z_sqrd/xn_sqrd); 00643 k = sqrt( ( electron_energy ) ); 00644 00645 K = (k/Z); 00646 00647 dim_rcsvV = (((n * 2) - 1) + 1); 00648 00649 for( i=0; i<dim_rcsvV; ++i ) 00650 { 00651 rcsvV[i] = 0.; 00652 } 00653 00654 /* rcsvV contains all results for quantum indices below n, l */ 00655 result = PHYSICAL_CONSTANT_TWO * (xn_sqrd/z_sqrd) * bh( K, n, l, rcsvV ); 00656 00657 ASSERT( result != 0. ); 00658 return result; 00659 } 00660 00661 /*****************************************************************************/ 00662 /*H_photo_cs_log10 returns hydrogenic photoionization cross section in cm-2 00663 * this routine is called by H_photo_cs when n is large */ 00664 /*****************************************************************************/ 00665 double H_photo_cs_log10( 00666 /* photon energy relative to threshold energy */ 00667 double rel_photon_energy, 00668 /* principal quantum number, 1 for ground */ 00669 long int n, 00670 /* angular momentum, 0 for s */ 00671 long int l, 00672 /* charge, 1 for H+, 2 for He++, etc */ 00673 long int iz 00674 ) 00675 { 00676 DEBUG_ENTRY( "H_photo_cs_log10()" ); 00677 00678 long int dim_rcsvV_mxq; 00679 00680 mxq *rcsvV_mxq = NULL; 00681 00682 double electron_energy; 00683 double t1; 00684 double result = 0.; 00685 double xn_sqrd = (double)(n*n); 00686 double z_sqrd = (double)(iz*iz); 00687 double Z = (double)iz; 00688 double K = 0.; /* K = k / Z */ 00689 double k = 0.; /* k^2 = ejected-electron-energy (Ryd) */ 00690 00691 /* expressions blow up at precisely threshold */ 00692 if( rel_photon_energy < 1.+FLT_EPSILON ) 00693 { 00694 /* below or very close to threshold, return zero */ 00695 fprintf( ioQQQ,"PROBLEM IN HYDRO_BAUMAN: rel_photon_energy, n, l, iz: %e\t%li\t%li\t%li\n", 00696 rel_photon_energy, 00697 n, 00698 l, 00699 iz ); 00700 cdEXIT(EXIT_FAILURE); 00701 } 00702 if( n < 1 || l >= n ) 00703 { 00704 fprintf(ioQQQ," The quantum numbers are impossible.\n"); 00705 cdEXIT(EXIT_FAILURE); 00706 } 00707 00708 /* k^2 is the ejected photoelectron energy in ryd */ 00709 /*electron_energy = SDIV( (photon_energy/ryd) - (z_sqrd/xn_sqrd) );*/ 00710 electron_energy = (rel_photon_energy-1.) * (z_sqrd/xn_sqrd); 00711 00712 k = sqrt( ( electron_energy ) ); 00713 /* k^2 is the ejected photoelectron energy in ryd */ 00714 /*k = sqrt( ( (photon_energy/ryd) - (z_sqrd/xn_sqrd) ) );*/ 00715 00716 K = (k/Z); 00717 00718 dim_rcsvV_mxq = (((n * 2) - 1) + 1); 00719 00720 /* create space */ 00721 rcsvV_mxq = (mxq*)CALLOC( (size_t)dim_rcsvV_mxq, sizeof(mxq) ); 00722 00723 t1 = bh_log( K, n, l, rcsvV_mxq ); 00724 00725 ASSERT( t1 > 0. ); 00726 00727 t1 = MAX2( t1, 1.0e-250 ); 00728 00729 result = PHYSICAL_CONSTANT_TWO * (xn_sqrd/z_sqrd) * t1; 00730 00731 free( rcsvV_mxq ); 00732 if( result <= 0. ) 00733 { 00734 fprintf( ioQQQ, "PROBLEM: Hydro_Bauman...t1\t%e\n", t1 ); 00735 } 00736 ASSERT( result > 0. ); 00737 return result; 00738 } 00739 00740 STATIC double bh( 00741 /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */ 00742 double K, 00743 /* principal quantum number */ 00744 long int n, 00745 /* angular momentum quantum number */ 00746 long int l, 00747 /* Temporary storage for intermediate */ 00748 /* results of the recursive routine */ 00749 double *rcsvV 00750 ) 00751 { 00752 DEBUG_ENTRY( "bh()" ); 00753 00754 long int lp = 0; /* l' */ 00755 double sigma = 0.; /* Sum in Brocklehurst eq. 3.13 */ 00756 00757 ASSERT( n > 0 ); 00758 ASSERT( l >= 0 ); 00759 ASSERT( n > l ); 00760 00761 if( l > 0 ) /* no lp=(l-1) for l=0 */ 00762 { 00763 for( lp = l - 1; lp <= l + 1; lp = lp + 2 ) 00764 { 00765 sigma += bhintegrand( K, n, l, lp, rcsvV ); 00766 } 00767 } 00768 else 00769 { 00770 lp = l + 1; 00771 sigma = bhintegrand( K, n, l, lp, rcsvV ); 00772 } 00773 ASSERT( sigma != 0. ); 00774 return sigma; 00775 } 00776 00777 STATIC double bh_log( 00778 /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */ 00779 double K, 00780 /* principal quantum number */ 00781 long int n, 00782 /* angular momentum quantum number */ 00783 long int l, 00784 /* Temporary storage for intermediate */ 00785 mxq *rcsvV_mxq 00786 /* results of the recursive routine */ 00787 ) 00788 { 00789 DEBUG_ENTRY( "bh_log()" ); 00790 00791 long int lp = 0; /* l' */ 00792 double sigma = 0.; /* Sum in Brocklehurst eq. 3.13 */ 00793 00794 ASSERT( n > 0 ); 00795 ASSERT( l >= 0 ); 00796 ASSERT( n > l ); 00797 00798 if( l > 0 ) /* no lp=(l-1) for l=0 */ 00799 { 00800 for( lp = l - 1; lp <= l + 1; lp = lp + 2 ) 00801 { 00802 sigma += bhintegrand_log( K, n, l, lp, rcsvV_mxq ); 00803 } 00804 } 00805 else 00806 { 00807 lp = l + 1; 00808 sigma = bhintegrand_log( K, n, l, lp, rcsvV_mxq ); 00809 } 00810 ASSERT( sigma != 0. ); 00811 return sigma; 00812 } 00813 00814 /********************************************************************************/ 00815 /* Here we calculate the integrand */ 00816 /* (as a function of K, so */ 00817 /* we need a dK^2 -> 2K dK ) */ 00818 /* for equation 3.14 of reference */ 00819 /* */ 00820 /* M. Brocklehurst Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 00821 /* */ 00822 /* namely: */ 00823 /* */ 00824 /* max(l,l') (1 + n^2 K^2)^2 Theta(n,l; K, l') exp( -K^2 y ) d(K^2) */ 00825 /* */ 00826 /* Note: the "y" is included in the code that called */ 00827 /* this function and we include here the n^2 from eq 3.13. */ 00828 /********************************************************************************/ 00829 00830 STATIC double bhintegrand( 00831 /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */ 00832 double K, 00833 long int n, 00834 long int l, 00835 long int lp, 00836 /* Temporary storage for intermediate */ 00837 /* results of the recursive routine */ 00838 double *rcsvV 00839 ) 00840 { 00841 DEBUG_ENTRY( "bhintegrand()" ); 00842 00843 double Two_L_Plus_One = (double)(2*l + 1); 00844 double lg = (double)(l > lp ? l : lp); 00845 00846 double n2 = (double)(n*n); 00847 00848 double Ksqrd = (K*K); 00849 00850 /**********************************************/ 00851 /* */ 00852 /* l> */ 00853 /* ------ Theta(nl,Kl') */ 00854 /* 2l+2 */ 00855 /* */ 00856 /* */ 00857 /* Theta(nl,Kl') = */ 00858 /* (1+n^2K^2) * | g(nl,Kl')|^2 */ 00859 /* */ 00860 /**********************************************/ 00861 00862 double d2 = 1. + n2*Ksqrd; 00863 double d5 = bhg( K, n, l, lp, rcsvV ); 00864 double Theta = d2 * d5 * d5; 00865 double d7 = (lg/Two_L_Plus_One) * Theta; 00866 00867 ASSERT( Two_L_Plus_One != 0. ); 00868 ASSERT( Theta != 0. ); 00869 ASSERT( Ksqrd != 0. ); 00870 ASSERT( d2 != 0. ); 00871 ASSERT( d5 != 0. ); 00872 ASSERT( d7 != 0. ); 00873 ASSERT( lp >= 0 ); 00874 ASSERT( lg != 0. ); 00875 ASSERT( n2 != 0. ); 00876 ASSERT( n > 0 ); 00877 ASSERT( l >= 0 ); 00878 ASSERT( K != 0. ); 00879 return d7; 00880 } 00881 00882 /************************************************************************************************/ 00883 /* The photoionization cross-section (see also Burgess 1965) */ 00884 /* is given by; */ 00885 /* _ _ l'=l+1 */ 00886 /* |4 PI alpha a_o^2 | n^2 --- max(l,l') */ 00887 /* a(Z';n,l,k)=|---------------- | --- > ---------- Theta(n,l; K, l') */ 00888 /* | 3 | Z^2 --- (2l + 1) */ 00889 /* _ _ l'=l-1 */ 00890 /* */ 00891 /* */ 00892 /* where Theta(n,l; K, l') is defined */ 00893 /* */ 00894 /* | |^2 */ 00895 /* Theta(n,l; K, l') = (1 + n^2 K^2) | g(n,l; K,l') | */ 00896 /* | | */ 00897 /* */ 00898 /* */ 00899 /* ---- Not sure if this is K or k */ 00900 /* OO / but think it is K */ 00901 /* where - v */ 00902 /* K^2 | */ 00903 /* g(n,l; K,l') = ----- | R_nl(r) F_(K,l) r dr */ 00904 /* n^2 | */ 00905 /* - */ 00906 /* 0 */ 00907 /************************************************************************************************/ 00908 00909 STATIC double bhintegrand_log( 00910 double K, /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */ 00911 long int n, 00912 long int l, 00913 long int lp, 00914 /* Temporary storage for intermediate */ 00915 /* results of the recursive routine */ 00916 mxq *rcsvV_mxq 00917 ) 00918 { 00919 DEBUG_ENTRY( "bhintegrand_log()" ); 00920 00921 double d2 = 0.; 00922 double d5 = 0.; 00923 double d7 = 0.; 00924 double Theta = 0.; 00925 double n2 = (double)(n*n); 00926 double Ksqrd = (K*K); 00927 double Two_L_Plus_One = (double)(2*l + 1); 00928 double lg = (double)(l > lp ? l : lp); 00929 00930 ASSERT( Ksqrd != 0. ); 00931 ASSERT( K != 0. ); 00932 ASSERT( lg != 0. ); 00933 ASSERT( n2 != 0. ); 00934 ASSERT( Two_L_Plus_One != 0. ); 00935 00936 ASSERT( n > 0); 00937 ASSERT( l >= 0); 00938 ASSERT( lp >= 0); 00939 00940 /**********************************************/ 00941 /* */ 00942 /* l> */ 00943 /* ------ Theta(nl,Kl') */ 00944 /* 2l+2 */ 00945 /* */ 00946 /* */ 00947 /* Theta(nl,Kl') = */ 00948 /* (1+n^2K^2) * | g(nl,Kl')|^2 */ 00949 /* */ 00950 /**********************************************/ 00951 00952 d2 = ( 1. + n2 * (Ksqrd) ); 00953 00954 ASSERT( d2 != 0. ); 00955 00956 d5 = bhg_log( K, n, l, lp, rcsvV_mxq ); 00957 00958 d5 = MAX2( d5, 1.0E-150 ); 00959 ASSERT( d5 != 0. ); 00960 00961 Theta = d2 * d5 * d5; 00962 ASSERT( Theta != 0. ); 00963 00964 d7 = (lg/Two_L_Plus_One) * Theta; 00965 00966 ASSERT( d7 != 0. ); 00967 return d7; 00968 } 00969 00970 /****************************************************************************************/ 00971 /* *** bhG *** */ 00972 /* Using various recursion relations */ 00973 /* (for l'=l+1) */ 00974 /* equation: (3.23) */ 00975 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2+l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */ 00976 /* - 4n^2 (n^2-l^2)[1+(l+1)^2 K^2] G(n,l; K,l+1) */ 00977 /* */ 00978 /* and (for l'=l-1) */ 00979 /* equation: (3.24) */ 00980 /* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l; K,l-1) */ 00981 /* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l; K,l+1) */ 00982 /* */ 00983 /* the starting point for the recursion relations are; */ 00984 /* equation: (3.18) */ 00985 /* | pi |(1/2) 8n */ 00986 /* G(n,n-1; 0,n) = | -- | ------- (4n)^n exp(-2n) */ 00987 /* | 2 | (2n-1)! */ 00988 /* */ 00989 /* equation: (3.20) */ 00990 /* exp(2n-2/K tan^(-1)(n K) */ 00991 /* G(n,n-1; K,n) = ----------------------------------------- * G(n,n-1; 0,n) */ 00992 /* sqrt(1 - exp(-2 pi K)) * (1+(n K)^2)^(n+2) */ 00993 /* */ 00994 /* equation: (3.20) */ 00995 /* G(n,n-2; K,n-1) = (2n-2)(1+(n K)^2) n G(n,n-1; K,n) */ 00996 /* */ 00997 /* equation: (3.21) */ 00998 /* (1+(n K)^2) */ 00999 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 01000 /* 2n */ 01001 /* */ 01002 /* equation: (3.22) */ 01003 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+n^2 K^2)) G(n,n-1; K,n-2) */ 01004 /****************************************************************************************/ 01005 01006 STATIC double bhG( 01007 double K, 01008 long int n, 01009 long int l, 01010 long int lp, 01011 /* Temporary storage for intermediate */ 01012 /* results of the recursive routine */ 01013 double *rcsvV 01014 ) 01015 { 01016 DEBUG_ENTRY( "bhG()" ); 01017 01018 double n1 = (double)n; 01019 double n2 = (double)(n * n); 01020 double Ksqrd = K * K; 01021 01022 double ld1 = factorial( 2*n - 1 ); 01023 double ld2 = powi((double)(4*n), n); 01024 double ld3 = exp(-(double)(2 * n)); 01025 01026 /****************************************************************************** 01027 * ********G0******* * 01028 * * 01029 * | pi |(1/2) 8n * 01030 * G0 = | -- | ------- (4n)^n exp(-2n) * 01031 * | 2 | (2n-1)! * 01032 ******************************************************************************/ 01033 01034 double G0 = SQRTPIBY2 * (8. * n1 * ld2 * ld3) / ld1; 01035 01036 double d1 = sqrt( 1. - exp(( -2. * PI )/ K )); 01037 double d2 = powi(( 1. + (n2 * Ksqrd)), ( n + 2 )); 01038 double d3 = atan( n1 * K ); 01039 double d4 = ((2. / K) * d3); 01040 double d5 = (double)( 2 * n ); 01041 double d6 = exp( d5 - d4 ); 01042 double GK = ( d6 /( d1 * d2 ) ) * G0; 01043 01044 /* l=l'-1 or l=l'+1 */ 01045 ASSERT( (l == lp - 1) || (l == lp + 1) ); 01046 ASSERT( K != 0. ); 01047 ASSERT( Ksqrd != 0. ); 01048 ASSERT( n1 != 0. ); 01049 ASSERT( n2 != 0. ); 01050 ASSERT( ((2*n) - 1) < 1755 ); 01051 ASSERT( ((2*n) - 1) >= 0 ); 01052 ASSERT( ld1 != 0. ); 01053 ASSERT( (1.0 / ld1) != 0. ); 01054 ASSERT( ld3 != 0. ); 01055 01056 ASSERT( K != 0. ); 01057 ASSERT( d1 != 0. ); 01058 ASSERT( d2 != 0. ); 01059 ASSERT( d3 != 0. ); 01060 ASSERT( d4 != 0. ); 01061 ASSERT( d5 != 0. ); 01062 ASSERT( d6 != 0. ); 01063 01064 ASSERT( G0 != 0. ); 01065 ASSERT( GK != 0. ); 01066 01067 /******************************************************************************/ 01068 /* *****GK***** */ 01069 /* */ 01070 /* exp(2n-2/K tan^(-1)(n K) */ 01071 /* G(n,n-1; K,n) = ----------------------------------------- * G0 */ 01072 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K))^(n+2) */ 01073 /******************************************************************************/ 01074 01075 /* GENERAL CASE: l = l'-1 */ 01076 if( l == lp - 1 ) 01077 { 01078 double result = bhGm( l, K, n, l, lp, rcsvV, GK ); 01079 /* Here the m in bhGm() refers */ 01080 /* to the minus sign(-) in l=l'-1 */ 01081 return result; 01082 } 01083 01084 /* GENERAL CASE: l = l'+1 */ 01085 else if( l == lp + 1 ) 01086 { 01087 double result = bhGp( l, K, n, l, lp, rcsvV, GK ); 01088 /* Here the p in bhGp() refers */ 01089 /* to the plus sign(+) in l=l'+1 */ 01090 return result; 01091 } 01092 else 01093 { 01094 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" ); 01095 cdEXIT(EXIT_FAILURE); 01096 } 01097 } 01098 01099 /*************log version********************************/ 01100 STATIC mx bhG_mx( 01101 double K, 01102 long int n, 01103 long int l, 01104 long int lp, 01105 /* Temporary storage for intermediate */ 01106 /* results of the recursive routine */ 01107 mxq *rcsvV_mxq 01108 ) 01109 { 01110 DEBUG_ENTRY( "bhG_mx()" ); 01111 01112 double log10_GK = 0.; 01113 double log10_G0 = 0.; 01114 01115 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.; 01116 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0.; 01117 01118 double n1 = (double)n; 01119 double n2 = n1 * n1; 01120 double Ksqrd = K * K; 01121 01122 mx GK_mx = {0.0,0}; 01123 01124 /* l=l'-1 or l=l'+1 */ 01125 ASSERT( (l == lp - 1) || (l == lp + 1) ); 01126 ASSERT( K != 0. ); 01127 ASSERT( n1 != 0. ); 01128 ASSERT( n2 != 0. ); 01129 ASSERT( Ksqrd != 0. ); 01130 ASSERT( ((2*n) - 1) >= 0 ); 01131 01132 /******************************/ 01133 /* n */ 01134 /* --- */ 01135 /* log( n! ) = > log(j) */ 01136 /* --- */ 01137 /* j=1 */ 01138 /******************************/ 01139 01140 /*************************************************************/ 01141 /* | pi |(1/2) 8n */ 01142 /* G(n,n-1; 0,n) = | -- | ------- (4n)^n exp(-2n) */ 01143 /* | 2 | (2n-1)! */ 01144 /*************************************************************/ 01145 01146 /******************************/ 01147 /* */ 01148 /* */ 01149 /* log10( (2n-1)! ) */ 01150 /* */ 01151 /* */ 01152 /******************************/ 01153 01154 ld1 = lfactorial( 2*n - 1 ); 01155 ASSERT( ld1 >= 0. ); 01156 01157 /**********************************************/ 01158 /* (4n)^n */ 01159 /**********************************************/ 01160 /* log10( 4n^n ) = n log10( 4n ) */ 01161 /**********************************************/ 01162 01163 ld2 = n1 * log10( 4. * n1 ); 01164 ASSERT( ld2 >= 0. ); 01165 01166 /**********************************************/ 01167 /* exp(-2n) */ 01168 /**********************************************/ 01169 /* log10( exp( -2n ) ) = (-2n) * log10(e) */ 01170 /**********************************************/ 01171 ld3 = (-(2. * n1)) * (LOG10_E); 01172 ASSERT( ld3 <= 0. ); 01173 01174 /******************************************************************************/ 01175 /* ********G0******* */ 01176 /* */ 01177 /* | pi |(1/2) 8n */ 01178 /* G0 = | -- | ------- (4n)^n exp(-2n) */ 01179 /* | 2 | (2n-1)! */ 01180 /******************************************************************************/ 01181 01182 log10_G0 = log10(SQRTPIBY2 * 8. * n1) + ( (ld2 + ld3) - ld1); 01183 01184 /******************************************************************************/ 01185 /* *****GK***** */ 01186 /* */ 01187 /* exp(2n- (2/K) tan^(-1)(n K) ) */ 01188 /* G(n,n-1; K,n) = ----------------------------------------- * G0 */ 01189 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K))^(n+2) */ 01190 /******************************************************************************/ 01191 01192 ASSERT( K != 0. ); 01193 01194 /**********************************************/ 01195 /* sqrt(1 - exp(-2 pi/ K)) */ 01196 /**********************************************/ 01197 /* log10(sqrt(1 - exp(-2 pi/ K))) = */ 01198 /* (1/2) log10(1 - exp(-2 pi/ K)) */ 01199 /**********************************************/ 01200 01201 d1 = (1. - exp(-(2. * PI )/ K )); 01202 ld4 = (1./2.) * log10( d1 ); 01203 ASSERT( K != 0. ); 01204 ASSERT( d1 != 0. ); 01205 01206 /**************************************/ 01207 /* (1+(n K)^2)^(n+2) */ 01208 /**************************************/ 01209 /* log10( (1+(n K)^2)^(n+2) ) = */ 01210 /* (n+2) log10( (1 + (n K)^2 ) ) */ 01211 /**************************************/ 01212 01213 d2 = ( 1. + (n2 * Ksqrd)); 01214 ld5 = (n1 + 2.) * log10( d2 ); 01215 ASSERT( d2 != 0. ); 01216 01217 ASSERT( ld5 >= 0. ); 01218 01219 /**********************************************/ 01220 /* exp(2n- (2/K)*tan^(-1)(n K) ) */ 01221 /**********************************************/ 01222 /* log10( exp(2n- (2/K) tan^(-1)(n K) ) = */ 01223 /* (2n- (2/K)*tan^(-1)(n K) ) * Log10(e) */ 01224 /**********************************************/ 01225 01226 /* tan^(-1)(n K) ) */ 01227 d3 = atan( n1 * K ); 01228 ASSERT( d3 != 0. ); 01229 01230 /* (2/K)*tan^(-1)(n K) ) */ 01231 d4 = (2. / K) * d3; 01232 ASSERT( d4 != 0. ); 01233 01234 /* 2n */ 01235 d5 = (double) ( 2. * n1 ); 01236 ASSERT( d5 != 0. ); 01237 01238 /* (2n-2/K tan^(-1)(n K)) */ 01239 d6 = d5 - d4; 01240 ASSERT( d6 != 0. ); 01241 01242 /* log10( exp(2n- (2/K) tan^(-1)(n K) ) */ 01243 ld6 = LOG10_E * d6; 01244 ASSERT( ld6 != 0. ); 01245 01246 /******************************************************************************/ 01247 /* *****GK***** */ 01248 /* */ 01249 /* exp(2n- (2/K) tan^(-1)(n K) ) */ 01250 /* G(n,n-1; K,n) = ----------------------------------------- * G0 */ 01251 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K))^(n+2) */ 01252 /******************************************************************************/ 01253 01254 log10_GK = (ld6 -(ld4 + ld5)) + log10_G0; 01255 ASSERT( log10_GK != 0. ); 01256 01257 GK_mx = mxify_log10( log10_GK ); 01258 01259 /* GENERAL CASE: l = l'-1 */ 01260 if( l == lp - 1 ) 01261 { 01262 mx result_mx = bhGm_mx( l, K, n, l, lp, rcsvV_mxq , GK_mx ); 01263 /* Here the m in bhGm() refers */ 01264 /* to the minus sign(-) in l=l'-1 */ 01265 return result_mx; 01266 } 01267 /* GENERAL CASE: l = l'+1 */ 01268 else if( l == lp + 1 ) 01269 { 01270 mx result_mx = bhGp_mx( l, K, n, l, lp, rcsvV_mxq , GK_mx ); 01271 /* Here the p in bhGp() refers */ 01272 /* to the plus sign(+) in l=l'+1 */ 01273 return result_mx; 01274 } 01275 else 01276 { 01277 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" ); 01278 cdEXIT(EXIT_FAILURE); 01279 } 01280 } 01281 01282 /************************************************************************************************/ 01283 /* *** bhGp.c *** */ 01284 /* */ 01285 /* Here we calculate G(n,l; K,l') with the recursive formula */ 01286 /* equation: (3.24) */ 01287 /* */ 01288 /* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */ 01289 /* */ 01290 /* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l+1; K,l) */ 01291 /* */ 01292 /* Under the transformation l -> l + 1 this gives */ 01293 /* */ 01294 /* G(n,l+1-1; K,l+1-2) = [ 4n^2-4(l+1)^2 + (l+1)(2(l+1)+1)(1+(n K)^2) ] G(n,l+1; K,l+1-1) */ 01295 /* */ 01296 /* - 4n^2 (n^2-((l+1)+1)^2)[ 1+((l+1)K)^2 ] G(n,l+1+1; K,l+1) */ 01297 /* */ 01298 /* or */ 01299 /* */ 01300 /* G(n,l; K,l-1) = [ 4n^2-4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l) */ 01301 /* */ 01302 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01303 /* */ 01304 /* from the reference */ 01305 /* M. Brocklehurst */ 01306 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 01307 /* */ 01308 /* */ 01309 /* * This is valid for the case l=l'+1 * */ 01310 /* * CASE: l = l'+1 * */ 01311 /* * Here the p in bhGp() refers * */ 01312 /* * to the Plus sign(+) in l=l'+1 * */ 01313 /************************************************************************************************/ 01314 01315 STATIC double bhGp( 01316 long int q, 01317 double K, 01318 long int n, 01319 long int l, 01320 long int lp, 01321 /* Temporary storage for intermediate */ 01322 /* results of the recursive routine */ 01323 double *rcsvV, 01324 double GK 01325 ) 01326 { 01327 DEBUG_ENTRY( "bhGp()" ); 01328 01329 /* static long int rcsv_Level = 1; 01330 printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */ 01331 01332 ASSERT( l == lp + 1 ); 01333 01334 long int rindx = 2*q; 01335 01336 if( rcsvV[rindx] == 0. ) 01337 { 01338 /* SPECIAL CASE: n = l+1 = l'+2 */ 01339 if( q == n - 1 ) 01340 { 01341 double Ksqrd = K * K; 01342 double n2 = (double)(n*n); 01343 01344 double dd1 = (double)(2 * n); 01345 double dd2 = ( 1. + ( n2 * Ksqrd)); 01346 01347 /* (1+(n K)^2) */ 01348 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 01349 /* 2n */ 01350 double G1 = ((dd2 * GK) / dd1); 01351 01352 ASSERT( l == lp + 1 ); 01353 ASSERT( Ksqrd != 0. ); 01354 ASSERT( dd1 != 0. ); 01355 ASSERT( dd2 != 0. ); 01356 ASSERT( G1 != 0. ); 01357 01358 rcsvV[rindx] = G1; 01359 return G1; 01360 } 01361 /* SPECIAL CASE: n = l+2 = l'+3 */ 01362 else if( q == (n - 2) ) 01363 { 01364 double Ksqrd = (K*K); 01365 double n2 = (double)(n*n); 01366 01367 /* */ 01368 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */ 01369 /* */ 01370 double dd1 = (double) (2 * n); 01371 double dd2 = ( 1. + ( n2 * Ksqrd)); 01372 double G1 = ((dd2 * GK) / dd1); 01373 01374 /* */ 01375 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */ 01376 /* */ 01377 double dd3 = (double)((2 * n) - 1); 01378 double dd4 = (double)(n - 1); 01379 double dd5 = (4. + (dd4 * dd2)); 01380 double G2 = (dd3 * dd5 * G1); 01381 01382 ASSERT( l == lp + 1 ); 01383 ASSERT( Ksqrd != 0. ); 01384 ASSERT( n2 != 0. ); 01385 ASSERT( dd1 != 0. ); 01386 ASSERT( dd2 != 0. ); 01387 ASSERT( dd3 != 0. ); 01388 ASSERT( dd4 != 0. ); 01389 ASSERT( dd5 != 0. ); 01390 ASSERT( G1 != 0. ); 01391 ASSERT( G2 != 0. ); 01392 01393 rcsvV[rindx] = G2; 01394 return G2; 01395 } 01396 /* The GENERAL CASE n > l + 2 */ 01397 else 01398 { 01399 /******************************************************************************/ 01400 /* G(n,l; K,l-1) = [ 4n^2-4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l) */ 01401 /* */ 01402 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01403 /* */ 01404 /* FROM Eq. 3.24 */ 01405 /* */ 01406 /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */ 01407 /* */ 01408 /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */ 01409 /******************************************************************************/ 01410 01411 long int lp1 = (q + 1); 01412 long int lp2 = (q + 2); 01413 01414 double Ksqrd = (K*K); 01415 double n2 = (double)(n * n); 01416 double lp1s = (double)(lp1 * lp1); 01417 double lp2s = (double)(lp2 * lp2); 01418 01419 double d1 = (4. * n2); 01420 double d2 = (4. * lp1s); 01421 double d3 = (double)((lp1)*((2 * q) + 3)); 01422 double d4 = (1. + (n2 * Ksqrd)); 01423 double d5 = (d1 - d2 + (d3 * d4)); 01424 double d5_1 = d5 * bhGp( (q+1), K, n, l, lp, rcsvV, GK ); 01425 01426 01427 /* G(n,l; K,l-1) = [ 4n^2-4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l) */ 01428 /* */ 01429 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01430 01431 double d6 = (n2 - lp2s); 01432 double d7 = (1. + (lp1s * Ksqrd)); 01433 double d8 = (d1 * d6 * d7); 01434 double d8_1 = d8 * bhGp( (q+2), K, n, l, lp, rcsvV, GK ); 01435 double d9 = (d5_1 - d8_1); 01436 01437 ASSERT( l == lp + 1 ); 01438 ASSERT( Ksqrd != 0. ); 01439 ASSERT( n2 != 0. ); 01440 01441 ASSERT( lp1s != 0. ); 01442 ASSERT( lp2s != 0. ); 01443 01444 ASSERT( d1 != 0. ); 01445 ASSERT( d2 != 0. ); 01446 ASSERT( d3 != 0. ); 01447 ASSERT( d4 != 0. ); 01448 ASSERT( d5 != 0. ); 01449 ASSERT( d6 != 0. ); 01450 ASSERT( d7 != 0. ); 01451 ASSERT( d8 != 0. ); 01452 ASSERT( d9 != 0. ); 01453 01454 rcsvV[rindx] = d9; 01455 return d9; 01456 } 01457 } 01458 else 01459 { 01460 ASSERT( rcsvV[rindx] != 0. ); 01461 return rcsvV[rindx]; 01462 } 01463 } 01464 01465 /***********************log version*******************************/ 01466 STATIC mx bhGp_mx( 01467 long int q, 01468 double K, 01469 long int n, 01470 long int l, 01471 long int lp, 01472 /* Temporary storage for intermediate */ 01473 /* results of the recursive routine */ 01474 mxq *rcsvV_mxq, 01475 const mx& GK_mx 01476 ) 01477 { 01478 DEBUG_ENTRY( "bhGp_mx()" ); 01479 01480 /* static long int rcsv_Level = 1; */ 01481 /* printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */ 01482 01483 ASSERT( l == lp + 1 ); 01484 01485 long int rindx = 2*q; 01486 01487 if( rcsvV_mxq[rindx].q == 0 ) 01488 { 01489 /* SPECIAL CASE: n = l+1 = l'+2 */ 01490 if( q == n - 1 ) 01491 { 01492 /******************************************************/ 01493 /* (1+(n K)^2) */ 01494 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 01495 /* 2n */ 01496 /******************************************************/ 01497 01498 double Ksqrd = (K * K); 01499 double n2 = (double)(n*n); 01500 01501 double dd1 = (double) (2 * n); 01502 double dd2 = ( 1. + ( n2 * Ksqrd)); 01503 double dd3 = dd2/dd1; 01504 01505 mx dd3_mx = mxify( dd3 ); 01506 mx G1_mx = mult_mx( dd3_mx, GK_mx); 01507 01508 normalize_mx( G1_mx ); 01509 01510 ASSERT( l == lp + 1 ); 01511 ASSERT( Ksqrd != 0. ); 01512 ASSERT( n2 != 0. ); 01513 ASSERT( dd1 != 0. ); 01514 ASSERT( dd2 != 0. ); 01515 01516 rcsvV_mxq[rindx].q = 1; 01517 rcsvV_mxq[rindx].mx = G1_mx; 01518 return G1_mx; 01519 } 01520 /* SPECIAL CASE: n = l+2 = l'+3 */ 01521 else if( q == (n - 2) ) 01522 { 01523 /****************************************************************/ 01524 /* */ 01525 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/ 01526 /* */ 01527 /****************************************************************/ 01528 /* (1+(n K)^2) */ 01529 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 01530 /* 2n */ 01531 /****************************************************************/ 01532 01533 double Ksqrd = (K*K); 01534 double n2 = (double)(n*n); 01535 double dd1 = (double)(2 * n); 01536 double dd2 = ( 1. + ( n2 * Ksqrd) ); 01537 double dd3 = (dd2/dd1); 01538 double dd4 = (double) ((2 * n) - 1); 01539 double dd5 = (double) (n - 1); 01540 double dd6 = (4. + (dd5 * dd2)); 01541 double dd7 = dd4 * dd6; 01542 01543 /****************************************************************/ 01544 /* */ 01545 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/ 01546 /* */ 01547 /****************************************************************/ 01548 01549 mx dd3_mx = mxify( dd3 ); 01550 mx dd7_mx = mxify( dd7 ); 01551 mx G1_mx = mult_mx( dd3_mx, GK_mx ); 01552 mx G2_mx = mult_mx( dd7_mx, G1_mx ); 01553 01554 normalize_mx( G2_mx ); 01555 01556 ASSERT( l == lp + 1 ); 01557 ASSERT( Ksqrd != 0. ); 01558 ASSERT( n2 != 0. ); 01559 ASSERT( dd1 != 0. ); 01560 ASSERT( dd2 != 0. ); 01561 ASSERT( dd3 != 0. ); 01562 ASSERT( dd4 != 0. ); 01563 ASSERT( dd5 != 0. ); 01564 ASSERT( dd6 != 0. ); 01565 ASSERT( dd7 != 0. ); 01566 01567 rcsvV_mxq[rindx].q = 1; 01568 rcsvV_mxq[rindx].mx = G2_mx; 01569 return G2_mx; 01570 } 01571 /* The GENERAL CASE n > l + 2*/ 01572 else 01573 { 01574 /**************************************************************************************/ 01575 /* G(n,l; K,l-1) = [ 4n^2-4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l) */ 01576 /* */ 01577 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01578 /* */ 01579 /* FROM Eq. 3.24 */ 01580 /* */ 01581 /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */ 01582 /* */ 01583 /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */ 01584 /**************************************************************************************/ 01585 01586 long int lp1 = (q + 1); 01587 long int lp2 = (q + 2); 01588 01589 double Ksqrd = (K * K); 01590 double n2 = (double)(n * n); 01591 double lp1s = (double)(lp1 * lp1); 01592 double lp2s = (double)(lp2 * lp2); 01593 01594 double d1 = (4. * n2); 01595 double d2 = (4. * lp1s); 01596 double d3 = (double)((lp1)*((2 * q) + 3)); 01597 double d4 = (1. + (n2 * Ksqrd)); 01598 /* [ 4n^2 - 4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] */ 01599 double d5 = d1 - d2 + (d3 * d4); 01600 01601 /* (n^2-(l+2)^2) */ 01602 double d6 = (n2 - lp2s); 01603 01604 /* [ 1+((l+1)K)^2 ] */ 01605 double d7 = (1. + (lp1s * Ksqrd)); 01606 01607 /* { 4n^2 (n^2-(l+1)^2)[ 1+((l+1) K)^2 ] } */ 01608 double d8 = (d1 * d6 * d7); 01609 01610 /**************************************************************************************/ 01611 /* G(n,l; K,l-1) = [ 4n^2 - 4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l) */ 01612 /* */ 01613 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01614 /**************************************************************************************/ 01615 01616 mx d5_mx=mxify( d5 ); 01617 mx d8_mx=mxify( d8 ); 01618 01619 mx t0_mx = bhGp_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx ); 01620 mx t1_mx = bhGp_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx ); 01621 01622 mx d9_mx = mult_mx( d5_mx, t0_mx ); 01623 mx d10_mx = mult_mx( d8_mx, t1_mx ); 01624 01625 mx result_mx = sub_mx( d9_mx, d10_mx ); 01626 normalize_mx( result_mx ); 01627 01628 ASSERT( d1 != 0. ); 01629 ASSERT( d2 != 0. ); 01630 ASSERT( d3 != 0. ); 01631 ASSERT( d4 != 0. ); 01632 ASSERT( d5 != 0. ); 01633 ASSERT( d6 != 0. ); 01634 ASSERT( d7 != 0. ); 01635 ASSERT( d8 != 0. ); 01636 01637 ASSERT( l == lp + 1 ); 01638 ASSERT( Ksqrd != 0. ); 01639 ASSERT( n2 != 0. ); 01640 ASSERT( lp1s != 0. ); 01641 ASSERT( lp2s != 0. ); 01642 01643 rcsvV_mxq[rindx].q = 1; 01644 rcsvV_mxq[rindx].mx = result_mx; 01645 return result_mx; 01646 } 01647 } 01648 else 01649 { 01650 ASSERT( rcsvV_mxq[rindx].q != 0 ); 01651 rcsvV_mxq[rindx].q = 1; 01652 return rcsvV_mxq[rindx].mx; 01653 } 01654 } 01655 01656 /************************************************************************************************/ 01657 /* *** bhGm.c *** */ 01658 /* */ 01659 /* Here we calculate G(n,l; K,l') with the recursive formula */ 01660 /* equation: (3.23) */ 01661 /* */ 01662 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */ 01663 /* */ 01664 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */ 01665 /* */ 01666 /* Under the transformation l -> l + 2 this gives */ 01667 /* */ 01668 /* G(n,l+2-2; K,l+2-1) = [ 4n^2-4(l+2)^2 + (l+2)(2(l+2)-1)(1+(n K)^2) ] G(n,l+2-1; K,l+2) */ 01669 /* */ 01670 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+2+1)^2 K^2 ] G(n,l+2; K,l+2+1) */ 01671 /* */ 01672 /* */ 01673 /* or */ 01674 /* */ 01675 /* G(n,l; K,l+1) = [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l+2) */ 01676 /* */ 01677 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */ 01678 /* */ 01679 /* */ 01680 /* from the reference */ 01681 /* M. Brocklehurst */ 01682 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 01683 /* */ 01684 /* */ 01685 /* * This is valid for the case l=l'-1 * */ 01686 /* * CASE: l = l'-1 * */ 01687 /* * Here the p in bhGm() refers * */ 01688 /* * to the Minus sign(-) in l=l'-1 * */ 01689 /************************************************************************************************/ 01690 01691 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000 01692 #pragma optimize("", off) 01693 #endif 01694 STATIC double bhGm( 01695 long int q, 01696 double K, 01697 long int n, 01698 long int l, 01699 long int lp, 01700 double *rcsvV, 01701 double GK 01702 ) 01703 { 01704 DEBUG_ENTRY( "bhGm()" ); 01705 01706 ASSERT( l == lp - 1 ); 01707 ASSERT( l < n ); 01708 01709 long int rindx = 2*q+1; 01710 01711 if( rcsvV[rindx] == 0. ) 01712 { 01713 /* CASE: l = n - 1 */ 01714 if( q == n - 1 ) 01715 { 01716 ASSERT( l == lp - 1 ); 01717 rcsvV[rindx] = GK; 01718 return GK; 01719 } 01720 /* CASE: l = n - 2 */ 01721 else if( q == n - 2 ) 01722 { 01723 double dd1 = 0.; 01724 double dd2 = 0.; 01725 01726 double G2 = 0.; 01727 01728 double Ksqrd = 0.; 01729 double n1 = 0.; 01730 double n2 = 0.; 01731 01732 ASSERT(l == lp - 1); 01733 01734 /* K^2 */ 01735 Ksqrd = K * K; 01736 ASSERT( Ksqrd != 0. ); 01737 01738 /* n */ 01739 n1 = (double)n; 01740 ASSERT( n1 != 0. ); 01741 01742 /* n^2 */ 01743 n2 = (double)(n*n); 01744 ASSERT( n2 != 0. ); 01745 01746 /* equation: (3.20) */ 01747 /* G(n,n-2; K,n-1) = */ 01748 /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */ 01749 dd1 = (double) ((2 * n) - 1); 01750 ASSERT( dd1 != 0. ); 01751 01752 dd2 = (1. + (n2 * Ksqrd)); 01753 ASSERT( dd2 != 0. ); 01754 01755 G2 = dd1 * dd2 * n1 * GK; 01756 ASSERT( G2 != 0. ); 01757 01758 rcsvV[rindx] = G2; 01759 ASSERT( G2 != 0. ); 01760 return G2; 01761 } 01762 else 01763 { 01764 long int lp2 = (q + 2); 01765 long int lp3 = (q + 3); 01766 01767 double lp2s = (double)(lp2 * lp2); 01768 double lp3s = (double)(lp3 * lp3); 01769 01770 /******************************************************************************/ 01771 /* G(n,l; K,l+1) = [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l+2) */ 01772 /* */ 01773 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */ 01774 /* */ 01775 /* */ 01776 /* FROM Eq. 3.23 */ 01777 /* */ 01778 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2 + (l+2)(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */ 01779 /* */ 01780 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */ 01781 /******************************************************************************/ 01782 01783 double Ksqrd = (K*K); 01784 double n2 = (double)(n*n); 01785 double d1 = (4. * n2); 01786 double d2 = (4. * lp2s); 01787 double d3 = (double)(lp2)*((2*q)+3); 01788 double d4 = (1. + (n2 * Ksqrd)); 01789 /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */ 01790 double d5 = d1 - d2 + (d3 * d4); 01791 01792 /******************************************************************************/ 01793 /* G(n,l; K,l+1) = [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l+2) */ 01794 /* */ 01795 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */ 01796 /******************************************************************************/ 01797 01798 double d6 = (n2 - lp2s); 01799 /* [ 1+((l+3)K)^2 ] */ 01800 double d7 = (1. + (lp3s * Ksqrd)); 01801 /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */ 01802 double d8 = d1 * d6 * d7; 01803 double d9 = d5 * bhGm( (q+1), K, n, l, lp, rcsvV, GK ); 01804 double d10 = d8 * bhGm( (q+2), K, n, l, lp, rcsvV, GK ); 01805 double d11 = d9 - d10; 01806 01807 ASSERT( l == lp - 1 ); 01808 ASSERT( lp2s != 0. ); 01809 ASSERT( Ksqrd != 0. ); 01810 ASSERT( n2 != 0. ); 01811 ASSERT( d1 != 0. ); 01812 ASSERT( d2 != 0. ); 01813 ASSERT( d3 != 0. ); 01814 ASSERT( d4 != 0. ); 01815 ASSERT( d5 != 0. ); 01816 ASSERT( d6 != 0. ); 01817 ASSERT( d7 != 0. ); 01818 ASSERT( d8 != 0. ); 01819 ASSERT( d9 != 0. ); 01820 ASSERT( d10 != 0. ); 01821 ASSERT( lp3s != 0. ); 01822 01823 rcsvV[rindx] = d11; 01824 return d11; 01825 } 01826 } 01827 else 01828 { 01829 ASSERT( rcsvV[rindx] != 0. ); 01830 return rcsvV[rindx]; 01831 } 01832 } 01833 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000 01834 #pragma optimize("", on) 01835 #endif 01836 01837 /************************log version***********************************/ 01838 STATIC mx bhGm_mx( 01839 long int q, 01840 double K, 01841 long int n, 01842 long int l, 01843 long int lp, 01844 mxq *rcsvV_mxq, 01845 const mx& GK_mx 01846 ) 01847 { 01848 DEBUG_ENTRY( "bhGm_mx()" ); 01849 01850 /*static long int rcsv_Level = 1; */ 01851 /*printf( "bhGm(): recursion level:\t%li\n",rcsv_Level++ ); */ 01852 01853 ASSERT( l == lp - 1 ); 01854 ASSERT( l < n ); 01855 01856 long int rindx = 2*q+1; 01857 01858 if( rcsvV_mxq[rindx].q == 0 ) 01859 { 01860 /* CASE: l = n - 1 */ 01861 if( q == n - 1 ) 01862 { 01863 mx result_mx = GK_mx; 01864 normalize_mx( result_mx ); 01865 01866 rcsvV_mxq[rindx].q = 1; 01867 rcsvV_mxq[rindx].mx = result_mx; 01868 01869 ASSERT(l == lp - 1); 01870 return result_mx; 01871 } 01872 /* CASE: l = n - 2 */ 01873 else if( q == n - 2 ) 01874 { 01875 double Ksqrd = (K * K); 01876 double n1 = (double)n; 01877 double n2 = (double) (n*n); 01878 double dd1 = (double)((2 * n) - 1); 01879 double dd2 = (1. + (n2 * Ksqrd)); 01880 /*(2n-1)(1+(n K)^2) n*/ 01881 double dd3 = (dd1*dd2*n1); 01882 01883 /******************************************************/ 01884 /* G(n,n-2; K,n-1) = */ 01885 /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */ 01886 /******************************************************/ 01887 01888 mx dd3_mx = mxify( dd3 ); 01889 mx G2_mx = mult_mx( dd3_mx, GK_mx ); 01890 01891 normalize_mx( G2_mx ); 01892 01893 ASSERT( l == lp - 1); 01894 ASSERT( n1 != 0. ); 01895 ASSERT( n2 != 0. ); 01896 ASSERT( dd1 != 0. ); 01897 ASSERT( dd2 != 0. ); 01898 ASSERT( dd3 != 0. ); 01899 ASSERT( Ksqrd != 0. ); 01900 01901 rcsvV_mxq[rindx].q = 1; 01902 rcsvV_mxq[rindx].mx = G2_mx; 01903 return G2_mx; 01904 } 01905 else 01906 { 01907 /******************************************************************************/ 01908 /* G(n,l; K,l+1) = [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l+2) */ 01909 /* */ 01910 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */ 01911 /* */ 01912 /* */ 01913 /* FROM Eq. 3.23 */ 01914 /* */ 01915 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2 + (l+2)(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */ 01916 /* */ 01917 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */ 01918 /******************************************************************************/ 01919 01920 long int lp2 = (q + 2); 01921 long int lp3 = (q + 3); 01922 01923 double lp2s = (double)(lp2 * lp2); 01924 double lp3s = (double)(lp3 * lp3); 01925 double n2 = (double)(n*n); 01926 double Ksqrd = (K * K); 01927 01928 /******************************************************/ 01929 /* [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] */ 01930 /******************************************************/ 01931 01932 double d1 = (4. * n2); 01933 double d2 = (4. * lp2s); 01934 double d3 = (double)(lp2)*((2*q)+3); 01935 double d4 = (1. + (n2 * Ksqrd)); 01936 /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */ 01937 double d5 = d1 - d2 + (d3 * d4); 01938 01939 mx d5_mx=mxify(d5); 01940 01941 /******************************************************/ 01942 /* 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] */ 01943 /******************************************************/ 01944 01945 double d6 = (n2 - lp2s); 01946 double d7 = (1. + (lp3s * Ksqrd)); 01947 /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */ 01948 double d8 = d1 * d6 * d7; 01949 01950 mx d8_mx = mxify(d8); 01951 01952 /******************************************************************************/ 01953 /* G(n,l; K,l+1) = [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] G(n,l+1; K,l+2) */ 01954 /* */ 01955 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */ 01956 /******************************************************************************/ 01957 01958 mx d9_mx = bhGm_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx ); 01959 mx d10_mx = bhGm_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx ); 01960 mx d11_mx = mult_mx( d5_mx, d9_mx ); 01961 mx d12_mx = mult_mx( d8_mx, d10_mx); 01962 mx result_mx = sub_mx( d11_mx , d12_mx ); 01963 rcsvV_mxq[rindx].q = 1; 01964 rcsvV_mxq[rindx].mx = result_mx; 01965 01966 ASSERT(l == lp - 1); 01967 ASSERT(n2 != 0. ); 01968 ASSERT(lp2s != 0.); 01969 ASSERT( lp3s != 0.); 01970 ASSERT(Ksqrd != 0.); 01971 01972 ASSERT(d1 != 0.); 01973 ASSERT(d2 != 0.); 01974 ASSERT(d3 != 0.); 01975 ASSERT(d4 != 0.); 01976 ASSERT(d5 != 0.); 01977 ASSERT(d6 != 0.); 01978 ASSERT(d7 != 0.); 01979 ASSERT(d8 != 0.); 01980 return result_mx; 01981 } 01982 } 01983 else 01984 { 01985 ASSERT( rcsvV_mxq[rindx].q != 0 ); 01986 return rcsvV_mxq[rindx].mx; 01987 } 01988 } 01989 01990 /****************************************************************************************/ 01991 /* */ 01992 /* bhg.c */ 01993 /* */ 01994 /* From reference; */ 01995 /* M. Brocklehurst */ 01996 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 01997 /* */ 01998 /* */ 01999 /* We wish to compute the following function, */ 02000 /* */ 02001 /* equation: (3.17) */ 02002 /* - s=l' - (1/2) */ 02003 /* | (n+l)! ----- | */ 02004 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */ 02005 /* | (n-l-1)! | | | */ 02006 /* - s=0 - */ 02007 /* */ 02008 /* Using various recursion relations (for l'=l+1) */ 02009 /* */ 02010 /* equation: (3.23) */ 02011 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2+l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */ 02012 /* */ 02013 /* - 4n^2 (n^2-l^2)[1+(l+1)^2 K^2] G(n,l; K,l+1) */ 02014 /* */ 02015 /* and (for l'=l-1) */ 02016 /* */ 02017 /* equation: (3.24) */ 02018 /* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l; K,l-1) */ 02019 /* */ 02020 /* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l; K,l+1) */ 02021 /* */ 02022 /* */ 02023 /* the starting point for the recursion relations are: */ 02024 /* */ 02025 /* */ 02026 /* equation (3.18): */ 02027 /* */ 02028 /* | pi |(1/2) 8n */ 02029 /* G(n,n-1; 0,n) = | -- | ------- (4n)^2 exp(-2n) */ 02030 /* | 2 | (2n-1)! */ 02031 /* */ 02032 /* equation (3.20): */ 02033 /* */ 02034 /* exp(2n-2/K tan^(-1)(n K) */ 02035 /* G(n,n-1; K,n) = --------------------------------------- */ 02036 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K)^(n+2) */ 02037 /* */ 02038 /* */ 02039 /* */ 02040 /* equation (3.20): */ 02041 /* G(n,n-2; K,n-1) = (2n-2)(1+(n K)^2) n G(n,n-1; K,n) */ 02042 /* */ 02043 /* */ 02044 /* equation (3.21): */ 02045 /* */ 02046 /* (1+(n K)^2) */ 02047 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 02048 /* 2n */ 02049 /****************************************************************************************/ 02050 02051 STATIC double bhg( 02052 double K, 02053 long int n, 02054 long int l, 02055 long int lp, 02056 /* Temporary storage for intermediate */ 02057 /* results of the recursive routine */ 02058 double *rcsvV 02059 ) 02060 { 02061 DEBUG_ENTRY( "bhg()" ); 02062 02063 double ld1 = factorial( n + l ); 02064 double ld2 = factorial( n - l - 1 ); 02065 double ld3 = (ld1 / ld2); 02066 02067 double partprod = local_product( K , lp ); 02068 02069 /**************************************************************************************/ 02070 /* equation: (3.17) */ 02071 /* - s=l' - (1/2) */ 02072 /* | (n+l)! ----- | */ 02073 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */ 02074 /* | (n-l-1)! | | | */ 02075 /* - s=0 - */ 02076 /**************************************************************************************/ 02077 02078 /**********************************************/ 02079 /* - s=l' - (1/2) */ 02080 /* | (n+l)! ----- | */ 02081 /* | -------- | | (1 + s^2 K^2) | */ 02082 /* | (n-l-1)! | | | */ 02083 /* - s=0 - */ 02084 /**********************************************/ 02085 02086 double d2 = sqrt( ld3 * partprod ); 02087 double d3 = powi( (2 * n) , (l - n) ); 02088 double d4 = bhG( K, n, l, lp, rcsvV ); 02089 double d5 = (d2 * d3); 02090 double d6 = (d5 * d4); 02091 02092 ASSERT(K != 0.); 02093 ASSERT( (n+l) >= 1 ); 02094 ASSERT( ((n-l)-1) >= 0 ); 02095 02096 ASSERT( partprod != 0. ); 02097 02098 ASSERT( ld1 != 0. ); 02099 ASSERT( ld2 != 0. ); 02100 ASSERT( ld3 != 0. ); 02101 02102 ASSERT( d2 != 0. ); 02103 ASSERT( d3 != 0. ); 02104 ASSERT( d4 != 0. ); 02105 ASSERT( d5 != 0. ); 02106 ASSERT( d6 != 0. ); 02107 return d6; 02108 } 02109 02110 /********************log version**************************/ 02111 STATIC double bhg_log( 02112 double K, 02113 long int n, 02114 long int l, 02115 long int lp, 02116 /* Temporary storage for intermediate */ 02117 /* results of the recursive routine */ 02118 mxq *rcsvV_mxq 02119 ) 02120 { 02121 /**************************************************************************************/ 02122 /* equation: (3.17) */ 02123 /* - s=l' - (1/2) */ 02124 /* | (n+l)! ----- | */ 02125 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */ 02126 /* | (n-l-1)! | | | */ 02127 /* - s=0 - */ 02128 /**************************************************************************************/ 02129 02130 DEBUG_ENTRY( "bhg_log()" ); 02131 02132 double d1 = (double)(2*n); 02133 double d2 = (double)(l-n); 02134 double Ksqrd = (K*K); 02135 02136 /**************************************************************************************/ 02137 /* */ 02138 /* | (n+l)! | */ 02139 /* log10 | -------- | = log10((n+1)!) - log10((n-l-1)!) */ 02140 /* | (n-l-1)! | */ 02141 /* */ 02142 /**************************************************************************************/ 02143 02144 double ld1 = lfactorial( n + l ); 02145 double ld2 = lfactorial( n - l - 1 ); 02146 02147 /**********************************************************************/ 02148 /* | s=l' | s=l' */ 02149 /* | ----- | --- */ 02150 /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */ 02151 /* | | | | --- */ 02152 /* | s=0 | s=0 */ 02153 /**********************************************************************/ 02154 02155 double ld3 = log10_prodxx( lp, Ksqrd ); 02156 02157 /**********************************************/ 02158 /* - s=l' - (1/2) */ 02159 /* | (n+l)! ----- | */ 02160 /* | -------- | | (1 + s^2 K^2) | */ 02161 /* | (n-l-1)! | | | */ 02162 /* - s=0 - */ 02163 /**********************************************/ 02164 02165 /***********************************************************************/ 02166 /* */ 02167 /* | - s=l' - (1/2) | */ 02168 /* | | (n+l)! ----- | | */ 02169 /* log10| | -------- | | (1 + s^2 K^2) | | == */ 02170 /* | | (n-l-1)! | | | | */ 02171 /* | - s=0 - | */ 02172 /* */ 02173 /* | | s=l' | | */ 02174 /* | | (n+l)! | | ----- | | */ 02175 /* (1/2)* |log10 | -------- | + log10 | | | (1 + s^2 K^2) | | */ 02176 /* | | (n-l-1)! | | | | | | */ 02177 /* | | s=0 | | */ 02178 /* */ 02179 /***********************************************************************/ 02180 02181 double ld4 = (1./2.) * ( ld3 + ld1 - ld2 ); 02182 02183 /**********************************************/ 02184 /* (2n)^(l-n) */ 02185 /**********************************************/ 02186 /* log10( 2n^(L-n) ) = (L-n) log10( 2n ) */ 02187 /**********************************************/ 02188 02189 double ld5 = d2 * log10( d1 ); 02190 02191 /**************************************************************************************/ 02192 /* equation: (3.17) */ 02193 /* - s=l' - (1/2) */ 02194 /* | (n+l)! ----- | */ 02195 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) * G(n,l; K,l') */ 02196 /* | (n-l-1)! | | | */ 02197 /* - s=0 - */ 02198 /**************************************************************************************/ 02199 02200 /****************************************************/ 02201 /* */ 02202 /* - s=l' - (1/2) */ 02203 /* | (n+l)! ----- | */ 02204 /* | -------- | | (1 + s^2 K^2) | * (2n)^(L-n) */ 02205 /* | (n-l-1)! | | | */ 02206 /* - s=0 - */ 02207 /****************************************************/ 02208 02209 double ld6 = (ld5+ld4); 02210 02211 mx d6_mx = mxify_log10( ld6 ); 02212 mx dd1_mx = bhG_mx( K, n, l, lp, rcsvV_mxq ); 02213 mx dd2_mx = mult_mx( d6_mx, dd1_mx ); 02214 normalize_mx( dd2_mx ); 02215 double result = unmxify( dd2_mx ); 02216 02217 ASSERT( result != 0. ); 02218 02219 ASSERT( Ksqrd != 0. ); 02220 ASSERT( ld3 >= 0. ); 02221 02222 ASSERT( d1 > 0. ); 02223 ASSERT( d2 < 0. ); 02224 return result; 02225 } 02226 02227 inline double local_product( double K , long int lp ) 02228 { 02229 long int s = 0; 02230 02231 double Ksqrd =(K*K); 02232 double partprod = 1.; 02233 02234 for( s = 0; s <= lp; s = s + 1 ) 02235 { 02236 double s2 = (double)(s*s); 02237 02238 /**************************/ 02239 /* s=l' */ 02240 /* ----- */ 02241 /* | | (1 + s^2 K^2) */ 02242 /* | | */ 02243 /* s=0 */ 02244 /**************************/ 02245 02246 partprod *= ( 1. + ( s2 * Ksqrd ) ); 02247 } 02248 return partprod; 02249 } 02250 02251 /************************************************************************/ 02252 /* Find the Einstein A's for hydrogen for a */ 02253 /* transition n,l --> n',l' */ 02254 /* */ 02255 /* In the following, we will assume n > n' */ 02256 /************************************************************************/ 02257 /*******************************************************************************/ 02258 /* */ 02259 /* Einstein A() for the transition from the */ 02260 /* initial state n,l to the finial state n',l' */ 02261 /* is given by oscillator f() */ 02262 /* */ 02263 /* hbar w max(l,l') | | 2 */ 02264 /* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */ 02265 /* 3 R_oo ( 2l + 1 ) | | */ 02266 /* */ 02267 /* */ 02268 /* E(n,l;n',l') max(l,l') | | 2 */ 02269 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 02270 /* 3 R_oo ( 2l + 1 ) | | */ 02271 /* */ 02272 /* */ 02273 /* See for example Gordan Drake's */ 02274 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 02275 /* */ 02276 /* Here R_oo is the infinite mass Rydberg length */ 02277 /* */ 02278 /* */ 02279 /* h c */ 02280 /* R_oo --- = 13.605698 eV */ 02281 /* {e} */ 02282 /* */ 02283 /* */ 02284 /* R_oo = 2.179874e-11 ergs */ 02285 /* */ 02286 /* w = omega */ 02287 /* = frequency of transition from n,l to n',l' */ 02288 /* */ 02289 /* */ 02290 /* */ 02291 /* here g_k are statistical weights obtained from */ 02292 /* the appropriate angular momentum quantum numbers */ 02293 /* */ 02294 /* */ 02295 /* - - 2 */ 02296 /* 64 pi^4 (e a_o)^2 max(l,l') | | */ 02297 /* A(n,l;n',l') = ------------------- ----------- v^3 | R(n,l;n',l') | */ 02298 /* 3 h c^3 2*l + 1 | | */ 02299 /* - - */ 02300 /* */ 02301 /* */ 02302 /* pi 3.141592654 */ 02303 /* plank_hbar 6.5821220 eV sec */ 02304 /* R_oo 2.179874e-11 ergs */ 02305 /* plank_h 6.6260755e-34 J sec */ 02306 /* e_charge 1.60217733e-19 C */ 02307 /* a_o 0.529177249e-10 m */ 02308 /* vel_light_c 299792458L m sec^-1 */ 02309 /* */ 02310 /* */ 02311 /* */ 02312 /* 64 pi^4 (e a_o)^2 64 pi^4 (a_o)^2 e^2 1 1 */ 02313 /* ----------------- = ----------------- -------- ---- = 7.5197711e-38 ----- */ 02314 /* 3 h c^3 3 c^2 hbar c 2 pi sec */ 02315 /* */ 02316 /* */ 02317 /* e^2 1 */ 02318 /* using ---------- = alpha = ---- */ 02319 /* hbar c 137 */ 02320 /*******************************************************************************/ 02321 02322 double H_Einstein_A(/* IN THE FOLLOWING WE HAVE n > n' */ 02323 /* principal quantum number, 1 for ground, upper level */ 02324 long int n, 02325 /* angular momentum, 0 for s */ 02326 long int l, 02327 /* principal quantum number, 1 for ground, lower level */ 02328 long int np, 02329 /* angular momentum, 0 for s */ 02330 long int lp, 02331 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 02332 long int iz 02333 ) 02334 { 02335 DEBUG_ENTRY( "H_Einstein_A()" ); 02336 02337 double result; 02338 if( n > 60 || np > 60 ) 02339 { 02340 result = H_Einstein_A_log10(n,l,np,lp,iz ); 02341 } 02342 else 02343 { 02344 result = H_Einstein_A_lin(n,l,np,lp,iz ); 02345 } 02346 return result; 02347 } 02348 02349 /************************************************************************/ 02350 /* Calculates the Einstein A's for hydrogen */ 02351 /* for the transition n,l --> n',l' */ 02352 /* units of sec^(-1) */ 02353 /* */ 02354 /* In the following, we have n > n' */ 02355 /************************************************************************/ 02356 STATIC double H_Einstein_A_lin(/* IN THE FOLLOWING WE HAVE n > n' */ 02357 /* principal quantum number, 1 for ground, upper level */ 02358 long int n, 02359 /* angular momentum, 0 for s */ 02360 long int l, 02361 /* principal quantum number, 1 for ground, lower level */ 02362 long int np, 02363 /* angular momentum, 0 for s */ 02364 long int lp, 02365 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 02366 long int iz 02367 ) 02368 { 02369 DEBUG_ENTRY( "H_Einstein_A_lin()" ); 02370 02371 /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */ 02372 double d1 = hv( n, np, iz ); 02373 double d2 = d1 / HPLANCK; /* v = hv / h */ 02374 double d3 = pow3(d2); 02375 double lg = (double)(l > lp ? l : lp); 02376 double Two_L_Plus_One = (double)(2*l + 1); 02377 double d6 = lg / Two_L_Plus_One; 02378 double d7 = hri( n, l, np, lp , iz ); 02379 double d8 = d7 * d7; 02380 double result = CONST_ONE * d3 * d6 * d8; 02381 02382 /* validate the incoming data */ 02383 if( n >= 70 ) 02384 { 02385 fprintf(ioQQQ,"Principal Quantum Number `n' too large.\n"); 02386 cdEXIT(EXIT_FAILURE); 02387 } 02388 if( iz < 1 ) 02389 { 02390 fprintf(ioQQQ," The charge is impossible.\n"); 02391 cdEXIT(EXIT_FAILURE); 02392 } 02393 if( n < 1 || np < 1 || l >= n || lp >= np ) 02394 { 02395 fprintf(ioQQQ," The quantum numbers are impossible.\n"); 02396 cdEXIT(EXIT_FAILURE); 02397 } 02398 if( n <= np ) 02399 { 02400 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n"); 02401 cdEXIT(EXIT_FAILURE); 02402 } 02403 return result; 02404 } 02405 02406 /**********************log version****************************/ 02407 double H_Einstein_A_log10(/* returns Einstein A in units of (sec)^-1 */ 02408 long int n, 02409 long int l, 02410 long int np, 02411 long int lp, 02412 long int iz 02413 ) 02414 { 02415 DEBUG_ENTRY( "H_Einstein_A_log10()" ); 02416 02417 /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */ 02418 double d1 = hv( n, np , iz ); 02419 double d2 = d1 / HPLANCK; /* v = hv / h */ 02420 double d3 = pow3(d2); 02421 double lg = (double)(l > lp ? l : lp); 02422 double Two_L_Plus_One = (double)(2*l + 1); 02423 double d6 = lg / Two_L_Plus_One; 02424 double d7 = hri_log10( n, l, np, lp , iz ); 02425 double d8 = d7 * d7; 02426 double result = CONST_ONE * d3 * d6 * d8; 02427 02428 /* validate the incoming data */ 02429 if( iz < 1 ) 02430 { 02431 fprintf(ioQQQ," The charge is impossible.\n"); 02432 cdEXIT(EXIT_FAILURE); 02433 } 02434 if( n < 1 || np < 1 || l >= n || lp >= np ) 02435 { 02436 fprintf(ioQQQ," The quantum numbers are impossible.\n"); 02437 cdEXIT(EXIT_FAILURE); 02438 } 02439 if( n <= np ) 02440 { 02441 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n"); 02442 cdEXIT(EXIT_FAILURE); 02443 } 02444 return result; 02445 } 02446 02447 /********************************************************************************/ 02448 /* hv calculates photon energy for n -> n' transitions for H and H-like ions */ 02449 /* simplest case of no "l" or "m" dependence */ 02450 /* epsilon_0 = 1 in vacu */ 02451 /* */ 02452 /* */ 02453 /* R_h */ 02454 /* Energy(n,Z) = - ------- */ 02455 /* n^2 */ 02456 /* */ 02457 /* */ 02458 /* */ 02459 /* Friedrich -- Theoretical Atomic Physics pg. 60 eq. 2.8 */ 02460 /* */ 02461 /* u */ 02462 /* R_h = --- R_oo where */ 02463 /* m_e */ 02464 /* */ 02465 /* h c */ 02466 /* R_oo --- = 2.179874e-11 ergs */ 02467 /* e */ 02468 /* */ 02469 /* (Harmin Lecture Notes for course phy-651 Spring 1994) */ 02470 /* where m_e (m_p) is the mass of and electron (proton) */ 02471 /* and u is the reduced electron mass for neutral hydrogen */ 02472 /* */ 02473 /* */ 02474 /* m_e m_p m_e */ 02475 /* u = --------- = ----------- */ 02476 /* m_e + m_p 1 + m_e/m_p */ 02477 /* */ 02478 /* m_e */ 02479 /* Now ----- = 0.000544617013 */ 02480 /* m_p */ 02481 /* u */ 02482 /* so that --- = 0.999455679 */ 02483 /* m_e */ 02484 /* */ 02485 /* */ 02486 /* returns energy of photon in ergs */ 02487 /* */ 02488 /* hv (n,n',Z) is for transitions n -> n' */ 02489 /* */ 02490 /* 1 erg = 1e-07 J */ 02491 /********************************************************************************/ 02492 /********************************************************************************/ 02493 /* WARNING: hv() use the electron reduced mass for hydrogen instead of */ 02494 /* the reduced mass associated with the apropriate ion */ 02495 /********************************************************************************/ 02496 02497 inline double hv( long int n, long int nprime, long int iz ) 02498 { 02499 DEBUG_ENTRY( "hv()" ); 02500 02501 double n1 = (double)n; 02502 double n2 = n1*n1; 02503 double np1 = (double)nprime; 02504 double np2 = np1*np1; 02505 double rmr = 1./(1. + ELECTRON_MASS/PROTON_MASS); /* 0.999455679 */ 02506 double izsqrd = (double)(iz*iz); 02507 02508 double d1 = 1. / n2; 02509 double d2 = 1. / np2; 02510 double d3 = izsqrd * rmr * EN1RYD; 02511 double d4 = d2 - d1; 02512 double result = d3 * d4; 02513 02514 ASSERT( n > 0 ); 02515 ASSERT( nprime > 0 ); 02516 ASSERT( n > nprime ); 02517 ASSERT( iz > 0 ); 02518 ASSERT( result > 0. ); 02519 02520 if( n <= nprime ) 02521 { 02522 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n"); 02523 cdEXIT(EXIT_FAILURE); 02524 } 02525 02526 return result; 02527 } 02528 02529 /************************************************************************/ 02530 /* hri() */ 02531 /* Calculate the hydrogen radial wavefunction integral */ 02532 /* for the dipole transition l'=l-1 or l'=l+1 */ 02533 /* for the higher energy state n,l to the lower energy state n',l' */ 02534 /* no "m" dependence */ 02535 /************************************************************************/ 02536 /* here we have a transition */ 02537 /* from the higher energy state n,l */ 02538 /* to the lower energy state n',l' */ 02539 /* with a dipole selection rule on l and l' */ 02540 /************************************************************************/ 02541 /* */ 02542 /* hri() test n,l,n',l' for domain errors and */ 02543 /* swaps n,l <--> n',l' for the case l'=l+1 */ 02544 /* */ 02545 /* It then calls hrii() */ 02546 /* */ 02547 /* Dec. 6, 1999 */ 02548 /* Robert Paul Bauman */ 02549 /************************************************************************/ 02550 02551 /************************************************************************/ 02552 /* This routine, hri(), calculates the hydrogen radial integral, */ 02553 /* for the transition n,l --> n',l' */ 02554 /* It is, of course, dimensionless. */ 02555 /* */ 02556 /* In the following, we have n > n' */ 02557 /************************************************************************/ 02558 02559 inline double hri( 02560 /* principal quantum number, 1 for ground, upper level */ 02561 long int n, 02562 /* angular momentum, 0 for s */ 02563 long int l, 02564 /* principal quantum number, 1 for ground, lower level */ 02565 long int np, 02566 /* angular momentum, 0 for s */ 02567 long int lp, 02568 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 02569 long int iz 02570 ) 02571 { 02572 DEBUG_ENTRY( "hri" ); 02573 02574 long int a; 02575 long int b; 02576 long int c; 02577 long int d; 02578 double ld1 = 0.; 02579 double Z = (double)iz; 02580 02581 /**********************************************************************/ 02582 /* from higher energy -> lower energy */ 02583 /* Selection Rule for l and l' */ 02584 /* dipole process only */ 02585 /**********************************************************************/ 02586 02587 ASSERT( n > 0 ); 02588 ASSERT( np > 0 ); 02589 ASSERT( l >= 0 ); 02590 ASSERT( lp >= 0 ); 02591 ASSERT( n > l ); 02592 ASSERT( np > lp ); 02593 ASSERT( n > np || ( n == np && l == lp + 1 )); 02594 ASSERT( iz > 0 ); 02595 ASSERT( lp == l + 1 || lp == l - 1 ); 02596 02597 if( l == lp + 1 ) 02598 { 02599 /* Keep variable the same */ 02600 a = n; 02601 b = l; 02602 c = np; 02603 d = lp; 02604 } 02605 else if( l == lp - 1 ) 02606 { 02607 /* swap n,l with n',l' */ 02608 a = np; 02609 b = lp; 02610 c = n; 02611 d = l; 02612 } 02613 else 02614 { 02615 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" ); 02616 cdEXIT(EXIT_FAILURE); 02617 } 02618 02619 /**********************************************/ 02620 /* Take care of the Z-dependence here. */ 02621 /**********************************************/ 02622 ld1 = hrii(a, b, c, d ) / Z; 02623 02624 return ld1; 02625 } 02626 02627 /************************************************************************/ 02628 /* hri_log10() */ 02629 /* Calculate the hydrogen radial wavefunction integral */ 02630 /* for the dipole transition l'=l-1 or l'=l+1 */ 02631 /* for the higher energy state n,l to the lower energy state n',l' */ 02632 /* no "m" dependence */ 02633 /************************************************************************/ 02634 /* here we have a transition */ 02635 /* from the higher energy state n,l */ 02636 /* to the lower energy state n',l' */ 02637 /* with a dipole selection rule on l and l' */ 02638 /************************************************************************/ 02639 /* */ 02640 /* hri_log10() test n,l,n',l' for domain errors and */ 02641 /* swaps n,l <--> n',l' for the case l'=l+1 */ 02642 /* */ 02643 /* It then calls hrii_log() */ 02644 /* */ 02645 /* Dec. 6, 1999 */ 02646 /* Robert Paul Bauman */ 02647 /************************************************************************/ 02648 02649 inline double hri_log10( long int n, long int l, long int np, long int lp , long int iz ) 02650 { 02651 /**********************************************************************/ 02652 /* from higher energy -> lower energy */ 02653 /* Selection Rule for l and l' */ 02654 /* dipole process only */ 02655 /**********************************************************************/ 02656 02657 DEBUG_ENTRY( "hri_log10()" ); 02658 02659 long int a; 02660 long int b; 02661 long int c; 02662 long int d; 02663 double ld1 = 0.; 02664 double Z = (double)iz; 02665 02666 ASSERT( n > 0); 02667 ASSERT( np > 0); 02668 ASSERT( l >= 0); 02669 ASSERT( lp >= 0 ); 02670 ASSERT( n > l ); 02671 ASSERT( np > lp ); 02672 ASSERT( n > np || ( n == np && l == lp + 1 )); 02673 ASSERT( iz > 0 ); 02674 ASSERT( lp == l + 1 || lp == l - 1 ); 02675 02676 if( l == lp + 1) 02677 { 02678 /* Keep variable the same */ 02679 a = n; 02680 b = l; 02681 c = np; 02682 d = lp; 02683 } 02684 else if( l == lp - 1 ) 02685 { 02686 /* swap n,l with n',l' */ 02687 a = np; 02688 b = lp; 02689 c = n; 02690 d = l; 02691 } 02692 else 02693 { 02694 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" ); 02695 cdEXIT(EXIT_FAILURE); 02696 } 02697 02698 /**********************************************/ 02699 /* Take care of the Z-dependence here. */ 02700 /**********************************************/ 02701 ld1 = hrii_log(a, b, c, d ) / Z; 02702 02703 return ld1; 02704 } 02705 02706 STATIC double hrii( long int n, long int l, long int np, long int lp) 02707 { 02708 /******************************************************************************/ 02709 /* this routine hrii() is internal to the parent routine hri() */ 02710 /* this internal routine only considers the case l=l'+1 */ 02711 /* the case l=l-1 is done in the parent routine hri() */ 02712 /* by the transformation n <--> n' and l <--> l' */ 02713 /* THUS WE TEST FOR */ 02714 /* l=l'-1 */ 02715 /******************************************************************************/ 02716 02717 DEBUG_ENTRY( "hrii()" ); 02718 02719 long int a = 0, b = 0, c = 0; 02720 long int i1 = 0, i2 = 0, i3 = 0, i4 = 0; 02721 02722 char A='a'; 02723 02724 double y = 0.; 02725 double fsf = 0.; 02726 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.; 02727 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.; 02728 double d00 = 0., d01 = 0.; 02729 02730 ASSERT( l == lp + 1 ); 02731 02732 if( n == np ) /* SPECIAL CASE 1 */ 02733 { 02734 /**********************************************************/ 02735 /* if lp= l + 1 then it has higher energy */ 02736 /* i.e. no photon */ 02737 /* this is the second time we check this, oh well */ 02738 /**********************************************************/ 02739 02740 if( lp != (l - 1) ) 02741 { 02742 printf( "BadMagic: Energy requirements not met.\n\n" ); 02743 cdEXIT(EXIT_FAILURE); 02744 } 02745 02746 d2 = 3. / 2.; 02747 i1 = n * n; 02748 i2 = l * l; 02749 d5 = (double)(i1 - i2); 02750 d6 = sqrt(d5); 02751 d7 = (double)n * d6; 02752 d8 = d2 * d7; 02753 return d8; 02754 } 02755 else if( l == np && lp == (l - 1) ) /* A Pair of Easy Special Cases */ 02756 { 02757 if( l == (n - 1) ) 02758 { 02759 /**********************************************************************/ 02760 /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */ 02761 /* */ 02762 /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */ 02763 /* [(2n-1) - 1/(2n-1)]/4 */ 02764 /**********************************************************************/ 02765 02766 d1 = (double)( 2*n - 2 ); 02767 d2 = (double)( 2*n - 1 ); 02768 d3 = d1 * d2; 02769 d4 = sqrt( d3 ); 02770 02771 d5 = (double)(4 * n * (n - 1)); 02772 i1 = (2*n - 1); 02773 d6 = (double)(i1 * i1); 02774 d7 = d5/ d6; 02775 d8 = powi( d7, n ); 02776 02777 d9 = 1./d2; 02778 d10 = d2 - d9; 02779 d11 = d10 / 4.; 02780 02781 /* Wrap it all up */ 02782 02783 d12 = d4 * d8 * d11; 02784 return d12; 02785 02786 } 02787 else 02788 { 02789 /******************************************************************************/ 02790 /* R(n,l;n',l') = R(n,l;l,l-1) */ 02791 /* */ 02792 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */ 02793 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */ 02794 /******************************************************************************/ 02795 02796 d2 = 1.; 02797 for( i1 = -l; i1 <= l; i1 = i1 + 1 ) /* from n-l to n+l INCLUSIVE */ 02798 { 02799 d1 = (double)(n - i1); 02800 d2 = d2 * d1; 02801 } 02802 i2 = (2*l - 1); 02803 d3 = factorial( i2 ); 02804 d4 = d2/d3; 02805 d4 = sqrt( d4 ); 02806 02807 02808 d5 = (double)( 4. * n * l ); 02809 i3 = (n - l); 02810 d6 = (double)( i3 * i3 ); 02811 d7 = d5 / d6; 02812 d8 = powi( d7, l+1 ); 02813 02814 02815 i4 = n + l; 02816 d9 = (double)( i3 ) / (double)( i4 ); 02817 d10 = powi( d9 , i4 ); 02818 02819 d11 = d9 * d9; 02820 d12 = 1. - d11; 02821 d13 = d12 / 4.; 02822 02823 /* Wrap it all up */ 02824 d14 = d4 * d8 * d10 * d13; 02825 return d14; 02826 } 02827 } 02828 02829 /*******************************************************************************************/ 02830 /* THE GENERAL CASE */ 02831 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 02832 /* REF: D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 02833 /* For F(a,b;c;x) we have from eq.4 */ 02834 /* */ 02835 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 02836 /* */ 02837 /* a (1-x) (a + bx - c) */ 02838 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 02839 /* (a-c) (a-c) */ 02840 /* */ 02841 /* */ 02842 /* A similiar recusion relation holds for b with a <--> b. */ 02843 /* */ 02844 /* */ 02845 /* we have initial conditions */ 02846 /* */ 02847 /* */ 02848 /* F(0) = 1 with a = -1 */ 02849 /* */ 02850 /* b */ 02851 /* F(-1) = 1 - (---) x with a = -1 */ 02852 /* c */ 02853 /*******************************************************************************************/ 02854 02855 if( lp == l - 1 ) /* use recursion over "b" */ 02856 { 02857 A='b'; 02858 } 02859 else if( lp == l + 1 ) /* use recursion over "a" */ 02860 { 02861 A='a'; 02862 } 02863 else 02864 { 02865 printf(" BadMagic: Don't know what to do here.\n\n"); 02866 cdEXIT(EXIT_FAILURE); 02867 } 02868 02869 /********************************************************************/ 02870 /* Calculate the whole shootin match */ 02871 /* - - (1/2) */ 02872 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 02873 /* ----------- * | ---------------- | */ 02874 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 02875 /* - - */ 02876 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 02877 /* */ 02878 /* This is used in the calculation of hydrogen */ 02879 /* radial wave function integral for dipole transition case */ 02880 /********************************************************************/ 02881 02882 fsf = fsff( n, l, np ); 02883 02884 /**************************************************************************************/ 02885 /* Use a -> a' + 1 */ 02886 /* _ _ */ 02887 /* (a' + 1) (1 - x) | | */ 02888 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */ 02889 /* (a' + 1 -c) | | */ 02890 /* - - */ 02891 /* */ 02892 /* For the first F() in the solution of the radial integral */ 02893 /* */ 02894 /* a = ( -n + l + 1 ) */ 02895 /* */ 02896 /* a = -n + l + 1 */ 02897 /* max(a) = max(-n) + max(l) + 1 */ 02898 /* = -n + max(n-1) + 1 */ 02899 /* = -n + n-1 + 1 */ 02900 /* = 0 */ 02901 /* */ 02902 /* similiarly */ 02903 /* */ 02904 /* min(a) = min(-n) + min(l) + 1 */ 02905 /* = min(-n) + 0 + 1 */ 02906 /* = -n + 1 */ 02907 /* */ 02908 /* a -> a' + 1 implies */ 02909 /* */ 02910 /* max(a') = -1 */ 02911 /* min(a') = -n */ 02912 /**************************************************************************************/ 02913 02914 /* a plus */ 02915 a = (-n + l + 1); 02916 02917 /* for the first 2_F_1 we use b = (-n' + l) */ 02918 b = (-np + l); 02919 02920 /* c is simple */ 02921 c = 2 * l; 02922 02923 /* -4 nn' */ 02924 /* where Y = -------- . */ 02925 /* (n-n')^2 */ 02926 d2 = (double)(n - np); 02927 d3 = d2 * d2; 02928 d4 = 1. / d3; 02929 d5 = (double)(n * np); 02930 d6 = d5 * 4.; 02931 d7 = - d6; 02932 y = d7 * d4; 02933 02934 d00 = F21( a, b, c, y, A ); 02935 02936 /**************************************************************/ 02937 /* For the second F() in the solution of the radial integral */ 02938 /* */ 02939 /* a = ( -n + l - 1 ) */ 02940 /* */ 02941 /* a = -n + l + 1 */ 02942 /* max(a) = max(-n) + max(l) - 1 */ 02943 /* = -n + (n - 1) - 1 */ 02944 /* = -2 */ 02945 /* */ 02946 /* similiarly */ 02947 /* */ 02948 /* min(a) = min(-n) + min(l) - 1 */ 02949 /* = (-n) + 0 - 1 */ 02950 /* = -n - 1 */ 02951 /* */ 02952 /* a -> a' + 1 implies */ 02953 /* */ 02954 /* max(a') = -3 */ 02955 /* */ 02956 /* min(a') = -n - 2 */ 02957 /**************************************************************/ 02958 02959 /* a minus */ 02960 a = (-n + l - 1); 02961 02962 /* for the first 2_F_1 we use b = (-n' + l) */ 02963 /* and does not change */ 02964 b = (-np + l); 02965 02966 /* c is simple */ 02967 c = 2 * l; 02968 02969 /**************************************************************/ 02970 /* -4 nn' */ 02971 /* where Y = -------- . */ 02972 /* (n-n')^2 */ 02973 /**************************************************************/ 02974 02975 /**************************************************************/ 02976 /* These are already calculated a few lines up */ 02977 /* */ 02978 /* d2 = (double) (n - np); */ 02979 /* d3 = d2 * d2; */ 02980 /* d4 = 1/ d3; */ 02981 /* d5 = (double) (n * np); */ 02982 /* d6 = d5 * 4.0; */ 02983 /* d7 = - d6; */ 02984 /* y = d7 * d4; */ 02985 /**************************************************************/ 02986 02987 d01 = F21(a, b, c, y, A ); 02988 02989 /* Calculate */ 02990 /* */ 02991 /* (n-n')^2 */ 02992 /* -------- */ 02993 /* (n+n')^2 */ 02994 02995 i1 = (n - np); 02996 d1 = pow2( (double)i1 ); 02997 i2 = (n + np); 02998 d2 = pow2( (double)i2 ); 02999 d3 = d1 / d2; 03000 03001 d4 = d01 * d3; 03002 d5 = d00 - d4; 03003 d6 = fsf * d5; 03004 03005 ASSERT( d6 != 0. ); 03006 return d6; 03007 } 03008 03009 03010 STATIC double hrii_log( long int n, long int l, long int np, long int lp) 03011 { 03012 /******************************************************************************/ 03013 /* this routine hrii_log() is internal to the parent routine hri_log10() */ 03014 /* this internal routine only considers the case l=l'+1 */ 03015 /* the case l=l-1 is done in the parent routine hri_log10() */ 03016 /* by the transformation n <--> n' and l <--> l' */ 03017 /* THUS WE TEST FOR */ 03018 /* l=l'-1 */ 03019 /******************************************************************************/ 03020 /**************************************************************************************/ 03021 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */ 03022 /* */ 03023 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */ 03024 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */ 03025 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */ 03026 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03027 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */ 03028 /**************************************************************************************/ 03029 03030 DEBUG_ENTRY( "hrii_log()" ); 03031 03032 char A='a'; 03033 03034 double y = 0.; 03035 double log10_fsf = 0.; 03036 03037 ASSERT( l == lp + 1 ); 03038 03039 if( n == np ) /* SPECIAL CASE 1 */ 03040 { 03041 /**********************************************************/ 03042 /* if lp= l + 1 then it has higher energy */ 03043 /* i.e. no photon */ 03044 /* this is the second time we check this, oh well */ 03045 /**********************************************************/ 03046 03047 if( lp != (l - 1) ) 03048 { 03049 printf( "BadMagic: l'= l+1 for n'= n.\n\n" ); 03050 cdEXIT(EXIT_FAILURE); 03051 } 03052 else 03053 { 03054 /**********************************************************/ 03055 /* 3 */ 03056 /* R(nl:n'=n,l'=l+1) = --- n sqrt( n^2 - l^2 ) */ 03057 /* 2 */ 03058 /**********************************************************/ 03059 03060 long int i1 = n * n; 03061 long int i2 = l * l; 03062 03063 double d1 = 3. / 2.; 03064 double d2 = (double)n; 03065 double d3 = (double)(i1 - i2); 03066 double d4 = sqrt(d3); 03067 double result = d1 * d2 * d4; 03068 03069 ASSERT( d3 >= 0. ); 03070 return result; 03071 } 03072 } 03073 else if( l == np && lp == l - 1 ) /* A Pair of Easy Special Cases */ 03074 { 03075 if( l == n - 1 ) 03076 { 03077 /**********************************************************************/ 03078 /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */ 03079 /* */ 03080 /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */ 03081 /* [(2n-1) - 1/(2n-1)]/4 */ 03082 /**********************************************************************/ 03083 03084 double d1 = (double)((2*n-2)*(2*n-1)); 03085 double d2 = sqrt( d1 ); 03086 double d3 = (double)(4*n*(n-1)); 03087 double d4 = (double)(2*n-1); 03088 double d5 = d4*d4; 03089 double d7 = d3/d5; 03090 double d8 = powi(d7,n); 03091 double d9 = 1./d4; 03092 double d10 = d4 - d9; 03093 double d11 = 0.25*d10; 03094 double result = (d2 * d8 * d11); /* Wrap it all up */ 03095 03096 ASSERT( d1 >= 0. ); 03097 ASSERT( d3 >= 0. ); 03098 return result; 03099 } 03100 else 03101 { 03102 double result = 0.; 03103 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0., ld7 = 0.; 03104 03105 /******************************************************************************/ 03106 /* R(n,l;n',l') = R(n,l;l,l-1) */ 03107 /* */ 03108 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */ 03109 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */ 03110 /******************************************************************************/ 03111 /**************************************/ 03112 /* [(n-l) ... (n+l)] */ 03113 /**************************************/ 03114 /* log10[(n-l) ... (n+l)] = */ 03115 /* */ 03116 /* n+l */ 03117 /* --- */ 03118 /* > log10(j) */ 03119 /* --- */ 03120 /* j=n-l */ 03121 /**************************************/ 03122 03123 ld1 = 0.; 03124 for( long int i1 = (n-l); i1 <= (n+l); i1++ ) /* from n-l to n+l INCLUSIVE */ 03125 { 03126 double d1 = (double)(i1); 03127 ld1 += log10( d1 ); 03128 } 03129 03130 /**************************************/ 03131 /* (2l-1)! */ 03132 /**************************************/ 03133 /* log10[ (2n-1)! ] */ 03134 /**************************************/ 03135 03136 ld2 = lfactorial( 2*l - 1 ); 03137 03138 ASSERT( ((2*l)+1) >= 0); 03139 03140 /**********************************************/ 03141 /* log10( [(n-l) ... (n+l)/(2l-1)!]^(1/2) ) = */ 03142 /* (1/2) log10[(n-l) ... (n+l)] - */ 03143 /* (1/2) log10[ (2n-1)! ] */ 03144 /**********************************************/ 03145 03146 ld3 = 0.5 * (ld1 - ld2); 03147 03148 /**********************************************/ 03149 /* [4nl/(n-l)^2]^(l+1) */ 03150 /**********************************************/ 03151 /* log10( [4nl/(n-l)^2]^(l+1) ) = */ 03152 /* (l+1) * log10( [4nl/(n-l)^2] ) */ 03153 /* */ 03154 /* = (l+1)*[ log10(4nl) - 2 log10(n-l) ] */ 03155 /* */ 03156 /**********************************************/ 03157 03158 double d1 = (double)(l+1); 03159 double d2 = (double)(4*n*l); 03160 double d3 = (double)(n-l); 03161 double d4 = log10(d2); 03162 double d5 = log10(d3); 03163 03164 ld4 = d1 * (d4 - 2.*d5); 03165 03166 /**********************************************/ 03167 /* [(n-l)/(n+l)]^(n+l) */ 03168 /**********************************************/ 03169 /* log10( [ (n-l)/(n+l) ]^(n+l) ) = */ 03170 /* */ 03171 /* (n+l) * [ log10(n-l) - log10(n+l) ] */ 03172 /* */ 03173 /**********************************************/ 03174 03175 d1 = (double)(n-l); 03176 d2 = (double)(n+l); 03177 d3 = log10( d1 ); 03178 d4 = log10( d2 ); 03179 03180 ld5 = d2 * (d3 - d4); 03181 03182 /**********************************************/ 03183 /* {1-[(n-l)/(n+l)]^2}/4 */ 03184 /**********************************************/ 03185 /* log10[ {1-[(n-l)/(n+l)]^2}/4 ] */ 03186 /**********************************************/ 03187 03188 d1 = (double)(n-l); 03189 d2 = (double)(n+l); 03190 d3 = d1/d2; 03191 d4 = d3*d3; 03192 d5 = 1. - d4; 03193 double d6 = 0.25*d5; 03194 03195 ld6 = log10(d6); 03196 03197 /******************************************************************************/ 03198 /* R(n,l;n',l') = R(n,l;l,l-1) */ 03199 /* */ 03200 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */ 03201 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */ 03202 /******************************************************************************/ 03203 03204 ld7 = ld3 + ld4 + ld5 + ld6; 03205 03206 result = pow( 10., ld7 ); 03207 03208 ASSERT( result > 0. ); 03209 return result; 03210 } 03211 } 03212 else 03213 { 03214 double result = 0.; 03215 long int a = 0, b = 0, c = 0; 03216 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.; 03217 mx d00={0.0,0}, d01={0.0,0}, d02={0.0,0}, d03={0.0,0}; 03218 03219 if( lp == l - 1 ) /* use recursion over "b" */ 03220 { 03221 A='b'; 03222 } 03223 else if( lp == l + 1 ) /* use recursion over "a" */ 03224 { 03225 A='a'; 03226 } 03227 else 03228 { 03229 printf(" BadMagic: Don't know what to do here.\n\n"); 03230 cdEXIT(EXIT_FAILURE); 03231 } 03232 03233 /**************************************************************************************/ 03234 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */ 03235 /* */ 03236 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */ 03237 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */ 03238 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */ 03239 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03240 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */ 03241 /**************************************************************************************/ 03242 03243 /****************************************************************************************************/ 03244 /* Calculate the whole shootin match */ 03245 /* - - (1/2) */ 03246 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 03247 /* fsff() = ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2)(n+n')^(-n-n') */ 03248 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03249 /* - - */ 03250 /* This is used in the calculation of hydrogen radial wave function integral for dipole transitions */ 03251 /****************************************************************************************************/ 03252 03253 log10_fsf = log10_fsff( n, l, np ); 03254 03255 /******************************************************************************************/ 03256 /* 2_F_1( a, b; c; y ) */ 03257 /* */ 03258 /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */ 03259 /* */ 03260 /* */ 03261 /* Use a -> a' + 1 */ 03262 /* _ _ */ 03263 /* (a' + 1) (1 - x) | | */ 03264 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */ 03265 /* (a' + 1 - c) | | */ 03266 /* - - */ 03267 /* */ 03268 /* For the first F() in the solution of the radial integral */ 03269 /* */ 03270 /* a = ( -n + l + 1 ) */ 03271 /* */ 03272 /* a = -n + l + 1 */ 03273 /* max(a) = max(-n) + max(l) + 1 */ 03274 /* = -n + max(n-1) + 1 */ 03275 /* = -n + n-1 + 1 */ 03276 /* = 0 */ 03277 /* */ 03278 /* similiarly */ 03279 /* */ 03280 /* min(a) = min(-n) + min(l) + 1 */ 03281 /* = min(-n) + 0 + 1 */ 03282 /* = -n + 1 */ 03283 /* */ 03284 /* a -> a' + 1 implies */ 03285 /* */ 03286 /* max(a') = -1 */ 03287 /* min(a') = -n */ 03288 /******************************************************************************************/ 03289 03290 /* a plus */ 03291 a = (-n + l + 1); 03292 03293 /* for the first 2_F_1 we use b = (-n' + l) */ 03294 b = (-np + l); 03295 03296 /* c is simple */ 03297 c = 2 * l; 03298 03299 /**********************************************************************************/ 03300 /* 2_F_1( a, b; c; y ) */ 03301 /* */ 03302 /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */ 03303 /* */ 03304 /* -4 nn' */ 03305 /* where Y = -------- . */ 03306 /* (n-n')^2 */ 03307 /* */ 03308 /**********************************************************************************/ 03309 03310 d2 = (double)(n - np); 03311 d3 = d2 * d2; 03312 03313 d4 = 1. / d3; 03314 d5 = (double)(n * np); 03315 d6 = d5 * 4.; 03316 03317 d7 = -d6; 03318 y = d7 * d4; 03319 03320 /**************************************************************************************************/ 03321 /* THE GENERAL CASE */ 03322 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 03323 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 03324 /* */ 03325 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 03326 /* */ 03327 /* a (1-x) (a + bx - c) */ 03328 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 03329 /* (a - c) (a - c) */ 03330 /* */ 03331 /* */ 03332 /* A similiar recusion relation holds for b with a <--> b. */ 03333 /* */ 03334 /* */ 03335 /* we have initial conditions */ 03336 /* */ 03337 /* */ 03338 /* F(0) = 1 with a = -1 */ 03339 /* */ 03340 /* b */ 03341 /* F(-1) = 1 - (---) x with a = -1 */ 03342 /* c */ 03343 /**************************************************************************************************/ 03344 03345 /* 2_F_1( long int a, long int b, long int c, (double) y, (string) "a" or "b") */ 03346 /* F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03347 d00 = F21_mx( a, b, c, y, A ); 03348 03349 /**************************************************************/ 03350 /* For the second F() in the solution of the radial integral */ 03351 /* */ 03352 /* a = ( -n + l - 1 ) */ 03353 /* */ 03354 /* a = -n + l + 1 */ 03355 /* max(a) = max(-n) + max(l) - 1 */ 03356 /* = -n + (n - 1) - 1 */ 03357 /* = -2 */ 03358 /* */ 03359 /* similiarly */ 03360 /* */ 03361 /* min(a) = min(-n) + min(l) - 1 */ 03362 /* = (-n) + 0 - 1 */ 03363 /* = -n - 1 */ 03364 /* */ 03365 /* a -> a' + 1 implies */ 03366 /* */ 03367 /* max(a') = -3 */ 03368 /* */ 03369 /* min(a') = -n - 2 */ 03370 /**************************************************************/ 03371 03372 /* a minus */ 03373 a = (-n + l - 1); 03374 03375 /* for the first 2_F_1 we use b = (-n' + l) */ 03376 /* and does not change */ 03377 b = (-np + l); 03378 03379 /* c is simple */ 03380 c = 2 * l; 03381 03382 /**************************************************************/ 03383 /* -4 nn' */ 03384 /* where Y = -------- . */ 03385 /* (n-n')^2 */ 03386 /**************************************************************/ 03387 03388 /**************************************************************/ 03389 /* These are already calculated a few lines up */ 03390 /* */ 03391 /* d2 = (double) (n - np); */ 03392 /* d3 = d2 * d2; */ 03393 /* d4 = 1/ d3; */ 03394 /* d5 = (double) (n * np); */ 03395 /* d6 = d5 * 4.0; */ 03396 /* d7 = - d6; */ 03397 /* y = d7 * d4; */ 03398 /**************************************************************/ 03399 03400 d01 = F21_mx(a, b, c, y, A ); 03401 03402 /**************************************************************************************/ 03403 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */ 03404 /* */ 03405 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */ 03406 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */ 03407 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */ 03408 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03409 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */ 03410 /* */ 03411 /* = fsf * ( F(a,b,c;y) - d3 * F(a',b',c';y) ) */ 03412 /* */ 03413 /* where d3 = (n-n')^2 (n+n')^2 */ 03414 /* */ 03415 /**************************************************************************************/ 03416 03417 /**************************************************************/ 03418 /* Calculate */ 03419 /* */ 03420 /* (n-n')^2 */ 03421 /* -------- */ 03422 /* (n+n')^2 */ 03423 /**************************************************************/ 03424 03425 d1 = (double)((n - np)*(n -np)); 03426 d2 = (double)((n + np)*(n + np)); 03427 d3 = d1 / d2; 03428 03429 d02.x = d01.x; 03430 d02.m = d01.m * d3; 03431 03432 while ( fabs(d02.m) > 1.0e+25 ) 03433 { 03434 d02.m /= 1.0e+25; 03435 d02.x += 25; 03436 } 03437 03438 d03.x = d00.x; 03439 d03.m = d00.m * (1. - (d02.m/d00.m) * powi( 10. , (d02.x - d00.x) ) ); 03440 03441 result = pow( 10., (log10_fsf + d03.x) ) * d03.m; 03442 03443 ASSERT( result != 0. ); 03444 03445 /* we don't care about the correct sign of result... */ 03446 return fabs(result); 03447 } 03448 } 03449 03450 STATIC double fsff( long int n, long int l, long int np ) 03451 { 03452 /****************************************************************/ 03453 /* Calculate the whole shootin match */ 03454 /* - - (1/2) */ 03455 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 03456 /* ----------- * | ---------------- | */ 03457 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03458 /* - - */ 03459 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03460 /* */ 03461 /****************************************************************/ 03462 03463 DEBUG_ENTRY( "fsff()" ); 03464 03465 long int i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0; 03466 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.; 03467 double sigma = 1.; 03468 03469 /**************************************************************** 03470 * Calculate the whole shootin match * 03471 * (-1)^(n'-1) | (n+l)! (n'+l-1)! | * 03472 * ----------- * | ---------------- | * 03473 * [4 (2l-1)!] | (n-l-1)! (n'-l)! | * 03474 * - - * 03475 * * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * 03476 * * 03477 ****************************************************************/ 03478 03479 /* Calculate (-1)^(n'-l) */ 03480 if( is_odd(np - l) ) 03481 { 03482 sigma *= -1.; 03483 } 03484 ASSERT( sigma != 0. ); 03485 03486 /*********************/ 03487 /* Calculate (2l-1)! */ 03488 /*********************/ 03489 i1 = (2*l - 1); 03490 if( i1 < 0 ) 03491 { 03492 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03493 cdEXIT(EXIT_FAILURE); 03494 } 03495 03496 /****************************************************************/ 03497 /* Calculate the whole shootin match */ 03498 /* - - (1/2) */ 03499 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 03500 /* ----------- * | ---------------- | */ 03501 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03502 /* - - */ 03503 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03504 /* */ 03505 /****************************************************************/ 03506 03507 d0 = factorial( i1 ); 03508 d1 = 4. * d0; 03509 d2 = 1. / d1; 03510 03511 /**********************************************************************/ 03512 /* We want the (negitive) of this */ 03513 /* since we really are interested in */ 03514 /* [(2l-1)!]^-1 */ 03515 /**********************************************************************/ 03516 03517 sigma = sigma * d2; 03518 ASSERT( sigma != 0. ); 03519 03520 /**********************************************************************/ 03521 /* Calculate (4 n n')^(l+1) */ 03522 /* powi( m , n) calcs m^n */ 03523 /* returns long double with m,n ints */ 03524 /**********************************************************************/ 03525 03526 i0 = 4 * n * np; 03527 i1 = l + 1; 03528 d2 = powi( (double)i0 , i1 ); 03529 sigma = sigma * d2; 03530 ASSERT( sigma != 0. ); 03531 03532 /* Calculate (n-n')^(n+n'-2l-2) */ 03533 i0 = n - np; 03534 i1 = n + np - 2*l - 2; 03535 d2 = powi( (double)i0 , i1 ); 03536 sigma = sigma * d2; 03537 ASSERT( sigma != 0. ); 03538 03539 /* Calculate (n+n')^(-n-n') */ 03540 i0 = n + np; 03541 i1 = -n - np; 03542 d2 = powi( (double)i0 , i1 ); 03543 sigma = sigma * d2; 03544 ASSERT( sigma != 0. ); 03545 03546 /**********************************************************************/ 03547 /* - - (1/2) */ 03548 /* | (n+l)! (n'+l-1)! | */ 03549 /* Calculate | ---------------- | */ 03550 /* | (n-l-1)! (n'-l)! | */ 03551 /* - - */ 03552 /**********************************************************************/ 03553 03554 i1 = n + l; 03555 if( i1 < 0 ) 03556 { 03557 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03558 cdEXIT(EXIT_FAILURE); 03559 } 03560 d1 = factorial( i1 ); 03561 03562 i2 = np + l - 1; 03563 if( i2 < 0 ) 03564 { 03565 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03566 cdEXIT(EXIT_FAILURE); 03567 } 03568 d2 = factorial( i2 ); 03569 03570 i3 = n - l - 1; 03571 if( i3 < 0 ) 03572 { 03573 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03574 cdEXIT(EXIT_FAILURE); 03575 } 03576 d3 = factorial( i3 ); 03577 03578 i4 = np - l; 03579 if( i4 < 0 ) 03580 { 03581 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03582 cdEXIT(EXIT_FAILURE); 03583 } 03584 d4 = factorial( i4 ); 03585 03586 ASSERT( d3 != 0. ); 03587 ASSERT( d4 != 0. ); 03588 03589 /* Do this a different way to prevent overflow */ 03590 /* d5 = (sqrt(d1 *d2)); */ 03591 d5 = sqrt(d1)*sqrt(d2); 03592 d5 /= sqrt(d3); 03593 d5 /= sqrt(d4); 03594 03595 sigma = sigma * d5; 03596 03597 ASSERT( sigma != 0. ); 03598 return sigma; 03599 } 03600 03601 /**************************log version*******************************/ 03602 STATIC double log10_fsff( long int n, long int l, long int np ) 03603 { 03604 /******************************************************************************************************/ 03605 /* Calculate the whole shootin match */ 03606 /* - - (1/2) */ 03607 /* 1 | (n+l)! (n'+l-1)! | */ 03608 /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03609 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03610 /* - - */ 03611 /******************************************************************************************************/ 03612 03613 DEBUG_ENTRY( "log10_fsff()" ); 03614 03615 double d0 = 0., d1 = 0.; 03616 double ld0 = 0., ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0.; 03617 double result = 0.; 03618 03619 /******************************************************************************************************/ 03620 /* Calculate the log10 of the whole shootin match */ 03621 /* - - (1/2) */ 03622 /* 1 | (n+l)! (n'+l-1)! | */ 03623 /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03624 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03625 /* - - */ 03626 /******************************************************************************************************/ 03627 03628 /**********************/ 03629 /* Calculate (2l-1)! */ 03630 /**********************/ 03631 03632 d0 = (double)(2*l - 1); 03633 ASSERT( d0 != 0. ); 03634 03635 /******************************************************************************************************/ 03636 /* Calculate the whole shootin match */ 03637 /* - - (1/2) */ 03638 /* 1 | (n+l)! (n'+l-1)! | */ 03639 /* ----------- * | ---------------- | * (4 n n')^(l+1) |(n-n')^(n+n'-2l-2)| (n+n')^(-n-n') */ 03640 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03641 /* - - */ 03642 /******************************************************************************************************/ 03643 03644 ld0 = lfactorial( 2*l - 1 ); 03645 ld1 = log10(4.); 03646 result = -(ld0 + ld1); 03647 ASSERT( result != 0. ); 03648 03649 /**********************************************************************/ 03650 /* Calculate (4 n n')^(l+1) */ 03651 /* powi( m , n) calcs m^n */ 03652 /* returns long double with m,n ints */ 03653 /**********************************************************************/ 03654 03655 d0 = (double)(4 * n * np); 03656 d1 = (double)(l + 1); 03657 result += d1 * log10(d0); 03658 ASSERT( d0 >= 0. ); 03659 ASSERT( d1 != 0. ); 03660 03661 /**********************************************************************/ 03662 /* Calculate |(n-n')^(n+n'-2l-2)| */ 03663 /* NOTE: Here we are interested only */ 03664 /* magnitude of (n-n')^(n+n'-2l-2) */ 03665 /**********************************************************************/ 03666 03667 d0 = (double)(n - np); 03668 d1 = (double)(n + np - 2*l - 2); 03669 result += d1 * log10(fabs(d0)); 03670 ASSERT( fabs(d0) > 0. ); 03671 ASSERT( d1 != 0. ); 03672 03673 /* Calculate (n+n')^(-n-n') */ 03674 d0 = (double)(n + np); 03675 d1 = (double)(-n - np); 03676 result += d1 * log10(d0); 03677 ASSERT( d0 > 0. ); 03678 ASSERT( d1 != 0. ); 03679 03680 /**********************************************************************/ 03681 /* - - (1/2) */ 03682 /* | (n+l)! (n'+l-1)! | */ 03683 /* Calculate | ---------------- | */ 03684 /* | (n-l-1)! (n'-l)! | */ 03685 /* - - */ 03686 /**********************************************************************/ 03687 03688 ASSERT( n+l > 0 ); 03689 ld0 = lfactorial( n + l ); 03690 03691 ASSERT( np+l-1 > 0 ); 03692 ld1 = lfactorial( np + l - 1 ); 03693 03694 ASSERT( n-l-1 >= 0 ); 03695 ld2 = lfactorial( n - l - 1 ); 03696 03697 ASSERT( np-l >= 0 ); 03698 ld3 = lfactorial( np - l ); 03699 03700 ld4 = 0.5*((ld0+ld1)-(ld2+ld3)); 03701 03702 result += ld4; 03703 ASSERT( result != 0. ); 03704 return result; 03705 } 03706 03707 /***************************************************************************/ 03708 /* Find the Oscillator Strength for hydrogen for any */ 03709 /* transition n,l --> n',l' */ 03710 /* returns a double */ 03711 /***************************************************************************/ 03712 03713 /************************************************************************/ 03714 /* Find the Oscillator Strength for hydrogen for any */ 03715 /* transition n,l --> n',l' */ 03716 /* returns a double */ 03717 /* */ 03718 /* Einstein A() for the transition from the */ 03719 /* initial state n,l to the finial state n',l' */ 03720 /* require the Oscillator Strength f() */ 03721 /* */ 03722 /* hbar w max(l,l') | | 2 */ 03723 /* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */ 03724 /* 3 R_oo ( 2l + 1 ) | | */ 03725 /* */ 03726 /* */ 03727 /* */ 03728 /* E(n,l;n',l') max(l,l') | | 2 */ 03729 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 03730 /* 3 R_oo ( 2l + 1 ) | | */ 03731 /* */ 03732 /* */ 03733 /* See for example Gordan Drake's */ 03734 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 03735 /* */ 03736 /* Here R_oo is the infinite mass Rydberg length */ 03737 /* */ 03738 /* */ 03739 /* h c */ 03740 /* R_oo --- = 13.605698 eV */ 03741 /* {e} */ 03742 /* */ 03743 /* */ 03744 /* R_oo = 2.179874e-11 ergs */ 03745 /* */ 03746 /* w = omega */ 03747 /* = frequency of transition from n,l to n',l' */ 03748 /* */ 03749 /* */ 03750 /* */ 03751 /* here g_k are statistical weights obtained from */ 03752 /* the appropriate angular momentum quantum numbers */ 03753 /************************************************************************/ 03754 03755 /********************************************************************************/ 03756 /* Calc the Oscillator Strength f(*) given by */ 03757 /* */ 03758 /* E(n,l;n',l') max(l,l') | | 2 */ 03759 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 03760 /* 3 R_oo ( 2l + 1 ) | | */ 03761 /* */ 03762 /* See for example Gordan Drake's */ 03763 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 03764 /********************************************************************************/ 03765 03766 /************************************************************************/ 03767 /* Calc the Oscillator Strength f(*) given by */ 03768 /* */ 03769 /* E(n,l;n',l') max(l,l') | | 2 */ 03770 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 03771 /* 3 R_oo ( 2l + 1 ) | | */ 03772 /* */ 03773 /* f(n,l;n',l') is dimensionless. */ 03774 /* */ 03775 /* See for example Gordan Drake's */ 03776 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 03777 /* */ 03778 /* In the following, we have n > n' */ 03779 /************************************************************************/ 03780 03781 inline double OscStr_f( 03782 /* IN THE FOLLOWING WE HAVE n > n' */ 03783 /* principal quantum number, 1 for ground, upper level */ 03784 long int n, 03785 /* angular momentum, 0 for s */ 03786 long int l, 03787 /* principal quantum number, 1 for ground, lower level */ 03788 long int np, 03789 /* angular momentum, 0 for s */ 03790 long int lp, 03791 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 03792 long int iz 03793 ) 03794 { 03795 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.; 03796 long int i1 = 0, i2 = 0; 03797 03798 if( l > lp ) 03799 i1 = l; 03800 else 03801 i1 = lp; 03802 03803 i2 = 2*lp + 1; 03804 d0 = 1. / 3.; 03805 d1 = (double)i1 / (double)i2; 03806 /* hv() returns energy in ergs */ 03807 d2 = hv( n, np, iz ); 03808 d3 = d2 / EN1RYD; 03809 d4 = hri( n, l, np, lp ,iz ); 03810 d5 = d4 * d4; 03811 03812 d6 = d0 * d1 * d3 * d5; 03813 03814 return d6; 03815 } 03816 03817 /************************log version ***************************/ 03818 inline double OscStr_f_log10( long int n , long int l , long int np , long int lp , long int iz ) 03819 { 03820 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.; 03821 long int i1 = 0, i2 = 0; 03822 03823 if( l > lp ) 03824 i1 = l; 03825 else 03826 i1 = lp; 03827 03828 i2 = 2*lp + 1; 03829 d0 = 1. / 3.; 03830 d1 = (double)i1 / (double)i2; 03831 /* hv() returns energy in ergs */ 03832 d2 = hv( n, np, iz ); 03833 d3 = d2 / EN1RYD; 03834 d4 = hri_log10( n, l, np, lp ,iz ); 03835 d5 = d4 * d4; 03836 03837 d6 = d0 * d1 * d3 * d5; 03838 03839 return d6; 03840 } 03841 03842 STATIC double F21( long int a , long int b, long int c, double y, char A ) 03843 { 03844 DEBUG_ENTRY( "F21()" ); 03845 03846 double d1 = 0.; 03847 long int i1 = 0; 03848 double *yV; 03849 03850 /**************************************************************/ 03851 /* A must be either 'a' or 'b' */ 03852 /* and is use to determine which */ 03853 /* variable recursion will be over */ 03854 /**************************************************************/ 03855 03856 ASSERT( A == 'a' || A == 'b' ); 03857 03858 if( A == 'b' ) 03859 { 03860 /* if the recursion is over "b" */ 03861 /* then make it over "a" by switching these around */ 03862 i1 = a; 03863 a = b; 03864 b = i1; 03865 A = 'a'; 03866 } 03867 03868 /**************************************************************************************/ 03869 /* malloc space for the (dynamic) 1-d array */ 03870 /* 2_F_1 works via recursion and needs room to store intermediate results */ 03871 /* Here the + 5 is a safety margin */ 03872 /**************************************************************************************/ 03873 03874 /*Must use CALLOC*/ 03875 yV = (double*)CALLOC( sizeof(double), (size_t)(-a + 5) ); 03876 03877 /**********************************************************************************************/ 03878 /* begin sanity check, check order, and that none negative */ 03879 /* THE GENERAL CASE */ 03880 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 03881 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 03882 /* */ 03883 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 03884 /* */ 03885 /* a (1-x) (a + bx - c) */ 03886 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 03887 /* (a-c) (a-c) */ 03888 /* */ 03889 /* */ 03890 /* A similiar recusion relation holds for b with a <--> b. */ 03891 /* */ 03892 /* */ 03893 /* we have initial conditions */ 03894 /* */ 03895 /* */ 03896 /* F(0) = 1 with a = -1 */ 03897 /* */ 03898 /* b */ 03899 /* F(-1) = 1 - (---) x with a = -1 */ 03900 /* c */ 03901 /* */ 03902 /* For the first F() in the solution of the radial integral */ 03903 /* */ 03904 /* a = ( -n + l + 1 ) */ 03905 /* */ 03906 /* a = -n + l + 1 */ 03907 /* max(a) = max(-n) + max(l) + 1 */ 03908 /* = max(-n) + max(n - 1) + 1 */ 03909 /* = max(-n + n - 1) + 1 */ 03910 /* = max(-1) + 1 */ 03911 /* = 0 */ 03912 /* */ 03913 /* similiarly */ 03914 /* */ 03915 /* min(a) = min(-n) + min(l) + 1 */ 03916 /* = min(-n) + 0 + 1 */ 03917 /* = (-n) + 0 + 1 */ 03918 /* = -n + 1 */ 03919 /* */ 03920 /* a -> a' + 1 implies */ 03921 /* */ 03922 /* max(a') = -1 */ 03923 /* min(a') = -n */ 03924 /* */ 03925 /* For the second F() in the solution of the radial integral */ 03926 /* */ 03927 /* a = ( -n + l - 1 ) */ 03928 /* */ 03929 /* a = -n + l + 1 */ 03930 /* max(a) = max(-n) + max(l) - 1 */ 03931 /* = -n + (n - 1) - 1 */ 03932 /* = -2 */ 03933 /* */ 03934 /* similiarly */ 03935 /* */ 03936 /* min(a) = min(-n) + min(l) - 1 */ 03937 /* = (-n) + 0 - 1 */ 03938 /* = -n - 1 */ 03939 /* */ 03940 /* a -> a' + 1 implies */ 03941 /* */ 03942 /* max(a') = -3 */ 03943 /* min(a') = -n - 2 */ 03944 /**********************************************************************************************/ 03945 03946 ASSERT( a <= 0 ); 03947 ASSERT( b <= 0 ); 03948 ASSERT( c >= 0 ); 03949 03950 d1 = F21i(a, b, c, y, yV ); 03951 free( yV ); 03952 return d1; 03953 } 03954 03955 STATIC mx F21_mx( long int a , long int b, long int c, double y, char A ) 03956 { 03957 DEBUG_ENTRY( "F21_mx()" ); 03958 03959 mx result_mx = {0.0,0}; 03960 mxq *yV = NULL; 03961 03962 /**************************************************************/ 03963 /* A must be either 'a' or 'b' */ 03964 /* and is use to determine which */ 03965 /* variable recursion will be over */ 03966 /**************************************************************/ 03967 03968 ASSERT( A == 'a' || A == 'b' ); 03969 03970 if( A == 'b' ) 03971 { 03972 /* if the recursion is over "b" */ 03973 /* then make it over "a" by switching these around */ 03974 long int i1 = a; 03975 a = b; 03976 b = i1; 03977 A = 'a'; 03978 } 03979 03980 /**************************************************************************************/ 03981 /* malloc space for the (dynamic) 1-d array */ 03982 /* 2_F_1 works via recursion and needs room to store intermediate results */ 03983 /* Here the + 5 is a safety margin */ 03984 /**************************************************************************************/ 03985 03986 /*Must use CALLOC*/ 03987 yV = (mxq *)CALLOC( sizeof(mxq), (size_t)(-a + 5) ); 03988 03989 /**********************************************************************************************/ 03990 /* begin sanity check, check order, and that none negative */ 03991 /* THE GENERAL CASE */ 03992 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 03993 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 03994 /* */ 03995 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 03996 /* */ 03997 /* a (1-x) (a + bx - c) */ 03998 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 03999 /* (a-c) (a-c) */ 04000 /* */ 04001 /* */ 04002 /* A similiar recusion relation holds for b with a <--> b. */ 04003 /* */ 04004 /* */ 04005 /* we have initial conditions */ 04006 /* */ 04007 /* */ 04008 /* F(0) = 1 with a = -1 */ 04009 /* */ 04010 /* b */ 04011 /* F(-1) = 1 - (---) x with a = -1 */ 04012 /* c */ 04013 /* */ 04014 /* For the first F() in the solution of the radial integral */ 04015 /* */ 04016 /* a = ( -n + l + 1 ) */ 04017 /* */ 04018 /* a = -n + l + 1 */ 04019 /* max(a) = max(-n) + max(l) + 1 */ 04020 /* = max(-n) + max(n - 1) + 1 */ 04021 /* = max(-n + n - 1) + 1 */ 04022 /* = max(-1) + 1 */ 04023 /* = 0 */ 04024 /* */ 04025 /* similiarly */ 04026 /* */ 04027 /* min(a) = min(-n) + min(l) + 1 */ 04028 /* = min(-n) + 0 + 1 */ 04029 /* = (-n) + 0 + 1 */ 04030 /* = -n + 1 */ 04031 /* */ 04032 /* a -> a' + 1 implies */ 04033 /* */ 04034 /* max(a') = -1 */ 04035 /* min(a') = -n */ 04036 /* */ 04037 /* For the second F() in the solution of the radial integral */ 04038 /* */ 04039 /* a = ( -n + l - 1 ) */ 04040 /* */ 04041 /* a = -n + l + 1 */ 04042 /* max(a) = max(-n) + max(l) - 1 */ 04043 /* = -n + (n - 1) - 1 */ 04044 /* = -2 */ 04045 /* */ 04046 /* similiarly */ 04047 /* */ 04048 /* min(a) = min(-n) + min(l) - 1 */ 04049 /* = (-n) + 0 - 1 */ 04050 /* = -n - 1 */ 04051 /* */ 04052 /* a -> a' + 1 implies */ 04053 /* */ 04054 /* max(a') = -3 */ 04055 /* min(a') = -n - 2 */ 04056 /**********************************************************************************************/ 04057 04058 ASSERT( a <= 0 ); 04059 ASSERT( b <= 0 ); 04060 ASSERT( c >= 0 ); 04061 04062 result_mx = F21i_log(a, b, c, y, yV ); 04063 free( yV ); 04064 return result_mx; 04065 } 04066 04067 STATIC double F21i(long int a, long int b, long int c, double y, double *yV ) 04068 { 04069 DEBUG_ENTRY( "F21i()" ); 04070 04071 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.; 04072 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.; 04073 long int i1 = 0, i2 = 0; 04074 04075 if( a == 0 ) 04076 { 04077 return 1.; 04078 } 04079 else if( a == -1 ) 04080 { 04081 ASSERT( c != 0 ); 04082 d1 = (double)b; 04083 d2 = (double)c; 04084 d3 = y * (d1/d2); 04085 d4 = 1. - d3; 04086 return d4; 04087 } 04088 /* Check to see if y(-a) != 0 in a very round about way to avoid lclint:error:13 */ 04089 else if( yV[-a] != 0. ) 04090 { 04091 /* Return the stored result */ 04092 return yV[-a]; 04093 } 04094 else 04095 { 04096 /******************************************************************************************/ 04097 /* - - */ 04098 /* (a)(1 - y) | | (a + bx + c) */ 04099 /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */ 04100 /* (a - c) | | (a - c) */ 04101 /* - - */ 04102 /* */ 04103 /* */ 04104 /* */ 04105 /* */ 04106 /* */ 04107 /* with F(0) = 1 */ 04108 /* b */ 04109 /* and F(-1) = 1 - (---) y */ 04110 /* c */ 04111 /* */ 04112 /* */ 04113 /* */ 04114 /* Use a -> a' + 1 */ 04115 /* _ _ */ 04116 /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */ 04117 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */ 04118 /* (a' + 1 - c) | | (a' + 1 - c) */ 04119 /* - - */ 04120 /* */ 04121 /* For the first F() in the solution of the radial integral */ 04122 /* */ 04123 /* a = ( -n + l + 1 ) */ 04124 /* */ 04125 /* a = -n + l + 1 */ 04126 /* max(a) = max(-n) + max(l) + 1 */ 04127 /* = -n + max(n-1) + 1 */ 04128 /* = -n + n-1 + 1 */ 04129 /* = 0 */ 04130 /* */ 04131 /* similiarly */ 04132 /* */ 04133 /* min(a) = min(-n) + min(l) + 1 */ 04134 /* = min(-n) + 0 + 1 */ 04135 /* = -n + 1 */ 04136 /* */ 04137 /* a -> a' + 1 implies */ 04138 /* */ 04139 /* max(a') = -1 */ 04140 /* min(a') = -n */ 04141 /******************************************************************************************/ 04142 04143 i1= a + 1; 04144 i2= a + 1 - c; 04145 d0= (double)i2; 04146 ASSERT( i2 != 0 ); 04147 d1= 1. - y; 04148 d2= (double)i1 * d1; 04149 d3= d2 / d0; 04150 d4= (double)b * y; 04151 d5= d0 + d4; 04152 04153 d8= F21i( (long int)(a + 1), b, c, y, yV ); 04154 d9= F21i( (long int)(a + 2), b, c, y, yV ); 04155 04156 d10= d8 - d9; 04157 d11= d3 * d10; 04158 d12= d5 / d0; 04159 d13= d12 * d8; 04160 d14= d11 + d13; 04161 04162 /* Store the result for later use */ 04163 yV[-a] = d14; 04164 return d14; 04165 } 04166 } 04167 04168 STATIC mx F21i_log( long int a, long int b, long int c, double y, mxq *yV ) 04169 { 04170 DEBUG_ENTRY( "F21i_log()" ); 04171 04172 mx result_mx = {0.0,0}; 04173 04174 if( yV[-a].q != 0. ) 04175 { 04176 /* Return the stored result */ 04177 return yV[-a].mx; 04178 } 04179 else if( a == 0 ) 04180 { 04181 ASSERT( yV[-a].q == 0. ); 04182 04183 result_mx.m = 1.; 04184 result_mx.x = 0; 04185 04186 ASSERT( yV[-a].mx.m == 0. ); 04187 ASSERT( yV[-a].mx.x == 0 ); 04188 04189 yV[-a].q = 1; 04190 yV[-a].mx = result_mx; 04191 return result_mx; 04192 } 04193 else if( a == -1 ) 04194 { 04195 double d1 = (double)b; 04196 double d2 = (double)c; 04197 double d3 = y * (d1/d2); 04198 04199 ASSERT( yV[-a].q == 0. ); 04200 ASSERT( c != 0 ); 04201 ASSERT( y != 0. ); 04202 04203 result_mx.m = 1. - d3; 04204 result_mx.x = 0; 04205 04206 while ( fabs(result_mx.m) > 1.0e+25 ) 04207 { 04208 result_mx.m /= 1.0e+25; 04209 result_mx.x += 25; 04210 } 04211 04212 ASSERT( yV[-a].mx.m == 0. ); 04213 ASSERT( yV[-a].mx.x == 0 ); 04214 04215 yV[-a].q = 1; 04216 yV[-a].mx = result_mx; 04217 return result_mx; 04218 } 04219 else 04220 { 04221 /******************************************************************************************/ 04222 /* - - */ 04223 /* (a)(1 - y) | | (a + bx + c) */ 04224 /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */ 04225 /* (a - c) | | (a - c) */ 04226 /* - - */ 04227 /* */ 04228 /* */ 04229 /* with F(0) = 1 */ 04230 /* b */ 04231 /* and F(-1) = 1 - (---) y */ 04232 /* c */ 04233 /* */ 04234 /* */ 04235 /* */ 04236 /* Use a -> a' + 1 */ 04237 /* _ _ */ 04238 /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */ 04239 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */ 04240 /* (a' + 1 - c) | | (a' + 1 - c) */ 04241 /* - - */ 04242 /* */ 04243 /* For the first F() in the solution of the radial integral */ 04244 /* */ 04245 /* a = ( -n + l + 1 ) */ 04246 /* */ 04247 /* a = -n + l + 1 */ 04248 /* max(a) = max(-n) + max(l) + 1 */ 04249 /* = -n + max(n-1) + 1 */ 04250 /* = -n + n-1 + 1 */ 04251 /* = 0 */ 04252 /* */ 04253 /* similiarly */ 04254 /* */ 04255 /* min(a) = min(-n) + min(l) + 1 */ 04256 /* = min(-n) + 0 + 1 */ 04257 /* = -n + 1 */ 04258 /* */ 04259 /* a -> a' + 1 implies */ 04260 /* */ 04261 /* max(a') = -1 */ 04262 /* min(a') = -n */ 04263 /******************************************************************************************/ 04264 04265 mx d8 = {0.0,0}, d9 = {0.0,0}, d10 = {0.0,0}, d11 = {0.0,0}; 04266 04267 double db = (double)b; 04268 double d00 = (double)(a + 1 - c); 04269 double d0 = (double)(a + 1); 04270 double d1 = 1. - y; 04271 double d2 = d0 * d1; 04272 double d3 = d2 / d00; 04273 double d4 = db * y; 04274 double d5 = d00 + d4; 04275 double d6 = d5 / d00; 04276 04277 ASSERT( yV[-a].q == 0. ); 04278 04279 /******************************************************************************************/ 04280 /* _ _ */ 04281 /* (a' + 1) (1 - x) | | (a' + 1 - c) + b*x */ 04282 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ------------------ F(a'+1) */ 04283 /* (a' + 1 - c) | | (a' + 1 - c) */ 04284 /* - - */ 04285 /******************************************************************************************/ 04286 04287 d8= F21i_log( (a + 1), b, c, y, yV ); 04288 d9= F21i_log( (a + 2), b, c, y, yV ); 04289 04290 if( (d8.m) != 0. ) 04291 { 04292 d10.x = d8.x; 04293 d10.m = d8.m * (1. - (d9.m/d8.m) * powi( 10., (d9.x - d8.x))); 04294 } 04295 else 04296 { 04297 d10.m = -d9.m; 04298 d10.x = d9.x; 04299 } 04300 04301 d10.m *= d3; 04302 04303 d11.x = d8.x; 04304 d11.m = d6 * d8.m; 04305 04306 if( (d11.m) != 0. ) 04307 { 04308 result_mx.x = d11.x; 04309 result_mx.m = d11.m * (1. + (d10.m/d11.m) * powi( 10. , (d10.x - d11.x) ) ); 04310 } 04311 else 04312 { 04313 result_mx = d10; 04314 } 04315 04316 while ( fabs(result_mx.m) > 1.0e+25 ) 04317 { 04318 result_mx.m /= 1.0e+25; 04319 result_mx.x += 25; 04320 } 04321 04322 /* Store the result for later use */ 04323 yV[-a].q = 1; 04324 yV[-a].mx = result_mx; 04325 return result_mx; 04326 } 04327 } 04328 04329 /********************************************************************************/ 04330 04331 inline void normalize_mx( mx& target ) 04332 { 04333 while( fabs(target.m) > 1.0e+25 ) 04334 { 04335 target.m /= 1.0e+25; 04336 target.x += 25; 04337 } 04338 while( fabs(target.m) < 1.0e-25 ) 04339 { 04340 target.m *= 1.0e+25; 04341 target.x -= 25; 04342 } 04343 return; 04344 } 04345 04346 inline mx add_mx( const mx& a, const mx& b ) 04347 { 04348 mx result = {0.0,0}; 04349 04350 if( a.m != 0. ) 04351 { 04352 result.x = a.x; 04353 result.m = a.m * (1. + (b.m/a.m) * powi( 10. , (b.x - a.x) ) ); 04354 } 04355 else 04356 { 04357 result = b; 04358 } 04359 normalize_mx( result ); 04360 return result; 04361 } 04362 04363 inline mx sub_mx( const mx& a, const mx& b ) 04364 { 04365 mx result = {0.0,0}; 04366 mx minusb = b; 04367 minusb.m = -minusb.m; 04368 04369 result = add_mx( a, minusb ); 04370 normalize_mx( result ); 04371 04372 return result; 04373 } 04374 04375 inline mx mxify( double a ) 04376 { 04377 mx result_mx = {0.0,0}; 04378 04379 result_mx.x = 0; 04380 result_mx.m = a; 04381 normalize_mx( result_mx ); 04382 04383 return result_mx; 04384 } 04385 04386 inline double unmxify( const mx& a_mx ) 04387 { 04388 return a_mx.m * powi( 10., a_mx.x ); 04389 } 04390 04391 inline mx mxify_log10( double log10_a ) 04392 { 04393 mx result_mx = {0.0,0}; 04394 04395 while ( log10_a > 25.0 ) 04396 { 04397 log10_a -= 25.0; 04398 result_mx.x += 25; 04399 } 04400 04401 while ( log10_a < -25.0 ) 04402 { 04403 log10_a += 25.0; 04404 result_mx.x -= 25; 04405 } 04406 04407 result_mx.m = pow(10., log10_a); 04408 04409 return result_mx; 04410 } 04411 04412 inline mx mult_mx( const mx& a, const mx& b ) 04413 { 04414 mx result = {0.0,0}; 04415 04416 result.m = (a.m * b.m); 04417 result.x = (a.x + b.x); 04418 normalize_mx( result ); 04419 04420 return result; 04421 } 04422 04423 inline double log10_prodxx( long int lp, double Ksqrd ) 04424 { 04425 /**********************************************************************/ 04426 /* | s=l' | s=l' */ 04427 /* | ----- | --- */ 04428 /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */ 04429 /* | | | | --- */ 04430 /* | s=0 | s=0 */ 04431 /**********************************************************************/ 04432 04433 if( lp == 0 ) 04434 return 0.; 04435 04436 double partsum = 0.; 04437 for( long int s = 1; s <= lp; s++ ) 04438 { 04439 double s2 = pow2((double)s); 04440 partsum += log10( 1. + s2*Ksqrd ); 04441 04442 ASSERT( partsum >= 0. ); 04443 } 04444 return partsum; 04445 }
1.4.7