00001 /* This file is part of Cloudy and is copyright (C)1978-2013 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)max(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)max(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 printf( "This code should be inaccessible\n\n" ); 01281 cdEXIT(EXIT_FAILURE); 01282 } 01283 01284 /************************************************************************************************/ 01285 /* *** bhGp.c *** */ 01286 /* */ 01287 /* Here we calculate G(n,l; K,l') with the recursive formula */ 01288 /* equation: (3.24) */ 01289 /* */ 01290 /* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */ 01291 /* */ 01292 /* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l+1; K,l) */ 01293 /* */ 01294 /* Under the transformation l -> l + 1 this gives */ 01295 /* */ 01296 /* 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) */ 01297 /* */ 01298 /* - 4n^2 (n^2-((l+1)+1)^2)[ 1+((l+1)K)^2 ] G(n,l+1+1; K,l+1) */ 01299 /* */ 01300 /* or */ 01301 /* */ 01302 /* 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) */ 01303 /* */ 01304 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01305 /* */ 01306 /* from the reference */ 01307 /* M. Brocklehurst */ 01308 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 01309 /* */ 01310 /* */ 01311 /* * This is valid for the case l=l'+1 * */ 01312 /* * CASE: l = l'+1 * */ 01313 /* * Here the p in bhGp() refers * */ 01314 /* * to the Plus sign(+) in l=l'+1 * */ 01315 /************************************************************************************************/ 01316 01317 STATIC double bhGp( 01318 long int q, 01319 double K, 01320 long int n, 01321 long int l, 01322 long int lp, 01323 /* Temporary storage for intermediate */ 01324 /* results of the recursive routine */ 01325 double *rcsvV, 01326 double GK 01327 ) 01328 { 01329 DEBUG_ENTRY( "bhGp()" ); 01330 01331 /* static long int rcsv_Level = 1; 01332 printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */ 01333 01334 ASSERT( l == lp + 1 ); 01335 01336 long int rindx = 2*q; 01337 01338 if( rcsvV[rindx] == 0. ) 01339 { 01340 /* SPECIAL CASE: n = l+1 = l'+2 */ 01341 if( q == n - 1 ) 01342 { 01343 double Ksqrd = K * K; 01344 double n2 = (double)(n*n); 01345 01346 double dd1 = (double)(2 * n); 01347 double dd2 = ( 1. + ( n2 * Ksqrd)); 01348 01349 /* (1+(n K)^2) */ 01350 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 01351 /* 2n */ 01352 double G1 = ((dd2 * GK) / dd1); 01353 01354 ASSERT( l == lp + 1 ); 01355 ASSERT( Ksqrd != 0. ); 01356 ASSERT( dd1 != 0. ); 01357 ASSERT( dd2 != 0. ); 01358 ASSERT( G1 != 0. ); 01359 01360 rcsvV[rindx] = G1; 01361 return G1; 01362 } 01363 /* SPECIAL CASE: n = l+2 = l'+3 */ 01364 else if( q == (n - 2) ) 01365 { 01366 double Ksqrd = (K*K); 01367 double n2 = (double)(n*n); 01368 01369 /* */ 01370 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */ 01371 /* */ 01372 double dd1 = (double) (2 * n); 01373 double dd2 = ( 1. + ( n2 * Ksqrd)); 01374 double G1 = ((dd2 * GK) / dd1); 01375 01376 /* */ 01377 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */ 01378 /* */ 01379 double dd3 = (double)((2 * n) - 1); 01380 double dd4 = (double)(n - 1); 01381 double dd5 = (4. + (dd4 * dd2)); 01382 double G2 = (dd3 * dd5 * G1); 01383 01384 ASSERT( l == lp + 1 ); 01385 ASSERT( Ksqrd != 0. ); 01386 ASSERT( n2 != 0. ); 01387 ASSERT( dd1 != 0. ); 01388 ASSERT( dd2 != 0. ); 01389 ASSERT( dd3 != 0. ); 01390 ASSERT( dd4 != 0. ); 01391 ASSERT( dd5 != 0. ); 01392 ASSERT( G1 != 0. ); 01393 ASSERT( G2 != 0. ); 01394 01395 rcsvV[rindx] = G2; 01396 return G2; 01397 } 01398 /* The GENERAL CASE n > l + 2 */ 01399 else 01400 { 01401 /******************************************************************************/ 01402 /* 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) */ 01403 /* */ 01404 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01405 /* */ 01406 /* FROM Eq. 3.24 */ 01407 /* */ 01408 /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */ 01409 /* */ 01410 /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */ 01411 /******************************************************************************/ 01412 01413 long int lp1 = (q + 1); 01414 long int lp2 = (q + 2); 01415 01416 double Ksqrd = (K*K); 01417 double n2 = (double)(n * n); 01418 double lp1s = (double)(lp1 * lp1); 01419 double lp2s = (double)(lp2 * lp2); 01420 01421 double d1 = (4. * n2); 01422 double d2 = (4. * lp1s); 01423 double d3 = (double)((lp1)*((2 * q) + 3)); 01424 double d4 = (1. + (n2 * Ksqrd)); 01425 double d5 = (d1 - d2 + (d3 * d4)); 01426 double d5_1 = d5 * bhGp( (q+1), K, n, l, lp, rcsvV, GK ); 01427 01428 01429 /* 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) */ 01430 /* */ 01431 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01432 01433 double d6 = (n2 - lp2s); 01434 double d7 = (1. + (lp1s * Ksqrd)); 01435 double d8 = (d1 * d6 * d7); 01436 double d8_1 = d8 * bhGp( (q+2), K, n, l, lp, rcsvV, GK ); 01437 double d9 = (d5_1 - d8_1); 01438 01439 ASSERT( l == lp + 1 ); 01440 ASSERT( Ksqrd != 0. ); 01441 ASSERT( n2 != 0. ); 01442 01443 ASSERT( lp1s != 0. ); 01444 ASSERT( lp2s != 0. ); 01445 01446 ASSERT( d1 != 0. ); 01447 ASSERT( d2 != 0. ); 01448 ASSERT( d3 != 0. ); 01449 ASSERT( d4 != 0. ); 01450 ASSERT( d5 != 0. ); 01451 ASSERT( d6 != 0. ); 01452 ASSERT( d7 != 0. ); 01453 ASSERT( d8 != 0. ); 01454 ASSERT( d9 != 0. ); 01455 01456 rcsvV[rindx] = d9; 01457 return d9; 01458 } 01459 } 01460 else 01461 { 01462 ASSERT( rcsvV[rindx] != 0. ); 01463 return rcsvV[rindx]; 01464 } 01465 } 01466 01467 /***********************log version*******************************/ 01468 STATIC mx bhGp_mx( 01469 long int q, 01470 double K, 01471 long int n, 01472 long int l, 01473 long int lp, 01474 /* Temporary storage for intermediate */ 01475 /* results of the recursive routine */ 01476 mxq *rcsvV_mxq, 01477 const mx& GK_mx 01478 ) 01479 { 01480 DEBUG_ENTRY( "bhGp_mx()" ); 01481 01482 /* static long int rcsv_Level = 1; */ 01483 /* printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */ 01484 01485 ASSERT( l == lp + 1 ); 01486 01487 long int rindx = 2*q; 01488 01489 if( rcsvV_mxq[rindx].q == 0 ) 01490 { 01491 /* SPECIAL CASE: n = l+1 = l'+2 */ 01492 if( q == n - 1 ) 01493 { 01494 /******************************************************/ 01495 /* (1+(n K)^2) */ 01496 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 01497 /* 2n */ 01498 /******************************************************/ 01499 01500 double Ksqrd = (K * K); 01501 double n2 = (double)(n*n); 01502 01503 double dd1 = (double) (2 * n); 01504 double dd2 = ( 1. + ( n2 * Ksqrd)); 01505 double dd3 = dd2/dd1; 01506 01507 mx dd3_mx = mxify( dd3 ); 01508 mx G1_mx = mult_mx( dd3_mx, GK_mx); 01509 01510 normalize_mx( G1_mx ); 01511 01512 ASSERT( l == lp + 1 ); 01513 ASSERT( Ksqrd != 0. ); 01514 ASSERT( n2 != 0. ); 01515 ASSERT( dd1 != 0. ); 01516 ASSERT( dd2 != 0. ); 01517 01518 rcsvV_mxq[rindx].q = 1; 01519 rcsvV_mxq[rindx].mx = G1_mx; 01520 return G1_mx; 01521 } 01522 /* SPECIAL CASE: n = l+2 = l'+3 */ 01523 else if( q == (n - 2) ) 01524 { 01525 /****************************************************************/ 01526 /* */ 01527 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/ 01528 /* */ 01529 /****************************************************************/ 01530 /* (1+(n K)^2) */ 01531 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 01532 /* 2n */ 01533 /****************************************************************/ 01534 01535 double Ksqrd = (K*K); 01536 double n2 = (double)(n*n); 01537 double dd1 = (double)(2 * n); 01538 double dd2 = ( 1. + ( n2 * Ksqrd) ); 01539 double dd3 = (dd2/dd1); 01540 double dd4 = (double) ((2 * n) - 1); 01541 double dd5 = (double) (n - 1); 01542 double dd6 = (4. + (dd5 * dd2)); 01543 double dd7 = dd4 * dd6; 01544 01545 /****************************************************************/ 01546 /* */ 01547 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/ 01548 /* */ 01549 /****************************************************************/ 01550 01551 mx dd3_mx = mxify( dd3 ); 01552 mx dd7_mx = mxify( dd7 ); 01553 mx G1_mx = mult_mx( dd3_mx, GK_mx ); 01554 mx G2_mx = mult_mx( dd7_mx, G1_mx ); 01555 01556 normalize_mx( G2_mx ); 01557 01558 ASSERT( l == lp + 1 ); 01559 ASSERT( Ksqrd != 0. ); 01560 ASSERT( n2 != 0. ); 01561 ASSERT( dd1 != 0. ); 01562 ASSERT( dd2 != 0. ); 01563 ASSERT( dd3 != 0. ); 01564 ASSERT( dd4 != 0. ); 01565 ASSERT( dd5 != 0. ); 01566 ASSERT( dd6 != 0. ); 01567 ASSERT( dd7 != 0. ); 01568 01569 rcsvV_mxq[rindx].q = 1; 01570 rcsvV_mxq[rindx].mx = G2_mx; 01571 return G2_mx; 01572 } 01573 /* The GENERAL CASE n > l + 2*/ 01574 else 01575 { 01576 /**************************************************************************************/ 01577 /* 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) */ 01578 /* */ 01579 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01580 /* */ 01581 /* FROM Eq. 3.24 */ 01582 /* */ 01583 /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */ 01584 /* */ 01585 /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */ 01586 /**************************************************************************************/ 01587 01588 long int lp1 = (q + 1); 01589 long int lp2 = (q + 2); 01590 01591 double Ksqrd = (K * K); 01592 double n2 = (double)(n * n); 01593 double lp1s = (double)(lp1 * lp1); 01594 double lp2s = (double)(lp2 * lp2); 01595 01596 double d1 = (4. * n2); 01597 double d2 = (4. * lp1s); 01598 double d3 = (double)((lp1)*((2 * q) + 3)); 01599 double d4 = (1. + (n2 * Ksqrd)); 01600 /* [ 4n^2 - 4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] */ 01601 double d5 = d1 - d2 + (d3 * d4); 01602 01603 /* (n^2-(l+2)^2) */ 01604 double d6 = (n2 - lp2s); 01605 01606 /* [ 1+((l+1)K)^2 ] */ 01607 double d7 = (1. + (lp1s * Ksqrd)); 01608 01609 /* { 4n^2 (n^2-(l+1)^2)[ 1+((l+1) K)^2 ] } */ 01610 double d8 = (d1 * d6 * d7); 01611 01612 /**************************************************************************************/ 01613 /* 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) */ 01614 /* */ 01615 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */ 01616 /**************************************************************************************/ 01617 01618 mx d5_mx=mxify( d5 ); 01619 mx d8_mx=mxify( d8 ); 01620 01621 mx t0_mx = bhGp_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx ); 01622 mx t1_mx = bhGp_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx ); 01623 01624 mx d9_mx = mult_mx( d5_mx, t0_mx ); 01625 mx d10_mx = mult_mx( d8_mx, t1_mx ); 01626 01627 mx result_mx = sub_mx( d9_mx, d10_mx ); 01628 normalize_mx( result_mx ); 01629 01630 ASSERT( d1 != 0. ); 01631 ASSERT( d2 != 0. ); 01632 ASSERT( d3 != 0. ); 01633 ASSERT( d4 != 0. ); 01634 ASSERT( d5 != 0. ); 01635 ASSERT( d6 != 0. ); 01636 ASSERT( d7 != 0. ); 01637 ASSERT( d8 != 0. ); 01638 01639 ASSERT( l == lp + 1 ); 01640 ASSERT( Ksqrd != 0. ); 01641 ASSERT( n2 != 0. ); 01642 ASSERT( lp1s != 0. ); 01643 ASSERT( lp2s != 0. ); 01644 01645 rcsvV_mxq[rindx].q = 1; 01646 rcsvV_mxq[rindx].mx = result_mx; 01647 return result_mx; 01648 } 01649 } 01650 else 01651 { 01652 ASSERT( rcsvV_mxq[rindx].q != 0 ); 01653 rcsvV_mxq[rindx].q = 1; 01654 return rcsvV_mxq[rindx].mx; 01655 } 01656 } 01657 01658 /************************************************************************************************/ 01659 /* *** bhGm.c *** */ 01660 /* */ 01661 /* Here we calculate G(n,l; K,l') with the recursive formula */ 01662 /* equation: (3.23) */ 01663 /* */ 01664 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */ 01665 /* */ 01666 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */ 01667 /* */ 01668 /* Under the transformation l -> l + 2 this gives */ 01669 /* */ 01670 /* 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) */ 01671 /* */ 01672 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+2+1)^2 K^2 ] G(n,l+2; K,l+2+1) */ 01673 /* */ 01674 /* */ 01675 /* or */ 01676 /* */ 01677 /* 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) */ 01678 /* */ 01679 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */ 01680 /* */ 01681 /* */ 01682 /* from the reference */ 01683 /* M. Brocklehurst */ 01684 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 01685 /* */ 01686 /* */ 01687 /* * This is valid for the case l=l'-1 * */ 01688 /* * CASE: l = l'-1 * */ 01689 /* * Here the p in bhGm() refers * */ 01690 /* * to the Minus sign(-) in l=l'-1 * */ 01691 /************************************************************************************************/ 01692 01693 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000 01694 #pragma optimize("", off) 01695 #endif 01696 STATIC double bhGm( 01697 long int q, 01698 double K, 01699 long int n, 01700 long int l, 01701 long int lp, 01702 double *rcsvV, 01703 double GK 01704 ) 01705 { 01706 DEBUG_ENTRY( "bhGm()" ); 01707 01708 ASSERT( l == lp - 1 ); 01709 ASSERT( l < n ); 01710 01711 long int rindx = 2*q+1; 01712 01713 if( rcsvV[rindx] == 0. ) 01714 { 01715 /* CASE: l = n - 1 */ 01716 if( q == n - 1 ) 01717 { 01718 ASSERT( l == lp - 1 ); 01719 rcsvV[rindx] = GK; 01720 return GK; 01721 } 01722 /* CASE: l = n - 2 */ 01723 else if( q == n - 2 ) 01724 { 01725 double dd1 = 0.; 01726 double dd2 = 0.; 01727 01728 double G2 = 0.; 01729 01730 double Ksqrd = 0.; 01731 double n1 = 0.; 01732 double n2 = 0.; 01733 01734 ASSERT(l == lp - 1); 01735 01736 /* K^2 */ 01737 Ksqrd = K * K; 01738 ASSERT( Ksqrd != 0. ); 01739 01740 /* n */ 01741 n1 = (double)n; 01742 ASSERT( n1 != 0. ); 01743 01744 /* n^2 */ 01745 n2 = (double)(n*n); 01746 ASSERT( n2 != 0. ); 01747 01748 /* equation: (3.20) */ 01749 /* G(n,n-2; K,n-1) = */ 01750 /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */ 01751 dd1 = (double) ((2 * n) - 1); 01752 ASSERT( dd1 != 0. ); 01753 01754 dd2 = (1. + (n2 * Ksqrd)); 01755 ASSERT( dd2 != 0. ); 01756 01757 G2 = dd1 * dd2 * n1 * GK; 01758 ASSERT( G2 != 0. ); 01759 01760 rcsvV[rindx] = G2; 01761 ASSERT( G2 != 0. ); 01762 return G2; 01763 } 01764 else 01765 { 01766 long int lp2 = (q + 2); 01767 long int lp3 = (q + 3); 01768 01769 double lp2s = (double)(lp2 * lp2); 01770 double lp3s = (double)(lp3 * lp3); 01771 01772 /******************************************************************************/ 01773 /* 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) */ 01774 /* */ 01775 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */ 01776 /* */ 01777 /* */ 01778 /* FROM Eq. 3.23 */ 01779 /* */ 01780 /* 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) */ 01781 /* */ 01782 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */ 01783 /******************************************************************************/ 01784 01785 double Ksqrd = (K*K); 01786 double n2 = (double)(n*n); 01787 double d1 = (4. * n2); 01788 double d2 = (4. * lp2s); 01789 double d3 = (double)(lp2)*((2*q)+3); 01790 double d4 = (1. + (n2 * Ksqrd)); 01791 /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */ 01792 double d5 = d1 - d2 + (d3 * d4); 01793 01794 /******************************************************************************/ 01795 /* 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) */ 01796 /* */ 01797 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */ 01798 /******************************************************************************/ 01799 01800 double d6 = (n2 - lp2s); 01801 /* [ 1+((l+3)K)^2 ] */ 01802 double d7 = (1. + (lp3s * Ksqrd)); 01803 /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */ 01804 double d8 = d1 * d6 * d7; 01805 double d9 = d5 * bhGm( (q+1), K, n, l, lp, rcsvV, GK ); 01806 double d10 = d8 * bhGm( (q+2), K, n, l, lp, rcsvV, GK ); 01807 double d11 = d9 - d10; 01808 01809 ASSERT( l == lp - 1 ); 01810 ASSERT( lp2s != 0. ); 01811 ASSERT( Ksqrd != 0. ); 01812 ASSERT( n2 != 0. ); 01813 ASSERT( d1 != 0. ); 01814 ASSERT( d2 != 0. ); 01815 ASSERT( d3 != 0. ); 01816 ASSERT( d4 != 0. ); 01817 ASSERT( d5 != 0. ); 01818 ASSERT( d6 != 0. ); 01819 ASSERT( d7 != 0. ); 01820 ASSERT( d8 != 0. ); 01821 ASSERT( d9 != 0. ); 01822 ASSERT( d10 != 0. ); 01823 ASSERT( lp3s != 0. ); 01824 01825 rcsvV[rindx] = d11; 01826 return d11; 01827 } 01828 } 01829 else 01830 { 01831 ASSERT( rcsvV[rindx] != 0. ); 01832 return rcsvV[rindx]; 01833 } 01834 } 01835 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000 01836 #pragma optimize("", on) 01837 #endif 01838 01839 /************************log version***********************************/ 01840 STATIC mx bhGm_mx( 01841 long int q, 01842 double K, 01843 long int n, 01844 long int l, 01845 long int lp, 01846 mxq *rcsvV_mxq, 01847 const mx& GK_mx 01848 ) 01849 { 01850 DEBUG_ENTRY( "bhGm_mx()" ); 01851 01852 /*static long int rcsv_Level = 1; */ 01853 /*printf( "bhGm(): recursion level:\t%li\n",rcsv_Level++ ); */ 01854 01855 ASSERT( l == lp - 1 ); 01856 ASSERT( l < n ); 01857 01858 long int rindx = 2*q+1; 01859 01860 if( rcsvV_mxq[rindx].q == 0 ) 01861 { 01862 /* CASE: l = n - 1 */ 01863 if( q == n - 1 ) 01864 { 01865 mx result_mx = GK_mx; 01866 normalize_mx( result_mx ); 01867 01868 rcsvV_mxq[rindx].q = 1; 01869 rcsvV_mxq[rindx].mx = result_mx; 01870 01871 ASSERT(l == lp - 1); 01872 return result_mx; 01873 } 01874 /* CASE: l = n - 2 */ 01875 else if( q == n - 2 ) 01876 { 01877 double Ksqrd = (K * K); 01878 double n1 = (double)n; 01879 double n2 = (double) (n*n); 01880 double dd1 = (double)((2 * n) - 1); 01881 double dd2 = (1. + (n2 * Ksqrd)); 01882 /*(2n-1)(1+(n K)^2) n*/ 01883 double dd3 = (dd1*dd2*n1); 01884 01885 /******************************************************/ 01886 /* G(n,n-2; K,n-1) = */ 01887 /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */ 01888 /******************************************************/ 01889 01890 mx dd3_mx = mxify( dd3 ); 01891 mx G2_mx = mult_mx( dd3_mx, GK_mx ); 01892 01893 normalize_mx( G2_mx ); 01894 01895 ASSERT( l == lp - 1); 01896 ASSERT( n1 != 0. ); 01897 ASSERT( n2 != 0. ); 01898 ASSERT( dd1 != 0. ); 01899 ASSERT( dd2 != 0. ); 01900 ASSERT( dd3 != 0. ); 01901 ASSERT( Ksqrd != 0. ); 01902 01903 rcsvV_mxq[rindx].q = 1; 01904 rcsvV_mxq[rindx].mx = G2_mx; 01905 return G2_mx; 01906 } 01907 else 01908 { 01909 /******************************************************************************/ 01910 /* 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) */ 01911 /* */ 01912 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */ 01913 /* */ 01914 /* */ 01915 /* FROM Eq. 3.23 */ 01916 /* */ 01917 /* 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) */ 01918 /* */ 01919 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */ 01920 /******************************************************************************/ 01921 01922 long int lp2 = (q + 2); 01923 long int lp3 = (q + 3); 01924 01925 double lp2s = (double)(lp2 * lp2); 01926 double lp3s = (double)(lp3 * lp3); 01927 double n2 = (double)(n*n); 01928 double Ksqrd = (K * K); 01929 01930 /******************************************************/ 01931 /* [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] */ 01932 /******************************************************/ 01933 01934 double d1 = (4. * n2); 01935 double d2 = (4. * lp2s); 01936 double d3 = (double)(lp2)*((2*q)+3); 01937 double d4 = (1. + (n2 * Ksqrd)); 01938 /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */ 01939 double d5 = d1 - d2 + (d3 * d4); 01940 01941 mx d5_mx=mxify(d5); 01942 01943 /******************************************************/ 01944 /* 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] */ 01945 /******************************************************/ 01946 01947 double d6 = (n2 - lp2s); 01948 double d7 = (1. + (lp3s * Ksqrd)); 01949 /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */ 01950 double d8 = d1 * d6 * d7; 01951 01952 mx d8_mx = mxify(d8); 01953 01954 /******************************************************************************/ 01955 /* 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) */ 01956 /* */ 01957 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */ 01958 /******************************************************************************/ 01959 01960 mx d9_mx = bhGm_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx ); 01961 mx d10_mx = bhGm_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx ); 01962 mx d11_mx = mult_mx( d5_mx, d9_mx ); 01963 mx d12_mx = mult_mx( d8_mx, d10_mx); 01964 mx result_mx = sub_mx( d11_mx , d12_mx ); 01965 rcsvV_mxq[rindx].q = 1; 01966 rcsvV_mxq[rindx].mx = result_mx; 01967 01968 ASSERT(l == lp - 1); 01969 ASSERT(n2 != 0. ); 01970 ASSERT(lp2s != 0.); 01971 ASSERT( lp3s != 0.); 01972 ASSERT(Ksqrd != 0.); 01973 01974 ASSERT(d1 != 0.); 01975 ASSERT(d2 != 0.); 01976 ASSERT(d3 != 0.); 01977 ASSERT(d4 != 0.); 01978 ASSERT(d5 != 0.); 01979 ASSERT(d6 != 0.); 01980 ASSERT(d7 != 0.); 01981 ASSERT(d8 != 0.); 01982 return result_mx; 01983 } 01984 } 01985 else 01986 { 01987 ASSERT( rcsvV_mxq[rindx].q != 0 ); 01988 return rcsvV_mxq[rindx].mx; 01989 } 01990 } 01991 01992 /****************************************************************************************/ 01993 /* */ 01994 /* bhg.c */ 01995 /* */ 01996 /* From reference; */ 01997 /* M. Brocklehurst */ 01998 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */ 01999 /* */ 02000 /* */ 02001 /* We wish to compute the following function, */ 02002 /* */ 02003 /* equation: (3.17) */ 02004 /* - s=l' - (1/2) */ 02005 /* | (n+l)! ----- | */ 02006 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */ 02007 /* | (n-l-1)! | | | */ 02008 /* - s=0 - */ 02009 /* */ 02010 /* Using various recursion relations (for l'=l+1) */ 02011 /* */ 02012 /* equation: (3.23) */ 02013 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2+l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */ 02014 /* */ 02015 /* - 4n^2 (n^2-l^2)[1+(l+1)^2 K^2] G(n,l; K,l+1) */ 02016 /* */ 02017 /* and (for l'=l-1) */ 02018 /* */ 02019 /* equation: (3.24) */ 02020 /* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l; K,l-1) */ 02021 /* */ 02022 /* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l; K,l+1) */ 02023 /* */ 02024 /* */ 02025 /* the starting point for the recursion relations are: */ 02026 /* */ 02027 /* */ 02028 /* equation (3.18): */ 02029 /* */ 02030 /* | pi |(1/2) 8n */ 02031 /* G(n,n-1; 0,n) = | -- | ------- (4n)^2 exp(-2n) */ 02032 /* | 2 | (2n-1)! */ 02033 /* */ 02034 /* equation (3.20): */ 02035 /* */ 02036 /* exp(2n-2/K tan^(-1)(n K) */ 02037 /* G(n,n-1; K,n) = --------------------------------------- */ 02038 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K)^(n+2) */ 02039 /* */ 02040 /* */ 02041 /* */ 02042 /* equation (3.20): */ 02043 /* G(n,n-2; K,n-1) = (2n-2)(1+(n K)^2) n G(n,n-1; K,n) */ 02044 /* */ 02045 /* */ 02046 /* equation (3.21): */ 02047 /* */ 02048 /* (1+(n K)^2) */ 02049 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */ 02050 /* 2n */ 02051 /****************************************************************************************/ 02052 02053 STATIC double bhg( 02054 double K, 02055 long int n, 02056 long int l, 02057 long int lp, 02058 /* Temporary storage for intermediate */ 02059 /* results of the recursive routine */ 02060 double *rcsvV 02061 ) 02062 { 02063 DEBUG_ENTRY( "bhg()" ); 02064 02065 double ld1 = factorial( n + l ); 02066 double ld2 = factorial( n - l - 1 ); 02067 double ld3 = (ld1 / ld2); 02068 02069 double partprod = local_product( K , lp ); 02070 02071 /**************************************************************************************/ 02072 /* equation: (3.17) */ 02073 /* - s=l' - (1/2) */ 02074 /* | (n+l)! ----- | */ 02075 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */ 02076 /* | (n-l-1)! | | | */ 02077 /* - s=0 - */ 02078 /**************************************************************************************/ 02079 02080 /**********************************************/ 02081 /* - s=l' - (1/2) */ 02082 /* | (n+l)! ----- | */ 02083 /* | -------- | | (1 + s^2 K^2) | */ 02084 /* | (n-l-1)! | | | */ 02085 /* - s=0 - */ 02086 /**********************************************/ 02087 02088 double d2 = sqrt( ld3 * partprod ); 02089 double d3 = powi( (2 * n) , (l - n) ); 02090 double d4 = bhG( K, n, l, lp, rcsvV ); 02091 double d5 = (d2 * d3); 02092 double d6 = (d5 * d4); 02093 02094 ASSERT(K != 0.); 02095 ASSERT( (n+l) >= 1 ); 02096 ASSERT( ((n-l)-1) >= 0 ); 02097 02098 ASSERT( partprod != 0. ); 02099 02100 ASSERT( ld1 != 0. ); 02101 ASSERT( ld2 != 0. ); 02102 ASSERT( ld3 != 0. ); 02103 02104 ASSERT( d2 != 0. ); 02105 ASSERT( d3 != 0. ); 02106 ASSERT( d4 != 0. ); 02107 ASSERT( d5 != 0. ); 02108 ASSERT( d6 != 0. ); 02109 return d6; 02110 } 02111 02112 /********************log version**************************/ 02113 STATIC double bhg_log( 02114 double K, 02115 long int n, 02116 long int l, 02117 long int lp, 02118 /* Temporary storage for intermediate */ 02119 /* results of the recursive routine */ 02120 mxq *rcsvV_mxq 02121 ) 02122 { 02123 /**************************************************************************************/ 02124 /* equation: (3.17) */ 02125 /* - s=l' - (1/2) */ 02126 /* | (n+l)! ----- | */ 02127 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */ 02128 /* | (n-l-1)! | | | */ 02129 /* - s=0 - */ 02130 /**************************************************************************************/ 02131 02132 DEBUG_ENTRY( "bhg_log()" ); 02133 02134 double d1 = (double)(2*n); 02135 double d2 = (double)(l-n); 02136 double Ksqrd = (K*K); 02137 02138 /**************************************************************************************/ 02139 /* */ 02140 /* | (n+l)! | */ 02141 /* log10 | -------- | = log10((n+1)!) - log10((n-l-1)!) */ 02142 /* | (n-l-1)! | */ 02143 /* */ 02144 /**************************************************************************************/ 02145 02146 double ld1 = lfactorial( n + l ); 02147 double ld2 = lfactorial( n - l - 1 ); 02148 02149 /**********************************************************************/ 02150 /* | s=l' | s=l' */ 02151 /* | ----- | --- */ 02152 /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */ 02153 /* | | | | --- */ 02154 /* | s=0 | s=0 */ 02155 /**********************************************************************/ 02156 02157 double ld3 = log10_prodxx( lp, Ksqrd ); 02158 02159 /**********************************************/ 02160 /* - s=l' - (1/2) */ 02161 /* | (n+l)! ----- | */ 02162 /* | -------- | | (1 + s^2 K^2) | */ 02163 /* | (n-l-1)! | | | */ 02164 /* - s=0 - */ 02165 /**********************************************/ 02166 02167 /***********************************************************************/ 02168 /* */ 02169 /* | - s=l' - (1/2) | */ 02170 /* | | (n+l)! ----- | | */ 02171 /* log10| | -------- | | (1 + s^2 K^2) | | == */ 02172 /* | | (n-l-1)! | | | | */ 02173 /* | - s=0 - | */ 02174 /* */ 02175 /* | | s=l' | | */ 02176 /* | | (n+l)! | | ----- | | */ 02177 /* (1/2)* |log10 | -------- | + log10 | | | (1 + s^2 K^2) | | */ 02178 /* | | (n-l-1)! | | | | | | */ 02179 /* | | s=0 | | */ 02180 /* */ 02181 /***********************************************************************/ 02182 02183 double ld4 = (1./2.) * ( ld3 + ld1 - ld2 ); 02184 02185 /**********************************************/ 02186 /* (2n)^(l-n) */ 02187 /**********************************************/ 02188 /* log10( 2n^(L-n) ) = (L-n) log10( 2n ) */ 02189 /**********************************************/ 02190 02191 double ld5 = d2 * log10( d1 ); 02192 02193 /**************************************************************************************/ 02194 /* equation: (3.17) */ 02195 /* - s=l' - (1/2) */ 02196 /* | (n+l)! ----- | */ 02197 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) * G(n,l; K,l') */ 02198 /* | (n-l-1)! | | | */ 02199 /* - s=0 - */ 02200 /**************************************************************************************/ 02201 02202 /****************************************************/ 02203 /* */ 02204 /* - s=l' - (1/2) */ 02205 /* | (n+l)! ----- | */ 02206 /* | -------- | | (1 + s^2 K^2) | * (2n)^(L-n) */ 02207 /* | (n-l-1)! | | | */ 02208 /* - s=0 - */ 02209 /****************************************************/ 02210 02211 double ld6 = (ld5+ld4); 02212 02213 mx d6_mx = mxify_log10( ld6 ); 02214 mx dd1_mx = bhG_mx( K, n, l, lp, rcsvV_mxq ); 02215 mx dd2_mx = mult_mx( d6_mx, dd1_mx ); 02216 normalize_mx( dd2_mx ); 02217 double result = unmxify( dd2_mx ); 02218 02219 ASSERT( result != 0. ); 02220 02221 ASSERT( Ksqrd != 0. ); 02222 ASSERT( ld3 >= 0. ); 02223 02224 ASSERT( d1 > 0. ); 02225 ASSERT( d2 < 0. ); 02226 return result; 02227 } 02228 02229 inline double local_product( double K , long int lp ) 02230 { 02231 long int s = 0; 02232 02233 double Ksqrd =(K*K); 02234 double partprod = 1.; 02235 02236 for( s = 0; s <= lp; s = s + 1 ) 02237 { 02238 double s2 = (double)(s*s); 02239 02240 /**************************/ 02241 /* s=l' */ 02242 /* ----- */ 02243 /* | | (1 + s^2 K^2) */ 02244 /* | | */ 02245 /* s=0 */ 02246 /**************************/ 02247 02248 partprod *= ( 1. + ( s2 * Ksqrd ) ); 02249 } 02250 return partprod; 02251 } 02252 02253 /************************************************************************/ 02254 /* Find the Einstein A's for hydrogen for a */ 02255 /* transition n,l --> n',l' */ 02256 /* */ 02257 /* In the following, we will assume n > n' */ 02258 /************************************************************************/ 02259 /*******************************************************************************/ 02260 /* */ 02261 /* Einstein A() for the transition from the */ 02262 /* initial state n,l to the finial state n',l' */ 02263 /* is given by oscillator f() */ 02264 /* */ 02265 /* hbar w max(l,l') | | 2 */ 02266 /* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */ 02267 /* 3 R_oo ( 2l + 1 ) | | */ 02268 /* */ 02269 /* */ 02270 /* E(n,l;n',l') max(l,l') | | 2 */ 02271 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 02272 /* 3 R_oo ( 2l + 1 ) | | */ 02273 /* */ 02274 /* */ 02275 /* See for example Gordan Drake's */ 02276 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 02277 /* */ 02278 /* Here R_oo is the infinite mass Rydberg length */ 02279 /* */ 02280 /* */ 02281 /* h c */ 02282 /* R_oo --- = 13.605698 eV */ 02283 /* {e} */ 02284 /* */ 02285 /* */ 02286 /* R_oo = 2.179874e-11 ergs */ 02287 /* */ 02288 /* w = omega */ 02289 /* = frequency of transition from n,l to n',l' */ 02290 /* */ 02291 /* */ 02292 /* */ 02293 /* here g_k are statistical weights obtained from */ 02294 /* the appropriate angular momentum quantum numbers */ 02295 /* */ 02296 /* */ 02297 /* - - 2 */ 02298 /* 64 pi^4 (e a_o)^2 max(l,l') | | */ 02299 /* A(n,l;n',l') = ------------------- ----------- v^3 | R(n,l;n',l') | */ 02300 /* 3 h c^3 2*l + 1 | | */ 02301 /* - - */ 02302 /* */ 02303 /* */ 02304 /* pi 3.141592654 */ 02305 /* plank_hbar 6.5821220 eV sec */ 02306 /* R_oo 2.179874e-11 ergs */ 02307 /* plank_h 6.6260755e-34 J sec */ 02308 /* e_charge 1.60217733e-19 C */ 02309 /* a_o 0.529177249e-10 m */ 02310 /* vel_light_c 299792458L m sec^-1 */ 02311 /* */ 02312 /* */ 02313 /* */ 02314 /* 64 pi^4 (e a_o)^2 64 pi^4 (a_o)^2 e^2 1 1 */ 02315 /* ----------------- = ----------------- -------- ---- = 7.5197711e-38 ----- */ 02316 /* 3 h c^3 3 c^2 hbar c 2 pi sec */ 02317 /* */ 02318 /* */ 02319 /* e^2 1 */ 02320 /* using ---------- = alpha = ---- */ 02321 /* hbar c 137 */ 02322 /*******************************************************************************/ 02323 02324 double H_Einstein_A(/* IN THE FOLLOWING WE HAVE n > n' */ 02325 /* principal quantum number, 1 for ground, upper level */ 02326 long int n, 02327 /* angular momentum, 0 for s */ 02328 long int l, 02329 /* principal quantum number, 1 for ground, lower level */ 02330 long int np, 02331 /* angular momentum, 0 for s */ 02332 long int lp, 02333 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 02334 long int iz 02335 ) 02336 { 02337 DEBUG_ENTRY( "H_Einstein_A()" ); 02338 02339 double result; 02340 if( n > 60 || np > 60 ) 02341 { 02342 result = H_Einstein_A_log10(n,l,np,lp,iz ); 02343 } 02344 else 02345 { 02346 result = H_Einstein_A_lin(n,l,np,lp,iz ); 02347 } 02348 return result; 02349 } 02350 02351 /************************************************************************/ 02352 /* Calculates the Einstein A's for hydrogen */ 02353 /* for the transition n,l --> n',l' */ 02354 /* units of sec^(-1) */ 02355 /* */ 02356 /* In the following, we have n > n' */ 02357 /************************************************************************/ 02358 STATIC double H_Einstein_A_lin(/* IN THE FOLLOWING WE HAVE n > n' */ 02359 /* principal quantum number, 1 for ground, upper level */ 02360 long int n, 02361 /* angular momentum, 0 for s */ 02362 long int l, 02363 /* principal quantum number, 1 for ground, lower level */ 02364 long int np, 02365 /* angular momentum, 0 for s */ 02366 long int lp, 02367 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 02368 long int iz 02369 ) 02370 { 02371 DEBUG_ENTRY( "H_Einstein_A_lin()" ); 02372 02373 /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */ 02374 double d1 = hv( n, np, iz ); 02375 double d2 = d1 / HPLANCK; /* v = hv / h */ 02376 double d3 = pow3(d2); 02377 double lg = (double)(l > lp ? l : lp); 02378 double Two_L_Plus_One = (double)(2*l + 1); 02379 double d6 = lg / Two_L_Plus_One; 02380 double d7 = hri( n, l, np, lp , iz ); 02381 double d8 = d7 * d7; 02382 double result = CONST_ONE * d3 * d6 * d8; 02383 02384 /* validate the incoming data */ 02385 if( n >= 70 ) 02386 { 02387 fprintf(ioQQQ,"Principal Quantum Number `n' too large.\n"); 02388 cdEXIT(EXIT_FAILURE); 02389 } 02390 if( iz < 1 ) 02391 { 02392 fprintf(ioQQQ," The charge is impossible.\n"); 02393 cdEXIT(EXIT_FAILURE); 02394 } 02395 if( n < 1 || np < 1 || l >= n || lp >= np ) 02396 { 02397 fprintf(ioQQQ," The quantum numbers are impossible.\n"); 02398 cdEXIT(EXIT_FAILURE); 02399 } 02400 if( n <= np ) 02401 { 02402 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n"); 02403 cdEXIT(EXIT_FAILURE); 02404 } 02405 return result; 02406 } 02407 02408 /**********************log version****************************/ 02409 double H_Einstein_A_log10(/* returns Einstein A in units of (sec)^-1 */ 02410 long int n, 02411 long int l, 02412 long int np, 02413 long int lp, 02414 long int iz 02415 ) 02416 { 02417 DEBUG_ENTRY( "H_Einstein_A_log10()" ); 02418 02419 /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */ 02420 double d1 = hv( n, np , iz ); 02421 double d2 = d1 / HPLANCK; /* v = hv / h */ 02422 double d3 = pow3(d2); 02423 double lg = (double)(l > lp ? l : lp); 02424 double Two_L_Plus_One = (double)(2*l + 1); 02425 double d6 = lg / Two_L_Plus_One; 02426 double d7 = hri_log10( n, l, np, lp , iz ); 02427 double d8 = d7 * d7; 02428 double result = CONST_ONE * d3 * d6 * d8; 02429 02430 /* validate the incoming data */ 02431 if( iz < 1 ) 02432 { 02433 fprintf(ioQQQ," The charge is impossible.\n"); 02434 cdEXIT(EXIT_FAILURE); 02435 } 02436 if( n < 1 || np < 1 || l >= n || lp >= np ) 02437 { 02438 fprintf(ioQQQ," The quantum numbers are impossible.\n"); 02439 cdEXIT(EXIT_FAILURE); 02440 } 02441 if( n <= np ) 02442 { 02443 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n"); 02444 cdEXIT(EXIT_FAILURE); 02445 } 02446 return result; 02447 } 02448 02449 /********************************************************************************/ 02450 /* hv calculates photon energy for n -> n' transitions for H and H-like ions */ 02451 /* simplest case of no "l" or "m" dependence */ 02452 /* epsilon_0 = 1 in vacu */ 02453 /* */ 02454 /* */ 02455 /* R_h */ 02456 /* Energy(n,Z) = - ------- */ 02457 /* n^2 */ 02458 /* */ 02459 /* */ 02460 /* */ 02461 /* Friedrich -- Theoretical Atomic Physics pg. 60 eq. 2.8 */ 02462 /* */ 02463 /* u */ 02464 /* R_h = --- R_oo where */ 02465 /* m_e */ 02466 /* */ 02467 /* h c */ 02468 /* R_oo --- = 2.179874e-11 ergs */ 02469 /* e */ 02470 /* */ 02471 /* (Harmin Lecture Notes for course phy-651 Spring 1994) */ 02472 /* where m_e (m_p) is the mass of and electron (proton) */ 02473 /* and u is the reduced electron mass for neutral hydrogen */ 02474 /* */ 02475 /* */ 02476 /* m_e m_p m_e */ 02477 /* u = --------- = ----------- */ 02478 /* m_e + m_p 1 + m_e/m_p */ 02479 /* */ 02480 /* m_e */ 02481 /* Now ----- = 0.000544617013 */ 02482 /* m_p */ 02483 /* u */ 02484 /* so that --- = 0.999455679 */ 02485 /* m_e */ 02486 /* */ 02487 /* */ 02488 /* returns energy of photon in ergs */ 02489 /* */ 02490 /* hv (n,n',Z) is for transitions n -> n' */ 02491 /* */ 02492 /* 1 erg = 1e-07 J */ 02493 /********************************************************************************/ 02494 /********************************************************************************/ 02495 /* WARNING: hv() use the electron reduced mass for hydrogen instead of */ 02496 /* the reduced mass associated with the apropriate ion */ 02497 /********************************************************************************/ 02498 02499 inline double hv( long int n, long int nprime, long int iz ) 02500 { 02501 DEBUG_ENTRY( "hv()" ); 02502 02503 double n1 = (double)n; 02504 double n2 = n1*n1; 02505 double np1 = (double)nprime; 02506 double np2 = np1*np1; 02507 double rmr = 1./(1. + ELECTRON_MASS/PROTON_MASS); /* 0.999455679 */ 02508 double izsqrd = (double)(iz*iz); 02509 02510 double d1 = 1. / n2; 02511 double d2 = 1. / np2; 02512 double d3 = izsqrd * rmr * EN1RYD; 02513 double d4 = d2 - d1; 02514 double result = d3 * d4; 02515 02516 ASSERT( n > 0 ); 02517 ASSERT( nprime > 0 ); 02518 ASSERT( n > nprime ); 02519 ASSERT( iz > 0 ); 02520 ASSERT( result > 0. ); 02521 02522 if( n <= nprime ) 02523 { 02524 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n"); 02525 cdEXIT(EXIT_FAILURE); 02526 } 02527 02528 return result; 02529 } 02530 02531 /************************************************************************/ 02532 /* hri() */ 02533 /* Calculate the hydrogen radial wavefunction integral */ 02534 /* for the dipole transition l'=l-1 or l'=l+1 */ 02535 /* for the higher energy state n,l to the lower energy state n',l' */ 02536 /* no "m" dependence */ 02537 /************************************************************************/ 02538 /* here we have a transition */ 02539 /* from the higher energy state n,l */ 02540 /* to the lower energy state n',l' */ 02541 /* with a dipole selection rule on l and l' */ 02542 /************************************************************************/ 02543 /* */ 02544 /* hri() test n,l,n',l' for domain errors and */ 02545 /* swaps n,l <--> n',l' for the case l'=l+1 */ 02546 /* */ 02547 /* It then calls hrii() */ 02548 /* */ 02549 /* Dec. 6, 1999 */ 02550 /* Robert Paul Bauman */ 02551 /************************************************************************/ 02552 02553 /************************************************************************/ 02554 /* This routine, hri(), calculates the hydrogen radial integral, */ 02555 /* for the transition n,l --> n',l' */ 02556 /* It is, of course, dimensionless. */ 02557 /* */ 02558 /* In the following, we have n > n' */ 02559 /************************************************************************/ 02560 02561 inline double hri( 02562 /* principal quantum number, 1 for ground, upper level */ 02563 long int n, 02564 /* angular momentum, 0 for s */ 02565 long int l, 02566 /* principal quantum number, 1 for ground, lower level */ 02567 long int np, 02568 /* angular momentum, 0 for s */ 02569 long int lp, 02570 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 02571 long int iz 02572 ) 02573 { 02574 DEBUG_ENTRY( "hri" ); 02575 02576 long int a; 02577 long int b; 02578 long int c; 02579 long int d; 02580 double ld1 = 0.; 02581 double Z = (double)iz; 02582 02583 /**********************************************************************/ 02584 /* from higher energy -> lower energy */ 02585 /* Selection Rule for l and l' */ 02586 /* dipole process only */ 02587 /**********************************************************************/ 02588 02589 ASSERT( n > 0 ); 02590 ASSERT( np > 0 ); 02591 ASSERT( l >= 0 ); 02592 ASSERT( lp >= 0 ); 02593 ASSERT( n > l ); 02594 ASSERT( np > lp ); 02595 ASSERT( n > np || ( n == np && l == lp + 1 )); 02596 ASSERT( iz > 0 ); 02597 ASSERT( lp == l + 1 || lp == l - 1 ); 02598 02599 if( l == lp + 1 ) 02600 { 02601 /* Keep variable the same */ 02602 a = n; 02603 b = l; 02604 c = np; 02605 d = lp; 02606 } 02607 else if( l == lp - 1 ) 02608 { 02609 /* swap n,l with n',l' */ 02610 a = np; 02611 b = lp; 02612 c = n; 02613 d = l; 02614 } 02615 else 02616 { 02617 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" ); 02618 cdEXIT(EXIT_FAILURE); 02619 } 02620 02621 /**********************************************/ 02622 /* Take care of the Z-dependence here. */ 02623 /**********************************************/ 02624 ld1 = hrii(a, b, c, d ) / Z; 02625 02626 return ld1; 02627 } 02628 02629 /************************************************************************/ 02630 /* hri_log10() */ 02631 /* Calculate the hydrogen radial wavefunction integral */ 02632 /* for the dipole transition l'=l-1 or l'=l+1 */ 02633 /* for the higher energy state n,l to the lower energy state n',l' */ 02634 /* no "m" dependence */ 02635 /************************************************************************/ 02636 /* here we have a transition */ 02637 /* from the higher energy state n,l */ 02638 /* to the lower energy state n',l' */ 02639 /* with a dipole selection rule on l and l' */ 02640 /************************************************************************/ 02641 /* */ 02642 /* hri_log10() test n,l,n',l' for domain errors and */ 02643 /* swaps n,l <--> n',l' for the case l'=l+1 */ 02644 /* */ 02645 /* It then calls hrii_log() */ 02646 /* */ 02647 /* Dec. 6, 1999 */ 02648 /* Robert Paul Bauman */ 02649 /************************************************************************/ 02650 02651 inline double hri_log10( long int n, long int l, long int np, long int lp , long int iz ) 02652 { 02653 /**********************************************************************/ 02654 /* from higher energy -> lower energy */ 02655 /* Selection Rule for l and l' */ 02656 /* dipole process only */ 02657 /**********************************************************************/ 02658 02659 DEBUG_ENTRY( "hri_log10()" ); 02660 02661 long int a; 02662 long int b; 02663 long int c; 02664 long int d; 02665 double ld1 = 0.; 02666 double Z = (double)iz; 02667 02668 ASSERT( n > 0); 02669 ASSERT( np > 0); 02670 ASSERT( l >= 0); 02671 ASSERT( lp >= 0 ); 02672 ASSERT( n > l ); 02673 ASSERT( np > lp ); 02674 ASSERT( n > np || ( n == np && l == lp + 1 )); 02675 ASSERT( iz > 0 ); 02676 ASSERT( lp == l + 1 || lp == l - 1 ); 02677 02678 if( l == lp + 1) 02679 { 02680 /* Keep variable the same */ 02681 a = n; 02682 b = l; 02683 c = np; 02684 d = lp; 02685 } 02686 else if( l == lp - 1 ) 02687 { 02688 /* swap n,l with n',l' */ 02689 a = np; 02690 b = lp; 02691 c = n; 02692 d = l; 02693 } 02694 else 02695 { 02696 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" ); 02697 cdEXIT(EXIT_FAILURE); 02698 } 02699 02700 /**********************************************/ 02701 /* Take care of the Z-dependence here. */ 02702 /**********************************************/ 02703 ld1 = hrii_log(a, b, c, d ) / Z; 02704 02705 return ld1; 02706 } 02707 02708 STATIC double hrii( long int n, long int l, long int np, long int lp) 02709 { 02710 /******************************************************************************/ 02711 /* this routine hrii() is internal to the parent routine hri() */ 02712 /* this internal routine only considers the case l=l'+1 */ 02713 /* the case l=l-1 is done in the parent routine hri() */ 02714 /* by the transformation n <--> n' and l <--> l' */ 02715 /* THUS WE TEST FOR */ 02716 /* l=l'-1 */ 02717 /******************************************************************************/ 02718 02719 DEBUG_ENTRY( "hrii()" ); 02720 02721 long int a = 0, b = 0, c = 0; 02722 long int i1 = 0, i2 = 0, i3 = 0, i4 = 0; 02723 02724 char A='a'; 02725 02726 double y = 0.; 02727 double fsf = 0.; 02728 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.; 02729 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.; 02730 double d00 = 0., d01 = 0.; 02731 02732 ASSERT( l == lp + 1 ); 02733 02734 if( n == np ) /* SPECIAL CASE 1 */ 02735 { 02736 /**********************************************************/ 02737 /* if lp= l + 1 then it has higher energy */ 02738 /* i.e. no photon */ 02739 /* this is the second time we check this, oh well */ 02740 /**********************************************************/ 02741 02742 if( lp != (l - 1) ) 02743 { 02744 printf( "BadMagic: Energy requirements not met.\n\n" ); 02745 cdEXIT(EXIT_FAILURE); 02746 } 02747 02748 d2 = 3. / 2.; 02749 i1 = n * n; 02750 i2 = l * l; 02751 d5 = (double)(i1 - i2); 02752 d6 = sqrt(d5); 02753 d7 = (double)n * d6; 02754 d8 = d2 * d7; 02755 return d8; 02756 } 02757 else if( l == np && lp == (l - 1) ) /* A Pair of Easy Special Cases */ 02758 { 02759 if( l == (n - 1) ) 02760 { 02761 /**********************************************************************/ 02762 /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */ 02763 /* */ 02764 /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */ 02765 /* [(2n-1) - 1/(2n-1)]/4 */ 02766 /**********************************************************************/ 02767 02768 d1 = (double)( 2*n - 2 ); 02769 d2 = (double)( 2*n - 1 ); 02770 d3 = d1 * d2; 02771 d4 = sqrt( d3 ); 02772 02773 d5 = (double)(4 * n * (n - 1)); 02774 i1 = (2*n - 1); 02775 d6 = (double)(i1 * i1); 02776 d7 = d5/ d6; 02777 d8 = powi( d7, n ); 02778 02779 d9 = 1./d2; 02780 d10 = d2 - d9; 02781 d11 = d10 / 4.; 02782 02783 /* Wrap it all up */ 02784 02785 d12 = d4 * d8 * d11; 02786 return d12; 02787 02788 } 02789 else 02790 { 02791 /******************************************************************************/ 02792 /* R(n,l;n',l') = R(n,l;l,l-1) */ 02793 /* */ 02794 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */ 02795 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */ 02796 /******************************************************************************/ 02797 02798 d2 = 1.; 02799 for( i1 = -l; i1 <= l; i1 = i1 + 1 ) /* from n-l to n+l INCLUSIVE */ 02800 { 02801 d1 = (double)(n - i1); 02802 d2 = d2 * d1; 02803 } 02804 i2 = (2*l - 1); 02805 d3 = factorial( i2 ); 02806 d4 = d2/d3; 02807 d4 = sqrt( d4 ); 02808 02809 02810 d5 = (double)( 4. * n * l ); 02811 i3 = (n - l); 02812 d6 = (double)( i3 * i3 ); 02813 d7 = d5 / d6; 02814 d8 = powi( d7, l+1 ); 02815 02816 02817 i4 = n + l; 02818 d9 = (double)( i3 ) / (double)( i4 ); 02819 d10 = powi( d9 , i4 ); 02820 02821 d11 = d9 * d9; 02822 d12 = 1. - d11; 02823 d13 = d12 / 4.; 02824 02825 /* Wrap it all up */ 02826 d14 = d4 * d8 * d10 * d13; 02827 return d14; 02828 } 02829 } 02830 02831 /*******************************************************************************************/ 02832 /* THE GENERAL CASE */ 02833 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 02834 /* REF: D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 02835 /* For F(a,b;c;x) we have from eq.4 */ 02836 /* */ 02837 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 02838 /* */ 02839 /* a (1-x) (a + bx - c) */ 02840 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 02841 /* (a-c) (a-c) */ 02842 /* */ 02843 /* */ 02844 /* A similiar recusion relation holds for b with a <--> b. */ 02845 /* */ 02846 /* */ 02847 /* we have initial conditions */ 02848 /* */ 02849 /* */ 02850 /* F(0) = 1 with a = -1 */ 02851 /* */ 02852 /* b */ 02853 /* F(-1) = 1 - (---) x with a = -1 */ 02854 /* c */ 02855 /*******************************************************************************************/ 02856 02857 if( lp == l - 1 ) /* use recursion over "b" */ 02858 { 02859 A='b'; 02860 } 02861 else if( lp == l + 1 ) /* use recursion over "a" */ 02862 { 02863 A='a'; 02864 } 02865 else 02866 { 02867 printf(" BadMagic: Don't know what to do here.\n\n"); 02868 cdEXIT(EXIT_FAILURE); 02869 } 02870 02871 /********************************************************************/ 02872 /* Calculate the whole shootin match */ 02873 /* - - (1/2) */ 02874 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 02875 /* ----------- * | ---------------- | */ 02876 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 02877 /* - - */ 02878 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 02879 /* */ 02880 /* This is used in the calculation of hydrogen */ 02881 /* radial wave function integral for dipole transition case */ 02882 /********************************************************************/ 02883 02884 fsf = fsff( n, l, np ); 02885 02886 /**************************************************************************************/ 02887 /* Use a -> a' + 1 */ 02888 /* _ _ */ 02889 /* (a' + 1) (1 - x) | | */ 02890 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */ 02891 /* (a' + 1 -c) | | */ 02892 /* - - */ 02893 /* */ 02894 /* For the first F() in the solution of the radial integral */ 02895 /* */ 02896 /* a = ( -n + l + 1 ) */ 02897 /* */ 02898 /* a = -n + l + 1 */ 02899 /* max(a) = max(-n) + max(l) + 1 */ 02900 /* = -n + max(n-1) + 1 */ 02901 /* = -n + n-1 + 1 */ 02902 /* = 0 */ 02903 /* */ 02904 /* similiarly */ 02905 /* */ 02906 /* min(a) = min(-n) + min(l) + 1 */ 02907 /* = min(-n) + 0 + 1 */ 02908 /* = -n + 1 */ 02909 /* */ 02910 /* a -> a' + 1 implies */ 02911 /* */ 02912 /* max(a') = -1 */ 02913 /* min(a') = -n */ 02914 /**************************************************************************************/ 02915 02916 /* a plus */ 02917 a = (-n + l + 1); 02918 02919 /* for the first 2_F_1 we use b = (-n' + l) */ 02920 b = (-np + l); 02921 02922 /* c is simple */ 02923 c = 2 * l; 02924 02925 /* -4 nn' */ 02926 /* where Y = -------- . */ 02927 /* (n-n')^2 */ 02928 d2 = (double)(n - np); 02929 d3 = d2 * d2; 02930 d4 = 1. / d3; 02931 d5 = (double)(n * np); 02932 d6 = d5 * 4.; 02933 d7 = - d6; 02934 y = d7 * d4; 02935 02936 d00 = F21( a, b, c, y, A ); 02937 02938 /**************************************************************/ 02939 /* For the second F() in the solution of the radial integral */ 02940 /* */ 02941 /* a = ( -n + l - 1 ) */ 02942 /* */ 02943 /* a = -n + l + 1 */ 02944 /* max(a) = max(-n) + max(l) - 1 */ 02945 /* = -n + (n - 1) - 1 */ 02946 /* = -2 */ 02947 /* */ 02948 /* similiarly */ 02949 /* */ 02950 /* min(a) = min(-n) + min(l) - 1 */ 02951 /* = (-n) + 0 - 1 */ 02952 /* = -n - 1 */ 02953 /* */ 02954 /* a -> a' + 1 implies */ 02955 /* */ 02956 /* max(a') = -3 */ 02957 /* */ 02958 /* min(a') = -n - 2 */ 02959 /**************************************************************/ 02960 02961 /* a minus */ 02962 a = (-n + l - 1); 02963 02964 /* for the first 2_F_1 we use b = (-n' + l) */ 02965 /* and does not change */ 02966 b = (-np + l); 02967 02968 /* c is simple */ 02969 c = 2 * l; 02970 02971 /**************************************************************/ 02972 /* -4 nn' */ 02973 /* where Y = -------- . */ 02974 /* (n-n')^2 */ 02975 /**************************************************************/ 02976 02977 /**************************************************************/ 02978 /* These are already calculated a few lines up */ 02979 /* */ 02980 /* d2 = (double) (n - np); */ 02981 /* d3 = d2 * d2; */ 02982 /* d4 = 1/ d3; */ 02983 /* d5 = (double) (n * np); */ 02984 /* d6 = d5 * 4.0; */ 02985 /* d7 = - d6; */ 02986 /* y = d7 * d4; */ 02987 /**************************************************************/ 02988 02989 d01 = F21(a, b, c, y, A ); 02990 02991 /* Calculate */ 02992 /* */ 02993 /* (n-n')^2 */ 02994 /* -------- */ 02995 /* (n+n')^2 */ 02996 02997 i1 = (n - np); 02998 d1 = pow2( (double)i1 ); 02999 i2 = (n + np); 03000 d2 = pow2( (double)i2 ); 03001 d3 = d1 / d2; 03002 03003 d4 = d01 * d3; 03004 d5 = d00 - d4; 03005 d6 = fsf * d5; 03006 03007 ASSERT( d6 != 0. ); 03008 return d6; 03009 } 03010 03011 03012 STATIC double hrii_log( long int n, long int l, long int np, long int lp) 03013 { 03014 /******************************************************************************/ 03015 /* this routine hrii_log() is internal to the parent routine hri_log10() */ 03016 /* this internal routine only considers the case l=l'+1 */ 03017 /* the case l=l-1 is done in the parent routine hri_log10() */ 03018 /* by the transformation n <--> n' and l <--> l' */ 03019 /* THUS WE TEST FOR */ 03020 /* l=l'-1 */ 03021 /******************************************************************************/ 03022 /**************************************************************************************/ 03023 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */ 03024 /* */ 03025 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */ 03026 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */ 03027 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */ 03028 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03029 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */ 03030 /**************************************************************************************/ 03031 03032 DEBUG_ENTRY( "hrii_log()" ); 03033 03034 char A='a'; 03035 03036 double y = 0.; 03037 double log10_fsf = 0.; 03038 03039 ASSERT( l == lp + 1 ); 03040 03041 if( n == np ) /* SPECIAL CASE 1 */ 03042 { 03043 /**********************************************************/ 03044 /* if lp= l + 1 then it has higher energy */ 03045 /* i.e. no photon */ 03046 /* this is the second time we check this, oh well */ 03047 /**********************************************************/ 03048 03049 if( lp != (l - 1) ) 03050 { 03051 printf( "BadMagic: l'= l+1 for n'= n.\n\n" ); 03052 cdEXIT(EXIT_FAILURE); 03053 } 03054 else 03055 { 03056 /**********************************************************/ 03057 /* 3 */ 03058 /* R(nl:n'=n,l'=l+1) = --- n sqrt( n^2 - l^2 ) */ 03059 /* 2 */ 03060 /**********************************************************/ 03061 03062 long int i1 = n * n; 03063 long int i2 = l * l; 03064 03065 double d1 = 3. / 2.; 03066 double d2 = (double)n; 03067 double d3 = (double)(i1 - i2); 03068 double d4 = sqrt(d3); 03069 double result = d1 * d2 * d4; 03070 03071 ASSERT( d3 >= 0. ); 03072 return result; 03073 } 03074 } 03075 else if( l == np && lp == l - 1 ) /* A Pair of Easy Special Cases */ 03076 { 03077 if( l == n - 1 ) 03078 { 03079 /**********************************************************************/ 03080 /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */ 03081 /* */ 03082 /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */ 03083 /* [(2n-1) - 1/(2n-1)]/4 */ 03084 /**********************************************************************/ 03085 03086 double d1 = (double)((2*n-2)*(2*n-1)); 03087 double d2 = sqrt( d1 ); 03088 double d3 = (double)(4*n*(n-1)); 03089 double d4 = (double)(2*n-1); 03090 double d5 = d4*d4; 03091 double d7 = d3/d5; 03092 double d8 = powi(d7,n); 03093 double d9 = 1./d4; 03094 double d10 = d4 - d9; 03095 double d11 = 0.25*d10; 03096 double result = (d2 * d8 * d11); /* Wrap it all up */ 03097 03098 ASSERT( d1 >= 0. ); 03099 ASSERT( d3 >= 0. ); 03100 return result; 03101 } 03102 else 03103 { 03104 double result = 0.; 03105 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0., ld7 = 0.; 03106 03107 /******************************************************************************/ 03108 /* R(n,l;n',l') = R(n,l;l,l-1) */ 03109 /* */ 03110 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */ 03111 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */ 03112 /******************************************************************************/ 03113 /**************************************/ 03114 /* [(n-l) ... (n+l)] */ 03115 /**************************************/ 03116 /* log10[(n-l) ... (n+l)] = */ 03117 /* */ 03118 /* n+l */ 03119 /* --- */ 03120 /* > log10(j) */ 03121 /* --- */ 03122 /* j=n-l */ 03123 /**************************************/ 03124 03125 ld1 = 0.; 03126 for( long int i1 = (n-l); i1 <= (n+l); i1++ ) /* from n-l to n+l INCLUSIVE */ 03127 { 03128 double d1 = (double)(i1); 03129 ld1 += log10( d1 ); 03130 } 03131 03132 /**************************************/ 03133 /* (2l-1)! */ 03134 /**************************************/ 03135 /* log10[ (2n-1)! ] */ 03136 /**************************************/ 03137 03138 ld2 = lfactorial( 2*l - 1 ); 03139 03140 ASSERT( ((2*l)+1) >= 0); 03141 03142 /**********************************************/ 03143 /* log10( [(n-l) ... (n+l)/(2l-1)!]^(1/2) ) = */ 03144 /* (1/2) log10[(n-l) ... (n+l)] - */ 03145 /* (1/2) log10[ (2n-1)! ] */ 03146 /**********************************************/ 03147 03148 ld3 = 0.5 * (ld1 - ld2); 03149 03150 /**********************************************/ 03151 /* [4nl/(n-l)^2]^(l+1) */ 03152 /**********************************************/ 03153 /* log10( [4nl/(n-l)^2]^(l+1) ) = */ 03154 /* (l+1) * log10( [4nl/(n-l)^2] ) */ 03155 /* */ 03156 /* = (l+1)*[ log10(4nl) - 2 log10(n-l) ] */ 03157 /* */ 03158 /**********************************************/ 03159 03160 double d1 = (double)(l+1); 03161 double d2 = (double)(4*n*l); 03162 double d3 = (double)(n-l); 03163 double d4 = log10(d2); 03164 double d5 = log10(d3); 03165 03166 ld4 = d1 * (d4 - 2.*d5); 03167 03168 /**********************************************/ 03169 /* [(n-l)/(n+l)]^(n+l) */ 03170 /**********************************************/ 03171 /* log10( [ (n-l)/(n+l) ]^(n+l) ) = */ 03172 /* */ 03173 /* (n+l) * [ log10(n-l) - log10(n+l) ] */ 03174 /* */ 03175 /**********************************************/ 03176 03177 d1 = (double)(n-l); 03178 d2 = (double)(n+l); 03179 d3 = log10( d1 ); 03180 d4 = log10( d2 ); 03181 03182 ld5 = d2 * (d3 - d4); 03183 03184 /**********************************************/ 03185 /* {1-[(n-l)/(n+l)]^2}/4 */ 03186 /**********************************************/ 03187 /* log10[ {1-[(n-l)/(n+l)]^2}/4 ] */ 03188 /**********************************************/ 03189 03190 d1 = (double)(n-l); 03191 d2 = (double)(n+l); 03192 d3 = d1/d2; 03193 d4 = d3*d3; 03194 d5 = 1. - d4; 03195 double d6 = 0.25*d5; 03196 03197 ld6 = log10(d6); 03198 03199 /******************************************************************************/ 03200 /* R(n,l;n',l') = R(n,l;l,l-1) */ 03201 /* */ 03202 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */ 03203 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */ 03204 /******************************************************************************/ 03205 03206 ld7 = ld3 + ld4 + ld5 + ld6; 03207 03208 result = pow( 10., ld7 ); 03209 03210 ASSERT( result > 0. ); 03211 return result; 03212 } 03213 } 03214 else 03215 { 03216 double result = 0.; 03217 long int a = 0, b = 0, c = 0; 03218 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.; 03219 mx d00={0.0,0}, d01={0.0,0}, d02={0.0,0}, d03={0.0,0}; 03220 03221 if( lp == l - 1 ) /* use recursion over "b" */ 03222 { 03223 A='b'; 03224 } 03225 else if( lp == l + 1 ) /* use recursion over "a" */ 03226 { 03227 A='a'; 03228 } 03229 else 03230 { 03231 printf(" BadMagic: Don't know what to do here.\n\n"); 03232 cdEXIT(EXIT_FAILURE); 03233 } 03234 03235 /**************************************************************************************/ 03236 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */ 03237 /* */ 03238 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */ 03239 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */ 03240 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */ 03241 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03242 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */ 03243 /**************************************************************************************/ 03244 03245 /****************************************************************************************************/ 03246 /* Calculate the whole shootin match */ 03247 /* - - (1/2) */ 03248 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 03249 /* fsff() = ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2)(n+n')^(-n-n') */ 03250 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03251 /* - - */ 03252 /* This is used in the calculation of hydrogen radial wave function integral for dipole transitions */ 03253 /****************************************************************************************************/ 03254 03255 log10_fsf = log10_fsff( n, l, np ); 03256 03257 /******************************************************************************************/ 03258 /* 2_F_1( a, b; c; y ) */ 03259 /* */ 03260 /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */ 03261 /* */ 03262 /* */ 03263 /* Use a -> a' + 1 */ 03264 /* _ _ */ 03265 /* (a' + 1) (1 - x) | | */ 03266 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */ 03267 /* (a' + 1 - c) | | */ 03268 /* - - */ 03269 /* */ 03270 /* For the first F() in the solution of the radial integral */ 03271 /* */ 03272 /* a = ( -n + l + 1 ) */ 03273 /* */ 03274 /* a = -n + l + 1 */ 03275 /* max(a) = max(-n) + max(l) + 1 */ 03276 /* = -n + max(n-1) + 1 */ 03277 /* = -n + n-1 + 1 */ 03278 /* = 0 */ 03279 /* */ 03280 /* similiarly */ 03281 /* */ 03282 /* min(a) = min(-n) + min(l) + 1 */ 03283 /* = min(-n) + 0 + 1 */ 03284 /* = -n + 1 */ 03285 /* */ 03286 /* a -> a' + 1 implies */ 03287 /* */ 03288 /* max(a') = -1 */ 03289 /* min(a') = -n */ 03290 /******************************************************************************************/ 03291 03292 /* a plus */ 03293 a = (-n + l + 1); 03294 03295 /* for the first 2_F_1 we use b = (-n' + l) */ 03296 b = (-np + l); 03297 03298 /* c is simple */ 03299 c = 2 * l; 03300 03301 /**********************************************************************************/ 03302 /* 2_F_1( a, b; c; y ) */ 03303 /* */ 03304 /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */ 03305 /* */ 03306 /* -4 nn' */ 03307 /* where Y = -------- . */ 03308 /* (n-n')^2 */ 03309 /* */ 03310 /**********************************************************************************/ 03311 03312 d2 = (double)(n - np); 03313 d3 = d2 * d2; 03314 03315 d4 = 1. / d3; 03316 d5 = (double)(n * np); 03317 d6 = d5 * 4.; 03318 03319 d7 = -d6; 03320 y = d7 * d4; 03321 03322 /**************************************************************************************************/ 03323 /* THE GENERAL CASE */ 03324 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 03325 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 03326 /* */ 03327 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 03328 /* */ 03329 /* a (1-x) (a + bx - c) */ 03330 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 03331 /* (a - c) (a - c) */ 03332 /* */ 03333 /* */ 03334 /* A similiar recusion relation holds for b with a <--> b. */ 03335 /* */ 03336 /* */ 03337 /* we have initial conditions */ 03338 /* */ 03339 /* */ 03340 /* F(0) = 1 with a = -1 */ 03341 /* */ 03342 /* b */ 03343 /* F(-1) = 1 - (---) x with a = -1 */ 03344 /* c */ 03345 /**************************************************************************************************/ 03346 03347 /* 2_F_1( long int a, long int b, long int c, (double) y, (string) "a" or "b") */ 03348 /* F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03349 d00 = F21_mx( a, b, c, y, A ); 03350 03351 /**************************************************************/ 03352 /* For the second F() in the solution of the radial integral */ 03353 /* */ 03354 /* a = ( -n + l - 1 ) */ 03355 /* */ 03356 /* a = -n + l + 1 */ 03357 /* max(a) = max(-n) + max(l) - 1 */ 03358 /* = -n + (n - 1) - 1 */ 03359 /* = -2 */ 03360 /* */ 03361 /* similiarly */ 03362 /* */ 03363 /* min(a) = min(-n) + min(l) - 1 */ 03364 /* = (-n) + 0 - 1 */ 03365 /* = -n - 1 */ 03366 /* */ 03367 /* a -> a' + 1 implies */ 03368 /* */ 03369 /* max(a') = -3 */ 03370 /* */ 03371 /* min(a') = -n - 2 */ 03372 /**************************************************************/ 03373 03374 /* a minus */ 03375 a = (-n + l - 1); 03376 03377 /* for the first 2_F_1 we use b = (-n' + l) */ 03378 /* and does not change */ 03379 b = (-np + l); 03380 03381 /* c is simple */ 03382 c = 2 * l; 03383 03384 /**************************************************************/ 03385 /* -4 nn' */ 03386 /* where Y = -------- . */ 03387 /* (n-n')^2 */ 03388 /**************************************************************/ 03389 03390 /**************************************************************/ 03391 /* These are already calculated a few lines up */ 03392 /* */ 03393 /* d2 = (double) (n - np); */ 03394 /* d3 = d2 * d2; */ 03395 /* d4 = 1/ d3; */ 03396 /* d5 = (double) (n * np); */ 03397 /* d6 = d5 * 4.0; */ 03398 /* d7 = - d6; */ 03399 /* y = d7 * d4; */ 03400 /**************************************************************/ 03401 03402 d01 = F21_mx(a, b, c, y, A ); 03403 03404 /**************************************************************************************/ 03405 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */ 03406 /* */ 03407 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */ 03408 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */ 03409 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */ 03410 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */ 03411 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */ 03412 /* */ 03413 /* = fsf * ( F(a,b,c;y) - d3 * F(a',b',c';y) ) */ 03414 /* */ 03415 /* where d3 = (n-n')^2 (n+n')^2 */ 03416 /* */ 03417 /**************************************************************************************/ 03418 03419 /**************************************************************/ 03420 /* Calculate */ 03421 /* */ 03422 /* (n-n')^2 */ 03423 /* -------- */ 03424 /* (n+n')^2 */ 03425 /**************************************************************/ 03426 03427 d1 = (double)((n - np)*(n -np)); 03428 d2 = (double)((n + np)*(n + np)); 03429 d3 = d1 / d2; 03430 03431 d02.x = d01.x; 03432 d02.m = d01.m * d3; 03433 03434 while ( fabs(d02.m) > 1.0e+25 ) 03435 { 03436 d02.m /= 1.0e+25; 03437 d02.x += 25; 03438 } 03439 03440 d03.x = d00.x; 03441 d03.m = d00.m * (1. - (d02.m/d00.m) * powi( 10. , (d02.x - d00.x) ) ); 03442 03443 result = pow( 10., (log10_fsf + d03.x) ) * d03.m; 03444 03445 ASSERT( result != 0. ); 03446 03447 /* we don't care about the correct sign of result... */ 03448 return fabs(result); 03449 } 03450 /* Shouldn't get here */ 03451 printf(" This code should be inaccessible\n\n"); 03452 cdEXIT(EXIT_FAILURE); 03453 } 03454 03455 STATIC double fsff( long int n, long int l, long int np ) 03456 { 03457 /****************************************************************/ 03458 /* Calculate the whole shootin match */ 03459 /* - - (1/2) */ 03460 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 03461 /* ----------- * | ---------------- | */ 03462 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03463 /* - - */ 03464 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03465 /* */ 03466 /****************************************************************/ 03467 03468 DEBUG_ENTRY( "fsff()" ); 03469 03470 long int i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0; 03471 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.; 03472 double sigma = 1.; 03473 03474 /**************************************************************** 03475 * Calculate the whole shootin match * 03476 * (-1)^(n'-1) | (n+l)! (n'+l-1)! | * 03477 * ----------- * | ---------------- | * 03478 * [4 (2l-1)!] | (n-l-1)! (n'-l)! | * 03479 * - - * 03480 * * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * 03481 * * 03482 ****************************************************************/ 03483 03484 /* Calculate (-1)^(n'-l) */ 03485 if( is_odd(np - l) ) 03486 { 03487 sigma *= -1.; 03488 } 03489 ASSERT( sigma != 0. ); 03490 03491 /*********************/ 03492 /* Calculate (2l-1)! */ 03493 /*********************/ 03494 i1 = (2*l - 1); 03495 if( i1 < 0 ) 03496 { 03497 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03498 cdEXIT(EXIT_FAILURE); 03499 } 03500 03501 /****************************************************************/ 03502 /* Calculate the whole shootin match */ 03503 /* - - (1/2) */ 03504 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */ 03505 /* ----------- * | ---------------- | */ 03506 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03507 /* - - */ 03508 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03509 /* */ 03510 /****************************************************************/ 03511 03512 d0 = factorial( i1 ); 03513 d1 = 4. * d0; 03514 d2 = 1. / d1; 03515 03516 /**********************************************************************/ 03517 /* We want the (negitive) of this */ 03518 /* since we really are interested in */ 03519 /* [(2l-1)!]^-1 */ 03520 /**********************************************************************/ 03521 03522 sigma = sigma * d2; 03523 ASSERT( sigma != 0. ); 03524 03525 /**********************************************************************/ 03526 /* Calculate (4 n n')^(l+1) */ 03527 /* powi( m , n) calcs m^n */ 03528 /* returns long double with m,n ints */ 03529 /**********************************************************************/ 03530 03531 i0 = 4 * n * np; 03532 i1 = l + 1; 03533 d2 = powi( (double)i0 , i1 ); 03534 sigma = sigma * d2; 03535 ASSERT( sigma != 0. ); 03536 03537 /* Calculate (n-n')^(n+n'-2l-2) */ 03538 i0 = n - np; 03539 i1 = n + np - 2*l - 2; 03540 d2 = powi( (double)i0 , i1 ); 03541 sigma = sigma * d2; 03542 ASSERT( sigma != 0. ); 03543 03544 /* Calculate (n+n')^(-n-n') */ 03545 i0 = n + np; 03546 i1 = -n - np; 03547 d2 = powi( (double)i0 , i1 ); 03548 sigma = sigma * d2; 03549 ASSERT( sigma != 0. ); 03550 03551 /**********************************************************************/ 03552 /* - - (1/2) */ 03553 /* | (n+l)! (n'+l-1)! | */ 03554 /* Calculate | ---------------- | */ 03555 /* | (n-l-1)! (n'-l)! | */ 03556 /* - - */ 03557 /**********************************************************************/ 03558 03559 i1 = n + l; 03560 if( i1 < 0 ) 03561 { 03562 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03563 cdEXIT(EXIT_FAILURE); 03564 } 03565 d1 = factorial( i1 ); 03566 03567 i2 = np + l - 1; 03568 if( i2 < 0 ) 03569 { 03570 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03571 cdEXIT(EXIT_FAILURE); 03572 } 03573 d2 = factorial( i2 ); 03574 03575 i3 = n - l - 1; 03576 if( i3 < 0 ) 03577 { 03578 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03579 cdEXIT(EXIT_FAILURE); 03580 } 03581 d3 = factorial( i3 ); 03582 03583 i4 = np - l; 03584 if( i4 < 0 ) 03585 { 03586 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" ); 03587 cdEXIT(EXIT_FAILURE); 03588 } 03589 d4 = factorial( i4 ); 03590 03591 ASSERT( d3 != 0. ); 03592 ASSERT( d4 != 0. ); 03593 03594 /* Do this a different way to prevent overflow */ 03595 /* d5 = (sqrt(d1 *d2)); */ 03596 d5 = sqrt(d1)*sqrt(d2); 03597 d5 /= sqrt(d3); 03598 d5 /= sqrt(d4); 03599 03600 sigma = sigma * d5; 03601 03602 ASSERT( sigma != 0. ); 03603 return sigma; 03604 } 03605 03606 /**************************log version*******************************/ 03607 STATIC double log10_fsff( long int n, long int l, long int np ) 03608 { 03609 /******************************************************************************************************/ 03610 /* Calculate the whole shootin match */ 03611 /* - - (1/2) */ 03612 /* 1 | (n+l)! (n'+l-1)! | */ 03613 /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03614 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03615 /* - - */ 03616 /******************************************************************************************************/ 03617 03618 DEBUG_ENTRY( "log10_fsff()" ); 03619 03620 double d0 = 0., d1 = 0.; 03621 double ld0 = 0., ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0.; 03622 double result = 0.; 03623 03624 /******************************************************************************************************/ 03625 /* Calculate the log10 of the whole shootin match */ 03626 /* - - (1/2) */ 03627 /* 1 | (n+l)! (n'+l-1)! | */ 03628 /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */ 03629 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03630 /* - - */ 03631 /******************************************************************************************************/ 03632 03633 /**********************/ 03634 /* Calculate (2l-1)! */ 03635 /**********************/ 03636 03637 d0 = (double)(2*l - 1); 03638 ASSERT( d0 != 0. ); 03639 03640 /******************************************************************************************************/ 03641 /* Calculate the whole shootin match */ 03642 /* - - (1/2) */ 03643 /* 1 | (n+l)! (n'+l-1)! | */ 03644 /* ----------- * | ---------------- | * (4 n n')^(l+1) |(n-n')^(n+n'-2l-2)| (n+n')^(-n-n') */ 03645 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */ 03646 /* - - */ 03647 /******************************************************************************************************/ 03648 03649 ld0 = lfactorial( 2*l - 1 ); 03650 ld1 = log10(4.); 03651 result = -(ld0 + ld1); 03652 ASSERT( result != 0. ); 03653 03654 /**********************************************************************/ 03655 /* Calculate (4 n n')^(l+1) */ 03656 /* powi( m , n) calcs m^n */ 03657 /* returns long double with m,n ints */ 03658 /**********************************************************************/ 03659 03660 d0 = (double)(4 * n * np); 03661 d1 = (double)(l + 1); 03662 result += d1 * log10(d0); 03663 ASSERT( d0 >= 0. ); 03664 ASSERT( d1 != 0. ); 03665 03666 /**********************************************************************/ 03667 /* Calculate |(n-n')^(n+n'-2l-2)| */ 03668 /* NOTE: Here we are interested only */ 03669 /* magnitude of (n-n')^(n+n'-2l-2) */ 03670 /**********************************************************************/ 03671 03672 d0 = (double)(n - np); 03673 d1 = (double)(n + np - 2*l - 2); 03674 result += d1 * log10(fabs(d0)); 03675 ASSERT( fabs(d0) > 0. ); 03676 ASSERT( d1 != 0. ); 03677 03678 /* Calculate (n+n')^(-n-n') */ 03679 d0 = (double)(n + np); 03680 d1 = (double)(-n - np); 03681 result += d1 * log10(d0); 03682 ASSERT( d0 > 0. ); 03683 ASSERT( d1 != 0. ); 03684 03685 /**********************************************************************/ 03686 /* - - (1/2) */ 03687 /* | (n+l)! (n'+l-1)! | */ 03688 /* Calculate | ---------------- | */ 03689 /* | (n-l-1)! (n'-l)! | */ 03690 /* - - */ 03691 /**********************************************************************/ 03692 03693 ASSERT( n+l > 0 ); 03694 ld0 = lfactorial( n + l ); 03695 03696 ASSERT( np+l-1 > 0 ); 03697 ld1 = lfactorial( np + l - 1 ); 03698 03699 ASSERT( n-l-1 >= 0 ); 03700 ld2 = lfactorial( n - l - 1 ); 03701 03702 ASSERT( np-l >= 0 ); 03703 ld3 = lfactorial( np - l ); 03704 03705 ld4 = 0.5*((ld0+ld1)-(ld2+ld3)); 03706 03707 result += ld4; 03708 ASSERT( result != 0. ); 03709 return result; 03710 } 03711 03712 /***************************************************************************/ 03713 /* Find the Oscillator Strength for hydrogen for any */ 03714 /* transition n,l --> n',l' */ 03715 /* returns a double */ 03716 /***************************************************************************/ 03717 03718 /************************************************************************/ 03719 /* Find the Oscillator Strength for hydrogen for any */ 03720 /* transition n,l --> n',l' */ 03721 /* returns a double */ 03722 /* */ 03723 /* Einstein A() for the transition from the */ 03724 /* initial state n,l to the finial state n',l' */ 03725 /* require the Oscillator Strength f() */ 03726 /* */ 03727 /* hbar w max(l,l') | | 2 */ 03728 /* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */ 03729 /* 3 R_oo ( 2l + 1 ) | | */ 03730 /* */ 03731 /* */ 03732 /* */ 03733 /* E(n,l;n',l') max(l,l') | | 2 */ 03734 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 03735 /* 3 R_oo ( 2l + 1 ) | | */ 03736 /* */ 03737 /* */ 03738 /* See for example Gordan Drake's */ 03739 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 03740 /* */ 03741 /* Here R_oo is the infinite mass Rydberg length */ 03742 /* */ 03743 /* */ 03744 /* h c */ 03745 /* R_oo --- = 13.605698 eV */ 03746 /* {e} */ 03747 /* */ 03748 /* */ 03749 /* R_oo = 2.179874e-11 ergs */ 03750 /* */ 03751 /* w = omega */ 03752 /* = frequency of transition from n,l to n',l' */ 03753 /* */ 03754 /* */ 03755 /* */ 03756 /* here g_k are statistical weights obtained from */ 03757 /* the appropriate angular momentum quantum numbers */ 03758 /************************************************************************/ 03759 03760 /********************************************************************************/ 03761 /* Calc the Oscillator Strength f(*) given by */ 03762 /* */ 03763 /* E(n,l;n',l') max(l,l') | | 2 */ 03764 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 03765 /* 3 R_oo ( 2l + 1 ) | | */ 03766 /* */ 03767 /* See for example Gordan Drake's */ 03768 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 03769 /********************************************************************************/ 03770 03771 /************************************************************************/ 03772 /* Calc the Oscillator Strength f(*) given by */ 03773 /* */ 03774 /* E(n,l;n',l') max(l,l') | | 2 */ 03775 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */ 03776 /* 3 R_oo ( 2l + 1 ) | | */ 03777 /* */ 03778 /* f(n,l;n',l') is dimensionless. */ 03779 /* */ 03780 /* See for example Gordan Drake's */ 03781 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */ 03782 /* */ 03783 /* In the following, we have n > n' */ 03784 /************************************************************************/ 03785 03786 inline double OscStr_f( 03787 /* IN THE FOLLOWING WE HAVE n > n' */ 03788 /* principal quantum number, 1 for ground, upper level */ 03789 long int n, 03790 /* angular momentum, 0 for s */ 03791 long int l, 03792 /* principal quantum number, 1 for ground, lower level */ 03793 long int np, 03794 /* angular momentum, 0 for s */ 03795 long int lp, 03796 /* Nuclear charge, 1 for H+, 2 for He++, etc */ 03797 long int iz 03798 ) 03799 { 03800 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.; 03801 long int i1 = 0, i2 = 0; 03802 03803 if( l > lp ) 03804 i1 = l; 03805 else 03806 i1 = lp; 03807 03808 i2 = 2*lp + 1; 03809 d0 = 1. / 3.; 03810 d1 = (double)i1 / (double)i2; 03811 /* hv() returns energy in ergs */ 03812 d2 = hv( n, np, iz ); 03813 d3 = d2 / EN1RYD; 03814 d4 = hri( n, l, np, lp ,iz ); 03815 d5 = d4 * d4; 03816 03817 d6 = d0 * d1 * d3 * d5; 03818 03819 return d6; 03820 } 03821 03822 /************************log version ***************************/ 03823 inline double OscStr_f_log10( long int n , long int l , long int np , long int lp , long int iz ) 03824 { 03825 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.; 03826 long int i1 = 0, i2 = 0; 03827 03828 if( l > lp ) 03829 i1 = l; 03830 else 03831 i1 = lp; 03832 03833 i2 = 2*lp + 1; 03834 d0 = 1. / 3.; 03835 d1 = (double)i1 / (double)i2; 03836 /* hv() returns energy in ergs */ 03837 d2 = hv( n, np, iz ); 03838 d3 = d2 / EN1RYD; 03839 d4 = hri_log10( n, l, np, lp ,iz ); 03840 d5 = d4 * d4; 03841 03842 d6 = d0 * d1 * d3 * d5; 03843 03844 return d6; 03845 } 03846 03847 STATIC double F21( long int a , long int b, long int c, double y, char A ) 03848 { 03849 DEBUG_ENTRY( "F21()" ); 03850 03851 double d1 = 0.; 03852 long int i1 = 0; 03853 double *yV; 03854 03855 /**************************************************************/ 03856 /* A must be either 'a' or 'b' */ 03857 /* and is use to determine which */ 03858 /* variable recursion will be over */ 03859 /**************************************************************/ 03860 03861 ASSERT( A == 'a' || A == 'b' ); 03862 03863 if( A == 'b' ) 03864 { 03865 /* if the recursion is over "b" */ 03866 /* then make it over "a" by switching these around */ 03867 i1 = a; 03868 a = b; 03869 b = i1; 03870 A = 'a'; 03871 } 03872 03873 /**************************************************************************************/ 03874 /* malloc space for the (dynamic) 1-d array */ 03875 /* 2_F_1 works via recursion and needs room to store intermediate results */ 03876 /* Here the + 5 is a safety margin */ 03877 /**************************************************************************************/ 03878 03879 /*Must use CALLOC*/ 03880 yV = (double*)CALLOC( sizeof(double), (size_t)(-a + 5) ); 03881 03882 /**********************************************************************************************/ 03883 /* begin sanity check, check order, and that none negative */ 03884 /* THE GENERAL CASE */ 03885 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 03886 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 03887 /* */ 03888 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 03889 /* */ 03890 /* a (1-x) (a + bx - c) */ 03891 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 03892 /* (a-c) (a-c) */ 03893 /* */ 03894 /* */ 03895 /* A similiar recusion relation holds for b with a <--> b. */ 03896 /* */ 03897 /* */ 03898 /* we have initial conditions */ 03899 /* */ 03900 /* */ 03901 /* F(0) = 1 with a = -1 */ 03902 /* */ 03903 /* b */ 03904 /* F(-1) = 1 - (---) x with a = -1 */ 03905 /* c */ 03906 /* */ 03907 /* For the first F() in the solution of the radial integral */ 03908 /* */ 03909 /* a = ( -n + l + 1 ) */ 03910 /* */ 03911 /* a = -n + l + 1 */ 03912 /* max(a) = max(-n) + max(l) + 1 */ 03913 /* = max(-n) + max(n - 1) + 1 */ 03914 /* = max(-n + n - 1) + 1 */ 03915 /* = max(-1) + 1 */ 03916 /* = 0 */ 03917 /* */ 03918 /* similiarly */ 03919 /* */ 03920 /* min(a) = min(-n) + min(l) + 1 */ 03921 /* = min(-n) + 0 + 1 */ 03922 /* = (-n) + 0 + 1 */ 03923 /* = -n + 1 */ 03924 /* */ 03925 /* a -> a' + 1 implies */ 03926 /* */ 03927 /* max(a') = -1 */ 03928 /* min(a') = -n */ 03929 /* */ 03930 /* For the second F() in the solution of the radial integral */ 03931 /* */ 03932 /* a = ( -n + l - 1 ) */ 03933 /* */ 03934 /* a = -n + l + 1 */ 03935 /* max(a) = max(-n) + max(l) - 1 */ 03936 /* = -n + (n - 1) - 1 */ 03937 /* = -2 */ 03938 /* */ 03939 /* similiarly */ 03940 /* */ 03941 /* min(a) = min(-n) + min(l) - 1 */ 03942 /* = (-n) + 0 - 1 */ 03943 /* = -n - 1 */ 03944 /* */ 03945 /* a -> a' + 1 implies */ 03946 /* */ 03947 /* max(a') = -3 */ 03948 /* min(a') = -n - 2 */ 03949 /**********************************************************************************************/ 03950 03951 ASSERT( a <= 0 ); 03952 ASSERT( b <= 0 ); 03953 ASSERT( c >= 0 ); 03954 03955 d1 = F21i(a, b, c, y, yV ); 03956 free( yV ); 03957 return d1; 03958 } 03959 03960 STATIC mx F21_mx( long int a , long int b, long int c, double y, char A ) 03961 { 03962 DEBUG_ENTRY( "F21_mx()" ); 03963 03964 mx result_mx = {0.0,0}; 03965 mxq *yV = NULL; 03966 03967 /**************************************************************/ 03968 /* A must be either 'a' or 'b' */ 03969 /* and is use to determine which */ 03970 /* variable recursion will be over */ 03971 /**************************************************************/ 03972 03973 ASSERT( A == 'a' || A == 'b' ); 03974 03975 if( A == 'b' ) 03976 { 03977 /* if the recursion is over "b" */ 03978 /* then make it over "a" by switching these around */ 03979 long int i1 = a; 03980 a = b; 03981 b = i1; 03982 A = 'a'; 03983 } 03984 03985 /**************************************************************************************/ 03986 /* malloc space for the (dynamic) 1-d array */ 03987 /* 2_F_1 works via recursion and needs room to store intermediate results */ 03988 /* Here the + 5 is a safety margin */ 03989 /**************************************************************************************/ 03990 03991 /*Must use CALLOC*/ 03992 yV = (mxq *)CALLOC( sizeof(mxq), (size_t)(-a + 5) ); 03993 03994 /**********************************************************************************************/ 03995 /* begin sanity check, check order, and that none negative */ 03996 /* THE GENERAL CASE */ 03997 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */ 03998 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */ 03999 /* */ 04000 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */ 04001 /* */ 04002 /* a (1-x) (a + bx - c) */ 04003 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */ 04004 /* (a-c) (a-c) */ 04005 /* */ 04006 /* */ 04007 /* A similiar recusion relation holds for b with a <--> b. */ 04008 /* */ 04009 /* */ 04010 /* we have initial conditions */ 04011 /* */ 04012 /* */ 04013 /* F(0) = 1 with a = -1 */ 04014 /* */ 04015 /* b */ 04016 /* F(-1) = 1 - (---) x with a = -1 */ 04017 /* c */ 04018 /* */ 04019 /* For the first F() in the solution of the radial integral */ 04020 /* */ 04021 /* a = ( -n + l + 1 ) */ 04022 /* */ 04023 /* a = -n + l + 1 */ 04024 /* max(a) = max(-n) + max(l) + 1 */ 04025 /* = max(-n) + max(n - 1) + 1 */ 04026 /* = max(-n + n - 1) + 1 */ 04027 /* = max(-1) + 1 */ 04028 /* = 0 */ 04029 /* */ 04030 /* similiarly */ 04031 /* */ 04032 /* min(a) = min(-n) + min(l) + 1 */ 04033 /* = min(-n) + 0 + 1 */ 04034 /* = (-n) + 0 + 1 */ 04035 /* = -n + 1 */ 04036 /* */ 04037 /* a -> a' + 1 implies */ 04038 /* */ 04039 /* max(a') = -1 */ 04040 /* min(a') = -n */ 04041 /* */ 04042 /* For the second F() in the solution of the radial integral */ 04043 /* */ 04044 /* a = ( -n + l - 1 ) */ 04045 /* */ 04046 /* a = -n + l + 1 */ 04047 /* max(a) = max(-n) + max(l) - 1 */ 04048 /* = -n + (n - 1) - 1 */ 04049 /* = -2 */ 04050 /* */ 04051 /* similiarly */ 04052 /* */ 04053 /* min(a) = min(-n) + min(l) - 1 */ 04054 /* = (-n) + 0 - 1 */ 04055 /* = -n - 1 */ 04056 /* */ 04057 /* a -> a' + 1 implies */ 04058 /* */ 04059 /* max(a') = -3 */ 04060 /* min(a') = -n - 2 */ 04061 /**********************************************************************************************/ 04062 04063 ASSERT( a <= 0 ); 04064 ASSERT( b <= 0 ); 04065 ASSERT( c >= 0 ); 04066 04067 result_mx = F21i_log(a, b, c, y, yV ); 04068 free( yV ); 04069 return result_mx; 04070 } 04071 04072 STATIC double F21i(long int a, long int b, long int c, double y, double *yV ) 04073 { 04074 DEBUG_ENTRY( "F21i()" ); 04075 04076 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.; 04077 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.; 04078 long int i1 = 0, i2 = 0; 04079 04080 if( a == 0 ) 04081 { 04082 return 1.; 04083 } 04084 else if( a == -1 ) 04085 { 04086 ASSERT( c != 0 ); 04087 d1 = (double)b; 04088 d2 = (double)c; 04089 d3 = y * (d1/d2); 04090 d4 = 1. - d3; 04091 return d4; 04092 } 04093 /* Check to see if y(-a) != 0 in a very round about way to avoid lclint:error:13 */ 04094 else if( yV[-a] != 0. ) 04095 { 04096 /* Return the stored result */ 04097 return yV[-a]; 04098 } 04099 else 04100 { 04101 /******************************************************************************************/ 04102 /* - - */ 04103 /* (a)(1 - y) | | (a + bx + c) */ 04104 /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */ 04105 /* (a - c) | | (a - c) */ 04106 /* - - */ 04107 /* */ 04108 /* */ 04109 /* */ 04110 /* */ 04111 /* */ 04112 /* with F(0) = 1 */ 04113 /* b */ 04114 /* and F(-1) = 1 - (---) y */ 04115 /* c */ 04116 /* */ 04117 /* */ 04118 /* */ 04119 /* Use a -> a' + 1 */ 04120 /* _ _ */ 04121 /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */ 04122 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */ 04123 /* (a' + 1 - c) | | (a' + 1 - c) */ 04124 /* - - */ 04125 /* */ 04126 /* For the first F() in the solution of the radial integral */ 04127 /* */ 04128 /* a = ( -n + l + 1 ) */ 04129 /* */ 04130 /* a = -n + l + 1 */ 04131 /* max(a) = max(-n) + max(l) + 1 */ 04132 /* = -n + max(n-1) + 1 */ 04133 /* = -n + n-1 + 1 */ 04134 /* = 0 */ 04135 /* */ 04136 /* similiarly */ 04137 /* */ 04138 /* min(a) = min(-n) + min(l) + 1 */ 04139 /* = min(-n) + 0 + 1 */ 04140 /* = -n + 1 */ 04141 /* */ 04142 /* a -> a' + 1 implies */ 04143 /* */ 04144 /* max(a') = -1 */ 04145 /* min(a') = -n */ 04146 /******************************************************************************************/ 04147 04148 i1= a + 1; 04149 i2= a + 1 - c; 04150 d0= (double)i2; 04151 ASSERT( i2 != 0 ); 04152 d1= 1. - y; 04153 d2= (double)i1 * d1; 04154 d3= d2 / d0; 04155 d4= (double)b * y; 04156 d5= d0 + d4; 04157 04158 d8= F21i( (long int)(a + 1), b, c, y, yV ); 04159 d9= F21i( (long int)(a + 2), b, c, y, yV ); 04160 04161 d10= d8 - d9; 04162 d11= d3 * d10; 04163 d12= d5 / d0; 04164 d13= d12 * d8; 04165 d14= d11 + d13; 04166 04167 /* Store the result for later use */ 04168 yV[-a] = d14; 04169 return d14; 04170 } 04171 } 04172 04173 STATIC mx F21i_log( long int a, long int b, long int c, double y, mxq *yV ) 04174 { 04175 DEBUG_ENTRY( "F21i_log()" ); 04176 04177 mx result_mx = {0.0,0}; 04178 04179 if( yV[-a].q != 0. ) 04180 { 04181 /* Return the stored result */ 04182 return yV[-a].mx; 04183 } 04184 else if( a == 0 ) 04185 { 04186 ASSERT( yV[-a].q == 0. ); 04187 04188 result_mx.m = 1.; 04189 result_mx.x = 0; 04190 04191 ASSERT( yV[-a].mx.m == 0. ); 04192 ASSERT( yV[-a].mx.x == 0 ); 04193 04194 yV[-a].q = 1; 04195 yV[-a].mx = result_mx; 04196 return result_mx; 04197 } 04198 else if( a == -1 ) 04199 { 04200 double d1 = (double)b; 04201 double d2 = (double)c; 04202 double d3 = y * (d1/d2); 04203 04204 ASSERT( yV[-a].q == 0. ); 04205 ASSERT( c != 0 ); 04206 ASSERT( y != 0. ); 04207 04208 result_mx.m = 1. - d3; 04209 result_mx.x = 0; 04210 04211 while ( fabs(result_mx.m) > 1.0e+25 ) 04212 { 04213 result_mx.m /= 1.0e+25; 04214 result_mx.x += 25; 04215 } 04216 04217 ASSERT( yV[-a].mx.m == 0. ); 04218 ASSERT( yV[-a].mx.x == 0 ); 04219 04220 yV[-a].q = 1; 04221 yV[-a].mx = result_mx; 04222 return result_mx; 04223 } 04224 else 04225 { 04226 /******************************************************************************************/ 04227 /* - - */ 04228 /* (a)(1 - y) | | (a + bx + c) */ 04229 /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */ 04230 /* (a - c) | | (a - c) */ 04231 /* - - */ 04232 /* */ 04233 /* */ 04234 /* with F(0) = 1 */ 04235 /* b */ 04236 /* and F(-1) = 1 - (---) y */ 04237 /* c */ 04238 /* */ 04239 /* */ 04240 /* */ 04241 /* Use a -> a' + 1 */ 04242 /* _ _ */ 04243 /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */ 04244 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */ 04245 /* (a' + 1 - c) | | (a' + 1 - c) */ 04246 /* - - */ 04247 /* */ 04248 /* For the first F() in the solution of the radial integral */ 04249 /* */ 04250 /* a = ( -n + l + 1 ) */ 04251 /* */ 04252 /* a = -n + l + 1 */ 04253 /* max(a) = max(-n) + max(l) + 1 */ 04254 /* = -n + max(n-1) + 1 */ 04255 /* = -n + n-1 + 1 */ 04256 /* = 0 */ 04257 /* */ 04258 /* similiarly */ 04259 /* */ 04260 /* min(a) = min(-n) + min(l) + 1 */ 04261 /* = min(-n) + 0 + 1 */ 04262 /* = -n + 1 */ 04263 /* */ 04264 /* a -> a' + 1 implies */ 04265 /* */ 04266 /* max(a') = -1 */ 04267 /* min(a') = -n */ 04268 /******************************************************************************************/ 04269 04270 mx d8 = {0.0,0}, d9 = {0.0,0}, d10 = {0.0,0}, d11 = {0.0,0}; 04271 04272 double db = (double)b; 04273 double d00 = (double)(a + 1 - c); 04274 double d0 = (double)(a + 1); 04275 double d1 = 1. - y; 04276 double d2 = d0 * d1; 04277 double d3 = d2 / d00; 04278 double d4 = db * y; 04279 double d5 = d00 + d4; 04280 double d6 = d5 / d00; 04281 04282 ASSERT( yV[-a].q == 0. ); 04283 04284 /******************************************************************************************/ 04285 /* _ _ */ 04286 /* (a' + 1) (1 - x) | | (a' + 1 - c) + b*x */ 04287 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ------------------ F(a'+1) */ 04288 /* (a' + 1 - c) | | (a' + 1 - c) */ 04289 /* - - */ 04290 /******************************************************************************************/ 04291 04292 d8= F21i_log( (a + 1), b, c, y, yV ); 04293 d9= F21i_log( (a + 2), b, c, y, yV ); 04294 04295 if( (d8.m) != 0. ) 04296 { 04297 d10.x = d8.x; 04298 d10.m = d8.m * (1. - (d9.m/d8.m) * powi( 10., (d9.x - d8.x))); 04299 } 04300 else 04301 { 04302 d10.m = -d9.m; 04303 d10.x = d9.x; 04304 } 04305 04306 d10.m *= d3; 04307 04308 d11.x = d8.x; 04309 d11.m = d6 * d8.m; 04310 04311 if( (d11.m) != 0. ) 04312 { 04313 result_mx.x = d11.x; 04314 result_mx.m = d11.m * (1. + (d10.m/d11.m) * powi( 10. , (d10.x - d11.x) ) ); 04315 } 04316 else 04317 { 04318 result_mx = d10; 04319 } 04320 04321 while ( fabs(result_mx.m) > 1.0e+25 ) 04322 { 04323 result_mx.m /= 1.0e+25; 04324 result_mx.x += 25; 04325 } 04326 04327 /* Store the result for later use */ 04328 yV[-a].q = 1; 04329 yV[-a].mx = result_mx; 04330 return result_mx; 04331 } 04332 } 04333 04334 /********************************************************************************/ 04335 04336 inline void normalize_mx( mx& target ) 04337 { 04338 while( fabs(target.m) > 1.0e+25 ) 04339 { 04340 target.m /= 1.0e+25; 04341 target.x += 25; 04342 } 04343 while( fabs(target.m) < 1.0e-25 ) 04344 { 04345 target.m *= 1.0e+25; 04346 target.x -= 25; 04347 } 04348 return; 04349 } 04350 04351 inline mx add_mx( const mx& a, const mx& b ) 04352 { 04353 mx result = {0.0,0}; 04354 04355 if( a.m != 0. ) 04356 { 04357 result.x = a.x; 04358 result.m = a.m * (1. + (b.m/a.m) * powi( 10. , (b.x - a.x) ) ); 04359 } 04360 else 04361 { 04362 result = b; 04363 } 04364 normalize_mx( result ); 04365 return result; 04366 } 04367 04368 inline mx sub_mx( const mx& a, const mx& b ) 04369 { 04370 mx result = {0.0,0}; 04371 mx minusb = b; 04372 minusb.m = -minusb.m; 04373 04374 result = add_mx( a, minusb ); 04375 normalize_mx( result ); 04376 04377 return result; 04378 } 04379 04380 inline mx mxify( double a ) 04381 { 04382 mx result_mx = {0.0,0}; 04383 04384 result_mx.x = 0; 04385 result_mx.m = a; 04386 normalize_mx( result_mx ); 04387 04388 return result_mx; 04389 } 04390 04391 inline double unmxify( const mx& a_mx ) 04392 { 04393 return a_mx.m * powi( 10., a_mx.x ); 04394 } 04395 04396 inline mx mxify_log10( double log10_a ) 04397 { 04398 mx result_mx = {0.0,0}; 04399 04400 while ( log10_a > 25.0 ) 04401 { 04402 log10_a -= 25.0; 04403 result_mx.x += 25; 04404 } 04405 04406 while ( log10_a < -25.0 ) 04407 { 04408 log10_a += 25.0; 04409 result_mx.x -= 25; 04410 } 04411 04412 result_mx.m = pow(10., log10_a); 04413 04414 return result_mx; 04415 } 04416 04417 inline mx mult_mx( const mx& a, const mx& b ) 04418 { 04419 mx result = {0.0,0}; 04420 04421 result.m = (a.m * b.m); 04422 result.x = (a.x + b.x); 04423 normalize_mx( result ); 04424 04425 return result; 04426 } 04427 04428 inline double log10_prodxx( long int lp, double Ksqrd ) 04429 { 04430 /**********************************************************************/ 04431 /* | s=l' | s=l' */ 04432 /* | ----- | --- */ 04433 /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */ 04434 /* | | | | --- */ 04435 /* | s=0 | s=0 */ 04436 /**********************************************************************/ 04437 04438 if( lp == 0 ) 04439 return 0.; 04440 04441 double partsum = 0.; 04442 for( long int s = 1; s <= lp; s++ ) 04443 { 04444 double s2 = pow2((double)s); 04445 partsum += log10( 1. + s2*Ksqrd ); 04446 04447 ASSERT( partsum >= 0. ); 04448 } 04449 return partsum; 04450 }