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