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