cloudy  trunk
 All Data Structures 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-2008 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 
337  long int a,
338  long int b,
339  long int c,
340  double y,
341  char A
342 );
343 
345  long int a,
346  long int b,
347  long int c,
348  double y,
349  mxq *yV
350 );
351 
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 
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 
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 
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 
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 */
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 /*****************************************************************************/
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 
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)(l > lp ? l : lp);
845 
846  double n2 = (double)(n*n);
847 
848  double Ksqrd = (K*K);
849 
850  /**********************************************/
851  /* */
852  /* l> */
853  /* ------ Theta(nl,Kl') */
854  /* 2l+2 */
855  /* */
856  /* */
857  /* Theta(nl,Kl') = */
858  /* (1+n^2K^2) * | g(nl,Kl')|^2 */
859  /* */
860  /**********************************************/
861 
862  double d2 = 1. + n2*Ksqrd;
863  double d5 = bhg( K, n, l, lp, rcsvV );
864  double Theta = d2 * d5 * d5;
865  double d7 = (lg/Two_L_Plus_One) * Theta;
866 
867  ASSERT( Two_L_Plus_One != 0. );
868  ASSERT( Theta != 0. );
869  ASSERT( Ksqrd != 0. );
870  ASSERT( d2 != 0. );
871  ASSERT( d5 != 0. );
872  ASSERT( d7 != 0. );
873  ASSERT( lp >= 0 );
874  ASSERT( lg != 0. );
875  ASSERT( n2 != 0. );
876  ASSERT( n > 0 );
877  ASSERT( l >= 0 );
878  ASSERT( K != 0. );
879  return d7;
880 }
881 
882 /************************************************************************************************/
883 /* The photoionization cross-section (see also Burgess 1965) */
884 /* is given by; */
885 /* _ _ l'=l+1 */
886 /* |4 PI alpha a_o^2 | n^2 --- max(l,l') */
887 /* a(Z';n,l,k)=|---------------- | --- > ---------- Theta(n,l; K, l') */
888 /* | 3 | Z^2 --- (2l + 1) */
889 /* _ _ l'=l-1 */
890 /* */
891 /* */
892 /* where Theta(n,l; K, l') is defined */
893 /* */
894 /* | |^2 */
895 /* Theta(n,l; K, l') = (1 + n^2 K^2) | g(n,l; K,l') | */
896 /* | | */
897 /* */
898 /* */
899 /* ---- Not sure if this is K or k */
900 /* OO / but think it is K */
901 /* where - v */
902 /* K^2 | */
903 /* g(n,l; K,l') = ----- | R_nl(r) F_(K,l) r dr */
904 /* n^2 | */
905 /* - */
906 /* 0 */
907 /************************************************************************************************/
908 
910  double K, /* K = k / Z ::: k^2 = ejected-electron-energy (Ryd) */
911  long int n,
912  long int l,
913  long int lp,
914  /* Temporary storage for intermediate */
915  /* results of the recursive routine */
916  mxq *rcsvV_mxq
917 )
918 {
919  DEBUG_ENTRY( "bhintegrand_log()" );
920 
921  double d2 = 0.;
922  double d5 = 0.;
923  double d7 = 0.;
924  double Theta = 0.;
925  double n2 = (double)(n*n);
926  double Ksqrd = (K*K);
927  double Two_L_Plus_One = (double)(2*l + 1);
928  double lg = (double)(l > lp ? 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********************************/
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 }
1281 
1282 /************************************************************************************************/
1283 /* *** bhGp.c *** */
1284 /* */
1285 /* Here we calculate G(n,l; K,l') with the recursive formula */
1286 /* equation: (3.24) */
1287 /* */
1288 /* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */
1289 /* */
1290 /* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l+1; K,l) */
1291 /* */
1292 /* Under the transformation l -> l + 1 this gives */
1293 /* */
1294 /* 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) */
1295 /* */
1296 /* - 4n^2 (n^2-((l+1)+1)^2)[ 1+((l+1)K)^2 ] G(n,l+1+1; K,l+1) */
1297 /* */
1298 /* or */
1299 /* */
1300 /* 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) */
1301 /* */
1302 /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1303 /* */
1304 /* from the reference */
1305 /* M. Brocklehurst */
1306 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
1307 /* */
1308 /* */
1309 /* * This is valid for the case l=l'+1 * */
1310 /* * CASE: l = l'+1 * */
1311 /* * Here the p in bhGp() refers * */
1312 /* * to the Plus sign(+) in l=l'+1 * */
1313 /************************************************************************************************/
1314 
1315 STATIC double bhGp(
1316  long int q,
1317  double K,
1318  long int n,
1319  long int l,
1320  long int lp,
1321  /* Temporary storage for intermediate */
1322  /* results of the recursive routine */
1323  double *rcsvV,
1324  double GK
1325 )
1326 {
1327  DEBUG_ENTRY( "bhGp()" );
1328 
1329  /* static long int rcsv_Level = 1;
1330  printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */
1331 
1332  ASSERT( l == lp + 1 );
1333 
1334  long int rindx = 2*q;
1335 
1336  if( rcsvV[rindx] == 0. )
1337  {
1338  /* SPECIAL CASE: n = l+1 = l'+2 */
1339  if( q == n - 1 )
1340  {
1341  double Ksqrd = K * K;
1342  double n2 = (double)(n*n);
1343 
1344  double dd1 = (double)(2 * n);
1345  double dd2 = ( 1. + ( n2 * Ksqrd));
1346 
1347  /* (1+(n K)^2) */
1348  /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
1349  /* 2n */
1350  double G1 = ((dd2 * GK) / dd1);
1351 
1352  ASSERT( l == lp + 1 );
1353  ASSERT( Ksqrd != 0. );
1354  ASSERT( dd1 != 0. );
1355  ASSERT( dd2 != 0. );
1356  ASSERT( G1 != 0. );
1357 
1358  rcsvV[rindx] = G1;
1359  return G1;
1360  }
1361  /* SPECIAL CASE: n = l+2 = l'+3 */
1362  else if( q == (n - 2) )
1363  {
1364  double Ksqrd = (K*K);
1365  double n2 = (double)(n*n);
1366 
1367  /* */
1368  /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */
1369  /* */
1370  double dd1 = (double) (2 * n);
1371  double dd2 = ( 1. + ( n2 * Ksqrd));
1372  double G1 = ((dd2 * GK) / dd1);
1373 
1374  /* */
1375  /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2) */
1376  /* */
1377  double dd3 = (double)((2 * n) - 1);
1378  double dd4 = (double)(n - 1);
1379  double dd5 = (4. + (dd4 * dd2));
1380  double G2 = (dd3 * dd5 * G1);
1381 
1382  ASSERT( l == lp + 1 );
1383  ASSERT( Ksqrd != 0. );
1384  ASSERT( n2 != 0. );
1385  ASSERT( dd1 != 0. );
1386  ASSERT( dd2 != 0. );
1387  ASSERT( dd3 != 0. );
1388  ASSERT( dd4 != 0. );
1389  ASSERT( dd5 != 0. );
1390  ASSERT( G1 != 0. );
1391  ASSERT( G2 != 0. );
1392 
1393  rcsvV[rindx] = G2;
1394  return G2;
1395  }
1396  /* The GENERAL CASE n > l + 2 */
1397  else
1398  {
1399  /******************************************************************************/
1400  /* 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) */
1401  /* */
1402  /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1403  /* */
1404  /* FROM Eq. 3.24 */
1405  /* */
1406  /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */
1407  /* */
1408  /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */
1409  /******************************************************************************/
1410 
1411  long int lp1 = (q + 1);
1412  long int lp2 = (q + 2);
1413 
1414  double Ksqrd = (K*K);
1415  double n2 = (double)(n * n);
1416  double lp1s = (double)(lp1 * lp1);
1417  double lp2s = (double)(lp2 * lp2);
1418 
1419  double d1 = (4. * n2);
1420  double d2 = (4. * lp1s);
1421  double d3 = (double)((lp1)*((2 * q) + 3));
1422  double d4 = (1. + (n2 * Ksqrd));
1423  double d5 = (d1 - d2 + (d3 * d4));
1424  double d5_1 = d5 * bhGp( (q+1), K, n, l, lp, rcsvV, GK );
1425 
1426 
1427  /* 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) */
1428  /* */
1429  /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1430 
1431  double d6 = (n2 - lp2s);
1432  double d7 = (1. + (lp1s * Ksqrd));
1433  double d8 = (d1 * d6 * d7);
1434  double d8_1 = d8 * bhGp( (q+2), K, n, l, lp, rcsvV, GK );
1435  double d9 = (d5_1 - d8_1);
1436 
1437  ASSERT( l == lp + 1 );
1438  ASSERT( Ksqrd != 0. );
1439  ASSERT( n2 != 0. );
1440 
1441  ASSERT( lp1s != 0. );
1442  ASSERT( lp2s != 0. );
1443 
1444  ASSERT( d1 != 0. );
1445  ASSERT( d2 != 0. );
1446  ASSERT( d3 != 0. );
1447  ASSERT( d4 != 0. );
1448  ASSERT( d5 != 0. );
1449  ASSERT( d6 != 0. );
1450  ASSERT( d7 != 0. );
1451  ASSERT( d8 != 0. );
1452  ASSERT( d9 != 0. );
1453 
1454  rcsvV[rindx] = d9;
1455  return d9;
1456  }
1457  }
1458  else
1459  {
1460  ASSERT( rcsvV[rindx] != 0. );
1461  return rcsvV[rindx];
1462  }
1463 }
1464 
1465 /***********************log version*******************************/
1467  long int q,
1468  double K,
1469  long int n,
1470  long int l,
1471  long int lp,
1472  /* Temporary storage for intermediate */
1473  /* results of the recursive routine */
1474  mxq *rcsvV_mxq,
1475  const mx& GK_mx
1476 )
1477 {
1478  DEBUG_ENTRY( "bhGp_mx()" );
1479 
1480  /* static long int rcsv_Level = 1; */
1481  /* printf( "bhGp(): recursion level:\t%li\n",rcsv_Level++ ); */
1482 
1483  ASSERT( l == lp + 1 );
1484 
1485  long int rindx = 2*q;
1486 
1487  if( rcsvV_mxq[rindx].q == 0 )
1488  {
1489  /* SPECIAL CASE: n = l+1 = l'+2 */
1490  if( q == n - 1 )
1491  {
1492  /******************************************************/
1493  /* (1+(n K)^2) */
1494  /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
1495  /* 2n */
1496  /******************************************************/
1497 
1498  double Ksqrd = (K * K);
1499  double n2 = (double)(n*n);
1500 
1501  double dd1 = (double) (2 * n);
1502  double dd2 = ( 1. + ( n2 * Ksqrd));
1503  double dd3 = dd2/dd1;
1504 
1505  mx dd3_mx = mxify( dd3 );
1506  mx G1_mx = mult_mx( dd3_mx, GK_mx);
1507 
1508  normalize_mx( G1_mx );
1509 
1510  ASSERT( l == lp + 1 );
1511  ASSERT( Ksqrd != 0. );
1512  ASSERT( n2 != 0. );
1513  ASSERT( dd1 != 0. );
1514  ASSERT( dd2 != 0. );
1515 
1516  rcsvV_mxq[rindx].q = 1;
1517  rcsvV_mxq[rindx].mx = G1_mx;
1518  return G1_mx;
1519  }
1520  /* SPECIAL CASE: n = l+2 = l'+3 */
1521  else if( q == (n - 2) )
1522  {
1523  /****************************************************************/
1524  /* */
1525  /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/
1526  /* */
1527  /****************************************************************/
1528  /* (1+(n K)^2) */
1529  /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
1530  /* 2n */
1531  /****************************************************************/
1532 
1533  double Ksqrd = (K*K);
1534  double n2 = (double)(n*n);
1535  double dd1 = (double)(2 * n);
1536  double dd2 = ( 1. + ( n2 * Ksqrd) );
1537  double dd3 = (dd2/dd1);
1538  double dd4 = (double) ((2 * n) - 1);
1539  double dd5 = (double) (n - 1);
1540  double dd6 = (4. + (dd5 * dd2));
1541  double dd7 = dd4 * dd6;
1542 
1543  /****************************************************************/
1544  /* */
1545  /* G(n,n-2; K,n-3) = (2n-1) (4+(n-1)(1+(n K)^2)) G(n,n-1; K,n-2)*/
1546  /* */
1547  /****************************************************************/
1548 
1549  mx dd3_mx = mxify( dd3 );
1550  mx dd7_mx = mxify( dd7 );
1551  mx G1_mx = mult_mx( dd3_mx, GK_mx );
1552  mx G2_mx = mult_mx( dd7_mx, G1_mx );
1553 
1554  normalize_mx( G2_mx );
1555 
1556  ASSERT( l == lp + 1 );
1557  ASSERT( Ksqrd != 0. );
1558  ASSERT( n2 != 0. );
1559  ASSERT( dd1 != 0. );
1560  ASSERT( dd2 != 0. );
1561  ASSERT( dd3 != 0. );
1562  ASSERT( dd4 != 0. );
1563  ASSERT( dd5 != 0. );
1564  ASSERT( dd6 != 0. );
1565  ASSERT( dd7 != 0. );
1566 
1567  rcsvV_mxq[rindx].q = 1;
1568  rcsvV_mxq[rindx].mx = G2_mx;
1569  return G2_mx;
1570  }
1571  /* The GENERAL CASE n > l + 2*/
1572  else
1573  {
1574  /**************************************************************************************/
1575  /* 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) */
1576  /* */
1577  /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1578  /* */
1579  /* FROM Eq. 3.24 */
1580  /* */
1581  /* G(n,l-1; K,l-2) = [ 4n^2-4l+^2 + l(2l+1)(1+(n K)^2) ] G(n,l; K,l-1) */
1582  /* */
1583  /* - 4n^2 (n^2-(l+1)^2)[ 1+((lK)^2 ] G(n,l+1; K,l) */
1584  /**************************************************************************************/
1585 
1586  long int lp1 = (q + 1);
1587  long int lp2 = (q + 2);
1588 
1589  double Ksqrd = (K * K);
1590  double n2 = (double)(n * n);
1591  double lp1s = (double)(lp1 * lp1);
1592  double lp2s = (double)(lp2 * lp2);
1593 
1594  double d1 = (4. * n2);
1595  double d2 = (4. * lp1s);
1596  double d3 = (double)((lp1)*((2 * q) + 3));
1597  double d4 = (1. + (n2 * Ksqrd));
1598  /* [ 4n^2 - 4(l+1)^2 + (l+1)(2l+3)(1+(n K)^2) ] */
1599  double d5 = d1 - d2 + (d3 * d4);
1600 
1601  /* (n^2-(l+2)^2) */
1602  double d6 = (n2 - lp2s);
1603 
1604  /* [ 1+((l+1)K)^2 ] */
1605  double d7 = (1. + (lp1s * Ksqrd));
1606 
1607  /* { 4n^2 (n^2-(l+1)^2)[ 1+((l+1) K)^2 ] } */
1608  double d8 = (d1 * d6 * d7);
1609 
1610  /**************************************************************************************/
1611  /* 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) */
1612  /* */
1613  /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+1)K)^2 ] G(n,l+2; K,l+1) */
1614  /**************************************************************************************/
1615 
1616  mx d5_mx=mxify( d5 );
1617  mx d8_mx=mxify( d8 );
1618 
1619  mx t0_mx = bhGp_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1620  mx t1_mx = bhGp_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1621 
1622  mx d9_mx = mult_mx( d5_mx, t0_mx );
1623  mx d10_mx = mult_mx( d8_mx, t1_mx );
1624 
1625  mx result_mx = sub_mx( d9_mx, d10_mx );
1626  normalize_mx( result_mx );
1627 
1628  ASSERT( d1 != 0. );
1629  ASSERT( d2 != 0. );
1630  ASSERT( d3 != 0. );
1631  ASSERT( d4 != 0. );
1632  ASSERT( d5 != 0. );
1633  ASSERT( d6 != 0. );
1634  ASSERT( d7 != 0. );
1635  ASSERT( d8 != 0. );
1636 
1637  ASSERT( l == lp + 1 );
1638  ASSERT( Ksqrd != 0. );
1639  ASSERT( n2 != 0. );
1640  ASSERT( lp1s != 0. );
1641  ASSERT( lp2s != 0. );
1642 
1643  rcsvV_mxq[rindx].q = 1;
1644  rcsvV_mxq[rindx].mx = result_mx;
1645  return result_mx;
1646  }
1647  }
1648  else
1649  {
1650  ASSERT( rcsvV_mxq[rindx].q != 0 );
1651  rcsvV_mxq[rindx].q = 1;
1652  return rcsvV_mxq[rindx].mx;
1653  }
1654 }
1655 
1656 /************************************************************************************************/
1657 /* *** bhGm.c *** */
1658 /* */
1659 /* Here we calculate G(n,l; K,l') with the recursive formula */
1660 /* equation: (3.23) */
1661 /* */
1662 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */
1663 /* */
1664 /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */
1665 /* */
1666 /* Under the transformation l -> l + 2 this gives */
1667 /* */
1668 /* 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) */
1669 /* */
1670 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+2+1)^2 K^2 ] G(n,l+2; K,l+2+1) */
1671 /* */
1672 /* */
1673 /* or */
1674 /* */
1675 /* 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) */
1676 /* */
1677 /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */
1678 /* */
1679 /* */
1680 /* from the reference */
1681 /* M. Brocklehurst */
1682 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
1683 /* */
1684 /* */
1685 /* * This is valid for the case l=l'-1 * */
1686 /* * CASE: l = l'-1 * */
1687 /* * Here the p in bhGm() refers * */
1688 /* * to the Minus sign(-) in l=l'-1 * */
1689 /************************************************************************************************/
1690 
1691 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1692 #pragma optimize("", off)
1693 #endif
1694 STATIC double bhGm(
1695  long int q,
1696  double K,
1697  long int n,
1698  long int l,
1699  long int lp,
1700  double *rcsvV,
1701  double GK
1702 )
1703 {
1704  DEBUG_ENTRY( "bhGm()" );
1705 
1706  ASSERT( l == lp - 1 );
1707  ASSERT( l < n );
1708 
1709  long int rindx = 2*q+1;
1710 
1711  if( rcsvV[rindx] == 0. )
1712  {
1713  /* CASE: l = n - 1 */
1714  if( q == n - 1 )
1715  {
1716  ASSERT( l == lp - 1 );
1717  rcsvV[rindx] = GK;
1718  return GK;
1719  }
1720  /* CASE: l = n - 2 */
1721  else if( q == n - 2 )
1722  {
1723  double dd1 = 0.;
1724  double dd2 = 0.;
1725 
1726  double G2 = 0.;
1727 
1728  double Ksqrd = 0.;
1729  double n1 = 0.;
1730  double n2 = 0.;
1731 
1732  ASSERT(l == lp - 1);
1733 
1734  /* K^2 */
1735  Ksqrd = K * K;
1736  ASSERT( Ksqrd != 0. );
1737 
1738  /* n */
1739  n1 = (double)n;
1740  ASSERT( n1 != 0. );
1741 
1742  /* n^2 */
1743  n2 = (double)(n*n);
1744  ASSERT( n2 != 0. );
1745 
1746  /* equation: (3.20) */
1747  /* G(n,n-2; K,n-1) = */
1748  /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */
1749  dd1 = (double) ((2 * n) - 1);
1750  ASSERT( dd1 != 0. );
1751 
1752  dd2 = (1. + (n2 * Ksqrd));
1753  ASSERT( dd2 != 0. );
1754 
1755  G2 = dd1 * dd2 * n1 * GK;
1756  ASSERT( G2 != 0. );
1757 
1758  rcsvV[rindx] = G2;
1759  ASSERT( G2 != 0. );
1760  return G2;
1761  }
1762  else
1763  {
1764  long int lp2 = (q + 2);
1765  long int lp3 = (q + 3);
1766 
1767  double lp2s = (double)(lp2 * lp2);
1768  double lp3s = (double)(lp3 * lp3);
1769 
1770  /******************************************************************************/
1771  /* 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) */
1772  /* */
1773  /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */
1774  /* */
1775  /* */
1776  /* FROM Eq. 3.23 */
1777  /* */
1778  /* 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) */
1779  /* */
1780  /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */
1781  /******************************************************************************/
1782 
1783  double Ksqrd = (K*K);
1784  double n2 = (double)(n*n);
1785  double d1 = (4. * n2);
1786  double d2 = (4. * lp2s);
1787  double d3 = (double)(lp2)*((2*q)+3);
1788  double d4 = (1. + (n2 * Ksqrd));
1789  /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */
1790  double d5 = d1 - d2 + (d3 * d4);
1791 
1792  /******************************************************************************/
1793  /* 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) */
1794  /* */
1795  /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */
1796  /******************************************************************************/
1797 
1798  double d6 = (n2 - lp2s);
1799  /* [ 1+((l+3)K)^2 ] */
1800  double d7 = (1. + (lp3s * Ksqrd));
1801  /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */
1802  double d8 = d1 * d6 * d7;
1803  double d9 = d5 * bhGm( (q+1), K, n, l, lp, rcsvV, GK );
1804  double d10 = d8 * bhGm( (q+2), K, n, l, lp, rcsvV, GK );
1805  double d11 = d9 - d10;
1806 
1807  ASSERT( l == lp - 1 );
1808  ASSERT( lp2s != 0. );
1809  ASSERT( Ksqrd != 0. );
1810  ASSERT( n2 != 0. );
1811  ASSERT( d1 != 0. );
1812  ASSERT( d2 != 0. );
1813  ASSERT( d3 != 0. );
1814  ASSERT( d4 != 0. );
1815  ASSERT( d5 != 0. );
1816  ASSERT( d6 != 0. );
1817  ASSERT( d7 != 0. );
1818  ASSERT( d8 != 0. );
1819  ASSERT( d9 != 0. );
1820  ASSERT( d10 != 0. );
1821  ASSERT( lp3s != 0. );
1822 
1823  rcsvV[rindx] = d11;
1824  return d11;
1825  }
1826  }
1827  else
1828  {
1829  ASSERT( rcsvV[rindx] != 0. );
1830  return rcsvV[rindx];
1831  }
1832 }
1833 #if defined (__ICC) && defined(__amd64) && __INTEL_COMPILER < 1000
1834 #pragma optimize("", on)
1835 #endif
1836 
1837 /************************log version***********************************/
1839  long int q,
1840  double K,
1841  long int n,
1842  long int l,
1843  long int lp,
1844  mxq *rcsvV_mxq,
1845  const mx& GK_mx
1846 )
1847 {
1848  DEBUG_ENTRY( "bhGm_mx()" );
1849 
1850  /*static long int rcsv_Level = 1; */
1851  /*printf( "bhGm(): recursion level:\t%li\n",rcsv_Level++ ); */
1852 
1853  ASSERT( l == lp - 1 );
1854  ASSERT( l < n );
1855 
1856  long int rindx = 2*q+1;
1857 
1858  if( rcsvV_mxq[rindx].q == 0 )
1859  {
1860  /* CASE: l = n - 1 */
1861  if( q == n - 1 )
1862  {
1863  mx result_mx = GK_mx;
1864  normalize_mx( result_mx );
1865 
1866  rcsvV_mxq[rindx].q = 1;
1867  rcsvV_mxq[rindx].mx = result_mx;
1868 
1869  ASSERT(l == lp - 1);
1870  return result_mx;
1871  }
1872  /* CASE: l = n - 2 */
1873  else if( q == n - 2 )
1874  {
1875  double Ksqrd = (K * K);
1876  double n1 = (double)n;
1877  double n2 = (double) (n*n);
1878  double dd1 = (double)((2 * n) - 1);
1879  double dd2 = (1. + (n2 * Ksqrd));
1880  /*(2n-1)(1+(n K)^2) n*/
1881  double dd3 = (dd1*dd2*n1);
1882 
1883  /******************************************************/
1884  /* G(n,n-2; K,n-1) = */
1885  /* (2n-1)(1+(n K)^2) n G(n,n-1; K,n) */
1886  /******************************************************/
1887 
1888  mx dd3_mx = mxify( dd3 );
1889  mx G2_mx = mult_mx( dd3_mx, GK_mx );
1890 
1891  normalize_mx( G2_mx );
1892 
1893  ASSERT( l == lp - 1);
1894  ASSERT( n1 != 0. );
1895  ASSERT( n2 != 0. );
1896  ASSERT( dd1 != 0. );
1897  ASSERT( dd2 != 0. );
1898  ASSERT( dd3 != 0. );
1899  ASSERT( Ksqrd != 0. );
1900 
1901  rcsvV_mxq[rindx].q = 1;
1902  rcsvV_mxq[rindx].mx = G2_mx;
1903  return G2_mx;
1904  }
1905  else
1906  {
1907  /******************************************************************************/
1908  /* 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) */
1909  /* */
1910  /* - 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] G(n,l+2; K,l+3) */
1911  /* */
1912  /* */
1913  /* FROM Eq. 3.23 */
1914  /* */
1915  /* 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) */
1916  /* */
1917  /* - 4n^2 (n^2-l^2)[ 1 + (l+1)^2 K^2 ] G(n,l; K,l+1) */
1918  /******************************************************************************/
1919 
1920  long int lp2 = (q + 2);
1921  long int lp3 = (q + 3);
1922 
1923  double lp2s = (double)(lp2 * lp2);
1924  double lp3s = (double)(lp3 * lp3);
1925  double n2 = (double)(n*n);
1926  double Ksqrd = (K * K);
1927 
1928  /******************************************************/
1929  /* [ 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) ] */
1930  /******************************************************/
1931 
1932  double d1 = (4. * n2);
1933  double d2 = (4. * lp2s);
1934  double d3 = (double)(lp2)*((2*q)+3);
1935  double d4 = (1. + (n2 * Ksqrd));
1936  /* 4n^2-4(l+2)^2 + (l+2)(2l+3)(1+(n K)^2) */
1937  double d5 = d1 - d2 + (d3 * d4);
1938 
1939  mx d5_mx=mxify(d5);
1940 
1941  /******************************************************/
1942  /* 4n^2 (n^2-(l+2)^2)[ 1+((l+3)^2 K^2) ] */
1943  /******************************************************/
1944 
1945  double d6 = (n2 - lp2s);
1946  double d7 = (1. + (lp3s * Ksqrd));
1947  /* 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] */
1948  double d8 = d1 * d6 * d7;
1949 
1950  mx d8_mx = mxify(d8);
1951 
1952  /******************************************************************************/
1953  /* 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) */
1954  /* */
1955  /* - 4n^2 (n^2-(l+2)^2)[ 1 + (l+3)^2 K^2 ] G(n,l+2; K,l+3) */
1956  /******************************************************************************/
1957 
1958  mx d9_mx = bhGm_mx( (q+1), K, n, l, lp, rcsvV_mxq, GK_mx );
1959  mx d10_mx = bhGm_mx( (q+2), K, n, l, lp, rcsvV_mxq, GK_mx );
1960  mx d11_mx = mult_mx( d5_mx, d9_mx );
1961  mx d12_mx = mult_mx( d8_mx, d10_mx);
1962  mx result_mx = sub_mx( d11_mx , d12_mx );
1963  rcsvV_mxq[rindx].q = 1;
1964  rcsvV_mxq[rindx].mx = result_mx;
1965 
1966  ASSERT(l == lp - 1);
1967  ASSERT(n2 != 0. );
1968  ASSERT(lp2s != 0.);
1969  ASSERT( lp3s != 0.);
1970  ASSERT(Ksqrd != 0.);
1971 
1972  ASSERT(d1 != 0.);
1973  ASSERT(d2 != 0.);
1974  ASSERT(d3 != 0.);
1975  ASSERT(d4 != 0.);
1976  ASSERT(d5 != 0.);
1977  ASSERT(d6 != 0.);
1978  ASSERT(d7 != 0.);
1979  ASSERT(d8 != 0.);
1980  return result_mx;
1981  }
1982  }
1983  else
1984  {
1985  ASSERT( rcsvV_mxq[rindx].q != 0 );
1986  return rcsvV_mxq[rindx].mx;
1987  }
1988 }
1989 
1990 /****************************************************************************************/
1991 /* */
1992 /* bhg.c */
1993 /* */
1994 /* From reference; */
1995 /* M. Brocklehurst */
1996 /* Mon. Note. R. astr. Soc. (1971) 153, 471-490 */
1997 /* */
1998 /* */
1999 /* We wish to compute the following function, */
2000 /* */
2001 /* equation: (3.17) */
2002 /* - s=l' - (1/2) */
2003 /* | (n+l)! ----- | */
2004 /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */
2005 /* | (n-l-1)! | | | */
2006 /* - s=0 - */
2007 /* */
2008 /* Using various recursion relations (for l'=l+1) */
2009 /* */
2010 /* equation: (3.23) */
2011 /* G(n,l-2; K,l-1) = [ 4n^2-4l^2+l(2l-1)(1+(n K)^2) ] G(n,l-1; K,l) */
2012 /* */
2013 /* - 4n^2 (n^2-l^2)[1+(l+1)^2 K^2] G(n,l; K,l+1) */
2014 /* */
2015 /* and (for l'=l-1) */
2016 /* */
2017 /* equation: (3.24) */
2018 /* G(n,l-1; K,l-2) = [ 4n^2-4l^2 + l(2l-1)(1+(n K)^2) ] G(n,l; K,l-1) */
2019 /* */
2020 /* - 4n^2 (n^2-(l+1)^2)[ 1+(lK)^2 ] G(n,l; K,l+1) */
2021 /* */
2022 /* */
2023 /* the starting point for the recursion relations are: */
2024 /* */
2025 /* */
2026 /* equation (3.18): */
2027 /* */
2028 /* | pi |(1/2) 8n */
2029 /* G(n,n-1; 0,n) = | -- | ------- (4n)^2 exp(-2n) */
2030 /* | 2 | (2n-1)! */
2031 /* */
2032 /* equation (3.20): */
2033 /* */
2034 /* exp(2n-2/K tan^(-1)(n K) */
2035 /* G(n,n-1; K,n) = --------------------------------------- */
2036 /* sqrt(1 - exp(-2 pi/ K)) * (1+(n K)^(n+2) */
2037 /* */
2038 /* */
2039 /* */
2040 /* equation (3.20): */
2041 /* G(n,n-2; K,n-1) = (2n-2)(1+(n K)^2) n G(n,n-1; K,n) */
2042 /* */
2043 /* */
2044 /* equation (3.21): */
2045 /* */
2046 /* (1+(n K)^2) */
2047 /* G(n,n-1; K,n-2) = ----------- G(n,n-1; K,n) */
2048 /* 2n */
2049 /****************************************************************************************/
2050 
2051 STATIC double bhg(
2052  double K,
2053  long int n,
2054  long int l,
2055  long int lp,
2056  /* Temporary storage for intermediate */
2057  /* results of the recursive routine */
2058  double *rcsvV
2059 )
2060 {
2061  DEBUG_ENTRY( "bhg()" );
2062 
2063  double ld1 = factorial( n + l );
2064  double ld2 = factorial( n - l - 1 );
2065  double ld3 = (ld1 / ld2);
2066 
2067  double partprod = local_product( K , lp );
2068 
2069  /**************************************************************************************/
2070  /* equation: (3.17) */
2071  /* - s=l' - (1/2) */
2072  /* | (n+l)! ----- | */
2073  /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */
2074  /* | (n-l-1)! | | | */
2075  /* - s=0 - */
2076  /**************************************************************************************/
2077 
2078  /**********************************************/
2079  /* - s=l' - (1/2) */
2080  /* | (n+l)! ----- | */
2081  /* | -------- | | (1 + s^2 K^2) | */
2082  /* | (n-l-1)! | | | */
2083  /* - s=0 - */
2084  /**********************************************/
2085 
2086  double d2 = sqrt( ld3 * partprod );
2087  double d3 = powi( (2 * n) , (l - n) );
2088  double d4 = bhG( K, n, l, lp, rcsvV );
2089  double d5 = (d2 * d3);
2090  double d6 = (d5 * d4);
2091 
2092  ASSERT(K != 0.);
2093  ASSERT( (n+l) >= 1 );
2094  ASSERT( ((n-l)-1) >= 0 );
2095 
2096  ASSERT( partprod != 0. );
2097 
2098  ASSERT( ld1 != 0. );
2099  ASSERT( ld2 != 0. );
2100  ASSERT( ld3 != 0. );
2101 
2102  ASSERT( d2 != 0. );
2103  ASSERT( d3 != 0. );
2104  ASSERT( d4 != 0. );
2105  ASSERT( d5 != 0. );
2106  ASSERT( d6 != 0. );
2107  return d6;
2108 }
2109 
2110 /********************log version**************************/
2112  double K,
2113  long int n,
2114  long int l,
2115  long int lp,
2116  /* Temporary storage for intermediate */
2117  /* results of the recursive routine */
2118  mxq *rcsvV_mxq
2119 )
2120 {
2121  /**************************************************************************************/
2122  /* equation: (3.17) */
2123  /* - s=l' - (1/2) */
2124  /* | (n+l)! ----- | */
2125  /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) G(n,l; K,l') */
2126  /* | (n-l-1)! | | | */
2127  /* - s=0 - */
2128  /**************************************************************************************/
2129 
2130  DEBUG_ENTRY( "bhg_log()" );
2131 
2132  double d1 = (double)(2*n);
2133  double d2 = (double)(l-n);
2134  double Ksqrd = (K*K);
2135 
2136  /**************************************************************************************/
2137  /* */
2138  /* | (n+l)! | */
2139  /* log10 | -------- | = log10((n+1)!) - log10((n-l-1)!) */
2140  /* | (n-l-1)! | */
2141  /* */
2142  /**************************************************************************************/
2143 
2144  double ld1 = lfactorial( n + l );
2145  double ld2 = lfactorial( n - l - 1 );
2146 
2147  /**********************************************************************/
2148  /* | s=l' | s=l' */
2149  /* | ----- | --- */
2150  /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */
2151  /* | | | | --- */
2152  /* | s=0 | s=0 */
2153  /**********************************************************************/
2154 
2155  double ld3 = log10_prodxx( lp, Ksqrd );
2156 
2157  /**********************************************/
2158  /* - s=l' - (1/2) */
2159  /* | (n+l)! ----- | */
2160  /* | -------- | | (1 + s^2 K^2) | */
2161  /* | (n-l-1)! | | | */
2162  /* - s=0 - */
2163  /**********************************************/
2164 
2165  /***********************************************************************/
2166  /* */
2167  /* | - s=l' - (1/2) | */
2168  /* | | (n+l)! ----- | | */
2169  /* log10| | -------- | | (1 + s^2 K^2) | | == */
2170  /* | | (n-l-1)! | | | | */
2171  /* | - s=0 - | */
2172  /* */
2173  /* | | s=l' | | */
2174  /* | | (n+l)! | | ----- | | */
2175  /* (1/2)* |log10 | -------- | + log10 | | | (1 + s^2 K^2) | | */
2176  /* | | (n-l-1)! | | | | | | */
2177  /* | | s=0 | | */
2178  /* */
2179  /***********************************************************************/
2180 
2181  double ld4 = (1./2.) * ( ld3 + ld1 - ld2 );
2182 
2183  /**********************************************/
2184  /* (2n)^(l-n) */
2185  /**********************************************/
2186  /* log10( 2n^(L-n) ) = (L-n) log10( 2n ) */
2187  /**********************************************/
2188 
2189  double ld5 = d2 * log10( d1 );
2190 
2191  /**************************************************************************************/
2192  /* equation: (3.17) */
2193  /* - s=l' - (1/2) */
2194  /* | (n+l)! ----- | */
2195  /* g(nl, Kl) = | -------- | | (1 + s^2 K^2) | * (2n)^(l-n) * G(n,l; K,l') */
2196  /* | (n-l-1)! | | | */
2197  /* - s=0 - */
2198  /**************************************************************************************/
2199 
2200  /****************************************************/
2201  /* */
2202  /* - s=l' - (1/2) */
2203  /* | (n+l)! ----- | */
2204  /* | -------- | | (1 + s^2 K^2) | * (2n)^(L-n) */
2205  /* | (n-l-1)! | | | */
2206  /* - s=0 - */
2207  /****************************************************/
2208 
2209  double ld6 = (ld5+ld4);
2210 
2211  mx d6_mx = mxify_log10( ld6 );
2212  mx dd1_mx = bhG_mx( K, n, l, lp, rcsvV_mxq );
2213  mx dd2_mx = mult_mx( d6_mx, dd1_mx );
2214  normalize_mx( dd2_mx );
2215  double result = unmxify( dd2_mx );
2216 
2217  ASSERT( result != 0. );
2218 
2219  ASSERT( Ksqrd != 0. );
2220  ASSERT( ld3 >= 0. );
2221 
2222  ASSERT( d1 > 0. );
2223  ASSERT( d2 < 0. );
2224  return result;
2225 }
2226 
2227 inline double local_product( double K , long int lp )
2228 {
2229  long int s = 0;
2230 
2231  double Ksqrd =(K*K);
2232  double partprod = 1.;
2233 
2234  for( s = 0; s <= lp; s = s + 1 )
2235  {
2236  double s2 = (double)(s*s);
2237 
2238  /**************************/
2239  /* s=l' */
2240  /* ----- */
2241  /* | | (1 + s^2 K^2) */
2242  /* | | */
2243  /* s=0 */
2244  /**************************/
2245 
2246  partprod *= ( 1. + ( s2 * Ksqrd ) );
2247  }
2248  return partprod;
2249 }
2250 
2251 /************************************************************************/
2252 /* Find the Einstein A's for hydrogen for a */
2253 /* transition n,l --> n',l' */
2254 /* */
2255 /* In the following, we will assume n > n' */
2256 /************************************************************************/
2257 /*******************************************************************************/
2258 /* */
2259 /* Einstein A() for the transition from the */
2260 /* initial state n,l to the finial state n',l' */
2261 /* is given by oscillator f() */
2262 /* */
2263 /* hbar w max(l,l') | | 2 */
2264 /* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */
2265 /* 3 R_oo ( 2l + 1 ) | | */
2266 /* */
2267 /* */
2268 /* E(n,l;n',l') max(l,l') | | 2 */
2269 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
2270 /* 3 R_oo ( 2l + 1 ) | | */
2271 /* */
2272 /* */
2273 /* See for example Gordan Drake's */
2274 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */
2275 /* */
2276 /* Here R_oo is the infinite mass Rydberg length */
2277 /* */
2278 /* */
2279 /* h c */
2280 /* R_oo --- = 13.605698 eV */
2281 /* {e} */
2282 /* */
2283 /* */
2284 /* R_oo = 2.179874e-11 ergs */
2285 /* */
2286 /* w = omega */
2287 /* = frequency of transition from n,l to n',l' */
2288 /* */
2289 /* */
2290 /* */
2291 /* here g_k are statistical weights obtained from */
2292 /* the appropriate angular momentum quantum numbers */
2293 /* */
2294 /* */
2295 /* - - 2 */
2296 /* 64 pi^4 (e a_o)^2 max(l,l') | | */
2297 /* A(n,l;n',l') = ------------------- ----------- v^3 | R(n,l;n',l') | */
2298 /* 3 h c^3 2*l + 1 | | */
2299 /* - - */
2300 /* */
2301 /* */
2302 /* pi 3.141592654 */
2303 /* plank_hbar 6.5821220 eV sec */
2304 /* R_oo 2.179874e-11 ergs */
2305 /* plank_h 6.6260755e-34 J sec */
2306 /* e_charge 1.60217733e-19 C */
2307 /* a_o 0.529177249e-10 m */
2308 /* vel_light_c 299792458L m sec^-1 */
2309 /* */
2310 /* */
2311 /* */
2312 /* 64 pi^4 (e a_o)^2 64 pi^4 (a_o)^2 e^2 1 1 */
2313 /* ----------------- = ----------------- -------- ---- = 7.5197711e-38 ----- */
2314 /* 3 h c^3 3 c^2 hbar c 2 pi sec */
2315 /* */
2316 /* */
2317 /* e^2 1 */
2318 /* using ---------- = alpha = ---- */
2319 /* hbar c 137 */
2320 /*******************************************************************************/
2321 
2322 double H_Einstein_A(/* IN THE FOLLOWING WE HAVE n > n' */
2323  /* principal quantum number, 1 for ground, upper level */
2324  long int n,
2325  /* angular momentum, 0 for s */
2326  long int l,
2327  /* principal quantum number, 1 for ground, lower level */
2328  long int np,
2329  /* angular momentum, 0 for s */
2330  long int lp,
2331  /* Nuclear charge, 1 for H+, 2 for He++, etc */
2332  long int iz
2333 )
2334 {
2335  DEBUG_ENTRY( "H_Einstein_A()" );
2336 
2337  double result;
2338  if( n > 60 || np > 60 )
2339  {
2340  result = H_Einstein_A_log10(n,l,np,lp,iz );
2341  }
2342  else
2343  {
2344  result = H_Einstein_A_lin(n,l,np,lp,iz );
2345  }
2346  return result;
2347 }
2348 
2349 /************************************************************************/
2350 /* Calculates the Einstein A's for hydrogen */
2351 /* for the transition n,l --> n',l' */
2352 /* units of sec^(-1) */
2353 /* */
2354 /* In the following, we have n > n' */
2355 /************************************************************************/
2356 STATIC double H_Einstein_A_lin(/* IN THE FOLLOWING WE HAVE n > n' */
2357  /* principal quantum number, 1 for ground, upper level */
2358  long int n,
2359  /* angular momentum, 0 for s */
2360  long int l,
2361  /* principal quantum number, 1 for ground, lower level */
2362  long int np,
2363  /* angular momentum, 0 for s */
2364  long int lp,
2365  /* Nuclear charge, 1 for H+, 2 for He++, etc */
2366  long int iz
2367 )
2368 {
2369  DEBUG_ENTRY( "H_Einstein_A_lin()" );
2370 
2371  /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */
2372  double d1 = hv( n, np, iz );
2373  double d2 = d1 / HPLANCK; /* v = hv / h */
2374  double d3 = pow3(d2);
2375  double lg = (double)(l > lp ? l : lp);
2376  double Two_L_Plus_One = (double)(2*l + 1);
2377  double d6 = lg / Two_L_Plus_One;
2378  double d7 = hri( n, l, np, lp , iz );
2379  double d8 = d7 * d7;
2380  double result = CONST_ONE * d3 * d6 * d8;
2381 
2382  /* validate the incoming data */
2383  if( n >= 70 )
2384  {
2385  fprintf(ioQQQ,"Principal Quantum Number `n' too large.\n");
2386  cdEXIT(EXIT_FAILURE);
2387  }
2388  if( iz < 1 )
2389  {
2390  fprintf(ioQQQ," The charge is impossible.\n");
2391  cdEXIT(EXIT_FAILURE);
2392  }
2393  if( n < 1 || np < 1 || l >= n || lp >= np )
2394  {
2395  fprintf(ioQQQ," The quantum numbers are impossible.\n");
2396  cdEXIT(EXIT_FAILURE);
2397  }
2398  if( n <= np )
2399  {
2400  fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n");
2401  cdEXIT(EXIT_FAILURE);
2402  }
2403  return result;
2404 }
2405 
2406 /**********************log version****************************/
2407 double H_Einstein_A_log10(/* returns Einstein A in units of (sec)^-1 */
2408  long int n,
2409  long int l,
2410  long int np,
2411  long int lp,
2412  long int iz
2413 )
2414 {
2415  DEBUG_ENTRY( "H_Einstein_A_log10()" );
2416 
2417  /* hv calculates photon energy in ergs for n -> n' transitions for H and H-like ions */
2418  double d1 = hv( n, np , iz );
2419  double d2 = d1 / HPLANCK; /* v = hv / h */
2420  double d3 = pow3(d2);
2421  double lg = (double)(l > lp ? l : lp);
2422  double Two_L_Plus_One = (double)(2*l + 1);
2423  double d6 = lg / Two_L_Plus_One;
2424  double d7 = hri_log10( n, l, np, lp , iz );
2425  double d8 = d7 * d7;
2426  double result = CONST_ONE * d3 * d6 * d8;
2427 
2428  /* validate the incoming data */
2429  if( iz < 1 )
2430  {
2431  fprintf(ioQQQ," The charge is impossible.\n");
2432  cdEXIT(EXIT_FAILURE);
2433  }
2434  if( n < 1 || np < 1 || l >= n || lp >= np )
2435  {
2436  fprintf(ioQQQ," The quantum numbers are impossible.\n");
2437  cdEXIT(EXIT_FAILURE);
2438  }
2439  if( n <= np )
2440  {
2441  fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n");
2442  cdEXIT(EXIT_FAILURE);
2443  }
2444  return result;
2445 }
2446 
2447 /********************************************************************************/
2448 /* hv calculates photon energy for n -> n' transitions for H and H-like ions */
2449 /* simplest case of no "l" or "m" dependence */
2450 /* epsilon_0 = 1 in vacu */
2451 /* */
2452 /* */
2453 /* R_h */
2454 /* Energy(n,Z) = - ------- */
2455 /* n^2 */
2456 /* */
2457 /* */
2458 /* */
2459 /* Friedrich -- Theoretical Atomic Physics pg. 60 eq. 2.8 */
2460 /* */
2461 /* u */
2462 /* R_h = --- R_oo where */
2463 /* m_e */
2464 /* */
2465 /* h c */
2466 /* R_oo --- = 2.179874e-11 ergs */
2467 /* e */
2468 /* */
2469 /* (Harmin Lecture Notes for course phy-651 Spring 1994) */
2470 /* where m_e (m_p) is the mass of and electron (proton) */
2471 /* and u is the reduced electron mass for neutral hydrogen */
2472 /* */
2473 /* */
2474 /* m_e m_p m_e */
2475 /* u = --------- = ----------- */
2476 /* m_e + m_p 1 + m_e/m_p */
2477 /* */
2478 /* m_e */
2479 /* Now ----- = 0.000544617013 */
2480 /* m_p */
2481 /* u */
2482 /* so that --- = 0.999455679 */
2483 /* m_e */
2484 /* */
2485 /* */
2486 /* returns energy of photon in ergs */
2487 /* */
2488 /* hv (n,n',Z) is for transitions n -> n' */
2489 /* */
2490 /* 1 erg = 1e-07 J */
2491 /********************************************************************************/
2492 /********************************************************************************/
2493 /* WARNING: hv() use the electron reduced mass for hydrogen instead of */
2494 /* the reduced mass associated with the apropriate ion */
2495 /********************************************************************************/
2496 
2497 inline double hv( long int n, long int nprime, long int iz )
2498 {
2499  DEBUG_ENTRY( "hv()" );
2500 
2501  double n1 = (double)n;
2502  double n2 = n1*n1;
2503  double np1 = (double)nprime;
2504  double np2 = np1*np1;
2505  double rmr = 1./(1. + ELECTRON_MASS/PROTON_MASS); /* 0.999455679 */
2506  double izsqrd = (double)(iz*iz);
2507 
2508  double d1 = 1. / n2;
2509  double d2 = 1. / np2;
2510  double d3 = izsqrd * rmr * EN1RYD;
2511  double d4 = d2 - d1;
2512  double result = d3 * d4;
2513 
2514  ASSERT( n > 0 );
2515  ASSERT( nprime > 0 );
2516  ASSERT( n > nprime );
2517  ASSERT( iz > 0 );
2518  ASSERT( result > 0. );
2519 
2520  if( n <= nprime )
2521  {
2522  fprintf(ioQQQ," The principal quantum numbers are such that n <= n'.\n");
2523  cdEXIT(EXIT_FAILURE);
2524  }
2525 
2526  return result;
2527 }
2528 
2529 /************************************************************************/
2530 /* hri() */
2531 /* Calculate the hydrogen radial wavefunction integral */
2532 /* for the dipole transition l'=l-1 or l'=l+1 */
2533 /* for the higher energy state n,l to the lower energy state n',l' */
2534 /* no "m" dependence */
2535 /************************************************************************/
2536 /* here we have a transition */
2537 /* from the higher energy state n,l */
2538 /* to the lower energy state n',l' */
2539 /* with a dipole selection rule on l and l' */
2540 /************************************************************************/
2541 /* */
2542 /* hri() test n,l,n',l' for domain errors and */
2543 /* swaps n,l <--> n',l' for the case l'=l+1 */
2544 /* */
2545 /* It then calls hrii() */
2546 /* */
2547 /* Dec. 6, 1999 */
2548 /* Robert Paul Bauman */
2549 /************************************************************************/
2550 
2551 /************************************************************************/
2552 /* This routine, hri(), calculates the hydrogen radial integral, */
2553 /* for the transition n,l --> n',l' */
2554 /* It is, of course, dimensionless. */
2555 /* */
2556 /* In the following, we have n > n' */
2557 /************************************************************************/
2558 
2559 inline double hri(
2560  /* principal quantum number, 1 for ground, upper level */
2561  long int n,
2562  /* angular momentum, 0 for s */
2563  long int l,
2564  /* principal quantum number, 1 for ground, lower level */
2565  long int np,
2566  /* angular momentum, 0 for s */
2567  long int lp,
2568  /* Nuclear charge, 1 for H+, 2 for He++, etc */
2569  long int iz
2570 )
2571 {
2572  DEBUG_ENTRY( "hri" );
2573 
2574  long int a;
2575  long int b;
2576  long int c;
2577  long int d;
2578  double ld1 = 0.;
2579  double Z = (double)iz;
2580 
2581  /**********************************************************************/
2582  /* from higher energy -> lower energy */
2583  /* Selection Rule for l and l' */
2584  /* dipole process only */
2585  /**********************************************************************/
2586 
2587  ASSERT( n > 0 );
2588  ASSERT( np > 0 );
2589  ASSERT( l >= 0 );
2590  ASSERT( lp >= 0 );
2591  ASSERT( n > l );
2592  ASSERT( np > lp );
2593  ASSERT( n > np || ( n == np && l == lp + 1 ));
2594  ASSERT( iz > 0 );
2595  ASSERT( lp == l + 1 || lp == l - 1 );
2596 
2597  if( l == lp + 1 )
2598  {
2599  /* Keep variable the same */
2600  a = n;
2601  b = l;
2602  c = np;
2603  d = lp;
2604  }
2605  else if( l == lp - 1 )
2606  {
2607  /* swap n,l with n',l' */
2608  a = np;
2609  b = lp;
2610  c = n;
2611  d = l;
2612  }
2613  else
2614  {
2615  printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2616  cdEXIT(EXIT_FAILURE);
2617  }
2618 
2619  /**********************************************/
2620  /* Take care of the Z-dependence here. */
2621  /**********************************************/
2622  ld1 = hrii(a, b, c, d ) / Z;
2623 
2624  return ld1;
2625 }
2626 
2627 /************************************************************************/
2628 /* hri_log10() */
2629 /* Calculate the hydrogen radial wavefunction integral */
2630 /* for the dipole transition l'=l-1 or l'=l+1 */
2631 /* for the higher energy state n,l to the lower energy state n',l' */
2632 /* no "m" dependence */
2633 /************************************************************************/
2634 /* here we have a transition */
2635 /* from the higher energy state n,l */
2636 /* to the lower energy state n',l' */
2637 /* with a dipole selection rule on l and l' */
2638 /************************************************************************/
2639 /* */
2640 /* hri_log10() test n,l,n',l' for domain errors and */
2641 /* swaps n,l <--> n',l' for the case l'=l+1 */
2642 /* */
2643 /* It then calls hrii_log() */
2644 /* */
2645 /* Dec. 6, 1999 */
2646 /* Robert Paul Bauman */
2647 /************************************************************************/
2648 
2649 inline double hri_log10( long int n, long int l, long int np, long int lp , long int iz )
2650 {
2651  /**********************************************************************/
2652  /* from higher energy -> lower energy */
2653  /* Selection Rule for l and l' */
2654  /* dipole process only */
2655  /**********************************************************************/
2656 
2657  DEBUG_ENTRY( "hri_log10()" );
2658 
2659  long int a;
2660  long int b;
2661  long int c;
2662  long int d;
2663  double ld1 = 0.;
2664  double Z = (double)iz;
2665 
2666  ASSERT( n > 0);
2667  ASSERT( np > 0);
2668  ASSERT( l >= 0);
2669  ASSERT( lp >= 0 );
2670  ASSERT( n > l );
2671  ASSERT( np > lp );
2672  ASSERT( n > np || ( n == np && l == lp + 1 ));
2673  ASSERT( iz > 0 );
2674  ASSERT( lp == l + 1 || lp == l - 1 );
2675 
2676  if( l == lp + 1)
2677  {
2678  /* Keep variable the same */
2679  a = n;
2680  b = l;
2681  c = np;
2682  d = lp;
2683  }
2684  else if( l == lp - 1 )
2685  {
2686  /* swap n,l with n',l' */
2687  a = np;
2688  b = lp;
2689  c = n;
2690  d = l;
2691  }
2692  else
2693  {
2694  printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
2695  cdEXIT(EXIT_FAILURE);
2696  }
2697 
2698  /**********************************************/
2699  /* Take care of the Z-dependence here. */
2700  /**********************************************/
2701  ld1 = hrii_log(a, b, c, d ) / Z;
2702 
2703  return ld1;
2704 }
2705 
2706 STATIC double hrii( long int n, long int l, long int np, long int lp)
2707 {
2708  /******************************************************************************/
2709  /* this routine hrii() is internal to the parent routine hri() */
2710  /* this internal routine only considers the case l=l'+1 */
2711  /* the case l=l-1 is done in the parent routine hri() */
2712  /* by the transformation n <--> n' and l <--> l' */
2713  /* THUS WE TEST FOR */
2714  /* l=l'-1 */
2715  /******************************************************************************/
2716 
2717  DEBUG_ENTRY( "hrii()" );
2718 
2719  long int a = 0, b = 0, c = 0;
2720  long int i1 = 0, i2 = 0, i3 = 0, i4 = 0;
2721 
2722  char A='a';
2723 
2724  double y = 0.;
2725  double fsf = 0.;
2726  double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
2727  double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
2728  double d00 = 0., d01 = 0.;
2729 
2730  ASSERT( l == lp + 1 );
2731 
2732  if( n == np ) /* SPECIAL CASE 1 */
2733  {
2734  /**********************************************************/
2735  /* if lp= l + 1 then it has higher energy */
2736  /* i.e. no photon */
2737  /* this is the second time we check this, oh well */
2738  /**********************************************************/
2739 
2740  if( lp != (l - 1) )
2741  {
2742  printf( "BadMagic: Energy requirements not met.\n\n" );
2743  cdEXIT(EXIT_FAILURE);
2744  }
2745 
2746  d2 = 3. / 2.;
2747  i1 = n * n;
2748  i2 = l * l;
2749  d5 = (double)(i1 - i2);
2750  d6 = sqrt(d5);
2751  d7 = (double)n * d6;
2752  d8 = d2 * d7;
2753  return d8;
2754  }
2755  else if( l == np && lp == (l - 1) ) /* A Pair of Easy Special Cases */
2756  {
2757  if( l == (n - 1) )
2758  {
2759  /**********************************************************************/
2760  /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */
2761  /* */
2762  /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */
2763  /* [(2n-1) - 1/(2n-1)]/4 */
2764  /**********************************************************************/
2765 
2766  d1 = (double)( 2*n - 2 );
2767  d2 = (double)( 2*n - 1 );
2768  d3 = d1 * d2;
2769  d4 = sqrt( d3 );
2770 
2771  d5 = (double)(4 * n * (n - 1));
2772  i1 = (2*n - 1);
2773  d6 = (double)(i1 * i1);
2774  d7 = d5/ d6;
2775  d8 = powi( d7, n );
2776 
2777  d9 = 1./d2;
2778  d10 = d2 - d9;
2779  d11 = d10 / 4.;
2780 
2781  /* Wrap it all up */
2782 
2783  d12 = d4 * d8 * d11;
2784  return d12;
2785 
2786  }
2787  else
2788  {
2789  /******************************************************************************/
2790  /* R(n,l;n',l') = R(n,l;l,l-1) */
2791  /* */
2792  /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */
2793  /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */
2794  /******************************************************************************/
2795 
2796  d2 = 1.;
2797  for( i1 = -l; i1 <= l; i1 = i1 + 1 ) /* from n-l to n+l INCLUSIVE */
2798  {
2799  d1 = (double)(n - i1);
2800  d2 = d2 * d1;
2801  }
2802  i2 = (2*l - 1);
2803  d3 = factorial( i2 );
2804  d4 = d2/d3;
2805  d4 = sqrt( d4 );
2806 
2807 
2808  d5 = (double)( 4. * n * l );
2809  i3 = (n - l);
2810  d6 = (double)( i3 * i3 );
2811  d7 = d5 / d6;
2812  d8 = powi( d7, l+1 );
2813 
2814 
2815  i4 = n + l;
2816  d9 = (double)( i3 ) / (double)( i4 );
2817  d10 = powi( d9 , i4 );
2818 
2819  d11 = d9 * d9;
2820  d12 = 1. - d11;
2821  d13 = d12 / 4.;
2822 
2823  /* Wrap it all up */
2824  d14 = d4 * d8 * d10 * d13;
2825  return d14;
2826  }
2827  }
2828 
2829  /*******************************************************************************************/
2830  /* THE GENERAL CASE */
2831  /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
2832  /* REF: D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
2833  /* For F(a,b;c;x) we have from eq.4 */
2834  /* */
2835  /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
2836  /* */
2837  /* a (1-x) (a + bx - c) */
2838  /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
2839  /* (a-c) (a-c) */
2840  /* */
2841  /* */
2842  /* A similiar recusion relation holds for b with a <--> b. */
2843  /* */
2844  /* */
2845  /* we have initial conditions */
2846  /* */
2847  /* */
2848  /* F(0) = 1 with a = -1 */
2849  /* */
2850  /* b */
2851  /* F(-1) = 1 - (---) x with a = -1 */
2852  /* c */
2853  /*******************************************************************************************/
2854 
2855  if( lp == l - 1 ) /* use recursion over "b" */
2856  {
2857  A='b';
2858  }
2859  else if( lp == l + 1 ) /* use recursion over "a" */
2860  {
2861  A='a';
2862  }
2863  else
2864  {
2865  printf(" BadMagic: Don't know what to do here.\n\n");
2866  cdEXIT(EXIT_FAILURE);
2867  }
2868 
2869  /********************************************************************/
2870  /* Calculate the whole shootin match */
2871  /* - - (1/2) */
2872  /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
2873  /* ----------- * | ---------------- | */
2874  /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
2875  /* - - */
2876  /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
2877  /* */
2878  /* This is used in the calculation of hydrogen */
2879  /* radial wave function integral for dipole transition case */
2880  /********************************************************************/
2881 
2882  fsf = fsff( n, l, np );
2883 
2884  /**************************************************************************************/
2885  /* Use a -> a' + 1 */
2886  /* _ _ */
2887  /* (a' + 1) (1 - x) | | */
2888  /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */
2889  /* (a' + 1 -c) | | */
2890  /* - - */
2891  /* */
2892  /* For the first F() in the solution of the radial integral */
2893  /* */
2894  /* a = ( -n + l + 1 ) */
2895  /* */
2896  /* a = -n + l + 1 */
2897  /* max(a) = max(-n) + max(l) + 1 */
2898  /* = -n + max(n-1) + 1 */
2899  /* = -n + n-1 + 1 */
2900  /* = 0 */
2901  /* */
2902  /* similiarly */
2903  /* */
2904  /* min(a) = min(-n) + min(l) + 1 */
2905  /* = min(-n) + 0 + 1 */
2906  /* = -n + 1 */
2907  /* */
2908  /* a -> a' + 1 implies */
2909  /* */
2910  /* max(a') = -1 */
2911  /* min(a') = -n */
2912  /**************************************************************************************/
2913 
2914  /* a plus */
2915  a = (-n + l + 1);
2916 
2917  /* for the first 2_F_1 we use b = (-n' + l) */
2918  b = (-np + l);
2919 
2920  /* c is simple */
2921  c = 2 * l;
2922 
2923  /* -4 nn' */
2924  /* where Y = -------- . */
2925  /* (n-n')^2 */
2926  d2 = (double)(n - np);
2927  d3 = d2 * d2;
2928  d4 = 1. / d3;
2929  d5 = (double)(n * np);
2930  d6 = d5 * 4.;
2931  d7 = - d6;
2932  y = d7 * d4;
2933 
2934  d00 = F21( a, b, c, y, A );
2935 
2936  /**************************************************************/
2937  /* For the second F() in the solution of the radial integral */
2938  /* */
2939  /* a = ( -n + l - 1 ) */
2940  /* */
2941  /* a = -n + l + 1 */
2942  /* max(a) = max(-n) + max(l) - 1 */
2943  /* = -n + (n - 1) - 1 */
2944  /* = -2 */
2945  /* */
2946  /* similiarly */
2947  /* */
2948  /* min(a) = min(-n) + min(l) - 1 */
2949  /* = (-n) + 0 - 1 */
2950  /* = -n - 1 */
2951  /* */
2952  /* a -> a' + 1 implies */
2953  /* */
2954  /* max(a') = -3 */
2955  /* */
2956  /* min(a') = -n - 2 */
2957  /**************************************************************/
2958 
2959  /* a minus */
2960  a = (-n + l - 1);
2961 
2962  /* for the first 2_F_1 we use b = (-n' + l) */
2963  /* and does not change */
2964  b = (-np + l);
2965 
2966  /* c is simple */
2967  c = 2 * l;
2968 
2969  /**************************************************************/
2970  /* -4 nn' */
2971  /* where Y = -------- . */
2972  /* (n-n')^2 */
2973  /**************************************************************/
2974 
2975  /**************************************************************/
2976  /* These are already calculated a few lines up */
2977  /* */
2978  /* d2 = (double) (n - np); */
2979  /* d3 = d2 * d2; */
2980  /* d4 = 1/ d3; */
2981  /* d5 = (double) (n * np); */
2982  /* d6 = d5 * 4.0; */
2983  /* d7 = - d6; */
2984  /* y = d7 * d4; */
2985  /**************************************************************/
2986 
2987  d01 = F21(a, b, c, y, A );
2988 
2989  /* Calculate */
2990  /* */
2991  /* (n-n')^2 */
2992  /* -------- */
2993  /* (n+n')^2 */
2994 
2995  i1 = (n - np);
2996  d1 = pow2( (double)i1 );
2997  i2 = (n + np);
2998  d2 = pow2( (double)i2 );
2999  d3 = d1 / d2;
3000 
3001  d4 = d01 * d3;
3002  d5 = d00 - d4;
3003  d6 = fsf * d5;
3004 
3005  ASSERT( d6 != 0. );
3006  return d6;
3007 }
3008 
3009 
3010 STATIC double hrii_log( long int n, long int l, long int np, long int lp)
3011 {
3012  /******************************************************************************/
3013  /* this routine hrii_log() is internal to the parent routine hri_log10() */
3014  /* this internal routine only considers the case l=l'+1 */
3015  /* the case l=l-1 is done in the parent routine hri_log10() */
3016  /* by the transformation n <--> n' and l <--> l' */
3017  /* THUS WE TEST FOR */
3018  /* l=l'-1 */
3019  /******************************************************************************/
3020  /**************************************************************************************/
3021  /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */
3022  /* */
3023  /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */
3024  /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */
3025  /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */
3026  /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3027  /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */
3028  /**************************************************************************************/
3029 
3030  DEBUG_ENTRY( "hrii_log()" );
3031 
3032  char A='a';
3033 
3034  double y = 0.;
3035  double log10_fsf = 0.;
3036 
3037  ASSERT( l == lp + 1 );
3038 
3039  if( n == np ) /* SPECIAL CASE 1 */
3040  {
3041  /**********************************************************/
3042  /* if lp= l + 1 then it has higher energy */
3043  /* i.e. no photon */
3044  /* this is the second time we check this, oh well */
3045  /**********************************************************/
3046 
3047  if( lp != (l - 1) )
3048  {
3049  printf( "BadMagic: l'= l+1 for n'= n.\n\n" );
3050  cdEXIT(EXIT_FAILURE);
3051  }
3052  else
3053  {
3054  /**********************************************************/
3055  /* 3 */
3056  /* R(nl:n'=n,l'=l+1) = --- n sqrt( n^2 - l^2 ) */
3057  /* 2 */
3058  /**********************************************************/
3059 
3060  long int i1 = n * n;
3061  long int i2 = l * l;
3062 
3063  double d1 = 3. / 2.;
3064  double d2 = (double)n;
3065  double d3 = (double)(i1 - i2);
3066  double d4 = sqrt(d3);
3067  double result = d1 * d2 * d4;
3068 
3069  ASSERT( d3 >= 0. );
3070  return result;
3071  }
3072  }
3073  else if( l == np && lp == l - 1 ) /* A Pair of Easy Special Cases */
3074  {
3075  if( l == n - 1 )
3076  {
3077  /**********************************************************************/
3078  /* R(n,l;n',l') = R(n,n-l;n-1,n-2) */
3079  /* */
3080  /* = [(2n-2)(2n-1)]^(1/2) [4n(n-1)/(2n-1)^2]^n * */
3081  /* [(2n-1) - 1/(2n-1)]/4 */
3082  /**********************************************************************/
3083 
3084  double d1 = (double)((2*n-2)*(2*n-1));
3085  double d2 = sqrt( d1 );
3086  double d3 = (double)(4*n*(n-1));
3087  double d4 = (double)(2*n-1);
3088  double d5 = d4*d4;
3089  double d7 = d3/d5;
3090  double d8 = powi(d7,n);
3091  double d9 = 1./d4;
3092  double d10 = d4 - d9;
3093  double d11 = 0.25*d10;
3094  double result = (d2 * d8 * d11); /* Wrap it all up */
3095 
3096  ASSERT( d1 >= 0. );
3097  ASSERT( d3 >= 0. );
3098  return result;
3099  }
3100  else
3101  {
3102  double result = 0.;
3103  double ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0., ld5 = 0., ld6 = 0., ld7 = 0.;
3104 
3105  /******************************************************************************/
3106  /* R(n,l;n',l') = R(n,l;l,l-1) */
3107  /* */
3108  /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */
3109  /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */
3110  /******************************************************************************/
3111  /**************************************/
3112  /* [(n-l) ... (n+l)] */
3113  /**************************************/
3114  /* log10[(n-l) ... (n+l)] = */
3115  /* */
3116  /* n+l */
3117  /* --- */
3118  /* > log10(j) */
3119  /* --- */
3120  /* j=n-l */
3121  /**************************************/
3122 
3123  ld1 = 0.;
3124  for( long int i1 = (n-l); i1 <= (n+l); i1++ ) /* from n-l to n+l INCLUSIVE */
3125  {
3126  double d1 = (double)(i1);
3127  ld1 += log10( d1 );
3128  }
3129 
3130  /**************************************/
3131  /* (2l-1)! */
3132  /**************************************/
3133  /* log10[ (2n-1)! ] */
3134  /**************************************/
3135 
3136  ld2 = lfactorial( 2*l - 1 );
3137 
3138  ASSERT( ((2*l)+1) >= 0);
3139 
3140  /**********************************************/
3141  /* log10( [(n-l) ... (n+l)/(2l-1)!]^(1/2) ) = */
3142  /* (1/2) log10[(n-l) ... (n+l)] - */
3143  /* (1/2) log10[ (2n-1)! ] */
3144  /**********************************************/
3145 
3146  ld3 = 0.5 * (ld1 - ld2);
3147 
3148  /**********************************************/
3149  /* [4nl/(n-l)^2]^(l+1) */
3150  /**********************************************/
3151  /* log10( [4nl/(n-l)^2]^(l+1) ) = */
3152  /* (l+1) * log10( [4nl/(n-l)^2] ) */
3153  /* */
3154  /* = (l+1)*[ log10(4nl) - 2 log10(n-l) ] */
3155  /* */
3156  /**********************************************/
3157 
3158  double d1 = (double)(l+1);
3159  double d2 = (double)(4*n*l);
3160  double d3 = (double)(n-l);
3161  double d4 = log10(d2);
3162  double d5 = log10(d3);
3163 
3164  ld4 = d1 * (d4 - 2.*d5);
3165 
3166  /**********************************************/
3167  /* [(n-l)/(n+l)]^(n+l) */
3168  /**********************************************/
3169  /* log10( [ (n-l)/(n+l) ]^(n+l) ) = */
3170  /* */
3171  /* (n+l) * [ log10(n-l) - log10(n+l) ] */
3172  /* */
3173  /**********************************************/
3174 
3175  d1 = (double)(n-l);
3176  d2 = (double)(n+l);
3177  d3 = log10( d1 );
3178  d4 = log10( d2 );
3179 
3180  ld5 = d2 * (d3 - d4);
3181 
3182  /**********************************************/
3183  /* {1-[(n-l)/(n+l)]^2}/4 */
3184  /**********************************************/
3185  /* log10[ {1-[(n-l)/(n+l)]^2}/4 ] */
3186  /**********************************************/
3187 
3188  d1 = (double)(n-l);
3189  d2 = (double)(n+l);
3190  d3 = d1/d2;
3191  d4 = d3*d3;
3192  d5 = 1. - d4;
3193  double d6 = 0.25*d5;
3194 
3195  ld6 = log10(d6);
3196 
3197  /******************************************************************************/
3198  /* R(n,l;n',l') = R(n,l;l,l-1) */
3199  /* */
3200  /* = [(n-l) ... (n+l)/(2l-1)!]^(1/2) [4nl/(n-l)^2]^(l+1) * */
3201  /* [(n-l)/(n+l)]^(n+l) {1-[(n-l)/(n+l)]^2}/4 */
3202  /******************************************************************************/
3203 
3204  ld7 = ld3 + ld4 + ld5 + ld6;
3205 
3206  result = pow( 10., ld7 );
3207 
3208  ASSERT( result > 0. );
3209  return result;
3210  }
3211  }
3212  else
3213  {
3214  double result = 0.;
3215  long int a = 0, b = 0, c = 0;
3216  double d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0., d7 = 0.;
3217  mx d00={0.0,0}, d01={0.0,0}, d02={0.0,0}, d03={0.0,0};
3218 
3219  if( lp == l - 1 ) /* use recursion over "b" */
3220  {
3221  A='b';
3222  }
3223  else if( lp == l + 1 ) /* use recursion over "a" */
3224  {
3225  A='a';
3226  }
3227  else
3228  {
3229  printf(" BadMagic: Don't know what to do here.\n\n");
3230  cdEXIT(EXIT_FAILURE);
3231  }
3232 
3233  /**************************************************************************************/
3234  /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */
3235  /* */
3236  /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */
3237  /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */
3238  /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */
3239  /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3240  /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */
3241  /**************************************************************************************/
3242 
3243  /****************************************************************************************************/
3244  /* Calculate the whole shootin match */
3245  /* - - (1/2) */
3246  /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
3247  /* fsff() = ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2)(n+n')^(-n-n') */
3248  /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3249  /* - - */
3250  /* This is used in the calculation of hydrogen radial wave function integral for dipole transitions */
3251  /****************************************************************************************************/
3252 
3253  log10_fsf = log10_fsff( n, l, np );
3254 
3255  /******************************************************************************************/
3256  /* 2_F_1( a, b; c; y ) */
3257  /* */
3258  /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */
3259  /* */
3260  /* */
3261  /* Use a -> a' + 1 */
3262  /* _ _ */
3263  /* (a' + 1) (1 - x) | | */
3264  /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + (a' + 1 + bx -c) F(a'+1) */
3265  /* (a' + 1 - c) | | */
3266  /* - - */
3267  /* */
3268  /* For the first F() in the solution of the radial integral */
3269  /* */
3270  /* a = ( -n + l + 1 ) */
3271  /* */
3272  /* a = -n + l + 1 */
3273  /* max(a) = max(-n) + max(l) + 1 */
3274  /* = -n + max(n-1) + 1 */
3275  /* = -n + n-1 + 1 */
3276  /* = 0 */
3277  /* */
3278  /* similiarly */
3279  /* */
3280  /* min(a) = min(-n) + min(l) + 1 */
3281  /* = min(-n) + 0 + 1 */
3282  /* = -n + 1 */
3283  /* */
3284  /* a -> a' + 1 implies */
3285  /* */
3286  /* max(a') = -1 */
3287  /* min(a') = -n */
3288  /******************************************************************************************/
3289 
3290  /* a plus */
3291  a = (-n + l + 1);
3292 
3293  /* for the first 2_F_1 we use b = (-n' + l) */
3294  b = (-np + l);
3295 
3296  /* c is simple */
3297  c = 2 * l;
3298 
3299  /**********************************************************************************/
3300  /* 2_F_1( a, b; c; y ) */
3301  /* */
3302  /* F21_mx(-n+l+1, -n'+l; 2l; -4nn'/[n-n']^2) */
3303  /* */
3304  /* -4 nn' */
3305  /* where Y = -------- . */
3306  /* (n-n')^2 */
3307  /* */
3308  /**********************************************************************************/
3309 
3310  d2 = (double)(n - np);
3311  d3 = d2 * d2;
3312 
3313  d4 = 1. / d3;
3314  d5 = (double)(n * np);
3315  d6 = d5 * 4.;
3316 
3317  d7 = -d6;
3318  y = d7 * d4;
3319 
3320  /**************************************************************************************************/
3321  /* THE GENERAL CASE */
3322  /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
3323  /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
3324  /* */
3325  /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
3326  /* */
3327  /* a (1-x) (a + bx - c) */
3328  /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
3329  /* (a - c) (a - c) */
3330  /* */
3331  /* */
3332  /* A similiar recusion relation holds for b with a <--> b. */
3333  /* */
3334  /* */
3335  /* we have initial conditions */
3336  /* */
3337  /* */
3338  /* F(0) = 1 with a = -1 */
3339  /* */
3340  /* b */
3341  /* F(-1) = 1 - (---) x with a = -1 */
3342  /* c */
3343  /**************************************************************************************************/
3344 
3345  /* 2_F_1( long int a, long int b, long int c, (double) y, (string) "a" or "b") */
3346  /* F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3347  d00 = F21_mx( a, b, c, y, A );
3348 
3349  /**************************************************************/
3350  /* For the second F() in the solution of the radial integral */
3351  /* */
3352  /* a = ( -n + l - 1 ) */
3353  /* */
3354  /* a = -n + l + 1 */
3355  /* max(a) = max(-n) + max(l) - 1 */
3356  /* = -n + (n - 1) - 1 */
3357  /* = -2 */
3358  /* */
3359  /* similiarly */
3360  /* */
3361  /* min(a) = min(-n) + min(l) - 1 */
3362  /* = (-n) + 0 - 1 */
3363  /* = -n - 1 */
3364  /* */
3365  /* a -> a' + 1 implies */
3366  /* */
3367  /* max(a') = -3 */
3368  /* */
3369  /* min(a') = -n - 2 */
3370  /**************************************************************/
3371 
3372  /* a minus */
3373  a = (-n + l - 1);
3374 
3375  /* for the first 2_F_1 we use b = (-n' + l) */
3376  /* and does not change */
3377  b = (-np + l);
3378 
3379  /* c is simple */
3380  c = 2 * l;
3381 
3382  /**************************************************************/
3383  /* -4 nn' */
3384  /* where Y = -------- . */
3385  /* (n-n')^2 */
3386  /**************************************************************/
3387 
3388  /**************************************************************/
3389  /* These are already calculated a few lines up */
3390  /* */
3391  /* d2 = (double) (n - np); */
3392  /* d3 = d2 * d2; */
3393  /* d4 = 1/ d3; */
3394  /* d5 = (double) (n * np); */
3395  /* d6 = d5 * 4.0; */
3396  /* d7 = - d6; */
3397  /* y = d7 * d4; */
3398  /**************************************************************/
3399 
3400  d01 = F21_mx(a, b, c, y, A );
3401 
3402  /**************************************************************************************/
3403  /* THIS HAS THE GENERAL FORM GIVEN BY (GORDAN 1929): */
3404  /* */
3405  /* R(n,l;n',l') = (-1)^(n'-1) [4(2l-1)!]^(-1) * */
3406  /* [(n+l)! (n'+l-1)!/(n-l-1)! (n'-l)!]^(1/2) * */
3407  /* (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') * */
3408  /* { F(-n+l+1,-n'+l;2l;-4nn'/[n-n']^2) */
3409  /* - (n-n')^2 (n+n')^2 F(-n+l-1,-n'+l;2l; -4nn'/[n-n']^2 ) } */
3410  /* */
3411  /* = fsf * ( F(a,b,c;y) - d3 * F(a',b',c';y) ) */
3412  /* */
3413  /* where d3 = (n-n')^2 (n+n')^2 */
3414  /* */
3415  /**************************************************************************************/
3416 
3417  /**************************************************************/
3418  /* Calculate */
3419  /* */
3420  /* (n-n')^2 */
3421  /* -------- */
3422  /* (n+n')^2 */
3423  /**************************************************************/
3424 
3425  d1 = (double)((n - np)*(n -np));
3426  d2 = (double)((n + np)*(n + np));
3427  d3 = d1 / d2;
3428 
3429  d02.x = d01.x;
3430  d02.m = d01.m * d3;
3431 
3432  while ( fabs(d02.m) > 1.0e+25 )
3433  {
3434  d02.m /= 1.0e+25;
3435  d02.x += 25;
3436  }
3437 
3438  d03.x = d00.x;
3439  d03.m = d00.m * (1. - (d02.m/d00.m) * powi( 10. , (d02.x - d00.x) ) );
3440 
3441  result = pow( 10., (log10_fsf + d03.x) ) * d03.m;
3442 
3443  ASSERT( result != 0. );
3444 
3445  /* we don't care about the correct sign of result... */
3446  return fabs(result);
3447  }
3448 }
3449 
3450 STATIC double fsff( long int n, long int l, long int np )
3451 {
3452  /****************************************************************/
3453  /* Calculate the whole shootin match */
3454  /* - - (1/2) */
3455  /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
3456  /* ----------- * | ---------------- | */
3457  /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3458  /* - - */
3459  /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3460  /* */
3461  /****************************************************************/
3462 
3463  DEBUG_ENTRY( "fsff()" );
3464 
3465  long int i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
3466  double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
3467  double sigma = 1.;
3468 
3469  /****************************************************************
3470  * Calculate the whole shootin match *
3471  * (-1)^(n'-1) | (n+l)! (n'+l-1)! | *
3472  * ----------- * | ---------------- | *
3473  * [4 (2l-1)!] | (n-l-1)! (n'-l)! | *
3474  * - - *
3475  * * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') *
3476  * *
3477  ****************************************************************/
3478 
3479  /* Calculate (-1)^(n'-l) */
3480  if( is_odd(np - l) )
3481  {
3482  sigma *= -1.;
3483  }
3484  ASSERT( sigma != 0. );
3485 
3486  /*********************/
3487  /* Calculate (2l-1)! */
3488  /*********************/
3489  i1 = (2*l - 1);
3490  if( i1 < 0 )
3491  {
3492  printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3493  cdEXIT(EXIT_FAILURE);
3494  }
3495 
3496  /****************************************************************/
3497  /* Calculate the whole shootin match */
3498  /* - - (1/2) */
3499  /* (-1)^(n'-1) | (n+l)! (n'+l-1)! | */
3500  /* ----------- * | ---------------- | */
3501  /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3502  /* - - */
3503  /* * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3504  /* */
3505  /****************************************************************/
3506 
3507  d0 = factorial( i1 );
3508  d1 = 4. * d0;
3509  d2 = 1. / d1;
3510 
3511  /**********************************************************************/
3512  /* We want the (negitive) of this */
3513  /* since we really are interested in */
3514  /* [(2l-1)!]^-1 */
3515  /**********************************************************************/
3516 
3517  sigma = sigma * d2;
3518  ASSERT( sigma != 0. );
3519 
3520  /**********************************************************************/
3521  /* Calculate (4 n n')^(l+1) */
3522  /* powi( m , n) calcs m^n */
3523  /* returns long double with m,n ints */
3524  /**********************************************************************/
3525 
3526  i0 = 4 * n * np;
3527  i1 = l + 1;
3528  d2 = powi( (double)i0 , i1 );
3529  sigma = sigma * d2;
3530  ASSERT( sigma != 0. );
3531 
3532  /* Calculate (n-n')^(n+n'-2l-2) */
3533  i0 = n - np;
3534  i1 = n + np - 2*l - 2;
3535  d2 = powi( (double)i0 , i1 );
3536  sigma = sigma * d2;
3537  ASSERT( sigma != 0. );
3538 
3539  /* Calculate (n+n')^(-n-n') */
3540  i0 = n + np;
3541  i1 = -n - np;
3542  d2 = powi( (double)i0 , i1 );
3543  sigma = sigma * d2;
3544  ASSERT( sigma != 0. );
3545 
3546  /**********************************************************************/
3547  /* - - (1/2) */
3548  /* | (n+l)! (n'+l-1)! | */
3549  /* Calculate | ---------------- | */
3550  /* | (n-l-1)! (n'-l)! | */
3551  /* - - */
3552  /**********************************************************************/
3553 
3554  i1 = n + l;
3555  if( i1 < 0 )
3556  {
3557  printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3558  cdEXIT(EXIT_FAILURE);
3559  }
3560  d1 = factorial( i1 );
3561 
3562  i2 = np + l - 1;
3563  if( i2 < 0 )
3564  {
3565  printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3566  cdEXIT(EXIT_FAILURE);
3567  }
3568  d2 = factorial( i2 );
3569 
3570  i3 = n - l - 1;
3571  if( i3 < 0 )
3572  {
3573  printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3574  cdEXIT(EXIT_FAILURE);
3575  }
3576  d3 = factorial( i3 );
3577 
3578  i4 = np - l;
3579  if( i4 < 0 )
3580  {
3581  printf( "BadMagic: Relational error amongst n, l, n' and l'\n" );
3582  cdEXIT(EXIT_FAILURE);
3583  }
3584  d4 = factorial( i4 );
3585 
3586  ASSERT( d3 != 0. );
3587  ASSERT( d4 != 0. );
3588 
3589  /* Do this a different way to prevent overflow */
3590  /* d5 = (sqrt(d1 *d2)); */
3591  d5 = sqrt(d1)*sqrt(d2);
3592  d5 /= sqrt(d3);
3593  d5 /= sqrt(d4);
3594 
3595  sigma = sigma * d5;
3596 
3597  ASSERT( sigma != 0. );
3598  return sigma;
3599 }
3600 
3601 /**************************log version*******************************/
3602 STATIC double log10_fsff( long int n, long int l, long int np )
3603 {
3604  /******************************************************************************************************/
3605  /* Calculate the whole shootin match */
3606  /* - - (1/2) */
3607  /* 1 | (n+l)! (n'+l-1)! | */
3608  /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3609  /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3610  /* - - */
3611  /******************************************************************************************************/
3612 
3613  DEBUG_ENTRY( "log10_fsff()" );
3614 
3615  double d0 = 0., d1 = 0.;
3616  double ld0 = 0., ld1 = 0., ld2 = 0., ld3 = 0., ld4 = 0.;
3617  double result = 0.;
3618 
3619  /******************************************************************************************************/
3620  /* Calculate the log10 of the whole shootin match */
3621  /* - - (1/2) */
3622  /* 1 | (n+l)! (n'+l-1)! | */
3623  /* ----------- * | ---------------- | * (4 n n')^(l+1) (n-n')^(n+n'-2l-2) (n+n')^(-n-n') */
3624  /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3625  /* - - */
3626  /******************************************************************************************************/
3627 
3628  /**********************/
3629  /* Calculate (2l-1)! */
3630  /**********************/
3631 
3632  d0 = (double)(2*l - 1);
3633  ASSERT( d0 != 0. );
3634 
3635  /******************************************************************************************************/
3636  /* Calculate the whole shootin match */
3637  /* - - (1/2) */
3638  /* 1 | (n+l)! (n'+l-1)! | */
3639  /* ----------- * | ---------------- | * (4 n n')^(l+1) |(n-n')^(n+n'-2l-2)| (n+n')^(-n-n') */
3640  /* [4 (2l-1)!] | (n-l-1)! (n'-l)! | */
3641  /* - - */
3642  /******************************************************************************************************/
3643 
3644  ld0 = lfactorial( 2*l - 1 );
3645  ld1 = log10(4.);
3646  result = -(ld0 + ld1);
3647  ASSERT( result != 0. );
3648 
3649  /**********************************************************************/
3650  /* Calculate (4 n n')^(l+1) */
3651  /* powi( m , n) calcs m^n */
3652  /* returns long double with m,n ints */
3653  /**********************************************************************/
3654 
3655  d0 = (double)(4 * n * np);
3656  d1 = (double)(l + 1);
3657  result += d1 * log10(d0);
3658  ASSERT( d0 >= 0. );
3659  ASSERT( d1 != 0. );
3660 
3661  /**********************************************************************/
3662  /* Calculate |(n-n')^(n+n'-2l-2)| */
3663  /* NOTE: Here we are interested only */
3664  /* magnitude of (n-n')^(n+n'-2l-2) */
3665  /**********************************************************************/
3666 
3667  d0 = (double)(n - np);
3668  d1 = (double)(n + np - 2*l - 2);
3669  result += d1 * log10(fabs(d0));
3670  ASSERT( fabs(d0) > 0. );
3671  ASSERT( d1 != 0. );
3672 
3673  /* Calculate (n+n')^(-n-n') */
3674  d0 = (double)(n + np);
3675  d1 = (double)(-n - np);
3676  result += d1 * log10(d0);
3677  ASSERT( d0 > 0. );
3678  ASSERT( d1 != 0. );
3679 
3680  /**********************************************************************/
3681  /* - - (1/2) */
3682  /* | (n+l)! (n'+l-1)! | */
3683  /* Calculate | ---------------- | */
3684  /* | (n-l-1)! (n'-l)! | */
3685  /* - - */
3686  /**********************************************************************/
3687 
3688  ASSERT( n+l > 0 );
3689  ld0 = lfactorial( n + l );
3690 
3691  ASSERT( np+l-1 > 0 );
3692  ld1 = lfactorial( np + l - 1 );
3693 
3694  ASSERT( n-l-1 >= 0 );
3695  ld2 = lfactorial( n - l - 1 );
3696 
3697  ASSERT( np-l >= 0 );
3698  ld3 = lfactorial( np - l );
3699 
3700  ld4 = 0.5*((ld0+ld1)-(ld2+ld3));
3701 
3702  result += ld4;
3703  ASSERT( result != 0. );
3704  return result;
3705 }
3706 
3707 /***************************************************************************/
3708 /* Find the Oscillator Strength for hydrogen for any */
3709 /* transition n,l --> n',l' */
3710 /* returns a double */
3711 /***************************************************************************/
3712 
3713 /************************************************************************/
3714 /* Find the Oscillator Strength for hydrogen for any */
3715 /* transition n,l --> n',l' */
3716 /* returns a double */
3717 /* */
3718 /* Einstein A() for the transition from the */
3719 /* initial state n,l to the finial state n',l' */
3720 /* require the Oscillator Strength f() */
3721 /* */
3722 /* hbar w max(l,l') | | 2 */
3723 /* f(n,l;n',l') = - -------- ------------ | R(n,l;n',l') | */
3724 /* 3 R_oo ( 2l + 1 ) | | */
3725 /* */
3726 /* */
3727 /* */
3728 /* E(n,l;n',l') max(l,l') | | 2 */
3729 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
3730 /* 3 R_oo ( 2l + 1 ) | | */
3731 /* */
3732 /* */
3733 /* See for example Gordan Drake's */
3734 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */
3735 /* */
3736 /* Here R_oo is the infinite mass Rydberg length */
3737 /* */
3738 /* */
3739 /* h c */
3740 /* R_oo --- = 13.605698 eV */
3741 /* {e} */
3742 /* */
3743 /* */
3744 /* R_oo = 2.179874e-11 ergs */
3745 /* */
3746 /* w = omega */
3747 /* = frequency of transition from n,l to n',l' */
3748 /* */
3749 /* */
3750 /* */
3751 /* here g_k are statistical weights obtained from */
3752 /* the appropriate angular momentum quantum numbers */
3753 /************************************************************************/
3754 
3755 /********************************************************************************/
3756 /* Calc the Oscillator Strength f(*) given by */
3757 /* */
3758 /* E(n,l;n',l') max(l,l') | | 2 */
3759 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
3760 /* 3 R_oo ( 2l + 1 ) | | */
3761 /* */
3762 /* See for example Gordan Drake's */
3763 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */
3764 /********************************************************************************/
3765 
3766 /************************************************************************/
3767 /* Calc the Oscillator Strength f(*) given by */
3768 /* */
3769 /* E(n,l;n',l') max(l,l') | | 2 */
3770 /* f(n,l;n',l') = - ------------ ------------ | R(n,l;n',l') | */
3771 /* 3 R_oo ( 2l + 1 ) | | */
3772 /* */
3773 /* f(n,l;n',l') is dimensionless. */
3774 /* */
3775 /* See for example Gordan Drake's */
3776 /* Atomic, Molecular, & Optical Physics Handbook pg.638 */
3777 /* */
3778 /* In the following, we have n > n' */
3779 /************************************************************************/
3780 
3781 inline double OscStr_f(
3782  /* IN THE FOLLOWING WE HAVE n > n' */
3783  /* principal quantum number, 1 for ground, upper level */
3784  long int n,
3785  /* angular momentum, 0 for s */
3786  long int l,
3787  /* principal quantum number, 1 for ground, lower level */
3788  long int np,
3789  /* angular momentum, 0 for s */
3790  long int lp,
3791  /* Nuclear charge, 1 for H+, 2 for He++, etc */
3792  long int iz
3793 )
3794 {
3795  double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3796  long int i1 = 0, i2 = 0;
3797 
3798  if( l > lp )
3799  i1 = l;
3800  else
3801  i1 = lp;
3802 
3803  i2 = 2*lp + 1;
3804  d0 = 1. / 3.;
3805  d1 = (double)i1 / (double)i2;
3806  /* hv() returns energy in ergs */
3807  d2 = hv( n, np, iz );
3808  d3 = d2 / EN1RYD;
3809  d4 = hri( n, l, np, lp ,iz );
3810  d5 = d4 * d4;
3811 
3812  d6 = d0 * d1 * d3 * d5;
3813 
3814  return d6;
3815 }
3816 
3817 /************************log version ***************************/
3818 inline double OscStr_f_log10( long int n , long int l , long int np , long int lp , long int iz )
3819 {
3820  double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0., d6 = 0.;
3821  long int i1 = 0, i2 = 0;
3822 
3823  if( l > lp )
3824  i1 = l;
3825  else
3826  i1 = lp;
3827 
3828  i2 = 2*lp + 1;
3829  d0 = 1. / 3.;
3830  d1 = (double)i1 / (double)i2;
3831  /* hv() returns energy in ergs */
3832  d2 = hv( n, np, iz );
3833  d3 = d2 / EN1RYD;
3834  d4 = hri_log10( n, l, np, lp ,iz );
3835  d5 = d4 * d4;
3836 
3837  d6 = d0 * d1 * d3 * d5;
3838 
3839  return d6;
3840 }
3841 
3842 STATIC double F21( long int a , long int b, long int c, double y, char A )
3843 {
3844  DEBUG_ENTRY( "F21()" );
3845 
3846  double d1 = 0.;
3847  long int i1 = 0;
3848  double *yV;
3849 
3850  /**************************************************************/
3851  /* A must be either 'a' or 'b' */
3852  /* and is use to determine which */
3853  /* variable recursion will be over */
3854  /**************************************************************/
3855 
3856  ASSERT( A == 'a' || A == 'b' );
3857 
3858  if( A == 'b' )
3859  {
3860  /* if the recursion is over "b" */
3861  /* then make it over "a" by switching these around */
3862  i1 = a;
3863  a = b;
3864  b = i1;
3865  A = 'a';
3866  }
3867 
3868  /**************************************************************************************/
3869  /* malloc space for the (dynamic) 1-d array */
3870  /* 2_F_1 works via recursion and needs room to store intermediate results */
3871  /* Here the + 5 is a safety margin */
3872  /**************************************************************************************/
3873 
3874  /*Must use CALLOC*/
3875  yV = (double*)CALLOC( sizeof(double), (size_t)(-a + 5) );
3876 
3877  /**********************************************************************************************/
3878  /* begin sanity check, check order, and that none negative */
3879  /* THE GENERAL CASE */
3880  /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
3881  /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
3882  /* */
3883  /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
3884  /* */
3885  /* a (1-x) (a + bx - c) */
3886  /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
3887  /* (a-c) (a-c) */
3888  /* */
3889  /* */
3890  /* A similiar recusion relation holds for b with a <--> b. */
3891  /* */
3892  /* */
3893  /* we have initial conditions */
3894  /* */
3895  /* */
3896  /* F(0) = 1 with a = -1 */
3897  /* */
3898  /* b */
3899  /* F(-1) = 1 - (---) x with a = -1 */
3900  /* c */
3901  /* */
3902  /* For the first F() in the solution of the radial integral */
3903  /* */
3904  /* a = ( -n + l + 1 ) */
3905  /* */
3906  /* a = -n + l + 1 */
3907  /* max(a) = max(-n) + max(l) + 1 */
3908  /* = max(-n) + max(n - 1) + 1 */
3909  /* = max(-n + n - 1) + 1 */
3910  /* = max(-1) + 1 */
3911  /* = 0 */
3912  /* */
3913  /* similiarly */
3914  /* */
3915  /* min(a) = min(-n) + min(l) + 1 */
3916  /* = min(-n) + 0 + 1 */
3917  /* = (-n) + 0 + 1 */
3918  /* = -n + 1 */
3919  /* */
3920  /* a -> a' + 1 implies */
3921  /* */
3922  /* max(a') = -1 */
3923  /* min(a') = -n */
3924  /* */
3925  /* For the second F() in the solution of the radial integral */
3926  /* */
3927  /* a = ( -n + l - 1 ) */
3928  /* */
3929  /* a = -n + l + 1 */
3930  /* max(a) = max(-n) + max(l) - 1 */
3931  /* = -n + (n - 1) - 1 */
3932  /* = -2 */
3933  /* */
3934  /* similiarly */
3935  /* */
3936  /* min(a) = min(-n) + min(l) - 1 */
3937  /* = (-n) + 0 - 1 */
3938  /* = -n - 1 */
3939  /* */
3940  /* a -> a' + 1 implies */
3941  /* */
3942  /* max(a') = -3 */
3943  /* min(a') = -n - 2 */
3944  /**********************************************************************************************/
3945 
3946  ASSERT( a <= 0 );
3947  ASSERT( b <= 0 );
3948  ASSERT( c >= 0 );
3949 
3950  d1 = F21i(a, b, c, y, yV );
3951  free( yV );
3952  return d1;
3953 }
3954 
3955 STATIC mx F21_mx( long int a , long int b, long int c, double y, char A )
3956 {
3957  DEBUG_ENTRY( "F21_mx()" );
3958 
3959  mx result_mx = {0.0,0};
3960  mxq *yV = NULL;
3961 
3962  /**************************************************************/
3963  /* A must be either 'a' or 'b' */
3964  /* and is use to determine which */
3965  /* variable recursion will be over */
3966  /**************************************************************/
3967 
3968  ASSERT( A == 'a' || A == 'b' );
3969 
3970  if( A == 'b' )
3971  {
3972  /* if the recursion is over "b" */
3973  /* then make it over "a" by switching these around */
3974  long int i1 = a;
3975  a = b;
3976  b = i1;
3977  A = 'a';
3978  }
3979 
3980  /**************************************************************************************/
3981  /* malloc space for the (dynamic) 1-d array */
3982  /* 2_F_1 works via recursion and needs room to store intermediate results */
3983  /* Here the + 5 is a safety margin */
3984  /**************************************************************************************/
3985 
3986  /*Must use CALLOC*/
3987  yV = (mxq *)CALLOC( sizeof(mxq), (size_t)(-a + 5) );
3988 
3989  /**********************************************************************************************/
3990  /* begin sanity check, check order, and that none negative */
3991  /* THE GENERAL CASE */
3992  /* USE RECURSION RELATION FOR HYPERGEOMETRIC FUNCTIONS */
3993  /* for F(a,b;c;x) we have from eq.4 D. Hoang-Bing Astron. Astrophys. 238: 449-451 (1990) */
3994  /* */
3995  /* (a-c) F(a-1) = a (1-x) [ F(a) - F(a-1) ] + (a + bx - c) F(a) */
3996  /* */
3997  /* a (1-x) (a + bx - c) */
3998  /* F(a-1) = --------- [ F(a) - F(a-1) ] + -------------- F(a) */
3999  /* (a-c) (a-c) */
4000  /* */
4001  /* */
4002  /* A similiar recusion relation holds for b with a <--> b. */
4003  /* */
4004  /* */
4005  /* we have initial conditions */
4006  /* */
4007  /* */
4008  /* F(0) = 1 with a = -1 */
4009  /* */
4010  /* b */
4011  /* F(-1) = 1 - (---) x with a = -1 */
4012  /* c */
4013  /* */
4014  /* For the first F() in the solution of the radial integral */
4015  /* */
4016  /* a = ( -n + l + 1 ) */
4017  /* */
4018  /* a = -n + l + 1 */
4019  /* max(a) = max(-n) + max(l) + 1 */
4020  /* = max(-n) + max(n - 1) + 1 */
4021  /* = max(-n + n - 1) + 1 */
4022  /* = max(-1) + 1 */
4023  /* = 0 */
4024  /* */
4025  /* similiarly */
4026  /* */
4027  /* min(a) = min(-n) + min(l) + 1 */
4028  /* = min(-n) + 0 + 1 */
4029  /* = (-n) + 0 + 1 */
4030  /* = -n + 1 */
4031  /* */
4032  /* a -> a' + 1 implies */
4033  /* */
4034  /* max(a') = -1 */
4035  /* min(a') = -n */
4036  /* */
4037  /* For the second F() in the solution of the radial integral */
4038  /* */
4039  /* a = ( -n + l - 1 ) */
4040  /* */
4041  /* a = -n + l + 1 */
4042  /* max(a) = max(-n) + max(l) - 1 */
4043  /* = -n + (n - 1) - 1 */
4044  /* = -2 */
4045  /* */
4046  /* similiarly */
4047  /* */
4048  /* min(a) = min(-n) + min(l) - 1 */
4049  /* = (-n) + 0 - 1 */
4050  /* = -n - 1 */
4051  /* */
4052  /* a -> a' + 1 implies */
4053  /* */
4054  /* max(a') = -3 */
4055  /* min(a') = -n - 2 */
4056  /**********************************************************************************************/
4057 
4058  ASSERT( a <= 0 );
4059  ASSERT( b <= 0 );
4060  ASSERT( c >= 0 );
4061 
4062  result_mx = F21i_log(a, b, c, y, yV );
4063  free( yV );
4064  return result_mx;
4065 }
4066 
4067 STATIC double F21i(long int a, long int b, long int c, double y, double *yV )
4068 {
4069  DEBUG_ENTRY( "F21i()" );
4070 
4071  double d0 = 0., d1 = 0., d2 = 0., d3 = 0., d4 = 0., d5 = 0.;
4072  double d8 = 0., d9 = 0., d10 = 0., d11 = 0., d12 = 0., d13 = 0., d14 = 0.;
4073  long int i1 = 0, i2 = 0;
4074 
4075  if( a == 0 )
4076  {
4077  return 1.;
4078  }
4079  else if( a == -1 )
4080  {
4081  ASSERT( c != 0 );
4082  d1 = (double)b;
4083  d2 = (double)c;
4084  d3 = y * (d1/d2);
4085  d4 = 1. - d3;
4086  return d4;
4087  }
4088  /* Check to see if y(-a) != 0 in a very round about way to avoid lclint:error:13 */
4089  else if( yV[-a] != 0. )
4090  {
4091  /* Return the stored result */
4092  return yV[-a];
4093  }
4094  else
4095  {
4096  /******************************************************************************************/
4097  /* - - */
4098  /* (a)(1 - y) | | (a + bx + c) */
4099  /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */
4100  /* (a - c) | | (a - c) */
4101  /* - - */
4102  /* */
4103  /* */
4104  /* */
4105  /* */
4106  /* */
4107  /* with F(0) = 1 */
4108  /* b */
4109  /* and F(-1) = 1 - (---) y */
4110  /* c */
4111  /* */
4112  /* */
4113  /* */
4114  /* Use a -> a' + 1 */
4115  /* _ _ */
4116  /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */
4117  /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */
4118  /* (a' + 1 - c) | | (a' + 1 - c) */
4119  /* - - */
4120  /* */
4121  /* For the first F() in the solution of the radial integral */
4122  /* */
4123  /* a = ( -n + l + 1 ) */
4124  /* */
4125  /* a = -n + l + 1 */
4126  /* max(a) = max(-n) + max(l) + 1 */
4127  /* = -n + max(n-1) + 1 */
4128  /* = -n + n-1 + 1 */
4129  /* = 0 */
4130  /* */
4131  /* similiarly */
4132  /* */
4133  /* min(a) = min(-n) + min(l) + 1 */
4134  /* = min(-n) + 0 + 1 */
4135  /* = -n + 1 */
4136  /* */
4137  /* a -> a' + 1 implies */
4138  /* */
4139  /* max(a') = -1 */
4140  /* min(a') = -n */
4141  /******************************************************************************************/
4142 
4143  i1= a + 1;
4144  i2= a + 1 - c;
4145  d0= (double)i2;
4146  ASSERT( i2 != 0 );
4147  d1= 1. - y;
4148  d2= (double)i1 * d1;
4149  d3= d2 / d0;
4150  d4= (double)b * y;
4151  d5= d0 + d4;
4152 
4153  d8= F21i( (long int)(a + 1), b, c, y, yV );
4154  d9= F21i( (long int)(a + 2), b, c, y, yV );
4155 
4156  d10= d8 - d9;
4157  d11= d3 * d10;
4158  d12= d5 / d0;
4159  d13= d12 * d8;
4160  d14= d11 + d13;
4161 
4162  /* Store the result for later use */
4163  yV[-a] = d14;
4164  return d14;
4165  }
4166 }
4167 
4168 STATIC mx F21i_log( long int a, long int b, long int c, double y, mxq *yV )
4169 {
4170  DEBUG_ENTRY( "F21i_log()" );
4171 
4172  mx result_mx = {0.0,0};
4173 
4174  if( yV[-a].q != 0. )
4175  {
4176  /* Return the stored result */
4177  return yV[-a].mx;
4178  }
4179  else if( a == 0 )
4180  {
4181  ASSERT( yV[-a].q == 0. );
4182 
4183  result_mx.m = 1.;
4184  result_mx.x = 0;
4185 
4186  ASSERT( yV[-a].mx.m == 0. );
4187  ASSERT( yV[-a].mx.x == 0 );
4188 
4189  yV[-a].q = 1;
4190  yV[-a].mx = result_mx;
4191  return result_mx;
4192  }
4193  else if( a == -1 )
4194  {
4195  double d1 = (double)b;
4196  double d2 = (double)c;
4197  double d3 = y * (d1/d2);
4198 
4199  ASSERT( yV[-a].q == 0. );
4200  ASSERT( c != 0 );
4201  ASSERT( y != 0. );
4202 
4203  result_mx.m = 1. - d3;
4204  result_mx.x = 0;
4205 
4206  while ( fabs(result_mx.m) > 1.0e+25 )
4207  {
4208  result_mx.m /= 1.0e+25;
4209  result_mx.x += 25;
4210  }
4211 
4212  ASSERT( yV[-a].mx.m == 0. );
4213  ASSERT( yV[-a].mx.x == 0 );
4214 
4215  yV[-a].q = 1;
4216  yV[-a].mx = result_mx;
4217  return result_mx;
4218  }
4219  else
4220  {
4221  /******************************************************************************************/
4222  /* - - */
4223  /* (a)(1 - y) | | (a + bx + c) */
4224  /* F(a-1) = -------------- | F(a) - F(a+1) | + --------------- F(a+1) */
4225  /* (a - c) | | (a - c) */
4226  /* - - */
4227  /* */
4228  /* */
4229  /* with F(0) = 1 */
4230  /* b */
4231  /* and F(-1) = 1 - (---) y */
4232  /* c */
4233  /* */
4234  /* */
4235  /* */
4236  /* Use a -> a' + 1 */
4237  /* _ _ */
4238  /* (a' + 1) (1 - x) | | (a' + 1 + bx - c) */
4239  /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ----------------- F(a'+1) */
4240  /* (a' + 1 - c) | | (a' + 1 - c) */
4241  /* - - */
4242  /* */
4243  /* For the first F() in the solution of the radial integral */
4244  /* */
4245  /* a = ( -n + l + 1 ) */
4246  /* */
4247  /* a = -n + l + 1 */
4248  /* max(a) = max(-n) + max(l) + 1 */
4249  /* = -n + max(n-1) + 1 */
4250  /* = -n + n-1 + 1 */
4251  /* = 0 */
4252  /* */
4253  /* similiarly */
4254  /* */
4255  /* min(a) = min(-n) + min(l) + 1 */
4256  /* = min(-n) + 0 + 1 */
4257  /* = -n + 1 */
4258  /* */
4259  /* a -> a' + 1 implies */
4260  /* */
4261  /* max(a') = -1 */
4262  /* min(a') = -n */
4263  /******************************************************************************************/
4264 
4265  mx d8 = {0.0,0}, d9 = {0.0,0}, d10 = {0.0,0}, d11 = {0.0,0};
4266 
4267  double db = (double)b;
4268  double d00 = (double)(a + 1 - c);
4269  double d0 = (double)(a + 1);
4270  double d1 = 1. - y;
4271  double d2 = d0 * d1;
4272  double d3 = d2 / d00;
4273  double d4 = db * y;
4274  double d5 = d00 + d4;
4275  double d6 = d5 / d00;
4276 
4277  ASSERT( yV[-a].q == 0. );
4278 
4279  /******************************************************************************************/
4280  /* _ _ */
4281  /* (a' + 1) (1 - x) | | (a' + 1 - c) + b*x */
4282  /* F(a') = ----------------- | F(a'+1) - F(a'+2) | + ------------------ F(a'+1) */
4283  /* (a' + 1 - c) | | (a' + 1 - c) */
4284  /* - - */
4285  /******************************************************************************************/
4286 
4287  d8= F21i_log( (a + 1), b, c, y, yV );
4288  d9= F21i_log( (a + 2), b, c, y, yV );
4289 
4290  if( (d8.m) != 0. )
4291  {
4292  d10.x = d8.x;
4293  d10.m = d8.m * (1. - (d9.m/d8.m) * powi( 10., (d9.x - d8.x)));
4294  }
4295  else
4296  {
4297  d10.m = -d9.m;
4298  d10.x = d9.x;
4299  }
4300 
4301  d10.m *= d3;
4302 
4303  d11.x = d8.x;
4304  d11.m = d6 * d8.m;
4305 
4306  if( (d11.m) != 0. )
4307  {
4308  result_mx.x = d11.x;
4309  result_mx.m = d11.m * (1. + (d10.m/d11.m) * powi( 10. , (d10.x - d11.x) ) );
4310  }
4311  else
4312  {
4313  result_mx = d10;
4314  }
4315 
4316  while ( fabs(result_mx.m) > 1.0e+25 )
4317  {
4318  result_mx.m /= 1.0e+25;
4319  result_mx.x += 25;
4320  }
4321 
4322  /* Store the result for later use */
4323  yV[-a].q = 1;
4324  yV[-a].mx = result_mx;
4325  return result_mx;
4326  }
4327 }
4328 
4329 /********************************************************************************/
4330 
4331 inline void normalize_mx( mx& target )
4332 {
4333  while( fabs(target.m) > 1.0e+25 )
4334  {
4335  target.m /= 1.0e+25;
4336  target.x += 25;
4337  }
4338  while( fabs(target.m) < 1.0e-25 )
4339  {
4340  target.m *= 1.0e+25;
4341  target.x -= 25;
4342  }
4343  return;
4344 }
4345 
4346 inline mx add_mx( const mx& a, const mx& b )
4347 {
4348  mx result = {0.0,0};
4349 
4350  if( a.m != 0. )
4351  {
4352  result.x = a.x;
4353  result.m = a.m * (1. + (b.m/a.m) * powi( 10. , (b.x - a.x) ) );
4354  }
4355  else
4356  {
4357  result = b;
4358  }
4359  normalize_mx( result );
4360  return result;
4361 }
4362 
4363 inline mx sub_mx( const mx& a, const mx& b )
4364 {
4365  mx result = {0.0,0};
4366  mx minusb = b;
4367  minusb.m = -minusb.m;
4368 
4369  result = add_mx( a, minusb );
4370  normalize_mx( result );
4371 
4372  return result;
4373 }
4374 
4375 inline mx mxify( double a )
4376 {
4377  mx result_mx = {0.0,0};
4378 
4379  result_mx.x = 0;
4380  result_mx.m = a;
4381  normalize_mx( result_mx );
4382 
4383  return result_mx;
4384 }
4385 
4386 inline double unmxify( const mx& a_mx )
4387 {
4388  return a_mx.m * powi( 10., a_mx.x );
4389 }
4390 
4391 inline mx mxify_log10( double log10_a )
4392 {
4393  mx result_mx = {0.0,0};
4394 
4395  while ( log10_a > 25.0 )
4396  {
4397  log10_a -= 25.0;
4398  result_mx.x += 25;
4399  }
4400 
4401  while ( log10_a < -25.0 )
4402  {
4403  log10_a += 25.0;
4404  result_mx.x -= 25;
4405  }
4406 
4407  result_mx.m = pow(10., log10_a);
4408 
4409  return result_mx;
4410 }
4411 
4412 inline mx mult_mx( const mx& a, const mx& b )
4413 {
4414  mx result = {0.0,0};
4415 
4416  result.m = (a.m * b.m);
4417  result.x = (a.x + b.x);
4418  normalize_mx( result );
4419 
4420  return result;
4421 }
4422 
4423 inline double log10_prodxx( long int lp, double Ksqrd )
4424 {
4425  /**********************************************************************/
4426  /* | s=l' | s=l' */
4427  /* | ----- | --- */
4428  /* log10 | | | (1 + s^2 K^2) | = > log10((1 + s^2 K^2)) */
4429  /* | | | | --- */
4430  /* | s=0 | s=0 */
4431  /**********************************************************************/
4432 
4433  if( lp == 0 )
4434  return 0.;
4435 
4436  double partsum = 0.;
4437  for( long int s = 1; s <= lp; s++ )
4438  {
4439  double s2 = pow2((double)s);
4440  partsum += log10( 1. + s2*Ksqrd );
4441 
4442  ASSERT( partsum >= 0. );
4443  }
4444  return partsum;
4445 }

Generated for cloudy by doxygen 1.8.3.1