//  This file contains useful library functions for computing
//  various polynomials, Clebsch-Gordan coefficients, factorials etc.

#include <iostream> 
#include <cmath> 
using namespace std;
// useful functions 
inline int MAX(int a, int b)   {return (((a) > (b)) ? (a) : (b) );}
inline int MIN(int a, int b)   {return ( ((a) < (b)) ? (a) : (b) );}


double fac_ratio(int, int);

/*  Function to compute generalized Laguerre polynomials
    alpha cannot be smaller or equal than  -1.0
    The variable n has to be larger or equal 0
*/
double LaguerreGeneral( int n, double alpha, double x)
{
  double *glaguerre; 
  if ( alpha <= -1.0 ) {
    cout << "LAGUERRE_GENERAL - Fatal error!" << endl;
    cout << "The input value of ALPHA is=  "  << alpha  << endl;
    cout << "but ALPHA must be greater than -1." << endl;
  }
  glaguerre = new double[n+1];
  if ( n >= 0 ) {
    for (int  i = 0; i <= n; i++) glaguerre[i] = 0.0;
    glaguerre[0] = 1.0;
    if ( n > 0 ){
      glaguerre[1] = 1.0+alpha-x;
      // recursion relation for generalized Laguerre polynomials
      for (int  i = 2; i <= n; i++){
          glaguerre[i] = ((2.0*i-1.0+alpha-x)*glaguerre[i-1]+
	              (1.0-i-alpha)*glaguerre[i-2])/((float) i);
      }
    }
  }
  double GLaguerre = glaguerre[n];
  delete [] glaguerre;
  return GLaguerre;
}   // end function glaguerre
 
/* Associated Legendre polynomials
   Before calling make sure that 
   you  check whether m, l and x are ok, that is, perform the test
   if ((m < 0) || (m > l) ||  (abs(x) > 1.0)), then don't call this 
   function!
*/
double AssociateLegendrePolynomials(int l, int m, double x)
{
  double fact,pll,pmm,pmmp1,somx2, legendrepolynomials;
  //  calculate now pmm as starting point for iterations
  pmm=1.0; 
  if (m > 0) {
    somx2=sqrt((1.0-x)*(1.0+x));
    fact=1.0;
    for (int i=1;i <=  m; i++){
      pmm = -fact*somx2*pmm; 
      fact += 2.0;
    }
  }
  //  if l == m we do not need to use recursion relation
  if (l == m){
    legendrepolynomials = pmm;
  }
  //  recursive relation for associated Legendre polynomials
  else{
    pmmp1=x*(2*m+1)*pmm;
    //  analytical formula for the case l == m+1
    if (l == (m+1)){
      legendrepolynomials=pmmp1;
    }
    else{
      for ( int ll=m+2; ll <= l; ll++){
	pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
	pmm=pmmp1;
	pmmp1=pll;
      }
      legendrepolynomials = pll;
    }
  }
  return legendrepolynomials;
} // end function AssociateLegendrePolynomials

// Compute factorial as sum of logarithms, should be stored as array
// by the calling function  if it is used many times
double factorial (int m)
{
  double fac = 0.0;
  for (int i=1; i <= m; i++)fac += log((float) i);
  return fac;
}  // end function fac

// The double factorial  (2n+1)!!
double dfactorial (int m)
{
  // make sure that if ( m%2 /= 1) then this function is not called!
  double dfac = 0.0;
  for( int i = 3; i <= m; i+=2) dfac += log((float) i);
  return dfac;
} // end function dfac

/*
  The function  clebsch_gordan()                   
  returns the value of the Clebsch-Gordan coefficient                                  
  <j1/2, m1/2, j2/2, m2/2 | j3/2, (m1 + m2)/2> 
*/

double clebsch_gordan(int j1, int j2, int j3, int m1, int m2)
{
  /*   fac[n] = n! / sqrt{10^n}  */
   
  static double fac[150] = {
    .1000000000000000E+01,.3162277660168379E-01,.2000000000000000E-02,
    .1897366596101028E-03,.2400000000000000E-04,.3794733192202056E-05,
    .7200000000000002E-06,.1593787940724864E-06,.4032000000000001E-07,
    .1147527317321902E-07,.3628800000000002E-08,.1262280049054092E-08,
    .4790016000000003E-09,.1969156876524384E-09,.8717829120000006E-10,
    .4135229440701207E-10,.2092278988800002E-10,.1124782407870728E-10,
    .6402373705728006E-11,.3846755834917892E-11,.2432902008176642E-11,
    .1615637450665515E-11,.1124000727777609E-11,.8175125500367506E-12,
    .6204484017332402E-12,.4905075300220504E-12,.4032914611266062E-12,
    .3443362860754794E-12,.3048883446117143E-12,.2796010642932893E-12,
    .2652528598121915E-12,.2600289897927591E-12,.2631308369336940E-12,
    .2745906132211536E-12,.2952327990396046E-12,.3267628297331728E-12,
    .3719933267899019E-12,.4352480892045862E-12,.5230226174666021E-12,
    .6450376682011969E-12,.8159152832478993E-12,.1057861775849963E-11,
    .1405006117752883E-11,.1910498367185033E-11,.2658271574788455E-11,
    .3782786767026367E-11,.5502622159812102E-11,.8178384990311005E-11,
    .1241391559253610E-10,.1923556149721149E-10,.3041409320171346E-10,
    .4905068181788930E-10,.8065817517094409E-10,.1351836790901029E-09,
    .2308436973392420E-09,.4014955268976057E-09,.7109985878048654E-09,
    .1281573721857157E-08,.2350561331282885E-08,.4385545276195193E-08,
    .8320987112741415E-08,.1605109571087441E-07,.3146997326038803E-07,
    .6269557984667544E-07,.1268869321858846E-06,.2608136121621699E-06,
    .5443449390774448E-06,.1153317792981115E-05,.2480035542436839E-05,
    .5411367084667394E-05,.1197857166996993E-04,.2689449441079695E-04,
    .6123445837688631E-04,.1413574626231488E-03,.3307885441519399E-03,
    .7845339175584759E-03,.1885494701666058E-02,.4591092485552201E-02,
    .1132428117820634E-01,.2829031189597267E-01,.7156945704626409E-01,
    .1833212210859029E+00,.4753643337012861E+00,.1247684230710655E+01,
    .3314240134565367E+01,.8908465407274079E+01,.2422709538367284E+02,
    .6665313817722467E+02,.1854826422573992E+03,.5220273782040236E+03,
    .1485715964481768E+04,.4275404227490954E+04,.1243841405464136E+05,
    .3658035857041260E+05,.1087366156656748E+06,.3266626020337846E+06,
    .9916779348709544E+06,.3041882150138602E+07,.9426890448883294E+07,
    .2951234062064472E+08,.9332621544394462E+08,.2980746402685117E+09,
    .9614466715035176E+09,.3131572170660985E+10,.1029901674514568E+11,
    .3419676810361796E+11,.1146280563734714E+12,.3878597438312349E+12,
    .1324641819451836E+13,.4565884904381298E+13,.1588245541522752E+14,
    .5574945468249565E+14,.1974506857221085E+15,.7055650984616650E+15,
    .2543559733472202E+16,.9249958440832429E+16,.3393108684451918E+17,
    .1255404359589777E+18,.4684525849754318E+18,.1762838801735966E+19,
    .6689502913449167E+19,.2559641940120622E+20,.9875044200833661E+20,
    .3840998695345006E+21,.1506141741511150E+22,.5953547977784760E+22,
    .2372173242880062E+23,.9526867474051174E+23,.3856204823625829E+24,
    .1573076357315330E+25,.6466855489220516E+25,.2678949036508007E+26,
    .1118248651196012E+27,.4703162928493458E+27,.1992942746161532E+28,
    .8508021737644667E+28,.3659042881952573E+29,.1585214610157954E+30,
    .6917786472619537E+30,.3040758665204989E+31,.1346201247571762E+32,
    .6002457605114648E+32,.2695364137888182E+33,.1218859041294581E+34,
    .5550293832739345E+34,.2544977678223085E+35,.1174997204390920E+36,
    .5462031093002385E+36,.2556323917872885E+37,.1204487096628886E+38};

  double value;       /* The return value */
  int temp, change;   /* interchange parameters */
  int jz1, jz2, jz3,
    jt1, jt2, jt3,
    jt4, jt5, j4, j5; /* ang. mom. relations */

  /* Local variables */

  int     loop_min, loop_max, phase, loop;
  int FAC_LIM = 300;
  double                           factor;

  value = 0.0;              /* initialization */

  /*
  ** Interchange the angular momenta such that the smallest j-value 
  ** is placed in position 2. A parameter is set to remember the interchange 
  **    change  = -1 - j1 is placed in position 2.  
  **            =  0 - no change                       
  **            =  1 - j3 is placed in position 2.  
  **/

  if(   (j1 < j2 )  && (j1 < j3)) {

    /* interchange 1 and 2 */

    temp = j1; j1 = j2; j2 = temp;
    temp = m1; m1 = m2; m2 = temp;
    change = -1;
  }
  else if(   (j3 < j1) && (j3 < j2))  {

    /* interchange 2 and 3 */

    temp = j2; j2 = j3; j3 = temp;
    m2 = -(m1 + m2);
    change = 1;
  }
  else change = 0;

  /* Test of angular momentum relations  */

  jz1 = (j1 + j2 - j3)/2;
  jz2 = (j1 + j3 - j2)/2;
  jz3 = (j2 + j3 - j1)/2;

  if((jz1 < 0) || (jz2 < 0) || (jz3 < 0))  return  value ;

  /* Test of angular projection relations */

  if((j1 < sqrt(m1*m1)) || (j2 < sqrt(m2*m2)) || (j3 < sqrt((m1 + m2)*(m1+m2))))  return  value;

  /* Definition of loop parameters */

  jt1 = (j1 - j3 + m2)/2;
  jt2 = (j2 - j3 - m1)/2;
  jt3 = (j1 - m1)/2;
  jt4 = (j2 + m2)/2;
  jt5 = (j2 - m2)/2;
  j4  = j1/2;
  j5  = j3/2;

  /* Loop limits */

  loop_min = MAX( MAX(jt1, jt2) , 0);
  loop_max = MIN( MIN(jt3, jt4) , jz1);

  /* Loop test */

  if( loop_min > loop_max) return  value;

  /*  test for maximum factorials   */

  if ((jz1 > FAC_LIM) || (jz3 > FAC_LIM) || (jt4 > FAC_LIM) || (jt5 > FAC_LIM)) {
    cout << "Termination from VCC. - too large factorials" << endl;
  }
  for (loop = loop_min,phase = pow(-1.0,loop_min); loop <= loop_max;
       phase *= -1, loop++) {
    value +=  (fac_ratio(jt3-loop,j4) / fac [jt4-loop] )
      * (fac_ratio(loop-jt2,j5) / fac[loop-jt1] )
      * ((phase / fac[jz1-loop]) / fac[loop] );
  }
  factor =  fac_ratio(j4,(j1 + m1)/2) * fac_ratio(j4,jt3)
    * fac_ratio((j1 + j2 + j3)/2+1,jz2)
    * fac_ratio(j5,(j3 + m1 + m2)/2)
    * fac_ratio(j5,(j3 - m1 - m2)/2)
    * fac[jt4] * fac[jt5] * fac[jz1] * fac[jz3]
    * (j3 + 1);
  value *= sqrt(factor);

  /* Add the angular momentum interchange factor */

  if(change == -1) value *= pow(-1.0,jz1);
  else if (change == 1)
    value *= pow(-1.0,jt3) * sqrt(  (j2+1.0)
				 / (j3+1.0));
  return (value);
} /* End: function clebsch_gordan() */

/*
  The function               
  double fac_ratio()     
  calculates and returns the ratio (n! / m!).           
*/

double fac_ratio(int m, int n)
{
  int i, high;
  double value;

  value = 1.0;      /* initialization */
  if( n == m)    return (value);          /* Return for n equal m */
  high = MAX(n,m); 
  for (i = MIN(n,m)+1; i <= high; i++) value*= i;
  if (n < m) return (1.0/value);       /* Return for n less  m  */
  return (value);                  /* Return for n greater than m */
} /* End: function fac_ratio() */