cloudy trunk
Loading...
Searching...
No Matches
hydro_bauman.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/************************************************************************************************/
4/*H_photo_cs_lin returns hydrogenic photoionization cross section in cm-2 */
5/*H_Einstein_A_lin calculates Einstein A for any nlz */
6/*hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */
7/************************************************************************************************/
8/*************************** LOG version of h_bauman.c ****************************************/
9/* In this version, quantities that would normally cause a 64-bit floating point processor */
10/* to either underflow or overflow are evaluated using logs instead of floating point math. */
11/* This allows us to use an upper principal quantum number `n' greater than which the */
12/* other version begins to fail. The trade-off is, of course, lower accuracy */
13/* ( or is it precision ). We use LOG_10 for convenience. */
14/************************************************************************************************/
15#include "cddefines.h"
16#include "physconst.h"
17#include "thirdparty.h"
18#include "hydro_bauman.h"
19
20struct t_mx
21{
22 double m;
23 long int x;
24};
25
26typedef struct t_mx mx;
27
28struct t_mxq
29{
30 struct t_mx mx;
31 long int q;
32};
33
34typedef struct t_mxq mxq;
35
36/************************************************************************************************/
37/* these routines were written by Robert Bauman */
38/* The main reference for this section of code is */
39/* M. Brocklehurst */
40/* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
41/* */
42/* The recombination coefficient is obtained from the */
43/* photoionization cross-section (see Burgess 1965). */
44/* We have: */
45/* */
46/* - - l'=l+1 */
47/* | 2 pi^(1/2) alpha^4 a_o^2 c | 2 y^(1/2) --- */
48/* alpha(n,l,Z,Te)=|-----------------------------| --------- Z > I(n, l, l', t) */
49/* | 3 | n^2 --- */
50/* - - l'=l-1 */
51/* */
52/* where OO */
53/* - */
54/* | */
55/* 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) */
56/* | */
57/* - */
58/* 0 */
59/* */
60/* Here K = k / Z */
61/* */
62/* */
63/* and */
64/* */
65/* */
66/* y = Z^2 Rhc/(k Te)= 15.778/t */
67/* */
68/* where "t" is the scaled temperature, and "Te" is the electron Temperature */
69/* */
70/* t = Te/(10^4 Z^2) */
71/* Te in kelvin */
72/* */
73/* | |^2 */
74/* Theta(n,l; K, l') = (1 + n^2 K^2) | g(n,l; K,l') | */
75/* | | */
76/* */
77/* */
78/* ---- Not sure if this is K or k */
79/* OO / but think it is K */
80/* where - v */
81/* K^2 | */
82/* g(n,l; K,l') = ----- | R_nl(r) F_(K,l) r dr */
83/* n^2 | */
84/* - */
85/* 0 */
86/* */
87/* */
88/* - - */
89/* | 2 pi^(1/2) alpha^4 a_o^2 c | */
90/* |-----------------------------| */
91/* | 3 | */
92/* - - */
93/* */
94/* = 2 * (3.141592654)^1/2 * (7.29735308e-3)^4 */
95/* * (0.529177249e-10)^2 * (2.99792458e8) / 3 */
96/* */
97/* = 2.8129897e-21 */
98/* Mathematica gives 2.4764282710571237e-21 */
99/* */
100/* The photoionization cross-section (see also Burgess 1965) */
101/* is given by; */
102/* _ _ l'=l+1 */
103/* |4 PI alpha a_o^2 | n^2 --- max(l,l') */
104/* a(Z';n,l,k)=|---------------- | --- > --------- Theta(n,l; K, l') */
105/* | 3 | Z^2 --- (2l + 1) */
106/* _ _ l'=l-1 */
107/* */
108/* */
109/* where Theta(n,l; K, l') is defined above */
110/************************************************************************************************/
111/************************************************************************************************/
112/* For the transformation: */
113/* Z -> rZ = Z' */
114/* */
115/* k -> rk = k' */
116/* then */
117/* */
118/* K -> K = K' */
119/* */
120/* and the cross-sections satisfy; */
121/* 1 */
122/* a(Z'; n,l,k') = --- a(Z; n,l,k) */
123/* r^2 */
124/* */
125/* Similiarly, the recombination coefficient satisfies */
126/************************************************************************************************/
127
128/************************************************************************/
129/* IN THE FOLLOWING WE HAVE n > n' */
130/************************************************************************/
131/* returns hydrogenic photoionization cross section in cm-2 */
132/* this routine is called by H_photo_cs when n is small */
134 /* photon energy relative to threshold energy */
135 double rel_photon_energy,
136 /* principal quantum number, 1 for ground */
137 long int n,
138 /* angular momentum, 0 for s */
139 long int l,
140 /* charge, 1 for H+, 2 for He++, etc */
141 long int iz );
142
167double H_photo_cs_log10(
168 double photon_energy,
169 long int n,
170 long int l,
171 long int iz
172);
173
174/****************************************************************************/
175/* Calculates the Einstein A's for hydrogen */
176STATIC double H_Einstein_A_lin( /* IN THE FOLLOWING WE HAVE n > n' */
177 /* principal quantum number, 1 for ground, upper level */
178 long int n,
179 /* angular momentum, 0 for s */
180 long int l,
181 /* principal quantum number, 1 for ground, lower level */
182 long int np,
183 /* angular momentum, 0 for s */
184 long int lp,
185 /* Nuclear charge, 1 for H+, 2 for He++, etc */
186 long int iz
187);
188
202double H_Einstein_A_log10(
203 long int n,
204 long int l,
205 long int np,
206 long int lp,
207 long int iz
208);
209
230inline double OscStr_f(
231 long int n,
232 long int l,
233 long int np,
234 long int lp,
235 long int iz
236);
237
258inline double OscStr_f_log10(
259 long int n,
260 long int l,
261 long int np,
262 long int lp,
263 long int iz
264);
265
266/******************************************************************************
267 * F21()
268 * Calculates the Hyper_Spherical_Function 2_F_1(a,b,c;y)
269 * Here a,b, and c are (long int)
270 * y is of type (double)
271 * A is of type (char) and specifies whether the recursion is over
272 * a or b. It has values A='a' or A='b'.
273 ******************************************************************************/
274
275STATIC double F21(
276 long int a,
277 long int b,
278 long int c,
279 double y,
280 char A
281);
282
283STATIC double F21i(
284 long int a,
285 long int b,
286 long int c,
287 double y,
288 double *yV
289);
290
291/****************************************************************************************/
292/* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */
293/* In the following, we have n > n' */
294/****************************************************************************************/
295
296inline double hv(
297 /* returns energy in ergs */
298 /* principal quantum number, 1 for ground, upper level */
299 long int n,
300 /* principal quantum number, 1 for ground, lower level */
301 long int nprime,
302 long int iz
303);
304
305/********************************************************************************/
306/* In the following, we have n > n' */
307/********************************************************************************/
308
309STATIC double fsff(
310 /* principal quantum number, 1 for ground, upper level */
311 long int n,
312 /* angular momentum, 0 for s */
313 long int l,
314 /* principal quantum number, 1 for ground, lower level */
315 long int np
316);
317
318STATIC double log10_fsff(
319 /* principal quantum number, 1 for ground, upper level */
320 long int n,
321 /* angular momentum, 0 for s */
322 long int l,
323 /* principal quantum number, 1 for ground, lower level */
324 long int np
325);
326
327/********************************************************************************/
328/* F21_mx() */
329/* Calculates the Hyper_Spherical_Function 2_F_1(a,b,c;y) */
330/* Here a,b, and c are (long int) */
331/* y is of type (double) */
332/* A is of type (char) and specifies whether the recursion is over */
333/* a or b. It has values A='a' or A='b'. */
334/********************************************************************************/
335
337 long int a,
338 long int b,
339 long int c,
340 double y,
341 char A
342);
343
345 long int a,
346 long int b,
347 long int c,
348 double y,
349 mxq *yV
350);
351
365inline double hri(
366 long int n,
367 long int l,
368 long int np,
369 long int lp,
370 long int iz
371);
372
385inline double hri_log10(
386 long int n,
387 long int l,
388 long int np,
389 long int lp,
390 long int iz
391);
392
393/******************************************************************************/
394/* In the following, we have n > n' */
395/******************************************************************************/
396STATIC double hrii(
397 /* principal quantum number, 1 for ground, upper level */
398 long int n,
399 /* angular momentum, 0 for s */
400 long int l,
401 /* principal quantum number, 1 for ground, lower level */
402 long int np,
403 /* angular momentum, 0 for s */
404 long int lp
405);
406
407STATIC double hrii_log(
408 /* principal quantum number, 1 for ground, upper level */
409 long int n,
410 /* angular momentum, 0 for s */
411 long int l,
412 /* principal quantum number, 1 for ground, lower level */
413 long int np,
414 /* angular momentum, 0 for s */
415 long int lp
416);
417
418STATIC double bh(
419 double k,
420 long int n,
421 long int l,
422 double *rcsvV
423);
424
425STATIC double bh_log(
426 double k,
427 long int n,
428 long int l,
429 mxq *rcsvV_mxq
430);
431
432STATIC double bhintegrand(
433 double k,
434 long int n,
435 long int l,
436 long int lp,
437 double *rcsvV
438);
439
441 double k,
442 long int n,
443 long int l,
444 long int lp,
445 mxq *rcsvV_mxq
446);
447
448STATIC double bhG(
449 double K,
450 long int n,
451 long int l,
452 long int lp,
453 double *rcsvV
454);
455
457 double K,
458 long int n,
459 long int l,
460 long int lp,
461 mxq *rcsvV_mxq
462);
463
464STATIC double bhGp(
465 long int q,
466 double K,
467 long int n,
468 long int l,
469 long int lp,
470 double *rcsvV,
471 double GK
472);
473
475 long int q,
476 double K,
477 long int n,
478 long int l,
479 long int lp,
480 mxq *rcsvV_mxq,
481 const mx& GK_mx
482);
483
484STATIC double bhGm(
485 long int q,
486 double K,
487 long int n,
488 long int l,
489 long int lp,
490 double *rcsvV,
491 double GK
492);
493
495 long int q,
496 double K,
497 long int n,
498 long int l,
499 long int lp,
500 mxq *rcsvV_mxq,
501 const mx& GK_mx
502);
503
504STATIC double bhg(
505 double K,
506 long int n,
507 long int l,
508 long int lp,
509 double *rcsvV
510);
511
512STATIC double bhg_log(
513 double K,
514 long int n,
515 long int l,
516 long int lp,
517 mxq *rcsvV_mxq
518);
519
520inline void normalize_mx( mx& target );
521inline mx add_mx( const mx& a, const mx& b );
522inline mx sub_mx( const mx& a, const mx& b );
523inline mx mxify( double a );
524inline double unmxify( const mx& a_mx );
525inline mx mxify_log10( double log10_a );
526inline mx mult_mx( const mx& a, const mx& b );
527
528inline double local_product( double K , long int lp );
529inline double log10_prodxx( long int lp, double Ksqrd );
530
531/****************************************************************************************/
532/* 64 pi^4 (e a_o)^2 64 pi^4 (a_o)^2 e^2 1 1 */
533/* ----------------- = ----------------- -------- ---- = 7.5197711e-38 ----- */
534/* 3 h c^3 3 c^2 hbar c 2 pi sec */
535/****************************************************************************************/
536
538
539/************************************************************************/
540/* (4.0/3.0) * PI * FINE_STRUCTURE_CONSTANT * BOHRRADIUS * BOHRRADIUS */
541/* */
542/* */
543/* 4 PI alpha a_o^2 */
544/* ---------------- */
545/* 3 */
546/* */
547/* where alpha = Fine Structure Constant */
548/* a_o = Bohr Radius */
549/* */
550/* = 3.056708^-02 (au Length)^2 */
551/* = 8.56x10^-23 (meters)^2 */
552/* = 8.56x10^-19 (cm)^2 */
553/* = 8.56x10^+05 (barns) */
554/* = 0.856 (MB or megabarns) */
555/* */
556/* */
557/* 1 barn = 10^-28 (meter)^2 */
558/************************************************************************/
559
561
562/************************Start of program***************************/
563
565 /* incident photon energy relative to threshold */
566 double rel_photon_energy,
567 /* principal quantum number, 1 for ground */
568 long int n,
569 /* angular momentum, 0 for s */
570 long int l,
571 /* charge, 1 for H+, 2 for He++, etc */
572 long int iz )
573{
574 DEBUG_ENTRY( "H_photo_cs()" );
575
576 double result;
577 if( n<= 25 )
578 {
579 result = H_photo_cs_lin( rel_photon_energy , n , l , iz );
580 }
581 else
582 {
583 result = H_photo_cs_log10( rel_photon_energy , n , l , iz );
584 }
585 return result;
586}
587
588/************************************************************************/
589/* IN THE FOLLOWING WE HAVE n > n' */
590/************************************************************************/
591
592/* returns hydrogenic photoionization cross section in cm-2 */
593/* this routine is called by H_photo_cs when n is small */
595 /* photon energy relative to threshold energy */
596 double rel_photon_energy,
597 /* principal quantum number, 1 for ground */
598 long int n,
599 /* angular momentum, 0 for s */
600 long int l,
601 /* charge, 1 for H+, 2 for He++, etc */
602 long int iz )
603{
604 DEBUG_ENTRY( "H_photo_cs_lin()" );
605
606 long int dim_rcsvV;
607
608 /* >>chng 02 sep 15, make rcsvV always NPRE_FACGTORIAL+3 long */
609 double rcsvV[NPRE_FACTORIAL+3];
610 int i;
611
612 double electron_energy;
613 double result = 0.;
614 double xn_sqrd = (double)(n*n);
615 double z_sqrd = (double)(iz*iz);
616 double Z = (double)iz;
617 double K = 0.; /* K = k / Z */
618 double k = 0.; /* k^2 = ejected-electron-energy (Ryd) */
619
620 /* expressions blow up at precisely threshold */
621 if( rel_photon_energy < 1.+FLT_EPSILON )
622 {
623 /* below or very close to threshold, return zero */
624 return 0.;
625 }
626
627 if( n < 1 || l >= n )
628 {
629 fprintf(ioQQQ," The quantum numbers are impossible.\n");
631 }
632
633 if( ((2*n) - 1) >= NPRE_FACTORIAL )
634 {
635 fprintf(ioQQQ," This value of n is too large.\n");
637 }
638
639 /* k^2 is the ejected photoelectron energy in ryd */
640 /*electron_energy = SDIV( (photon_energy/ryd) - (z_sqrd/xn_sqrd) );*/
641
642 electron_energy = (rel_photon_energy-1.) * (z_sqrd/xn_sqrd);
643 k = sqrt( ( electron_energy ) );
644
645 K = (k/Z);
646
647 dim_rcsvV = (((n * 2) - 1) + 1);
648
649 for( i=0; i<dim_rcsvV; ++i )
650 {
651 rcsvV[i] = 0.;
652 }
653
654 /* rcsvV contains all results for quantum indices below n, l */
655 result = PHYSICAL_CONSTANT_TWO * (xn_sqrd/z_sqrd) * bh( K, n, l, rcsvV );
656
657 ASSERT( result != 0. );
658 return result;
659}
660
661/*****************************************************************************/
662/*H_photo_cs_log10 returns hydrogenic photoionization cross section in cm-2
663 * this routine is called by H_photo_cs when n is large */
664/*****************************************************************************/
666 /* photon energy relative to threshold energy */
667 double rel_photon_energy,
668 /* principal quantum number, 1 for ground */
669 long int n,
670 /* angular momentum, 0 for s */
671 long int l,
672 /* charge, 1 for H+, 2 for He++, etc */
673 long int iz
674)
675{
676 DEBUG_ENTRY( "H_photo_cs_log10()" );
677
678 long int dim_rcsvV_mxq;
679
680 mxq *rcsvV_mxq = NULL;
681
682 double electron_energy;
683 double t1;
684 double result = 0.;
685 double xn_sqrd = (double)(n*n);
686 double z_sqrd = (double)(iz*iz);
687 double Z = (double)iz;
688 double K = 0.; /* K = k / Z */
689 double k = 0.; /* k^2 = ejected-electron-energy (Ryd) */
690
691 /* expressions blow up at precisely threshold */
692 if( rel_photon_energy < 1.+FLT_EPSILON )
693 {
694 /* below or very close to threshold, return zero */
695 fprintf( ioQQQ,"PROBLEM IN HYDRO_BAUMAN: rel_photon_energy, n, l, iz: %e\t%li\t%li\t%li\n",
696 rel_photon_energy,
697 n,
698 l,
699 iz );
701 }
702 if( n < 1 || l >= n )
703 {
704 fprintf(ioQQQ," The quantum numbers are impossible.\n");
706 }
707
708 /* k^2 is the ejected photoelectron energy in ryd */
709 /*electron_energy = SDIV( (photon_energy/ryd) - (z_sqrd/xn_sqrd) );*/
710 electron_energy = (rel_photon_energy-1.) * (z_sqrd/xn_sqrd);
711
712 k = sqrt( ( electron_energy ) );
713 /* k^2 is the ejected photoelectron energy in ryd */
714 /*k = sqrt( ( (photon_energy/ryd) - (z_sqrd/xn_sqrd) ) );*/
715
716 K = (k/Z);
717
718 dim_rcsvV_mxq = (((n * 2) - 1) + 1);
719
720 /* create space */
721 rcsvV_mxq = (mxq*)CALLOC( (size_t)dim_rcsvV_mxq, sizeof(mxq) );
722
723 t1 = bh_log( K, n, l, rcsvV_mxq );
724
725 ASSERT( t1 > 0. );
726
727 t1 = MAX2( t1, 1.0e-250 );
728
729 result = PHYSICAL_CONSTANT_TWO * (xn_sqrd/z_sqrd) * t1;
730
731 free( rcsvV_mxq );
732 if( result <= 0. )
733 {
734 fprintf( ioQQQ, "PROBLEM: Hydro_Bauman...t1\t%e\n", t1 );
735 }
736 ASSERT( result > 0. );
737 return result;
738}
739
740STATIC double bh(
741 /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */
742 double K,
743 /* principal quantum number */
744 long int n,
745 /* angular momentum quantum number */
746 long int l,
747 /* Temporary storage for intermediate */
748 /* results of the recursive routine */
749 double *rcsvV
750)
751{
752 DEBUG_ENTRY( "bh()" );
753
754 long int lp = 0; /* l' */
755 double sigma = 0.; /* Sum in Brocklehurst eq. 3.13 */
756
757 ASSERT( n > 0 );
758 ASSERT( l >= 0 );
759 ASSERT( n > l );
760
761 if( l > 0 ) /* no lp=(l-1) for l=0 */
762 {
763 for( lp = l - 1; lp <= l + 1; lp = lp + 2 )
764 {
765 sigma += bhintegrand( K, n, l, lp, rcsvV );
766 }
767 }
768 else
769 {
770 lp = l + 1;
771 sigma = bhintegrand( K, n, l, lp, rcsvV );
772 }
773 ASSERT( sigma != 0. );
774 return sigma;
775}
776
778 /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */
779 double K,
780 /* principal quantum number */
781 long int n,
782 /* angular momentum quantum number */
783 long int l,
784 /* Temporary storage for intermediate */
785 mxq *rcsvV_mxq
786 /* results of the recursive routine */
787)
788{
789 DEBUG_ENTRY( "bh_log()" );
790
791 long int lp = 0; /* l' */
792 double sigma = 0.; /* Sum in Brocklehurst eq. 3.13 */
793
794 ASSERT( n > 0 );
795 ASSERT( l >= 0 );
796 ASSERT( n > l );
797
798 if( l > 0 ) /* no lp=(l-1) for l=0 */
799 {
800 for( lp = l - 1; lp <= l + 1; lp = lp + 2 )
801 {
802 sigma += bhintegrand_log( K, n, l, lp, rcsvV_mxq );
803 }
804 }
805 else
806 {
807 lp = l + 1;
808 sigma = bhintegrand_log( K, n, l, lp, rcsvV_mxq );
809 }
810 ASSERT( sigma != 0. );
811 return sigma;
812}
813
814/********************************************************************************/
815/* Here we calculate the integrand */
816/* (as a function of K, so */
817/* we need a dK^2 -> 2K dK ) */
818/* for equation 3.14 of reference */
819/* */
820/* M. Brocklehurst Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
821/* */
822/* namely: */
823/* */
824/* max(l,l') (1 + n^2 K^2)^2 Theta(n,l; K, l') exp( -K^2 y ) d(K^2) */
825/* */
826/* Note: the "y" is included in the code that called */
827/* this function and we include here the n^2 from eq 3.13. */
828/********************************************************************************/
829
831 /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */
832 double K,
833 long int n,
834 long int l,
835 long int lp,
836 /* Temporary storage for intermediate */
837 /* results of the recursive routine */
838 double *rcsvV
839)
840{
841 DEBUG_ENTRY( "bhintegrand()" );
842
843 double Two_L_Plus_One = (double)(2*l + 1);
844 double lg = (double)max(l,lp);
845
846 double n2 = (double)(n*n);
847
848 double Ksqrd = (K*K);
849
850 /**********************************************/
851 /* */
852 /* l> */
853 /* ------ Theta(nl,Kl') */
854 /* 2l+2 */
855 /* */
856 /* */
857 /* Theta(nl,Kl') = */
858 /* (1+n^2K^2) * | g(nl,Kl')|^2 */
859 /* */
860 /**********************************************/
861
862 double d2 = 1. + n2*Ksqrd;
863 double d5 = bhg( K, n, l, lp, rcsvV );
864 double Theta = d2 * d5 * d5;
865 double d7 = (lg/Two_L_Plus_One) * Theta;
866
867 ASSERT( Two_L_Plus_One != 0. );
868 ASSERT( Theta != 0. );
869 ASSERT( Ksqrd != 0. );
870 ASSERT( d2 != 0. );
871 ASSERT( d5 != 0. );
872 ASSERT( d7 != 0. );
873 ASSERT( lp >= 0 );
874 ASSERT( lg != 0. );
875 ASSERT( n2 != 0. );
876 ASSERT( n > 0 );
877 ASSERT( l >= 0 );
878 ASSERT( K != 0. );
879 return d7;
880}
881
882/************************************************************************************************/
883/* The photoionization cross-section (see also Burgess 1965) */
884/* is given by; */
885/* _ _ l'=l+1 */
886/* |4 PI alpha a_o^2 | n^2 --- max(l,l') */
887/* a(Z';n,l,k)=|---------------- | --- > ---------- Theta(n,l; K, l') */
888/* | 3 | Z^2 --- (2l + 1) */
889/* _ _ l'=l-1 */
890/* */
891/* */
892/* where Theta(n,l; K, l') is defined */
893/* */
894/* | |^2 */
895/* Theta(n,l; K, l') = (1 + n^2 K^2) | g(n,l; K,l') | */
896/* | | */
897/* */
898/* */
899/* ---- Not sure if this is K or k */
900/* OO / but think it is K */
901/* where - v */
902/* K^2 | */
903/* g(n,l; K,l') = ----- | R_nl(r) F_(K,l) r dr */
904/* n^2 | */
905/* - */
906/* 0 */
907/************************************************************************************************/
908
910 double K, /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */
911 long int n,
912 long int l,
913 long int lp,
914 /* Temporary storage for intermediate */
915 /* results of the recursive routine */
916 mxq *rcsvV_mxq
917)
918{
919 DEBUG_ENTRY( "bhintegrand_log()" );
920
921 double d2 = 0.;
922 double d5 = 0.;
923 double d7 = 0.;
924 double Theta = 0.;
925 double n2 = (double)(n*n);
926 double Ksqrd = (K*K);
927 double Two_L_Plus_One = (double)(2*l + 1);
928 double lg = (double)max(l,lp);
929
930 ASSERT( Ksqrd != 0. );
931 ASSERT( K != 0. );
932 ASSERT( lg != 0. );
933 ASSERT( n2 != 0. );
934 ASSERT( Two_L_Plus_One != 0. );
935
936 ASSERT( n > 0);
937 ASSERT( l >= 0);
938 ASSERT( lp >= 0);
939
940 /**********************************************/
941 /* */
942 /* l> */
943 /* ------ Theta(nl,Kl') */
944 /* 2l+2 */
945 /* */
946 /* */
947 /* Theta(nl,Kl') = */
948 /* (1+n^2K^2) * | g(nl,Kl')|^2 */
949 /* */
950 /**********************************************/
951
952 d2 = ( 1. + n2 * (Ksqrd) );
953
954 ASSERT( d2 != 0. );
955
956 d5 = bhg_log( K, n, l, lp, rcsvV_mxq );
957
958 d5 = MAX2( d5, 1.0E-150 );
959 ASSERT( d5 != 0. );
960
961 Theta = d2 * d5 * d5;
962 ASSERT( Theta != 0. );
963
964 d7 = (lg/Two_L_Plus_One) * Theta;
965
966 ASSERT( d7 != 0. );
967 return d7;
968}
969
970/****************************************************************************************/
971/* *** bhG *** */
972/* Using various recursion relations */
973/* (for l'=l+1) */
974/* equation: (3.23) */
975/* G(n,l-2; K,l-1) = [ 4n^2-4l^2+l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */
976/* - 4n^2 (n^2-l^2)[1+(l+1)^2 K^2] G(n,l; K,l+1) */
977/* */
978/* and (for l'=l-1) */
979/* equation: (3.24) */
980/* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l; K,l-1) */
981/* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l; K,l+1) */
982/* */
983/* the starting point for the recursion relations are; */
984/* equation: (3.18) */
985/* | pi |(1/2) 8n */
986/* G(n,n-1; 0,n) = | -- | ------- (4n)^n exp(-2n) */
987/* | 2 | (2n-1)! */
988/* */
989/* equation: (3.20) */
990/* exp(2n-2/K tan^(-1)(n K) */
991/* G(n,n-1; K,n) = ----------------------------------------- * G(n,n-1; 0,n) */
992/* sqrt(1 - exp(-2 pi K)) * (1+(n K)^2)^(n+2) */
993/* */
994/* equation: (3.20) */
995/* G(n,n-2; K,n-1) = (2n-2)(1+(n K)^2) n G(n,n-1; K,n) */
996/* */
997/* equation: (3.21) */
998/* (1+(n K)^2) */
999/* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
1000/* 2n */
1001/* */
1002/* equation: (3.22) */
1003/* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+n^2 K^2)) G(n,n-1; K,n-2) */
1004/****************************************************************************************/
1005
1007 double K,
1008 long int n,
1009 long int l,
1010 long int lp,
1011 /* Temporary storage for intermediate */
1012 /* results of the recursive routine */
1013 double *rcsvV
1014)
1015{
1016 DEBUG_ENTRY( "bhG()" );
1017
1018 double n1 = (double)n;
1019 double n2 = (double)(n * n);
1020 double Ksqrd = K * K;
1021
1022 double ld1 = factorial( 2*n - 1 );
1023 double ld2 = powi((double)(4*n), n);
1024 double ld3 = exp(-(double)(2 * n));
1025
1026 /******************************************************************************
1027 * ********G0******* *
1028 * *
1029 * | pi |(1/2) 8n *
1030 * G0 = | -- | ------- (4n)^n exp(-2n) *
1031 * | 2 | (2n-1)! *
1032 ******************************************************************************/
1033
1034 double G0 = SQRTPIBY2 * (8. * n1 * ld2 * ld3) / ld1;
1035
1036 double d1 = sqrt( 1. - exp(( -2. * PI )/ K ));
1037 double d2 = powi(( 1. + (n2 * Ksqrd)), ( n + 2 ));
1038 double d3 = atan( n1 * K );
1039 double d4 = ((2. / K) * d3);
1040 double d5 = (double)( 2 * n );
1041 double d6 = exp( d5 - d4 );
1042 double GK = ( d6 /( d1 * d2 ) ) * G0;
1043
1044 /* l=l'-1 or l=l'+1 */
1045 ASSERT( (l == lp - 1) || (l == lp + 1) );
1046 ASSERT( K != 0. );
1047 ASSERT( Ksqrd != 0. );
1048 ASSERT( n1 != 0. );
1049 ASSERT( n2 != 0. );
1050 ASSERT( ((2*n) - 1) < 1755 );
1051 ASSERT( ((2*n) - 1) >= 0 );
1052 ASSERT( ld1 != 0. );
1053 ASSERT( (1.0 / ld1) != 0. );
1054 ASSERT( ld3 != 0. );
1055
1056 ASSERT( K != 0. );
1057 ASSERT( d1 != 0. );
1058 ASSERT( d2 != 0. );
1059 ASSERT( d3 != 0. );
1060 ASSERT( d4 != 0. );
1061 ASSERT( d5 != 0. );
1062 ASSERT( d6 != 0. );
1063
1064 ASSERT( G0 != 0. );
1065 ASSERT( GK != 0. );
1066
1067 /******************************************************************************/
1068 /* *****GK***** */
1069 /* */
1070 /* exp(2n-2/K tan^(-1)(n K) */
1071 /* G(n,n-1; K,n) = ----------------------------------------- * G0 */
1072 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K))^(n+2) */
1073 /******************************************************************************/
1074
1075 /* GENERAL CASE: l = l'-1 */
1076 if( l == lp - 1 )
1077 {
1078 double result = bhGm( l, K, n, l, lp, rcsvV, GK );
1079 /* Here the m in bhGm() refers */
1080 /* to the minus sign(-) in l=l'-1 */
1081 return result;
1082 }
1083
1084 /* GENERAL CASE: l = l'+1 */
1085 else if( l == lp + 1 )
1086 {
1087 double result = bhGp( l, K, n, l, lp, rcsvV, GK );
1088 /* Here the p in bhGp() refers */
1089 /* to the plus sign(+) in l=l'+1 */
1090 return result;
1091 }
1092 else
1093 {
1094 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
1096 }
1097}
1098
1099/*************log version********************************/
1101 double K,
1102 long int n,
1103 long int l,
1104 long int lp,
1105 /* Temporary storage for intermediate */
1106 /* results of the recursive routine */
1107 mxq *rcsvV_mxq
1108)
1109{
1110 DEBUG_ENTRY( "bhG_mx()" );
1111
1112 double log10_GK = 0.;
1113 double log10_G0 = 0.;
1114
1115 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
1116 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0.;
1117
1118 double n1 = (double)n;
1119 double n2 = n1 * n1;
1120 double Ksqrd = K * K;
1121
1122 mx GK_mx = {0.0,0};
1123
1124 /* l=l'-1 or l=l'+1 */
1125 ASSERT( (l == lp - 1) || (l == lp + 1) );
1126 ASSERT( K != 0. );
1127 ASSERT( n1 != 0. );
1128 ASSERT( n2 != 0. );
1129 ASSERT( Ksqrd != 0. );
1130 ASSERT( ((2*n) - 1) >= 0 );
1131
1132 /******************************/
1133 /* n */
1134 /* --- */
1135 /* log( n! ) = > log(j) */
1136 /* --- */
1137 /* j=1 */
1138 /******************************/
1139
1140 /*************************************************************/
1141 /* | pi |(1/2) 8n */
1142 /* G(n,n-1; 0,n) = | -- | ------- (4n)^n exp(-2n) */
1143 /* | 2 | (2n-1)! */
1144 /*************************************************************/
1145
1146 /******************************/
1147 /* */
1148 /* */
1149 /* log10( (2n-1)! ) */
1150 /* */
1151 /* */
1152 /******************************/
1153
1154 ld1 = lfactorial( 2*n - 1 );
1155 ASSERT( ld1 >= 0. );
1156
1157 /**********************************************/
1158 /* (4n)^n */
1159 /**********************************************/
1160 /* log10( 4n^n ) = n log10( 4n ) */
1161 /**********************************************/
1162
1163 ld2 = n1 * log10( 4. * n1 );
1164 ASSERT( ld2 >= 0. );
1165
1166 /**********************************************/
1167 /* exp(-2n) */
1168 /**********************************************/
1169 /* log10( exp( -2n ) ) = (-2n) * log10(e) */
1170 /**********************************************/
1171 ld3 = (-(2. * n1)) * (LOG10_E);
1172 ASSERT( ld3 <= 0. );
1173
1174 /******************************************************************************/
1175 /* ********G0******* */
1176 /* */
1177 /* | pi |(1/2) 8n */
1178 /* G0 = | -- | ------- (4n)^n exp(-2n) */
1179 /* | 2 | (2n-1)! */
1180 /******************************************************************************/
1181
1182 log10_G0 = log10(SQRTPIBY2 * 8. * n1) + ( (ld2 + ld3) - ld1);
1183
1184 /******************************************************************************/
1185 /* *****GK***** */
1186 /* */
1187 /* exp(2n- (2/K) tan^(-1)(n K) ) */
1188 /* G(n,n-1; K,n) = ----------------------------------------- * G0 */
1189 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K))^(n+2) */
1190 /******************************************************************************/
1191
1192 ASSERT( K != 0. );
1193
1194 /**********************************************/
1195 /* sqrt(1 - exp(-2 pi/ K)) */
1196 /**********************************************/
1197 /* log10(sqrt(1 - exp(-2 pi/ K))) = */
1198 /* (1/2) log10(1 - exp(-2 pi/ K)) */
1199 /**********************************************/
1200
1201 d1 = (1. - exp(-(2. * PI )/ K ));
1202 ld4 = (1./2.) * log10( d1 );
1203 ASSERT( K != 0. );
1204 ASSERT( d1 != 0. );
1205
1206 /**************************************/
1207 /* (1+(n K)^2)^(n+2) */
1208 /**************************************/
1209 /* log10( (1+(n K)^2)^(n+2) ) = */
1210 /* (n+2) log10( (1 + (n K)^2 ) ) */
1211 /**************************************/
1212
1213 d2 = ( 1. + (n2 * Ksqrd));
1214 ld5 = (n1 + 2.) * log10( d2 );
1215 ASSERT( d2 != 0. );
1216
1217 ASSERT( ld5 >= 0. );
1218
1219 /**********************************************/
1220 /* exp(2n- (2/K)*tan^(-1)(n K) ) */
1221 /**********************************************/
1222 /* log10( exp(2n- (2/K) tan^(-1)(n K) ) = */
1223 /* (2n- (2/K)*tan^(-1)(n K) ) * Log10(e) */
1224 /**********************************************/
1225
1226 /* tan^(-1)(n K) ) */
1227 d3 = atan( n1 * K );
1228 ASSERT( d3 != 0. );
1229
1230 /* (2/K)*tan^(-1)(n K) ) */
1231 d4 = (2. / K) * d3;
1232 ASSERT( d4 != 0. );
1233
1234 /* 2n */
1235 d5 = (double) ( 2. * n1 );
1236 ASSERT( d5 != 0. );
1237
1238 /* (2n-2/K tan^(-1)(n K)) */
1239 d6 = d5 - d4;
1240 ASSERT( d6 != 0. );
1241
1242 /* log10( exp(2n- (2/K) tan^(-1)(n K) ) */
1243 ld6 = LOG10_E * d6;
1244 ASSERT( ld6 != 0. );
1245
1246 /******************************************************************************/
1247 /* *****GK***** */
1248 /* */
1249 /* exp(2n- (2/K) tan^(-1)(n K) ) */
1250 /* G(n,n-1; K,n) = ----------------------------------------- * G0 */
1251 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K))^(n+2) */
1252 /******************************************************************************/
1253
1254 log10_GK = (ld6 -(ld4 + ld5)) + log10_G0;
1255 ASSERT( log10_GK != 0. );
1256
1257 GK_mx = mxify_log10( log10_GK );
1258
1259 /* GENERAL CASE: l = l'-1 */
1260 if( l == lp - 1 )
1261 {
1262 mx result_mx = bhGm_mx( l, K, n, l, lp, rcsvV_mxq , GK_mx );
1263 /* Here the m in bhGm() refers */
1264 /* to the minus sign(-) in l=l'-1 */
1265 return result_mx;
1266 }
1267 /* GENERAL CASE: l = l'+1 */
1268 else if( l == lp + 1 )
1269 {
1270 mx result_mx = bhGp_mx( l, K, n, l, lp, rcsvV_mxq , GK_mx );
1271 /* Here the p in bhGp() refers */
1272 /* to the plus sign(+) in l=l'+1 */
1273 return result_mx;
1274 }
1275 else
1276 {
1277 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
1279 }
1280 printf( "This code should be inaccessible\n\n" );
1282}
1283
1284/************************************************************************************************/
1285/* *** bhGp.c *** */
1286/* */
1287/* Here we calculate G(n,l; K,l') with the recursive formula */
1288/* equation: (3.24) */
1289/* */
1290/* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */
1291/* */
1292/* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l+1; K,l) */
1293/* */
1294/* Under the transformation l -> l + 1 this gives */
1295/* */
1296/* 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) */
1297/* */
1298/* - 4n^2 (n^2-((l+1)+1)^2)[ 1+((l+1)K)^2 ] G(n,l+1+1; K,l+1) */
1299/* */
1300/* or */
1301/* */
1302/* 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) */
1303/* */
1304/* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1305/* */
1306/* from the reference */
1307/* M. Brocklehurst */
1308/* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
1309/* */
1310/* */
1311/* * This is valid for the case l=l'+1 * */
1312/* * CASE: l = l'+1 * */
1313/* * Here the p in bhGp() refers * */
1314/* * to the Plus sign(+) in l=l'+1 * */
1315/************************************************************************************************/
1316
1318 long int q,
1319 double K,
1320 long int n,
1321 long int l,
1322 long int lp,
1323 /* Temporary storage for intermediate */
1324 /* results of the recursive routine */
1325 double *rcsvV,
1326 double GK
1327)
1328{
1329 DEBUG_ENTRY( "bhGp()" );
1330
1331 /* static long int rcsv_Level = 1;
1332 printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */
1333
1334 ASSERT( l == lp + 1 );
1335
1336 long int rindx = 2*q;
1337
1338 if( rcsvV[rindx] == 0. )
1339 {
1340 /* SPECIAL CASE: n = l+1 = l'+2 */
1341 if( q == n - 1 )
1342 {
1343 double Ksqrd = K * K;
1344 double n2 = (double)(n*n);
1345
1346 double dd1 = (double)(2 * n);
1347 double dd2 = ( 1. + ( n2 * Ksqrd));
1348
1349 /* (1+(n K)^2) */
1350 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
1351 /* 2n */
1352 double G1 = ((dd2 * GK) / dd1);
1353
1354 ASSERT( l == lp + 1 );
1355 ASSERT( Ksqrd != 0. );
1356 ASSERT( dd1 != 0. );
1357 ASSERT( dd2 != 0. );
1358 ASSERT( G1 != 0. );
1359
1360 rcsvV[rindx] = G1;
1361 return G1;
1362 }
1363 /* SPECIAL CASE: n = l+2 = l'+3 */
1364 else if( q == (n - 2) )
1365 {
1366 double Ksqrd = (K*K);
1367 double n2 = (double)(n*n);
1368
1369 /* */
1370 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */
1371 /* */
1372 double dd1 = (double) (2 * n);
1373 double dd2 = ( 1. + ( n2 * Ksqrd));
1374 double G1 = ((dd2 * GK) / dd1);
1375
1376 /* */
1377 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */
1378 /* */
1379 double dd3 = (double)((2 * n) - 1);
1380 double dd4 = (double)(n - 1);
1381 double dd5 = (4. + (dd4 * dd2));
1382 double G2 = (dd3 * dd5 * G1);
1383
1384 ASSERT( l == lp + 1 );
1385 ASSERT( Ksqrd != 0. );
1386 ASSERT( n2 != 0. );
1387 ASSERT( dd1 != 0. );
1388 ASSERT( dd2 != 0. );
1389 ASSERT( dd3 != 0. );
1390 ASSERT( dd4 != 0. );
1391 ASSERT( dd5 != 0. );
1392 ASSERT( G1 != 0. );
1393 ASSERT( G2 != 0. );
1394
1395 rcsvV[rindx] = G2;
1396 return G2;
1397 }
1398 /* The GENERAL CASE n > l + 2 */
1399 else
1400 {
1401 /******************************************************************************/
1402 /* 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) */
1403 /* */
1404 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1405 /* */
1406 /* FROM Eq. 3.24 */
1407 /* */
1408 /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */
1409 /* */
1410 /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */
1411 /******************************************************************************/
1412
1413 long int lp1 = (q + 1);
1414 long int lp2 = (q + 2);
1415
1416 double Ksqrd = (K*K);
1417 double n2 = (double)(n * n);
1418 double lp1s = (double)(lp1 * lp1);
1419 double lp2s = (double)(lp2 * lp2);
1420
1421 double d1 = (4. * n2);
1422 double d2 = (4. * lp1s);
1423 double d3 = (double)((lp1)*((2 * q) + 3));
1424 double d4 = (1. + (n2 * Ksqrd));
1425 double d5 = (d1 - d2 + (d3 * d4));
1426 double d5_1 = d5 * bhGp( (q+1), K, n, l, lp, rcsvV, GK );
1427
1428
1429 /* 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) */
1430 /* */
1431 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1432
1433 double d6 = (n2 - lp2s);
1434 double d7 = (1. + (lp1s * Ksqrd));
1435 double d8 = (d1 * d6 * d7);
1436 double d8_1 = d8 * bhGp( (q+2), K, n, l, lp, rcsvV, GK );
1437 double d9 = (d5_1 - d8_1);
1438
1439 ASSERT( l == lp + 1 );
1440 ASSERT( Ksqrd != 0. );
1441 ASSERT( n2 != 0. );
1442
1443 ASSERT( lp1s != 0. );
1444 ASSERT( lp2s != 0. );
1445
1446 ASSERT( d1 != 0. );
1447 ASSERT( d2 != 0. );
1448 ASSERT( d3 != 0. );
1449 ASSERT( d4 != 0. );
1450 ASSERT( d5 != 0. );
1451 ASSERT( d6 != 0. );
1452 ASSERT( d7 != 0. );
1453 ASSERT( d8 != 0. );
1454 ASSERT( d9 != 0. );
1455
1456 rcsvV[rindx] = d9;
1457 return d9;
1458 }
1459 }
1460 else
1461 {
1462 ASSERT( rcsvV[rindx] != 0. );
1463 return rcsvV[rindx];
1464 }
1465}
1466
1467/***********************log version*******************************/
1469 long int q,
1470 double K,
1471 long int n,
1472 long int l,
1473 long int lp,
1474 /* Temporary storage for intermediate */
1475 /* results of the recursive routine */
1476 mxq *rcsvV_mxq,
1477 const mx& GK_mx
1478)
1479{
1480 DEBUG_ENTRY( "bhGp_mx()" );
1481
1482 /* static long int rcsv_Level = 1; */
1483 /* printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */
1484
1485 ASSERT( l == lp + 1 );
1486
1487 long int rindx = 2*q;
1488
1489 if( rcsvV_mxq[rindx].q == 0 )
1490 {
1491 /* SPECIAL CASE: n = l+1 = l'+2 */
1492 if( q == n - 1 )
1493 {
1494 /******************************************************/
1495 /* (1+(n K)^2) */
1496 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
1497 /* 2n */
1498 /******************************************************/
1499
1500 double Ksqrd = (K * K);
1501 double n2 = (double)(n*n);
1502
1503 double dd1 = (double) (2 * n);
1504 double dd2 = ( 1. + ( n2 * Ksqrd));
1505 double dd3 = dd2/dd1;
1506
1507 mx dd3_mx = mxify( dd3 );
1508 mx G1_mx = mult_mx( dd3_mx, GK_mx);
1509
1510 normalize_mx( G1_mx );
1511
1512 ASSERT( l == lp + 1 );
1513 ASSERT( Ksqrd != 0. );
1514 ASSERT( n2 != 0. );
1515 ASSERT( dd1 != 0. );
1516 ASSERT( dd2 != 0. );
1517
1518 rcsvV_mxq[rindx].q = 1;
1519 rcsvV_mxq[rindx].mx = G1_mx;
1520 return G1_mx;
1521 }
1522 /* SPECIAL CASE: n = l+2 = l'+3 */
1523 else if( q == (n - 2) )
1524 {
1525 /****************************************************************/
1526 /* */
1527 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/
1528 /* */
1529 /****************************************************************/
1530 /* (1+(n K)^2) */
1531 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
1532 /* 2n */
1533 /****************************************************************/
1534
1535 double Ksqrd = (K*K);
1536 double n2 = (double)(n*n);
1537 double dd1 = (double)(2 * n);
1538 double dd2 = ( 1. + ( n2 * Ksqrd) );
1539 double dd3 = (dd2/dd1);
1540 double dd4 = (double) ((2 * n) - 1);
1541 double dd5 = (double) (n - 1);
1542 double dd6 = (4. + (dd5 * dd2));
1543 double dd7 = dd4 * dd6;
1544
1545 /****************************************************************/
1546 /* */
1547 /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/
1548 /* */
1549 /****************************************************************/
1550
1551 mx dd3_mx = mxify( dd3 );
1552 mx dd7_mx = mxify( dd7 );
1553 mx G1_mx = mult_mx( dd3_mx, GK_mx );
1554 mx G2_mx = mult_mx( dd7_mx, G1_mx );
1555
1556 normalize_mx( G2_mx );
1557
1558 ASSERT( l == lp + 1 );
1559 ASSERT( Ksqrd != 0. );
1560 ASSERT( n2 != 0. );
1561 ASSERT( dd1 != 0. );
1562 ASSERT( dd2 != 0. );
1563 ASSERT( dd3 != 0. );
1564 ASSERT( dd4 != 0. );
1565 ASSERT( dd5 != 0. );
1566 ASSERT( dd6 != 0. );
1567 ASSERT( dd7 != 0. );
1568
1569 rcsvV_mxq[rindx].q = 1;
1570 rcsvV_mxq[rindx].mx = G2_mx;
1571 return G2_mx;
1572 }
1573 /* The GENERAL CASE n > l + 2*/
1574 else
1575 {
1576 /**************************************************************************************/
1577 /* 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) */
1578 /* */
1579 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1580 /* */
1581 /* FROM Eq. 3.24 */
1582 /* */
1583 /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */
1584 /* */
1585 /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */
1586 /**************************************************************************************/
1587
1588 long int lp1 = (q + 1);
1589 long int lp2 = (q + 2);
1590
1591 double Ksqrd = (K * K);
1592 double n2 = (double)(n * n);
1593 double lp1s = (double)(lp1 * lp1);
1594 double lp2s = (double)(lp2 * lp2);
1595
1596 double d1 = (4. * n2);
1597 double d2 = (4. * lp1s);
1598 double d3 = (double)((lp1)*((2 * q) + 3));
1599 double d4 = (1. + (n2 * Ksqrd));
1600 /* [ 4n^2 - 4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] */
1601 double d5 = d1 - d2 + (d3 * d4);
1602
1603 /* (n^2-(l+2)^2) */
1604 double d6 = (n2 - lp2s);
1605
1606 /* [ 1+((l+1)K)^2 ] */
1607 double d7 = (1. + (lp1s * Ksqrd));
1608
1609 /* { 4n^2 (n^2-(l+1)^2)[ 1+((l+1) K)^2 ] } */
1610 double d8 = (d1 * d6 * d7);
1611
1612 /**************************************************************************************/
1613 /* 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) */
1614 /* */
1615 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1616 /**************************************************************************************/
1617
1618 mx d5_mx=mxify( d5 );
1619 mx d8_mx=mxify( d8 );
1620
1621 mx t0_mx = bhGp_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1622 mx t1_mx = bhGp_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1623
1624 mx d9_mx = mult_mx( d5_mx, t0_mx );
1625 mx d10_mx = mult_mx( d8_mx, t1_mx );
1626
1627 mx result_mx = sub_mx( d9_mx, d10_mx );
1628 normalize_mx( result_mx );
1629
1630 ASSERT( d1 != 0. );
1631 ASSERT( d2 != 0. );
1632 ASSERT( d3 != 0. );
1633 ASSERT( d4 != 0. );
1634 ASSERT( d5 != 0. );
1635 ASSERT( d6 != 0. );
1636 ASSERT( d7 != 0. );
1637 ASSERT( d8 != 0. );
1638
1639 ASSERT( l == lp + 1 );
1640 ASSERT( Ksqrd != 0. );
1641 ASSERT( n2 != 0. );
1642 ASSERT( lp1s != 0. );
1643 ASSERT( lp2s != 0. );
1644
1645 rcsvV_mxq[rindx].q = 1;
1646 rcsvV_mxq[rindx].mx = result_mx;
1647 return result_mx;
1648 }
1649 }
1650 else
1651 {
1652 ASSERT( rcsvV_mxq[rindx].q != 0 );
1653 rcsvV_mxq[rindx].q = 1;
1654 return rcsvV_mxq[rindx].mx;
1655 }
1656}
1657
1658/************************************************************************************************/
1659/* *** bhGm.c *** */
1660/* */
1661/* Here we calculate G(n,l; K,l') with the recursive formula */
1662/* equation: (3.23) */
1663/* */
1664/* G(n,l-2; K,l-1) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */
1665/* */
1666/* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */
1667/* */
1668/* Under the transformation l -> l + 2 this gives */
1669/* */
1670/* 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) */
1671/* */
1672/* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+2+1)^2 K^2 ] G(n,l+2; K,l+2+1) */
1673/* */
1674/* */
1675/* or */
1676/* */
1677/* 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) */
1678/* */
1679/* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */
1680/* */
1681/* */
1682/* from the reference */
1683/* M. Brocklehurst */
1684/* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
1685/* */
1686/* */
1687/* * This is valid for the case l=l'-1 * */
1688/* * CASE: l = l'-1 * */
1689/* * Here the p in bhGm() refers * */
1690/* * to the Minus sign(-) in l=l'-1 * */
1691/************************************************************************************************/
1692
1693#if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1694#pragma optimize("", off)
1695#endif
1697 long int q,
1698 double K,
1699 long int n,
1700 long int l,
1701 long int lp,
1702 double *rcsvV,
1703 double GK
1704)
1705{
1706 DEBUG_ENTRY( "bhGm()" );
1707
1708 ASSERT( l == lp - 1 );
1709 ASSERT( l < n );
1710
1711 long int rindx = 2*q+1;
1712
1713 if( rcsvV[rindx] == 0. )
1714 {
1715 /* CASE: l = n - 1 */
1716 if( q == n - 1 )
1717 {
1718 ASSERT( l == lp - 1 );
1719 rcsvV[rindx] = GK;
1720 return GK;
1721 }
1722 /* CASE: l = n - 2 */
1723 else if( q == n - 2 )
1724 {
1725 double dd1 = 0.;
1726 double dd2 = 0.;
1727
1728 double G2 = 0.;
1729
1730 double Ksqrd = 0.;
1731 double n1 = 0.;
1732 double n2 = 0.;
1733
1734 ASSERT(l == lp - 1);
1735
1736 /* K^2 */
1737 Ksqrd = K * K;
1738 ASSERT( Ksqrd != 0. );
1739
1740 /* n */
1741 n1 = (double)n;
1742 ASSERT( n1 != 0. );
1743
1744 /* n^2 */
1745 n2 = (double)(n*n);
1746 ASSERT( n2 != 0. );
1747
1748 /* equation: (3.20) */
1749 /* G(n,n-2; K,n-1) = */
1750 /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */
1751 dd1 = (double) ((2 * n) - 1);
1752 ASSERT( dd1 != 0. );
1753
1754 dd2 = (1. + (n2 * Ksqrd));
1755 ASSERT( dd2 != 0. );
1756
1757 G2 = dd1 * dd2 * n1 * GK;
1758 ASSERT( G2 != 0. );
1759
1760 rcsvV[rindx] = G2;
1761 ASSERT( G2 != 0. );
1762 return G2;
1763 }
1764 else
1765 {
1766 long int lp2 = (q + 2);
1767 long int lp3 = (q + 3);
1768
1769 double lp2s = (double)(lp2 * lp2);
1770 double lp3s = (double)(lp3 * lp3);
1771
1772 /******************************************************************************/
1773 /* 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) */
1774 /* */
1775 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */
1776 /* */
1777 /* */
1778 /* FROM Eq. 3.23 */
1779 /* */
1780 /* 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) */
1781 /* */
1782 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */
1783 /******************************************************************************/
1784
1785 double Ksqrd = (K*K);
1786 double n2 = (double)(n*n);
1787 double d1 = (4. * n2);
1788 double d2 = (4. * lp2s);
1789 double d3 = (double)(lp2)*((2*q)+3);
1790 double d4 = (1. + (n2 * Ksqrd));
1791 /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */
1792 double d5 = d1 - d2 + (d3 * d4);
1793
1794 /******************************************************************************/
1795 /* 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) */
1796 /* */
1797 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */
1798 /******************************************************************************/
1799
1800 double d6 = (n2 - lp2s);
1801 /* [ 1+((l+3)K)^2 ] */
1802 double d7 = (1. + (lp3s * Ksqrd));
1803 /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */
1804 double d8 = d1 * d6 * d7;
1805 double d9 = d5 * bhGm( (q+1), K, n, l, lp, rcsvV, GK );
1806 double d10 = d8 * bhGm( (q+2), K, n, l, lp, rcsvV, GK );
1807 double d11 = d9 - d10;
1808
1809 ASSERT( l == lp - 1 );
1810 ASSERT( lp2s != 0. );
1811 ASSERT( Ksqrd != 0. );
1812 ASSERT( n2 != 0. );
1813 ASSERT( d1 != 0. );
1814 ASSERT( d2 != 0. );
1815 ASSERT( d3 != 0. );
1816 ASSERT( d4 != 0. );
1817 ASSERT( d5 != 0. );
1818 ASSERT( d6 != 0. );
1819 ASSERT( d7 != 0. );
1820 ASSERT( d8 != 0. );
1821 ASSERT( d9 != 0. );
1822 ASSERT( d10 != 0. );
1823 ASSERT( lp3s != 0. );
1824
1825 rcsvV[rindx] = d11;
1826 return d11;
1827 }
1828 }
1829 else
1830 {
1831 ASSERT( rcsvV[rindx] != 0. );
1832 return rcsvV[rindx];
1833 }
1834}
1835#if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1836#pragma optimize("", on)
1837#endif
1838
1839/************************log version***********************************/
1841 long int q,
1842 double K,
1843 long int n,
1844 long int l,
1845 long int lp,
1846 mxq *rcsvV_mxq,
1847 const mx& GK_mx
1848)
1849{
1850 DEBUG_ENTRY( "bhGm_mx()" );
1851
1852 /*static long int rcsv_Level = 1; */
1853 /*printf( "bhGm(): recursion level:\t%li\n",rcsv_Level++ ); */
1854
1855 ASSERT( l == lp - 1 );
1856 ASSERT( l < n );
1857
1858 long int rindx = 2*q+1;
1859
1860 if( rcsvV_mxq[rindx].q == 0 )
1861 {
1862 /* CASE: l = n - 1 */
1863 if( q == n - 1 )
1864 {
1865 mx result_mx = GK_mx;
1866 normalize_mx( result_mx );
1867
1868 rcsvV_mxq[rindx].q = 1;
1869 rcsvV_mxq[rindx].mx = result_mx;
1870
1871 ASSERT(l == lp - 1);
1872 return result_mx;
1873 }
1874 /* CASE: l = n - 2 */
1875 else if( q == n - 2 )
1876 {
1877 double Ksqrd = (K * K);
1878 double n1 = (double)n;
1879 double n2 = (double) (n*n);
1880 double dd1 = (double)((2 * n) - 1);
1881 double dd2 = (1. + (n2 * Ksqrd));
1882 /*(2n-1)(1+(n K)^2) n*/
1883 double dd3 = (dd1*dd2*n1);
1884
1885 /******************************************************/
1886 /* G(n,n-2; K,n-1) = */
1887 /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */
1888 /******************************************************/
1889
1890 mx dd3_mx = mxify( dd3 );
1891 mx G2_mx = mult_mx( dd3_mx, GK_mx );
1892
1893 normalize_mx( G2_mx );
1894
1895 ASSERT( l == lp - 1);
1896 ASSERT( n1 != 0. );
1897 ASSERT( n2 != 0. );
1898 ASSERT( dd1 != 0. );
1899 ASSERT( dd2 != 0. );
1900 ASSERT( dd3 != 0. );
1901 ASSERT( Ksqrd != 0. );
1902
1903 rcsvV_mxq[rindx].q = 1;
1904 rcsvV_mxq[rindx].mx = G2_mx;
1905 return G2_mx;
1906 }
1907 else
1908 {
1909 /******************************************************************************/
1910 /* 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) */
1911 /* */
1912 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */
1913 /* */
1914 /* */
1915 /* FROM Eq. 3.23 */
1916 /* */
1917 /* 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) */
1918 /* */
1919 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */
1920 /******************************************************************************/
1921
1922 long int lp2 = (q + 2);
1923 long int lp3 = (q + 3);
1924
1925 double lp2s = (double)(lp2 * lp2);
1926 double lp3s = (double)(lp3 * lp3);
1927 double n2 = (double)(n*n);
1928 double Ksqrd = (K * K);
1929
1930 /******************************************************/
1931 /* [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] */
1932 /******************************************************/
1933
1934 double d1 = (4. * n2);
1935 double d2 = (4. * lp2s);
1936 double d3 = (double)(lp2)*((2*q)+3);
1937 double d4 = (1. + (n2 * Ksqrd));
1938 /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */
1939 double d5 = d1 - d2 + (d3 * d4);
1940
1941 mx d5_mx=mxify(d5);
1942
1943 /******************************************************/
1944 /* 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] */
1945 /******************************************************/
1946
1947 double d6 = (n2 - lp2s);
1948 double d7 = (1. + (lp3s * Ksqrd));
1949 /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */
1950 double d8 = d1 * d6 * d7;
1951
1952 mx d8_mx = mxify(d8);
1953
1954 /******************************************************************************/
1955 /* 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) */
1956 /* */
1957 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */
1958 /******************************************************************************/
1959
1960 mx d9_mx = bhGm_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1961 mx d10_mx = bhGm_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1962 mx d11_mx = mult_mx( d5_mx, d9_mx );
1963 mx d12_mx = mult_mx( d8_mx, d10_mx);
1964 mx result_mx = sub_mx( d11_mx , d12_mx );
1965 rcsvV_mxq[rindx].q = 1;
1966 rcsvV_mxq[rindx].mx = result_mx;
1967
1968 ASSERT(l == lp - 1);
1969 ASSERT(n2 != 0. );
1970 ASSERT(lp2s != 0.);
1971 ASSERT( lp3s != 0.);
1972 ASSERT(Ksqrd != 0.);
1973
1974 ASSERT(d1 != 0.);
1975 ASSERT(d2 != 0.);
1976 ASSERT(d3 != 0.);
1977 ASSERT(d4 != 0.);
1978 ASSERT(d5 != 0.);
1979 ASSERT(d6 != 0.);
1980 ASSERT(d7 != 0.);
1981 ASSERT(d8 != 0.);
1982 return result_mx;
1983 }
1984 }
1985 else
1986 {
1987 ASSERT( rcsvV_mxq[rindx].q != 0 );
1988 return rcsvV_mxq[rindx].mx;
1989 }
1990}
1991
1992/****************************************************************************************/
1993/* */
1994/* bhg.c */
1995/* */
1996/* From reference; */
1997/* M. Brocklehurst */
1998/* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
1999/* */
2000/* */
2001/* We wish to compute the following function, */
2002/* */
2003/* equation: (3.17) */
2004/* - s=l' - (1/2) */
2005/* | (n+l)! ----- | */
2006/* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */
2007/* | (n-l-1)! | | | */
2008/* - s=0 - */
2009/* */
2010/* Using various recursion relations (for l'=l+1) */
2011/* */
2012/* equation: (3.23) */
2013/* G(n,l-2; K,l-1) = [ 4n^2-4l^2+l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */
2014/* */
2015/* - 4n^2 (n^2-l^2)[1+(l+1)^2 K^2] G(n,l; K,l+1) */
2016/* */
2017/* and (for l'=l-1) */
2018/* */
2019/* equation: (3.24) */
2020/* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l; K,l-1) */
2021/* */
2022/* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l; K,l+1) */
2023/* */
2024/* */
2025/* the starting point for the recursion relations are: */
2026/* */
2027/* */
2028/* equation (3.18): */
2029/* */
2030/* | pi |(1/2) 8n */
2031/* G(n,n-1; 0,n) = | -- | ------- (4n)^2 exp(-2n) */
2032/* | 2 | (2n-1)! */
2033/* */
2034/* equation (3.20): */
2035/* */
2036/* exp(2n-2/K tan^(-1)(n K) */
2037/* G(n,n-1; K,n) = --------------------------------------- */
2038/* sqrt(1 - exp(-2 pi/ K)) * (1+(n K)^(n+2) */
2039/* */
2040/* */
2041/* */
2042/* equation (3.20): */
2043/* G(n,n-2; K,n-1) = (2n-2)(1+(n K)^2) n G(n,n-1; K,n) */
2044/* */
2045/* */
2046/* equation (3.21): */
2047/* */
2048/* (1+(n K)^2) */
2049/* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
2050/* 2n */
2051/****************************************************************************************/
2052
2054 double K,
2055 long int n,
2056 long int l,
2057 long int lp,
2058 /* Temporary storage for intermediate */
2059 /* results of the recursive routine */
2060 double *rcsvV
2061)
2062{
2063 DEBUG_ENTRY( "bhg()" );
2064
2065 double ld1 = factorial( n + l );
2066 double ld2 = factorial( n - l - 1 );
2067 double ld3 = (ld1 / ld2);
2068
2069 double partprod = local_product( K , lp );
2070
2071 /**************************************************************************************/
2072 /* equation: (3.17) */
2073 /* - s=l' - (1/2) */
2074 /* | (n+l)! ----- | */
2075 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */
2076 /* | (n-l-1)! | | | */
2077 /* - s=0 - */
2078 /**************************************************************************************/
2079
2080 /**********************************************/
2081 /* - s=l' - (1/2) */
2082 /* | (n+l)! ----- | */
2083 /* | -------- | | (1 + s^2 K^2) | */
2084 /* | (n-l-1)! | | | */
2085 /* - s=0 - */
2086 /**********************************************/
2087
2088 double d2 = sqrt( ld3 * partprod );
2089 double d3 = powi( (2 * n) , (l - n) );
2090 double d4 = bhG( K, n, l, lp, rcsvV );
2091 double d5 = (d2 * d3);
2092 double d6 = (d5 * d4);
2093
2094 ASSERT(K != 0.);
2095 ASSERT( (n+l) >= 1 );
2096 ASSERT( ((n-l)-1) >= 0 );
2097
2098 ASSERT( partprod != 0. );
2099
2100 ASSERT( ld1 != 0. );
2101 ASSERT( ld2 != 0. );
2102 ASSERT( ld3 != 0. );
2103
2104 ASSERT( d2 != 0. );
2105 ASSERT( d3 != 0. );
2106 ASSERT( d4 != 0. );
2107 ASSERT( d5 != 0. );
2108 ASSERT( d6 != 0. );
2109 return d6;
2110}
2111
2112/********************log version**************************/
2114 double K,
2115 long int n,
2116 long int l,
2117 long int lp,
2118 /* Temporary storage for intermediate */
2119 /* results of the recursive routine */
2120 mxq *rcsvV_mxq
2121)
2122{
2123 /**************************************************************************************/
2124 /* equation: (3.17) */
2125 /* - s=l' - (1/2) */
2126 /* | (n+l)! ----- | */
2127 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */
2128 /* | (n-l-1)! | | | */
2129 /* - s=0 - */
2130 /**************************************************************************************/
2131
2132 DEBUG_ENTRY( "bhg_log()" );
2133
2134 double d1 = (double)(2*n);
2135 double d2 = (double)(l-n);
2136 double Ksqrd = (K*K);
2137
2138 /**************************************************************************************/
2139 /* */
2140 /* | (n+l)! | */
2141 /* log10 | -------- | = log10((n+1)!) - log10((n-l-1)!) */
2142 /* | (n-l-1)! | */
2143 /* */
2144 /**************************************************************************************/
2145
2146 double ld1 = lfactorial( n + l );
2147 double ld2 = lfactorial( n - l - 1 );
2148
2149 /**********************************************************************/
2150 /* | s=l' | s=l' */
2151 /* | ----- | --- */
2152 /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */
2153 /* | | | | --- */
2154 /* | s=0 | s=0 */
2155 /**********************************************************************/
2156
2157 double ld3 = log10_prodxx( lp, Ksqrd );
2158
2159 /**********************************************/
2160 /* - s=l' - (1/2) */
2161 /* | (n+l)! ----- | */
2162 /* | -------- | | (1 + s^2 K^2) | */
2163 /* | (n-l-1)! | | | */
2164 /* - s=0 - */
2165 /**********************************************/
2166
2167 /***********************************************************************/
2168 /* */
2169 /* | - s=l' - (1/2) | */
2170 /* | | (n+l)! ----- | | */
2171 /* log10| | -------- | | (1 + s^2 K^2) | | == */
2172 /* | | (n-l-1)! | | | | */
2173 /* | - s=0 - | */
2174 /* */
2175 /* | | s=l' | | */
2176 /* | | (n+l)! | | ----- | | */
2177 /* (1/2)* |log10 | -------- | + log10 | | | (1 + s^2 K^2) | | */
2178 /* | | (n-l-1)! | | | | | | */
2179 /* | | s=0 | | */
2180 /* */
2181 /***********************************************************************/
2182
2183 double ld4 = (1./2.) * ( ld3 + ld1 - ld2 );
2184
2185 /**********************************************/
2186 /* (2n)^(l-n) */
2187 /**********************************************/
2188 /* log10( 2n^(L-n) ) = (L-n) log10( 2n ) */
2189 /**********************************************/
2190
2191 double ld5 = d2 * log10( d1 );
2192
2193 /**************************************************************************************/
2194 /* equation: (3.17) */
2195 /* - s=l' - (1/2) */
2196 /* | (n+l)! ----- | */
2197 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) * G(n,l; K,l') */
2198 /* | (n-l-1)! | | | */
2199 /* - s=0 - */
2200 /**************************************************************************************/
2201
2202 /****************************************************/
2203 /* */
2204 /* - s=l' - (1/2) */
2205 /* | (n+l)! ----- | */
2206 /* | -------- | | (1 + s^2 K^2) | * (2n)^(L-n) */
2207 /* | (n-l-1)! | | | */
2208 /* - s=0 - */
2209 /****************************************************/
2210
2211 double ld6 = (ld5+ld4);
2212
2213 mx d6_mx = mxify_log10( ld6 );
2214 mx dd1_mx = bhG_mx( K, n, l, lp, rcsvV_mxq );
2215 mx dd2_mx = mult_mx( d6_mx, dd1_mx );
2216 normalize_mx( dd2_mx );
2217 double result = unmxify( dd2_mx );
2218
2219 ASSERT( result != 0. );
2220
2221 ASSERT( Ksqrd != 0. );
2222 ASSERT( ld3 >= 0. );
2223
2224 ASSERT( d1 > 0. );
2225 ASSERT( d2 < 0. );
2226 return result;
2227}
2228
2229inline double local_product( double K , long int lp )
2230{
2231 long int s = 0;
2232
2233 double Ksqrd =(K*K);
2234 double partprod = 1.;
2235
2236 for( s = 0; s <= lp; s = s + 1 )
2237 {
2238 double s2 = (double)(s*s);
2239
2240 /**************************/
2241 /* s=l' */
2242 /* ----- */
2243 /* | | (1 + s^2 K^2) */
2244 /* | | */
2245 /* s=0 */
2246 /**************************/
2247
2248 partprod *= ( 1. + ( s2 * Ksqrd ) );
2249 }
2250 return partprod;
2251}
2252
2253/************************************************************************/
2254/* Find the Einstein A's for hydrogen for a */
2255/* transition n,l --> n',l' */
2256/* */
2257/* In the following, we will assume n > n' */
2258/************************************************************************/
2259/*******************************************************************************/
2260/* */
2261/* Einstein A() for the transition from the */
2262/* initial state n,l to the finial state n',l' */
2263/* is given by oscillator f() */
2264/* */
2265/* hbar w max(l,l') | | 2 */
2266/* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */
2267/* 3 R_oo ( 2l + 1 ) | | */
2268/* */
2269/* */
2270/* E(n,l;n',l') max(l,l') | | 2 */
2271/* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
2272/* 3 R_oo ( 2l + 1 ) | | */
2273/* */
2274/* */
2275/* See for example Gordan Drake's */
2276/* Atomic, Molecular, & Optical Physics Handbook pg.638 */
2277/* */
2278/* Here R_oo is the infinite mass Rydberg length */
2279/* */
2280/* */
2281/* h c */
2282/* R_oo --- = 13.605698 eV */
2283/* {e} */
2284/* */
2285/* */
2286/* R_oo = 2.179874e-11 ergs */
2287/* */
2288/* w = omega */
2289/* = frequency of transition from n,l to n',l' */
2290/* */
2291/* */
2292/* */
2293/* here g_k are statistical weights obtained from */
2294/* the appropriate angular momentum quantum numbers */
2295/* */
2296/* */
2297/* - - 2 */
2298/* 64 pi^4 (e a_o)^2 max(l,l') | | */
2299/* A(n,l;n',l') = ------------------- ----------- v^3 | R(n,l;n',l') | */
2300/* 3 h c^3 2*l + 1 | | */
2301/* - - */
2302/* */
2303/* */
2304/* pi 3.141592654 */
2305/* plank_hbar 6.5821220 eV sec */
2306/* R_oo 2.179874e-11 ergs */
2307/* plank_h 6.6260755e-34 J sec */
2308/* e_charge 1.60217733e-19 C */
2309/* a_o 0.529177249e-10 m */
2310/* vel_light_c 299792458L m sec^-1 */
2311/* */
2312/* */
2313/* */
2314/* 64 pi^4 (e a_o)^2 64 pi^4 (a_o)^2 e^2 1 1 */
2315/* ----------------- = ----------------- -------- ---- = 7.5197711e-38 ----- */
2316/* 3 h c^3 3 c^2 hbar c 2 pi sec */
2317/* */
2318/* */
2319/* e^2 1 */
2320/* using ---------- = alpha = ---- */
2321/* hbar c 137 */
2322/*******************************************************************************/
2323
2324double H_Einstein_A(/* IN THE FOLLOWING WE HAVE n > n' */
2325 /* principal quantum number, 1 for ground, upper level */
2326 long int n,
2327 /* angular momentum, 0 for s */
2328 long int l,
2329 /* principal quantum number, 1 for ground, lower level */
2330 long int np,
2331 /* angular momentum, 0 for s */
2332 long int lp,
2333 /* Nuclear charge, 1 for H+, 2 for He++, etc */
2334 long int iz
2335)
2336{
2337 DEBUG_ENTRY( "H_Einstein_A()" );
2338
2339 double result;
2340 if( n > 60 || np > 60 )
2341 {
2342 result = H_Einstein_A_log10(n,l,np,lp,iz );
2343 }
2344 else
2345 {
2346 result = H_Einstein_A_lin(n,l,np,lp,iz );
2347 }
2348 return result;
2349}
2350
2351/************************************************************************/
2352/* Calculates the Einstein A's for hydrogen */
2353/* for the transition n,l --> n',l' */
2354/* units of sec^(-1) */
2355/* */
2356/* In the following, we have n > n' */
2357/************************************************************************/
2358STATIC double H_Einstein_A_lin(/* IN THE FOLLOWING WE HAVE n > n' */
2359 /* principal quantum number, 1 for ground, upper level */
2360 long int n,
2361 /* angular momentum, 0 for s */
2362 long int l,
2363 /* principal quantum number, 1 for ground, lower level */
2364 long int np,
2365 /* angular momentum, 0 for s */
2366 long int lp,
2367 /* Nuclear charge, 1 for H+, 2 for He++, etc */
2368 long int iz
2369)
2370{
2371 DEBUG_ENTRY( "H_Einstein_A_lin()" );
2372
2373 /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */
2374 double d1 = hv( n, np, iz );
2375 double d2 = d1 / HPLANCK; /* v = hv / h */
2376 double d3 = pow3(d2);
2377 double lg = (double)(l > lp ? l : lp);
2378 double Two_L_Plus_One = (double)(2*l + 1);
2379 double d6 = lg / Two_L_Plus_One;
2380 double d7 = hri( n, l, np, lp , iz );
2381 double d8 = d7 * d7;
2382 double result = CONST_ONE * d3 * d6 * d8;
2383
2384 /* validate the incoming data */
2385 if( n >= 70 )
2386 {
2387 fprintf(ioQQQ,"Principal Quantum Number `n' too large.\n");
2389 }
2390 if( iz < 1 )
2391 {
2392 fprintf(ioQQQ," The charge is impossible.\n");
2394 }
2395 if( n < 1 || np < 1 || l >= n || lp >= np )
2396 {
2397 fprintf(ioQQQ," The quantum numbers are impossible.\n");
2399 }
2400 if( n <= np )
2401 {
2402 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n");
2404 }
2405 return result;
2406}
2407
2408/**********************log version****************************/
2409double H_Einstein_A_log10(/* returns Einstein A in units of (sec)^-1 */
2410 long int n,
2411 long int l,
2412 long int np,
2413 long int lp,
2414 long int iz
2415)
2416{
2417 DEBUG_ENTRY( "H_Einstein_A_log10()" );
2418
2419 /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */
2420 double d1 = hv( n, np , iz );
2421 double d2 = d1 / HPLANCK; /* v = hv / h */
2422 double d3 = pow3(d2);
2423 double lg = (double)(l > lp ? l : lp);
2424 double Two_L_Plus_One = (double)(2*l + 1);
2425 double d6 = lg / Two_L_Plus_One;
2426 double d7 = hri_log10( n, l, np, lp , iz );
2427 double d8 = d7 * d7;
2428 double result = CONST_ONE * d3 * d6 * d8;
2429
2430 /* validate the incoming data */
2431 if( iz < 1 )
2432 {
2433 fprintf(ioQQQ," The charge is impossible.\n");
2435 }
2436 if( n < 1 || np < 1 || l >= n || lp >= np )
2437 {
2438 fprintf(ioQQQ," The quantum numbers are impossible.\n");
2440 }
2441 if( n <= np )
2442 {
2443 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n");
2445 }
2446 return result;
2447}
2448
2449/********************************************************************************/
2450/* hv calculates photon energy for n -> n' transitions for H and H-like ions */
2451/* simplest case of no "l" or "m" dependence */
2452/* epsilon_0 = 1 in vacu */
2453/* */
2454/* */
2455/* R_h */
2456/* Energy(n,Z) = - ------- */
2457/* n^2 */
2458/* */
2459/* */
2460/* */
2461/* Friedrich -- Theoretical Atomic Physics pg. 60 eq. 2.8 */
2462/* */
2463/* u */
2464/* R_h = --- R_oo where */
2465/* m_e */
2466/* */
2467/* h c */
2468/* R_oo --- = 2.179874e-11 ergs */
2469/* e */
2470/* */
2471/* (Harmin Lecture Notes for course phy-651 Spring 1994) */
2472/* where m_e (m_p) is the mass of and electron (proton) */
2473/* and u is the reduced electron mass for neutral hydrogen */
2474/* */
2475/* */
2476/* m_e m_p m_e */
2477/* u = --------- = ----------- */
2478/* m_e + m_p 1 + m_e/m_p */
2479/* */
2480/* m_e */
2481/* Now ----- = 0.000544617013 */
2482/* m_p */
2483/* u */
2484/* so that --- = 0.999455679 */
2485/* m_e */
2486/* */
2487/* */
2488/* returns energy of photon in ergs */
2489/* */
2490/* hv (n,n',Z) is for transitions n -> n' */
2491/* */
2492/* 1 erg = 1e-07 J */
2493/********************************************************************************/
2494/********************************************************************************/
2495/* WARNING: hv() use the electron reduced mass for hydrogen instead of */
2496/* the reduced mass associated with the apropriate ion */
2497/********************************************************************************/
2498
2499inline double hv( long int n, long int nprime, long int iz )
2500{
2501 DEBUG_ENTRY( "hv()" );
2502
2503 double n1 = (double)n;
2504 double n2 = n1*n1;
2505 double np1 = (double)nprime;
2506 double np2 = np1*np1;
2507 double rmr = 1./(1. + ELECTRON_MASS/PROTON_MASS); /* 0.999455679 */
2508 double izsqrd = (double)(iz*iz);
2509
2510 double d1 = 1. / n2;
2511 double d2 = 1. / np2;
2512 double d3 = izsqrd * rmr * EN1RYD;
2513 double d4 = d2 - d1;
2514 double result = d3 * d4;
2515
2516 ASSERT( n > 0 );
2517 ASSERT( nprime > 0 );
2518 ASSERT( n > nprime );
2519 ASSERT( iz > 0 );
2520 ASSERT( result > 0. );
2521
2522 if( n <= nprime )
2523 {
2524 fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n");
2526 }
2527
2528 return result;
2529}
2530
2531/************************************************************************/
2532/* hri() */
2533/* Calculate the hydrogen radial wavefunction integral */
2534/* for the dipole transition l'=l-1 or l'=l+1 */
2535/* for the higher energy state n,l to the lower energy state n',l' */
2536/* no "m" dependence */
2537/************************************************************************/
2538/* here we have a transition */
2539/* from the higher energy state n,l */
2540/* to the lower energy state n',l' */
2541/* with a dipole selection rule on l and l' */
2542/************************************************************************/
2543/* */
2544/* hri() test n,l,n',l' for domain errors and */
2545/* swaps n,l <--> n',l' for the case l'=l+1 */
2546/* */
2547/* It then calls hrii() */
2548/* */
2549/* Dec. 6, 1999 */
2550/* Robert Paul Bauman */
2551/************************************************************************/
2552
2553/************************************************************************/
2554/* This routine, hri(), calculates the hydrogen radial integral, */
2555/* for the transition n,l --> n',l' */
2556/* It is, of course, dimensionless. */
2557/* */
2558/* In the following, we have n > n' */
2559/************************************************************************/
2560
2561inline double hri(
2562 /* principal quantum number, 1 for ground, upper level */
2563 long int n,
2564 /* angular momentum, 0 for s */
2565 long int l,
2566 /* principal quantum number, 1 for ground, lower level */
2567 long int np,
2568 /* angular momentum, 0 for s */
2569 long int lp,
2570 /* Nuclear charge, 1 for H+, 2 for He++, etc */
2571 long int iz
2572)
2573{
2574 DEBUG_ENTRY( "hri" );
2575
2576 long int a;
2577 long int b;
2578 long int c;
2579 long int d;
2580 double ld1 = 0.;
2581 double Z = (double)iz;
2582
2583 /**********************************************************************/
2584 /* from higher energy -> lower energy */
2585 /* Selection Rule for l and l' */
2586 /* dipole process only */
2587 /**********************************************************************/
2588
2589 ASSERT( n > 0 );
2590 ASSERT( np > 0 );
2591 ASSERT( l >= 0 );
2592 ASSERT( lp >= 0 );
2593 ASSERT( n > l );
2594 ASSERT( np > lp );
2595 ASSERT( n > np || ( n == np && l == lp + 1 ));
2596 ASSERT( iz > 0 );
2597 ASSERT( lp == l + 1 || lp == l - 1 );
2598
2599 if( l == lp + 1 )
2600 {
2601 /* Keep variable the same */
2602 a = n;
2603 b = l;
2604 c = np;
2605 d = lp;
2606 }
2607 else if( l == lp - 1 )
2608 {
2609 /* swap n,l with n',l' */
2610 a = np;
2611 b = lp;
2612 c = n;
2613 d = l;
2614 }
2615 else
2616 {
2617 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2619 }
2620
2621 /**********************************************/
2622 /* Take care of the Z-dependence here. */
2623 /**********************************************/
2624 ld1 = hrii(a, b, c, d ) / Z;
2625
2626 return ld1;
2627}
2628
2629/************************************************************************/
2630/* hri_log10() */
2631/* Calculate the hydrogen radial wavefunction integral */
2632/* for the dipole transition l'=l-1 or l'=l+1 */
2633/* for the higher energy state n,l to the lower energy state n',l' */
2634/* no "m" dependence */
2635/************************************************************************/
2636/* here we have a transition */
2637/* from the higher energy state n,l */
2638/* to the lower energy state n',l' */
2639/* with a dipole selection rule on l and l' */
2640/************************************************************************/
2641/* */
2642/* hri_log10() test n,l,n',l' for domain errors and */
2643/* swaps n,l <--> n',l' for the case l'=l+1 */
2644/* */
2645/* It then calls hrii_log() */
2646/* */
2647/* Dec. 6, 1999 */
2648/* Robert Paul Bauman */
2649/************************************************************************/
2650
2651inline double hri_log10( long int n, long int l, long int np, long int lp , long int iz )
2652{
2653 /**********************************************************************/
2654 /* from higher energy -> lower energy */
2655 /* Selection Rule for l and l' */
2656 /* dipole process only */
2657 /**********************************************************************/
2658
2659 DEBUG_ENTRY( "hri_log10()" );
2660
2661 long int a;
2662 long int b;
2663 long int c;
2664 long int d;
2665 double ld1 = 0.;
2666 double Z = (double)iz;
2667
2668 ASSERT( n > 0);
2669 ASSERT( np > 0);
2670 ASSERT( l >= 0);
2671 ASSERT( lp >= 0 );
2672 ASSERT( n > l );
2673 ASSERT( np > lp );
2674 ASSERT( n > np || ( n == np && l == lp + 1 ));
2675 ASSERT( iz > 0 );
2676 ASSERT( lp == l + 1 || lp == l - 1 );
2677
2678 if( l == lp + 1)
2679 {
2680 /* Keep variable the same */
2681 a = n;
2682 b = l;
2683 c = np;
2684 d = lp;
2685 }
2686 else if( l == lp - 1 )
2687 {
2688 /* swap n,l with n',l' */
2689 a = np;
2690 b = lp;
2691 c = n;
2692 d = l;
2693 }
2694 else
2695 {
2696 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2698 }
2699
2700 /**********************************************/
2701 /* Take care of the Z-dependence here. */
2702 /**********************************************/
2703 ld1 = hrii_log(a, b, c, d ) / Z;
2704
2705 return ld1;
2706}
2707
2708STATIC double hrii( long int n, long int l, long int np, long int lp)
2709{
2710 /******************************************************************************/
2711 /* this routine hrii() is internal to the parent routine hri() */
2712 /* this internal routine only considers the case l=l'+1 */
2713 /* the case l=l-1 is done in the parent routine hri() */
2714 /* by the transformation n <--> n' and l <--> l' */
2715 /* THUS WE TEST FOR */
2716 /* l=l'-1 */
2717 /******************************************************************************/
2718
2719 DEBUG_ENTRY( "hrii()" );
2720
2721 long int a = 0, b = 0, c = 0;
2722 long int i1 = 0, i2 = 0, i3 = 0, i4 = 0;
2723
2724 char A='a';
2725
2726 double y = 0.;
2727 double fsf = 0.;
2728 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
2729 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
2730 double d00 = 0., d01 = 0.;
2731
2732 ASSERT( l == lp + 1 );
2733
2734 if( n == np ) /* SPECIAL CASE 1 */
2735 {
2736 /**********************************************************/
2737 /* if lp= l + 1 then it has higher energy */
2738 /* i.e. no photon */
2739 /* this is the second time we check this, oh well */
2740 /**********************************************************/
2741
2742 if( lp != (l - 1) )
2743 {
2744 printf( "BadMagic: Energy requirements not met.\n\n" );
2746 }
2747
2748 d2 = 3. / 2.;
2749 i1 = n * n;
2750 i2 = l * l;
2751 d5 = (double)(i1 - i2);
2752 d6 = sqrt(d5);
2753 d7 = (double)n * d6;
2754 d8 = d2 * d7;
2755 return d8;
2756 }
2757 else if( l == np && lp == (l - 1) ) /* A Pair of Easy Special Cases */
2758 {
2759 if( l == (n - 1) )
2760 {
2761 /**********************************************************************/
2762 /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */
2763 /* */
2764 /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */
2765 /* [(2n-1) - 1/(2n-1)]/4 */
2766 /**********************************************************************/
2767
2768 d1 = (double)( 2*n - 2 );
2769 d2 = (double)( 2*n - 1 );
2770 d3 = d1 * d2;
2771 d4 = sqrt( d3 );
2772
2773 d5 = (double)(4 * n * (n - 1));
2774 i1 = (2*n - 1);
2775 d6 = (double)(i1 * i1);
2776 d7 = d5/ d6;
2777 d8 = powi( d7, n );
2778
2779 d9 = 1./d2;
2780 d10 = d2 - d9;
2781 d11 = d10 / 4.;
2782
2783 /* Wrap it all up */
2784
2785 d12 = d4 * d8 * d11;
2786 return d12;
2787
2788 }
2789 else
2790 {
2791 /******************************************************************************/
2792 /* R(n,l;n',l') = R(n,l;l,l-1) */
2793 /* */
2794 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */
2795 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */
2796 /******************************************************************************/
2797
2798 d2 = 1.;
2799 for( i1 = -l; i1 <= l; i1 = i1 + 1 ) /* from n-l to n+l INCLUSIVE */
2800 {
2801 d1 = (double)(n - i1);
2802 d2 = d2 * d1;
2803 }
2804 i2 = (2*l - 1);
2805 d3 = factorial( i2 );
2806 d4 = d2/d3;
2807 d4 = sqrt( d4 );
2808
2809
2810 d5 = (double)( 4. * n * l );
2811 i3 = (n - l);
2812 d6 = (double)( i3 * i3 );
2813 d7 = d5 / d6;
2814 d8 = powi( d7, l+1 );
2815
2816
2817 i4 = n + l;
2818 d9 = (double)( i3 ) / (double)( i4 );
2819 d10 = powi( d9 , i4 );
2820
2821 d11 = d9 * d9;
2822 d12 = 1. - d11;
2823 d13 = d12 / 4.;
2824
2825 /* Wrap it all up */
2826 d14 = d4 * d8 * d10 * d13;
2827 return d14;
2828 }
2829 }
2830
2831 /*******************************************************************************************/
2832 /* THE GENERAL CASE */
2833 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
2834 /* REF: D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
2835 /* For F(a,b;c;x) we have from eq.4 */
2836 /* */
2837 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
2838 /* */
2839 /* a (1-x) (a + bx - c) */
2840 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
2841 /* (a-c) (a-c) */
2842 /* */
2843 /* */
2844 /* A similiar recusion relation holds for b with a <--> b. */
2845 /* */
2846 /* */
2847 /* we have initial conditions */
2848 /* */
2849 /* */
2850 /* F(0) = 1 with a = -1 */
2851 /* */
2852 /* b */
2853 /* F(-1) = 1 - (---) x with a = -1 */
2854 /* c */
2855 /*******************************************************************************************/
2856
2857 if( lp == l - 1 ) /* use recursion over "b" */
2858 {
2859 A='b';
2860 }
2861 else if( lp == l + 1 ) /* use recursion over "a" */
2862 {
2863 A='a';
2864 }
2865 else
2866 {
2867 printf(" BadMagic: Don't know what to do here.\n\n");
2869 }
2870
2871 /********************************************************************/
2872 /* Calculate the whole shootin match */
2873 /* - - (1/2) */
2874 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
2875 /* ----------- * | ---------------- | */
2876 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
2877 /* - - */
2878 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
2879 /* */
2880 /* This is used in the calculation of hydrogen */
2881 /* radial wave function integral for dipole transition case */
2882 /********************************************************************/
2883
2884 fsf = fsff( n, l, np );
2885
2886 /**************************************************************************************/
2887 /* Use a -> a' + 1 */
2888 /* _ _ */
2889 /* (a' + 1) (1 - x) | | */
2890 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */
2891 /* (a' + 1 -c) | | */
2892 /* - - */
2893 /* */
2894 /* For the first F() in the solution of the radial integral */
2895 /* */
2896 /* a = ( -n + l + 1 ) */
2897 /* */
2898 /* a = -n + l + 1 */
2899 /* max(a) = max(-n) + max(l) + 1 */
2900 /* = -n + max(n-1) + 1 */
2901 /* = -n + n-1 + 1 */
2902 /* = 0 */
2903 /* */
2904 /* similiarly */
2905 /* */
2906 /* min(a) = min(-n) + min(l) + 1 */
2907 /* = min(-n) + 0 + 1 */
2908 /* = -n + 1 */
2909 /* */
2910 /* a -> a' + 1 implies */
2911 /* */
2912 /* max(a') = -1 */
2913 /* min(a') = -n */
2914 /**************************************************************************************/
2915
2916 /* a plus */
2917 a = (-n + l + 1);
2918
2919 /* for the first 2_F_1 we use b = (-n' + l) */
2920 b = (-np + l);
2921
2922 /* c is simple */
2923 c = 2 * l;
2924
2925 /* -4 nn' */
2926 /* where Y = -------- . */
2927 /* (n-n')^2 */
2928 d2 = (double)(n - np);
2929 d3 = d2 * d2;
2930 d4 = 1. / d3;
2931 d5 = (double)(n * np);
2932 d6 = d5 * 4.;
2933 d7 = - d6;
2934 y = d7 * d4;
2935
2936 d00 = F21( a, b, c, y, A );
2937
2938 /**************************************************************/
2939 /* For the second F() in the solution of the radial integral */
2940 /* */
2941 /* a = ( -n + l - 1 ) */
2942 /* */
2943 /* a = -n + l + 1 */
2944 /* max(a) = max(-n) + max(l) - 1 */
2945 /* = -n + (n - 1) - 1 */
2946 /* = -2 */
2947 /* */
2948 /* similiarly */
2949 /* */
2950 /* min(a) = min(-n) + min(l) - 1 */
2951 /* = (-n) + 0 - 1 */
2952 /* = -n - 1 */
2953 /* */
2954 /* a -> a' + 1 implies */
2955 /* */
2956 /* max(a') = -3 */
2957 /* */
2958 /* min(a') = -n - 2 */
2959 /**************************************************************/
2960
2961 /* a minus */
2962 a = (-n + l - 1);
2963
2964 /* for the first 2_F_1 we use b = (-n' + l) */
2965 /* and does not change */
2966 b = (-np + l);
2967
2968 /* c is simple */
2969 c = 2 * l;
2970
2971 /**************************************************************/
2972 /* -4 nn' */
2973 /* where Y = -------- . */
2974 /* (n-n')^2 */
2975 /**************************************************************/
2976
2977 /**************************************************************/
2978 /* These are already calculated a few lines up */
2979 /* */
2980 /* d2 = (double) (n - np); */
2981 /* d3 = d2 * d2; */
2982 /* d4 = 1/ d3; */
2983 /* d5 = (double) (n * np); */
2984 /* d6 = d5 * 4.0; */
2985 /* d7 = - d6; */
2986 /* y = d7 * d4; */
2987 /**************************************************************/
2988
2989 d01 = F21(a, b, c, y, A );
2990
2991 /* Calculate */
2992 /* */
2993 /* (n-n')^2 */
2994 /* -------- */
2995 /* (n+n')^2 */
2996
2997 i1 = (n - np);
2998 d1 = pow2( (double)i1 );
2999 i2 = (n + np);
3000 d2 = pow2( (double)i2 );
3001 d3 = d1 / d2;
3002
3003 d4 = d01 * d3;
3004 d5 = d00 - d4;
3005 d6 = fsf * d5;
3006
3007 ASSERT( d6 != 0. );
3008 return d6;
3009}
3010
3011
3012STATIC double hrii_log( long int n, long int l, long int np, long int lp)
3013{
3014 /******************************************************************************/
3015 /* this routine hrii_log() is internal to the parent routine hri_log10() */
3016 /* this internal routine only considers the case l=l'+1 */
3017 /* the case l=l-1 is done in the parent routine hri_log10() */
3018 /* by the transformation n <--> n' and l <--> l' */
3019 /* THUS WE TEST FOR */
3020 /* l=l'-1 */
3021 /******************************************************************************/
3022 /**************************************************************************************/
3023 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */
3024 /* */
3025 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */
3026 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */
3027 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */
3028 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3029 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */
3030 /**************************************************************************************/
3031
3032 DEBUG_ENTRY( "hrii_log()" );
3033
3034 char A='a';
3035
3036 double y = 0.;
3037 double log10_fsf = 0.;
3038
3039 ASSERT( l == lp + 1 );
3040
3041 if( n == np ) /* SPECIAL CASE 1 */
3042 {
3043 /**********************************************************/
3044 /* if lp= l + 1 then it has higher energy */
3045 /* i.e. no photon */
3046 /* this is the second time we check this, oh well */
3047 /**********************************************************/
3048
3049 if( lp != (l - 1) )
3050 {
3051 printf( "BadMagic: l'= l+1 for n'= n.\n\n" );
3053 }
3054 else
3055 {
3056 /**********************************************************/
3057 /* 3 */
3058 /* R(nl:n'=n,l'=l+1) = --- n sqrt( n^2 - l^2 ) */
3059 /* 2 */
3060 /**********************************************************/
3061
3062 long int i1 = n * n;
3063 long int i2 = l * l;
3064
3065 double d1 = 3. / 2.;
3066 double d2 = (double)n;
3067 double d3 = (double)(i1 - i2);
3068 double d4 = sqrt(d3);
3069 double result = d1 * d2 * d4;
3070
3071 ASSERT( d3 >= 0. );
3072 return result;
3073 }
3074 }
3075 else if( l == np && lp == l - 1 ) /* A Pair of Easy Special Cases */
3076 {
3077 if( l == n - 1 )
3078 {
3079 /**********************************************************************/
3080 /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */
3081 /* */
3082 /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */
3083 /* [(2n-1) - 1/(2n-1)]/4 */
3084 /**********************************************************************/
3085
3086 double d1 = (double)((2*n-2)*(2*n-1));
3087 double d2 = sqrt( d1 );
3088 double d3 = (double)(4*n*(n-1));
3089 double d4 = (double)(2*n-1);
3090 double d5 = d4*d4;
3091 double d7 = d3/d5;
3092 double d8 = powi(d7,n);
3093 double d9 = 1./d4;
3094 double d10 = d4 - d9;
3095 double d11 = 0.25*d10;
3096 double result = (d2 * d8 * d11); /* Wrap it all up */
3097
3098 ASSERT( d1 >= 0. );
3099 ASSERT( d3 >= 0. );
3100 return result;
3101 }
3102 else
3103 {
3104 double result = 0.;
3105 double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0., ld7 = 0.;
3106
3107 /******************************************************************************/
3108 /* R(n,l;n',l') = R(n,l;l,l-1) */
3109 /* */
3110 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */
3111 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */
3112 /******************************************************************************/
3113 /**************************************/
3114 /* [(n-l) ... (n+l)] */
3115 /**************************************/
3116 /* log10[(n-l) ... (n+l)] = */
3117 /* */
3118 /* n+l */
3119 /* --- */
3120 /* > log10(j) */
3121 /* --- */
3122 /* j=n-l */
3123 /**************************************/
3124
3125 ld1 = 0.;
3126 for( long int i1 = (n-l); i1 <= (n+l); i1++ ) /* from n-l to n+l INCLUSIVE */
3127 {
3128 double d1 = (double)(i1);
3129 ld1 += log10( d1 );
3130 }
3131
3132 /**************************************/
3133 /* (2l-1)! */
3134 /**************************************/
3135 /* log10[ (2n-1)! ] */
3136 /**************************************/
3137
3138 ld2 = lfactorial( 2*l - 1 );
3139
3140 ASSERT( ((2*l)+1) >= 0);
3141
3142 /**********************************************/
3143 /* log10( [(n-l) ... (n+l)/(2l-1)!]^(1/2) ) = */
3144 /* (1/2) log10[(n-l) ... (n+l)] - */
3145 /* (1/2) log10[ (2n-1)! ] */
3146 /**********************************************/
3147
3148 ld3 = 0.5 * (ld1 - ld2);
3149
3150 /**********************************************/
3151 /* [4nl/(n-l)^2]^(l+1) */
3152 /**********************************************/
3153 /* log10( [4nl/(n-l)^2]^(l+1) ) = */
3154 /* (l+1) * log10( [4nl/(n-l)^2] ) */
3155 /* */
3156 /* = (l+1)*[ log10(4nl) - 2 log10(n-l) ] */
3157 /* */
3158 /**********************************************/
3159
3160 double d1 = (double)(l+1);
3161 double d2 = (double)(4*n*l);
3162 double d3 = (double)(n-l);
3163 double d4 = log10(d2);
3164 double d5 = log10(d3);
3165
3166 ld4 = d1 * (d4 - 2.*d5);
3167
3168 /**********************************************/
3169 /* [(n-l)/(n+l)]^(n+l) */
3170 /**********************************************/
3171 /* log10( [ (n-l)/(n+l) ]^(n+l) ) = */
3172 /* */
3173 /* (n+l) * [ log10(n-l) - log10(n+l) ] */
3174 /* */
3175 /**********************************************/
3176
3177 d1 = (double)(n-l);
3178 d2 = (double)(n+l);
3179 d3 = log10( d1 );
3180 d4 = log10( d2 );
3181
3182 ld5 = d2 * (d3 - d4);
3183
3184 /**********************************************/
3185 /* {1-[(n-l)/(n+l)]^2}/4 */
3186 /**********************************************/
3187 /* log10[ {1-[(n-l)/(n+l)]^2}/4 ] */
3188 /**********************************************/
3189
3190 d1 = (double)(n-l);
3191 d2 = (double)(n+l);
3192 d3 = d1/d2;
3193 d4 = d3*d3;
3194 d5 = 1. - d4;
3195 double d6 = 0.25*d5;
3196
3197 ld6 = log10(d6);
3198
3199 /******************************************************************************/
3200 /* R(n,l;n',l') = R(n,l;l,l-1) */
3201 /* */
3202 /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */
3203 /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */
3204 /******************************************************************************/
3205
3206 ld7 = ld3 + ld4 + ld5 + ld6;
3207
3208 result = pow( 10., ld7 );
3209
3210 ASSERT( result > 0. );
3211 return result;
3212 }
3213 }
3214 else
3215 {
3216 double result = 0.;
3217 long int a = 0, b = 0, c = 0;
3218 double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
3219 mx d00={0.0,0}, d01={0.0,0}, d02={0.0,0}, d03={0.0,0};
3220
3221 if( lp == l - 1 ) /* use recursion over "b" */
3222 {
3223 A='b';
3224 }
3225 else if( lp == l + 1 ) /* use recursion over "a" */
3226 {
3227 A='a';
3228 }
3229 else
3230 {
3231 printf(" BadMagic: Don't know what to do here.\n\n");
3233 }
3234
3235 /**************************************************************************************/
3236 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */
3237 /* */
3238 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */
3239 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */
3240 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */
3241 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3242 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */
3243 /**************************************************************************************/
3244
3245 /****************************************************************************************************/
3246 /* Calculate the whole shootin match */
3247 /* - - (1/2) */
3248 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
3249 /* fsff() = ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2)(n+n')^(-n-n') */
3250 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3251 /* - - */
3252 /* This is used in the calculation of hydrogen radial wave function integral for dipole transitions */
3253 /****************************************************************************************************/
3254
3255 log10_fsf = log10_fsff( n, l, np );
3256
3257 /******************************************************************************************/
3258 /* 2_F_1( a, b; c; y ) */
3259 /* */
3260 /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */
3261 /* */
3262 /* */
3263 /* Use a -> a' + 1 */
3264 /* _ _ */
3265 /* (a' + 1) (1 - x) | | */
3266 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */
3267 /* (a' + 1 - c) | | */
3268 /* - - */
3269 /* */
3270 /* For the first F() in the solution of the radial integral */
3271 /* */
3272 /* a = ( -n + l + 1 ) */
3273 /* */
3274 /* a = -n + l + 1 */
3275 /* max(a) = max(-n) + max(l) + 1 */
3276 /* = -n + max(n-1) + 1 */
3277 /* = -n + n-1 + 1 */
3278 /* = 0 */
3279 /* */
3280 /* similiarly */
3281 /* */
3282 /* min(a) = min(-n) + min(l) + 1 */
3283 /* = min(-n) + 0 + 1 */
3284 /* = -n + 1 */
3285 /* */
3286 /* a -> a' + 1 implies */
3287 /* */
3288 /* max(a') = -1 */
3289 /* min(a') = -n */
3290 /******************************************************************************************/
3291
3292 /* a plus */
3293 a = (-n + l + 1);
3294
3295 /* for the first 2_F_1 we use b = (-n' + l) */
3296 b = (-np + l);
3297
3298 /* c is simple */
3299 c = 2 * l;
3300
3301 /**********************************************************************************/
3302 /* 2_F_1( a, b; c; y ) */
3303 /* */
3304 /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */
3305 /* */
3306 /* -4 nn' */
3307 /* where Y = -------- . */
3308 /* (n-n')^2 */
3309 /* */
3310 /**********************************************************************************/
3311
3312 d2 = (double)(n - np);
3313 d3 = d2 * d2;
3314
3315 d4 = 1. / d3;
3316 d5 = (double)(n * np);
3317 d6 = d5 * 4.;
3318
3319 d7 = -d6;
3320 y = d7 * d4;
3321
3322 /**************************************************************************************************/
3323 /* THE GENERAL CASE */
3324 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
3325 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
3326 /* */
3327 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
3328 /* */
3329 /* a (1-x) (a + bx - c) */
3330 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
3331 /* (a - c) (a - c) */
3332 /* */
3333 /* */
3334 /* A similiar recusion relation holds for b with a <--> b. */
3335 /* */
3336 /* */
3337 /* we have initial conditions */
3338 /* */
3339 /* */
3340 /* F(0) = 1 with a = -1 */
3341 /* */
3342 /* b */
3343 /* F(-1) = 1 - (---) x with a = -1 */
3344 /* c */
3345 /**************************************************************************************************/
3346
3347 /* 2_F_1( long int a, long int b, long int c, (double) y, (string) "a" or "b") */
3348 /* F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3349 d00 = F21_mx( a, b, c, y, A );
3350
3351 /**************************************************************/
3352 /* For the second F() in the solution of the radial integral */
3353 /* */
3354 /* a = ( -n + l - 1 ) */
3355 /* */
3356 /* a = -n + l + 1 */
3357 /* max(a) = max(-n) + max(l) - 1 */
3358 /* = -n + (n - 1) - 1 */
3359 /* = -2 */
3360 /* */
3361 /* similiarly */
3362 /* */
3363 /* min(a) = min(-n) + min(l) - 1 */
3364 /* = (-n) + 0 - 1 */
3365 /* = -n - 1 */
3366 /* */
3367 /* a -> a' + 1 implies */
3368 /* */
3369 /* max(a') = -3 */
3370 /* */
3371 /* min(a') = -n - 2 */
3372 /**************************************************************/
3373
3374 /* a minus */
3375 a = (-n + l - 1);
3376
3377 /* for the first 2_F_1 we use b = (-n' + l) */
3378 /* and does not change */
3379 b = (-np + l);
3380
3381 /* c is simple */
3382 c = 2 * l;
3383
3384 /**************************************************************/
3385 /* -4 nn' */
3386 /* where Y = -------- . */
3387 /* (n-n')^2 */
3388 /**************************************************************/
3389
3390 /**************************************************************/
3391 /* These are already calculated a few lines up */
3392 /* */
3393 /* d2 = (double) (n - np); */
3394 /* d3 = d2 * d2; */
3395 /* d4 = 1/ d3; */
3396 /* d5 = (double) (n * np); */
3397 /* d6 = d5 * 4.0; */
3398 /* d7 = - d6; */
3399 /* y = d7 * d4; */
3400 /**************************************************************/
3401
3402 d01 = F21_mx(a, b, c, y, A );
3403
3404 /**************************************************************************************/
3405 /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */
3406 /* */
3407 /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */
3408 /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */
3409 /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */
3410 /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3411 /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */
3412 /* */
3413 /* = fsf * ( F(a,b,c;y) - d3 * F(a',b',c';y) ) */
3414 /* */
3415 /* where d3 = (n-n')^2 (n+n')^2 */
3416 /* */
3417 /**************************************************************************************/
3418
3419 /**************************************************************/
3420 /* Calculate */
3421 /* */
3422 /* (n-n')^2 */
3423 /* -------- */
3424 /* (n+n')^2 */
3425 /**************************************************************/
3426
3427 d1 = (double)((n - np)*(n -np));
3428 d2 = (double)((n + np)*(n + np));
3429 d3 = d1 / d2;
3430
3431 d02.x = d01.x;
3432 d02.m = d01.m * d3;
3433
3434 while ( fabs(d02.m) > 1.0e+25 )
3435 {
3436 d02.m /= 1.0e+25;
3437 d02.x += 25;
3438 }
3439
3440 d03.x = d00.x;
3441 d03.m = d00.m * (1. - (d02.m/d00.m) * powi( 10. , (d02.x - d00.x) ) );
3442
3443 result = pow( 10., (log10_fsf + d03.x) ) * d03.m;
3444
3445 ASSERT( result != 0. );
3446
3447 /* we don't care about the correct sign of result... */
3448 return fabs(result);
3449 }
3450 /* Shouldn't get here */
3451 printf(" This code should be inaccessible\n\n");
3453}
3454
3455STATIC double fsff( long int n, long int l, long int np )
3456{
3457 /****************************************************************/
3458 /* Calculate the whole shootin match */
3459 /* - - (1/2) */
3460 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
3461 /* ----------- * | ---------------- | */
3462 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3463 /* - - */
3464 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3465 /* */
3466 /****************************************************************/
3467
3468 DEBUG_ENTRY( "fsff()" );
3469
3470 long int i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
3471 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
3472 double sigma = 1.;
3473
3474 /****************************************************************
3475 * Calculate the whole shootin match *
3476 * (-1)^(n'-1) | (n+l)! (n'+l-1)! | *
3477 * ----------- * | ---------------- | *
3478 * [4 (2l-1)!] | (n-l-1)! (n'-l)! | *
3479 * - - *
3480 * * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') *
3481 * *
3482 ****************************************************************/
3483
3484 /* Calculate (-1)^(n'-l) */
3485 if( is_odd(np - l) )
3486 {
3487 sigma *= -1.;
3488 }
3489 ASSERT( sigma != 0. );
3490
3491 /*********************/
3492 /* Calculate (2l-1)! */
3493 /*********************/
3494 i1 = (2*l - 1);
3495 if( i1 < 0 )
3496 {
3497 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3499 }
3500
3501 /****************************************************************/
3502 /* Calculate the whole shootin match */
3503 /* - - (1/2) */
3504 /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
3505 /* ----------- * | ---------------- | */
3506 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3507 /* - - */
3508 /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3509 /* */
3510 /****************************************************************/
3511
3512 d0 = factorial( i1 );
3513 d1 = 4. * d0;
3514 d2 = 1. / d1;
3515
3516 /**********************************************************************/
3517 /* We want the (negitive) of this */
3518 /* since we really are interested in */
3519 /* [(2l-1)!]^-1 */
3520 /**********************************************************************/
3521
3522 sigma = sigma * d2;
3523 ASSERT( sigma != 0. );
3524
3525 /**********************************************************************/
3526 /* Calculate (4 n n')^(l+1) */
3527 /* powi( m , n) calcs m^n */
3528 /* returns long double with m,n ints */
3529 /**********************************************************************/
3530
3531 i0 = 4 * n * np;
3532 i1 = l + 1;
3533 d2 = powi( (double)i0 , i1 );
3534 sigma = sigma * d2;
3535 ASSERT( sigma != 0. );
3536
3537 /* Calculate (n-n')^(n+n'-2l-2) */
3538 i0 = n - np;
3539 i1 = n + np - 2*l - 2;
3540 d2 = powi( (double)i0 , i1 );
3541 sigma = sigma * d2;
3542 ASSERT( sigma != 0. );
3543
3544 /* Calculate (n+n')^(-n-n') */
3545 i0 = n + np;
3546 i1 = -n - np;
3547 d2 = powi( (double)i0 , i1 );
3548 sigma = sigma * d2;
3549 ASSERT( sigma != 0. );
3550
3551 /**********************************************************************/
3552 /* - - (1/2) */
3553 /* | (n+l)! (n'+l-1)! | */
3554 /* Calculate | ---------------- | */
3555 /* | (n-l-1)! (n'-l)! | */
3556 /* - - */
3557 /**********************************************************************/
3558
3559 i1 = n + l;
3560 if( i1 < 0 )
3561 {
3562 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3564 }
3565 d1 = factorial( i1 );
3566
3567 i2 = np + l - 1;
3568 if( i2 < 0 )
3569 {
3570 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3572 }
3573 d2 = factorial( i2 );
3574
3575 i3 = n - l - 1;
3576 if( i3 < 0 )
3577 {
3578 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3580 }
3581 d3 = factorial( i3 );
3582
3583 i4 = np - l;
3584 if( i4 < 0 )
3585 {
3586 printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3588 }
3589 d4 = factorial( i4 );
3590
3591 ASSERT( d3 != 0. );
3592 ASSERT( d4 != 0. );
3593
3594 /* Do this a different way to prevent overflow */
3595 /* d5 = (sqrt(d1 *d2)); */
3596 d5 = sqrt(d1)*sqrt(d2);
3597 d5 /= sqrt(d3);
3598 d5 /= sqrt(d4);
3599
3600 sigma = sigma * d5;
3601
3602 ASSERT( sigma != 0. );
3603 return sigma;
3604}
3605
3606/**************************log version*******************************/
3607STATIC double log10_fsff( long int n, long int l, long int np )
3608{
3609 /******************************************************************************************************/
3610 /* Calculate the whole shootin match */
3611 /* - - (1/2) */
3612 /* 1 | (n+l)! (n'+l-1)! | */
3613 /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3614 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3615 /* - - */
3616 /******************************************************************************************************/
3617
3618 DEBUG_ENTRY( "log10_fsff()" );
3619
3620 double d0 = 0., d1 = 0.;
3621 double ld0 = 0., ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0.;
3622 double result = 0.;
3623
3624 /******************************************************************************************************/
3625 /* Calculate the log10 of the whole shootin match */
3626 /* - - (1/2) */
3627 /* 1 | (n+l)! (n'+l-1)! | */
3628 /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3629 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3630 /* - - */
3631 /******************************************************************************************************/
3632
3633 /**********************/
3634 /* Calculate (2l-1)! */
3635 /**********************/
3636
3637 d0 = (double)(2*l - 1);
3638 ASSERT( d0 != 0. );
3639
3640 /******************************************************************************************************/
3641 /* Calculate the whole shootin match */
3642 /* - - (1/2) */
3643 /* 1 | (n+l)! (n'+l-1)! | */
3644 /* ----------- * | ---------------- | * (4 n n')^(l+1) |(n-n')^(n+n'-2l-2)| (n+n')^(-n-n') */
3645 /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3646 /* - - */
3647 /******************************************************************************************************/
3648
3649 ld0 = lfactorial( 2*l - 1 );
3650 ld1 = log10(4.);
3651 result = -(ld0 + ld1);
3652 ASSERT( result != 0. );
3653
3654 /**********************************************************************/
3655 /* Calculate (4 n n')^(l+1) */
3656 /* powi( m , n) calcs m^n */
3657 /* returns long double with m,n ints */
3658 /**********************************************************************/
3659
3660 d0 = (double)(4 * n * np);
3661 d1 = (double)(l + 1);
3662 result += d1 * log10(d0);
3663 ASSERT( d0 >= 0. );
3664 ASSERT( d1 != 0. );
3665
3666 /**********************************************************************/
3667 /* Calculate |(n-n')^(n+n'-2l-2)| */
3668 /* NOTE: Here we are interested only */
3669 /* magnitude of (n-n')^(n+n'-2l-2) */
3670 /**********************************************************************/
3671
3672 d0 = (double)(n - np);
3673 d1 = (double)(n + np - 2*l - 2);
3674 result += d1 * log10(fabs(d0));
3675 ASSERT( fabs(d0) > 0. );
3676 ASSERT( d1 != 0. );
3677
3678 /* Calculate (n+n')^(-n-n') */
3679 d0 = (double)(n + np);
3680 d1 = (double)(-n - np);
3681 result += d1 * log10(d0);
3682 ASSERT( d0 > 0. );
3683 ASSERT( d1 != 0. );
3684
3685 /**********************************************************************/
3686 /* - - (1/2) */
3687 /* | (n+l)! (n'+l-1)! | */
3688 /* Calculate | ---------------- | */
3689 /* | (n-l-1)! (n'-l)! | */
3690 /* - - */
3691 /**********************************************************************/
3692
3693 ASSERT( n+l > 0 );
3694 ld0 = lfactorial( n + l );
3695
3696 ASSERT( np+l-1 > 0 );
3697 ld1 = lfactorial( np + l - 1 );
3698
3699 ASSERT( n-l-1 >= 0 );
3700 ld2 = lfactorial( n - l - 1 );
3701
3702 ASSERT( np-l >= 0 );
3703 ld3 = lfactorial( np - l );
3704
3705 ld4 = 0.5*((ld0+ld1)-(ld2+ld3));
3706
3707 result += ld4;
3708 ASSERT( result != 0. );
3709 return result;
3710}
3711
3712/***************************************************************************/
3713/* Find the Oscillator Strength for hydrogen for any */
3714/* transition n,l --> n',l' */
3715/* returns a double */
3716/***************************************************************************/
3717
3718/************************************************************************/
3719/* Find the Oscillator Strength for hydrogen for any */
3720/* transition n,l --> n',l' */
3721/* returns a double */
3722/* */
3723/* Einstein A() for the transition from the */
3724/* initial state n,l to the finial state n',l' */
3725/* require the Oscillator Strength f() */
3726/* */
3727/* hbar w max(l,l') | | 2 */
3728/* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */
3729/* 3 R_oo ( 2l + 1 ) | | */
3730/* */
3731/* */
3732/* */
3733/* E(n,l;n',l') max(l,l') | | 2 */
3734/* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
3735/* 3 R_oo ( 2l + 1 ) | | */
3736/* */
3737/* */
3738/* See for example Gordan Drake's */
3739/* Atomic, Molecular, & Optical Physics Handbook pg.638 */
3740/* */
3741/* Here R_oo is the infinite mass Rydberg length */
3742/* */
3743/* */
3744/* h c */
3745/* R_oo --- = 13.605698 eV */
3746/* {e} */
3747/* */
3748/* */
3749/* R_oo = 2.179874e-11 ergs */
3750/* */
3751/* w = omega */
3752/* = frequency of transition from n,l to n',l' */
3753/* */
3754/* */
3755/* */
3756/* here g_k are statistical weights obtained from */
3757/* the appropriate angular momentum quantum numbers */
3758/************************************************************************/
3759
3760/********************************************************************************/
3761/* Calc the Oscillator Strength f(*) given by */
3762/* */
3763/* E(n,l;n',l') max(l,l') | | 2 */
3764/* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
3765/* 3 R_oo ( 2l + 1 ) | | */
3766/* */
3767/* See for example Gordan Drake's */
3768/* Atomic, Molecular, & Optical Physics Handbook pg.638 */
3769/********************************************************************************/
3770
3771/************************************************************************/
3772/* Calc the Oscillator Strength f(*) given by */
3773/* */
3774/* E(n,l;n',l') max(l,l') | | 2 */
3775/* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
3776/* 3 R_oo ( 2l + 1 ) | | */
3777/* */
3778/* f(n,l;n',l') is dimensionless. */
3779/* */
3780/* See for example Gordan Drake's */
3781/* Atomic, Molecular, & Optical Physics Handbook pg.638 */
3782/* */
3783/* In the following, we have n > n' */
3784/************************************************************************/
3785
3786inline double OscStr_f(
3787 /* IN THE FOLLOWING WE HAVE n > n' */
3788 /* principal quantum number, 1 for ground, upper level */
3789 long int n,
3790 /* angular momentum, 0 for s */
3791 long int l,
3792 /* principal quantum number, 1 for ground, lower level */
3793 long int np,
3794 /* angular momentum, 0 for s */
3795 long int lp,
3796 /* Nuclear charge, 1 for H+, 2 for He++, etc */
3797 long int iz
3798)
3799{
3800 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3801 long int i1 = 0, i2 = 0;
3802
3803 if( l > lp )
3804 i1 = l;
3805 else
3806 i1 = lp;
3807
3808 i2 = 2*lp + 1;
3809 d0 = 1. / 3.;
3810 d1 = (double)i1 / (double)i2;
3811 /* hv() returns energy in ergs */
3812 d2 = hv( n, np, iz );
3813 d3 = d2 / EN1RYD;
3814 d4 = hri( n, l, np, lp ,iz );
3815 d5 = d4 * d4;
3816
3817 d6 = d0 * d1 * d3 * d5;
3818
3819 return d6;
3820}
3821
3822/************************log version ***************************/
3823inline double OscStr_f_log10( long int n , long int l , long int np , long int lp , long int iz )
3824{
3825 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3826 long int i1 = 0, i2 = 0;
3827
3828 if( l > lp )
3829 i1 = l;
3830 else
3831 i1 = lp;
3832
3833 i2 = 2*lp + 1;
3834 d0 = 1. / 3.;
3835 d1 = (double)i1 / (double)i2;
3836 /* hv() returns energy in ergs */
3837 d2 = hv( n, np, iz );
3838 d3 = d2 / EN1RYD;
3839 d4 = hri_log10( n, l, np, lp ,iz );
3840 d5 = d4 * d4;
3841
3842 d6 = d0 * d1 * d3 * d5;
3843
3844 return d6;
3845}
3846
3847STATIC double F21( long int a , long int b, long int c, double y, char A )
3848{
3849 DEBUG_ENTRY( "F21()" );
3850
3851 double d1 = 0.;
3852 long int i1 = 0;
3853 double *yV;
3854
3855 /**************************************************************/
3856 /* A must be either 'a' or 'b' */
3857 /* and is use to determine which */
3858 /* variable recursion will be over */
3859 /**************************************************************/
3860
3861 ASSERT( A == 'a' || A == 'b' );
3862
3863 if( A == 'b' )
3864 {
3865 /* if the recursion is over "b" */
3866 /* then make it over "a" by switching these around */
3867 i1 = a;
3868 a = b;
3869 b = i1;
3870 A = 'a';
3871 }
3872
3873 /**************************************************************************************/
3874 /* malloc space for the (dynamic) 1-d array */
3875 /* 2_F_1 works via recursion and needs room to store intermediate results */
3876 /* Here the + 5 is a safety margin */
3877 /**************************************************************************************/
3878
3879 /*Must use CALLOC*/
3880 yV = (double*)CALLOC( sizeof(double), (size_t)(-a + 5) );
3881
3882 /**********************************************************************************************/
3883 /* begin sanity check, check order, and that none negative */
3884 /* THE GENERAL CASE */
3885 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
3886 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
3887 /* */
3888 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
3889 /* */
3890 /* a (1-x) (a + bx - c) */
3891 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
3892 /* (a-c) (a-c) */
3893 /* */
3894 /* */
3895 /* A similiar recusion relation holds for b with a <--> b. */
3896 /* */
3897 /* */
3898 /* we have initial conditions */
3899 /* */
3900 /* */
3901 /* F(0) = 1 with a = -1 */
3902 /* */
3903 /* b */
3904 /* F(-1) = 1 - (---) x with a = -1 */
3905 /* c */
3906 /* */
3907 /* For the first F() in the solution of the radial integral */
3908 /* */
3909 /* a = ( -n + l + 1 ) */
3910 /* */
3911 /* a = -n + l + 1 */
3912 /* max(a) = max(-n) + max(l) + 1 */
3913 /* = max(-n) + max(n - 1) + 1 */
3914 /* = max(-n + n - 1) + 1 */
3915 /* = max(-1) + 1 */
3916 /* = 0 */
3917 /* */
3918 /* similiarly */
3919 /* */
3920 /* min(a) = min(-n) + min(l) + 1 */
3921 /* = min(-n) + 0 + 1 */
3922 /* = (-n) + 0 + 1 */
3923 /* = -n + 1 */
3924 /* */
3925 /* a -> a' + 1 implies */
3926 /* */
3927 /* max(a') = -1 */
3928 /* min(a') = -n */
3929 /* */
3930 /* For the second F() in the solution of the radial integral */
3931 /* */
3932 /* a = ( -n + l - 1 ) */
3933 /* */
3934 /* a = -n + l + 1 */
3935 /* max(a) = max(-n) + max(l) - 1 */
3936 /* = -n + (n - 1) - 1 */
3937 /* = -2 */
3938 /* */
3939 /* similiarly */
3940 /* */
3941 /* min(a) = min(-n) + min(l) - 1 */
3942 /* = (-n) + 0 - 1 */
3943 /* = -n - 1 */
3944 /* */
3945 /* a -> a' + 1 implies */
3946 /* */
3947 /* max(a') = -3 */
3948 /* min(a') = -n - 2 */
3949 /**********************************************************************************************/
3950
3951 ASSERT( a <= 0 );
3952 ASSERT( b <= 0 );
3953 ASSERT( c >= 0 );
3954
3955 d1 = F21i(a, b, c, y, yV );
3956 free( yV );
3957 return d1;
3958}
3959
3960STATIC mx F21_mx( long int a , long int b, long int c, double y, char A )
3961{
3962 DEBUG_ENTRY( "F21_mx()" );
3963
3964 mx result_mx = {0.0,0};
3965 mxq *yV = NULL;
3966
3967 /**************************************************************/
3968 /* A must be either 'a' or 'b' */
3969 /* and is use to determine which */
3970 /* variable recursion will be over */
3971 /**************************************************************/
3972
3973 ASSERT( A == 'a' || A == 'b' );
3974
3975 if( A == 'b' )
3976 {
3977 /* if the recursion is over "b" */
3978 /* then make it over "a" by switching these around */
3979 long int i1 = a;
3980 a = b;
3981 b = i1;
3982 A = 'a';
3983 }
3984
3985 /**************************************************************************************/
3986 /* malloc space for the (dynamic) 1-d array */
3987 /* 2_F_1 works via recursion and needs room to store intermediate results */
3988 /* Here the + 5 is a safety margin */
3989 /**************************************************************************************/
3990
3991 /*Must use CALLOC*/
3992 yV = (mxq *)CALLOC( sizeof(mxq), (size_t)(-a + 5) );
3993
3994 /**********************************************************************************************/
3995 /* begin sanity check, check order, and that none negative */
3996 /* THE GENERAL CASE */
3997 /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
3998 /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
3999 /* */
4000 /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
4001 /* */
4002 /* a (1-x) (a + bx - c) */
4003 /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
4004 /* (a-c) (a-c) */
4005 /* */
4006 /* */
4007 /* A similiar recusion relation holds for b with a <--> b. */
4008 /* */
4009 /* */
4010 /* we have initial conditions */
4011 /* */
4012 /* */
4013 /* F(0) = 1 with a = -1 */
4014 /* */
4015 /* b */
4016 /* F(-1) = 1 - (---) x with a = -1 */
4017 /* c */
4018 /* */
4019 /* For the first F() in the solution of the radial integral */
4020 /* */
4021 /* a = ( -n + l + 1 ) */
4022 /* */
4023 /* a = -n + l + 1 */
4024 /* max(a) = max(-n) + max(l) + 1 */
4025 /* = max(-n) + max(n - 1) + 1 */
4026 /* = max(-n + n - 1) + 1 */
4027 /* = max(-1) + 1 */
4028 /* = 0 */
4029 /* */
4030 /* similiarly */
4031 /* */
4032 /* min(a) = min(-n) + min(l) + 1 */
4033 /* = min(-n) + 0 + 1 */
4034 /* = (-n) + 0 + 1 */
4035 /* = -n + 1 */
4036 /* */
4037 /* a -> a' + 1 implies */
4038 /* */
4039 /* max(a') = -1 */
4040 /* min(a') = -n */
4041 /* */
4042 /* For the second F() in the solution of the radial integral */
4043 /* */
4044 /* a = ( -n + l - 1 ) */
4045 /* */
4046 /* a = -n + l + 1 */
4047 /* max(a) = max(-n) + max(l) - 1 */
4048 /* = -n + (n - 1) - 1 */
4049 /* = -2 */
4050 /* */
4051 /* similiarly */
4052 /* */
4053 /* min(a) = min(-n) + min(l) - 1 */
4054 /* = (-n) + 0 - 1 */
4055 /* = -n - 1 */
4056 /* */
4057 /* a -> a' + 1 implies */
4058 /* */
4059 /* max(a') = -3 */
4060 /* min(a') = -n - 2 */
4061 /**********************************************************************************************/
4062
4063 ASSERT( a <= 0 );
4064 ASSERT( b <= 0 );
4065 ASSERT( c >= 0 );
4066
4067 result_mx = F21i_log(a, b, c, y, yV );
4068 free( yV );
4069 return result_mx;
4070}
4071
4072STATIC double F21i(long int a, long int b, long int c, double y, double *yV )
4073{
4074 DEBUG_ENTRY( "F21i()" );
4075
4076 double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
4077 double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
4078 long int i1 = 0, i2 = 0;
4079
4080 if( a == 0 )
4081 {
4082 return 1.;
4083 }
4084 else if( a == -1 )
4085 {
4086 ASSERT( c != 0 );
4087 d1 = (double)b;
4088 d2 = (double)c;
4089 d3 = y * (d1/d2);
4090 d4 = 1. - d3;
4091 return d4;
4092 }
4093 /* Check to see if y(-a) != 0 in a very round about way to avoid lclint:error:13 */
4094 else if( yV[-a] != 0. )
4095 {
4096 /* Return the stored result */
4097 return yV[-a];
4098 }
4099 else
4100 {
4101 /******************************************************************************************/
4102 /* - - */
4103 /* (a)(1 - y) | | (a + bx + c) */
4104 /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */
4105 /* (a - c) | | (a - c) */
4106 /* - - */
4107 /* */
4108 /* */
4109 /* */
4110 /* */
4111 /* */
4112 /* with F(0) = 1 */
4113 /* b */
4114 /* and F(-1) = 1 - (---) y */
4115 /* c */
4116 /* */
4117 /* */
4118 /* */
4119 /* Use a -> a' + 1 */
4120 /* _ _ */
4121 /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */
4122 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */
4123 /* (a' + 1 - c) | | (a' + 1 - c) */
4124 /* - - */
4125 /* */
4126 /* For the first F() in the solution of the radial integral */
4127 /* */
4128 /* a = ( -n + l + 1 ) */
4129 /* */
4130 /* a = -n + l + 1 */
4131 /* max(a) = max(-n) + max(l) + 1 */
4132 /* = -n + max(n-1) + 1 */
4133 /* = -n + n-1 + 1 */
4134 /* = 0 */
4135 /* */
4136 /* similiarly */
4137 /* */
4138 /* min(a) = min(-n) + min(l) + 1 */
4139 /* = min(-n) + 0 + 1 */
4140 /* = -n + 1 */
4141 /* */
4142 /* a -> a' + 1 implies */
4143 /* */
4144 /* max(a') = -1 */
4145 /* min(a') = -n */
4146 /******************************************************************************************/
4147
4148 i1= a + 1;
4149 i2= a + 1 - c;
4150 d0= (double)i2;
4151 ASSERT( i2 != 0 );
4152 d1= 1. - y;
4153 d2= (double)i1 * d1;
4154 d3= d2 / d0;
4155 d4= (double)b * y;
4156 d5= d0 + d4;
4157
4158 d8= F21i( (long int)(a + 1), b, c, y, yV );
4159 d9= F21i( (long int)(a + 2), b, c, y, yV );
4160
4161 d10= d8 - d9;
4162 d11= d3 * d10;
4163 d12= d5 / d0;
4164 d13= d12 * d8;
4165 d14= d11 + d13;
4166
4167 /* Store the result for later use */
4168 yV[-a] = d14;
4169 return d14;
4170 }
4171}
4172
4173STATIC mx F21i_log( long int a, long int b, long int c, double y, mxq *yV )
4174{
4175 DEBUG_ENTRY( "F21i_log()" );
4176
4177 mx result_mx = {0.0,0};
4178
4179 if( yV[-a].q != 0. )
4180 {
4181 /* Return the stored result */
4182 return yV[-a].mx;
4183 }
4184 else if( a == 0 )
4185 {
4186 ASSERT( yV[-a].q == 0. );
4187
4188 result_mx.m = 1.;
4189 result_mx.x = 0;
4190
4191 ASSERT( yV[-a].mx.m == 0. );
4192 ASSERT( yV[-a].mx.x == 0 );
4193
4194 yV[-a].q = 1;
4195 yV[-a].mx = result_mx;
4196 return result_mx;
4197 }
4198 else if( a == -1 )
4199 {
4200 double d1 = (double)b;
4201 double d2 = (double)c;
4202 double d3 = y * (d1/d2);
4203
4204 ASSERT( yV[-a].q == 0. );
4205 ASSERT( c != 0 );
4206 ASSERT( y != 0. );
4207
4208 result_mx.m = 1. - d3;
4209 result_mx.x = 0;
4210
4211 while ( fabs(result_mx.m) > 1.0e+25 )
4212 {
4213 result_mx.m /= 1.0e+25;
4214 result_mx.x += 25;
4215 }
4216
4217 ASSERT( yV[-a].mx.m == 0. );
4218 ASSERT( yV[-a].mx.x == 0 );
4219
4220 yV[-a].q = 1;
4221 yV[-a].mx = result_mx;
4222 return result_mx;
4223 }
4224 else
4225 {
4226 /******************************************************************************************/
4227 /* - - */
4228 /* (a)(1 - y) | | (a + bx + c) */
4229 /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */
4230 /* (a - c) | | (a - c) */
4231 /* - - */
4232 /* */
4233 /* */
4234 /* with F(0) = 1 */
4235 /* b */
4236 /* and F(-1) = 1 - (---) y */
4237 /* c */
4238 /* */
4239 /* */
4240 /* */
4241 /* Use a -> a' + 1 */
4242 /* _ _ */
4243 /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */
4244 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */
4245 /* (a' + 1 - c) | | (a' + 1 - c) */
4246 /* - - */
4247 /* */
4248 /* For the first F() in the solution of the radial integral */
4249 /* */
4250 /* a = ( -n + l + 1 ) */
4251 /* */
4252 /* a = -n + l + 1 */
4253 /* max(a) = max(-n) + max(l) + 1 */
4254 /* = -n + max(n-1) + 1 */
4255 /* = -n + n-1 + 1 */
4256 /* = 0 */
4257 /* */
4258 /* similiarly */
4259 /* */
4260 /* min(a) = min(-n) + min(l) + 1 */
4261 /* = min(-n) + 0 + 1 */
4262 /* = -n + 1 */
4263 /* */
4264 /* a -> a' + 1 implies */
4265 /* */
4266 /* max(a') = -1 */
4267 /* min(a') = -n */
4268 /******************************************************************************************/
4269
4270 mx d8 = {0.0,0}, d9 = {0.0,0}, d10 = {0.0,0}, d11 = {0.0,0};
4271
4272 double db = (double)b;
4273 double d00 = (double)(a + 1 - c);
4274 double d0 = (double)(a + 1);
4275 double d1 = 1. - y;
4276 double d2 = d0 * d1;
4277 double d3 = d2 / d00;
4278 double d4 = db * y;
4279 double d5 = d00 + d4;
4280 double d6 = d5 / d00;
4281
4282 ASSERT( yV[-a].q == 0. );
4283
4284 /******************************************************************************************/
4285 /* _ _ */
4286 /* (a' + 1) (1 - x) | | (a' + 1 - c) + b*x */
4287 /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ------------------ F(a'+1) */
4288 /* (a' + 1 - c) | | (a' + 1 - c) */
4289 /* - - */
4290 /******************************************************************************************/
4291
4292 d8= F21i_log( (a + 1), b, c, y, yV );
4293 d9= F21i_log( (a + 2), b, c, y, yV );
4294
4295 if( (d8.m) != 0. )
4296 {
4297 d10.x = d8.x;
4298 d10.m = d8.m * (1. - (d9.m/d8.m) * powi( 10., (d9.x - d8.x)));
4299 }
4300 else
4301 {
4302 d10.m = -d9.m;
4303 d10.x = d9.x;
4304 }
4305
4306 d10.m *= d3;
4307
4308 d11.x = d8.x;
4309 d11.m = d6 * d8.m;
4310
4311 if( (d11.m) != 0. )
4312 {
4313 result_mx.x = d11.x;
4314 result_mx.m = d11.m * (1. + (d10.m/d11.m) * powi( 10. , (d10.x - d11.x) ) );
4315 }
4316 else
4317 {
4318 result_mx = d10;
4319 }
4320
4321 while ( fabs(result_mx.m) > 1.0e+25 )
4322 {
4323 result_mx.m /= 1.0e+25;
4324 result_mx.x += 25;
4325 }
4326
4327 /* Store the result for later use */
4328 yV[-a].q = 1;
4329 yV[-a].mx = result_mx;
4330 return result_mx;
4331 }
4332}
4333
4334/********************************************************************************/
4335
4336inline void normalize_mx( mx& target )
4337{
4338 while( fabs(target.m) > 1.0e+25 )
4339 {
4340 target.m /= 1.0e+25;
4341 target.x += 25;
4342 }
4343 while( fabs(target.m) < 1.0e-25 )
4344 {
4345 target.m *= 1.0e+25;
4346 target.x -= 25;
4347 }
4348 return;
4349}
4350
4351inline mx add_mx( const mx& a, const mx& b )
4352{
4353 mx result = {0.0,0};
4354
4355 if( a.m != 0. )
4356 {
4357 result.x = a.x;
4358 result.m = a.m * (1. + (b.m/a.m) * powi( 10. , (b.x - a.x) ) );
4359 }
4360 else
4361 {
4362 result = b;
4363 }
4364 normalize_mx( result );
4365 return result;
4366}
4367
4368inline mx sub_mx( const mx& a, const mx& b )
4369{
4370 mx result = {0.0,0};
4371 mx minusb = b;
4372 minusb.m = -minusb.m;
4373
4374 result = add_mx( a, minusb );
4375 normalize_mx( result );
4376
4377 return result;
4378}
4379
4380inline mx mxify( double a )
4381{
4382 mx result_mx = {0.0,0};
4383
4384 result_mx.x = 0;
4385 result_mx.m = a;
4386 normalize_mx( result_mx );
4387
4388 return result_mx;
4389}
4390
4391inline double unmxify( const mx& a_mx )
4392{
4393 return a_mx.m * powi( 10., a_mx.x );
4394}
4395
4396inline mx mxify_log10( double log10_a )
4397{
4398 mx result_mx = {0.0,0};
4399
4400 while ( log10_a > 25.0 )
4401 {
4402 log10_a -= 25.0;
4403 result_mx.x += 25;
4404 }
4405
4406 while ( log10_a < -25.0 )
4407 {
4408 log10_a += 25.0;
4409 result_mx.x -= 25;
4410 }
4411
4412 result_mx.m = pow(10., log10_a);
4413
4414 return result_mx;
4415}
4416
4417inline mx mult_mx( const mx& a, const mx& b )
4418{
4419 mx result = {0.0,0};
4420
4421 result.m = (a.m * b.m);
4422 result.x = (a.x + b.x);
4423 normalize_mx( result );
4424
4425 return result;
4426}
4427
4428inline double log10_prodxx( long int lp, double Ksqrd )
4429{
4430 /**********************************************************************/
4431 /* | s=l' | s=l' */
4432 /* | ----- | --- */
4433 /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */
4434 /* | | | | --- */
4435 /* | s=0 | s=0 */
4436 /**********************************************************************/
4437
4438 if( lp == 0 )
4439 return 0.;
4440
4441 double partsum = 0.;
4442 for( long int s = 1; s <= lp; s++ )
4443 {
4444 double s2 = pow2((double)s);
4445 partsum += log10( 1. + s2*Ksqrd );
4446
4447 ASSERT( partsum >= 0. );
4448 }
4449 return partsum;
4450}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
bool is_odd(int j)
Definition cddefines.h:714
#define STATIC
Definition cddefines.h:97
T pow2(T a)
Definition cddefines.h:931
#define CALLOC
Definition cddefines.h:510
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
#define MAX2
Definition cddefines.h:782
long max(int a, long b)
Definition cddefines.h:775
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
T pow3(T a)
Definition cddefines.h:938
double powi(double, long int)
Definition service.cpp:604
double hri(long int n, long int l, long int np, long int lp, long int iz)
mx mxify_log10(double log10_a)
mx add_mx(const mx &a, const mx &b)
mx mxify(double a)
double H_photo_cs_log10(double photon_energy, long int n, long int l, long int iz)
STATIC double H_Einstein_A_lin(long int n, long int l, long int np, long int lp, long int iz)
double OscStr_f(long int n, long int l, long int np, long int lp, long int iz)
STATIC mx F21i_log(long int a, long int b, long int c, double y, mxq *yV)
static const double PHYSICAL_CONSTANT_TWO
STATIC double fsff(long int n, long int l, long int np)
double local_product(double K, long int lp)
STATIC mx bhGm_mx(long int q, double K, long int n, long int l, long int lp, mxq *rcsvV_mxq, const mx &GK_mx)
double unmxify(const mx &a_mx)
struct t_mxq mxq
double hv(long int n, long int nprime, long int iz)
void normalize_mx(mx &target)
STATIC double log10_fsff(long int n, long int l, long int np)
STATIC double hrii(long int n, long int l, long int np, long int lp)
double OscStr_f_log10(long int n, long int l, long int np, long int lp, long int iz)
STATIC double bh(double k, long int n, long int l, double *rcsvV)
double H_Einstein_A_log10(long int n, long int l, long int np, long int lp, long int iz)
STATIC double bhintegrand_log(double k, long int n, long int l, long int lp, mxq *rcsvV_mxq)
mx sub_mx(const mx &a, const mx &b)
STATIC double bhg(double K, long int n, long int l, long int lp, double *rcsvV)
STATIC mx F21_mx(long int a, long int b, long int c, double y, char A)
double hri_log10(long int n, long int l, long int np, long int lp, long int iz)
STATIC double bhGp(long int q, double K, long int n, long int l, long int lp, double *rcsvV, double GK)
STATIC mx bhG_mx(double K, long int n, long int l, long int lp, mxq *rcsvV_mxq)
mx mult_mx(const mx &a, const mx &b)
STATIC double bh_log(double k, long int n, long int l, mxq *rcsvV_mxq)
double H_Einstein_A(long int n, long int l, long int np, long int lp, long int iz)
double log10_prodxx(long int lp, double Ksqrd)
STATIC double bhg_log(double K, long int n, long int l, long int lp, mxq *rcsvV_mxq)
double H_photo_cs(double rel_photon_energy, long int n, long int l, long int iz)
static const double CONST_ONE
STATIC double bhG(double K, long int n, long int l, long int lp, double *rcsvV)
STATIC double F21i(long int a, long int b, long int c, double y, double *yV)
struct t_mx mx
STATIC double bhintegrand(double k, long int n, long int l, long int lp, double *rcsvV)
STATIC double F21(long int a, long int b, long int c, double y, char A)
STATIC double H_photo_cs_lin(double rel_photon_energy, long int n, long int l, long int iz)
STATIC mx bhGp_mx(long int q, double K, long int n, long int l, long int lp, mxq *rcsvV_mxq, const mx &GK_mx)
STATIC double bhGm(long int q, double K, long int n, long int l, long int lp, double *rcsvV, double GK)
STATIC double hrii_log(long int n, long int l, long int np, long int lp)
UNUSED const double PI
Definition physconst.h:29
UNUSED const double BOHR_RADIUS_CM
Definition physconst.h:222
UNUSED const double LOG10_E
Definition physconst.h:56
UNUSED const double HPLANCK
Definition physconst.h:103
UNUSED const double FINE_STRUCTURE
Definition physconst.h:216
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double SQRTPIBY2
Definition physconst.h:47
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double PROTON_MASS
Definition physconst.h:94
double m
long int x
struct t_mx mx
long int q
double factorial(long n)
double lfactorial(long n)
static const int NPRE_FACTORIAL
Definition thirdparty.h:24