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