/* This file contains the functions  */
/* that are needed to describe the angular motion in 3 dimensions */
/* and the symbols to couple angular momenta. */
/* Contents : associated Legendre functions, spherical harmonics, */
/* Clebsch-Gordon coefficients, three-j and six-j symbols.  */
/* The Condon-Shortley phase convention is used. */
/* For definitions, see The nuclear shell model, K.Heyde, Springer 1990. */




#include <angular.h>


/* The function plgndr(l,m,x) returns the associated Legendre function 
   P(lower l, upper m) with l and m integers satisfying 0 <= m <= l , 
   while x lies in the range -1 <= x <= 1 .
   For definition and implementation : see 
   Press et al., Numerical recipes in C, Cambridge 1992. 
   !!! plgndr(l,m,x) has a DIFFERENT phase as in Joachain !!! */

double plgndr(int l, int m, double x)
{
     double fact, pll, pmm, pmmp1, somx2;
     int i, ll;

     if(l < 0)
        return 0.0;
        
     else
     {   
     pmm = 1.0;
     if (m > 0)
     {
	somx2 = sqrt((1.0 - x) * (1.0 + x));
	fact = 1.0;
	for (i = 1; i <= m; i++)
	{
	   pmm *= - fact*somx2;
	   fact += 2.0;
	}
     }
     if (l == m)
	return pmm;
     else
     {
	pmmp1 = x * (2.0 * m + 1.) * pmm;
	if (l == (m + 1))
	   return pmmp1;
	else
	{
	   for (ll = m + 2; ll <= l; ll++)
	   {
	      pll = (x * (2. * ll - 1.) * pmmp1 - (ll + m - 1.) * pmm)/(ll - m);
	      pmm = pmmp1;
	      pmmp1 = pll;
	   }
	   return pll;
	}
     }
     }
}


/* The function spherharm(l,m,th,phi) returns the spherical harmonic */
/* Y(l, m) ( th, phi) with l and m satisfying -l <= m <= l , */
/* while th is the polar and phi the azimuthal angle */
/* with 0 <= th <= PI and 0 <= phi <= 2 PI. */
/* For definition and implementation : see */
/* Press et al., Numerical recipes in C, Cambridge 1992. */

dcomplex spherharm(int l, int m, double th, double phi)
{
     double prefactor;
     dcomplex yp;
     int mpos;

     mpos = abs(m);
     prefactor = sqrt( ( (2.0 * l + 1.) * fact_int( l - mpos) ) /
     ( 4.0 * PI * fact_int( l + mpos) ) );
     yp.re = prefactor * plgndr(l,mpos,cos(th)) * cos(mpos * phi);
     yp.im = prefactor * plgndr(l,mpos,cos(th)) * sin(mpos * phi);

     if (m >= 0)
	return yp;
     else
	return ( CRmul( phase(mpos) , Conjg(yp) ) );

}


/* the function phase returns the phase (-1)^(t) of an integer t */

double phase(int t)
{
  if (abs(t % 2) == 1) return -1.;
  else return 1.;
}


/* the function fact_int(n) returns the factorial of an integer n <= 32 */

double fact_int(int n)
{
     static int ntop = 4;
     static double a[33] = {1.0,1.0,2.0,6.0,24.0};
     int j;

     while (ntop < n)
     {
	j = ntop++;
	a[ntop] = a[j] * ntop;
     }
     return a[n];
}


/* The function fact_double() calculates the factorial of a double. /*
/* When working with half-valued angular momenta (0.5, 1.5, ...), we */
/* have to make sure that the sum of such angular momenta which */
/* appears as argument of the fact_int() function is an integer. */

double fact_double(double a)
{
     int n;

     if ( (a != floor(a)) || (a < 0.) )
     {
	  fprintf(stderr, "\nFactorial not determined");
	  exit(1);
     }
     else
     {
	n = (int)a;
	return fact_int(n);
     }
}


/* The function threej() returns the Wigner 3-j symbol. */

double threej( double j1, double j2, double j3,
	       double m1, double m2, double m3 )
{
     int t, mint, maxt;
     double sum = 0., factor1, factor2, min, max, temp;

     if ( (m1 + m2 + m3) != 0 )
	return 0.;
     if ( (j3 < fabs(j1 - j2)) || (j3 > (j1 + j2)) )
	return 0.;
     if ( (fabs(m1) > j1) || (fabs(m2) > j2) || (fabs(m3) > j3) )
	return 0.;

     factor1 = sqrt( fact_double(j1 + j2 - j3) *
		     fact_double(j2 + j3 - j1) *
		     fact_double(j3 + j1 - j2) /
		     fact_double(j1 + j2 + j3 + 1) );

     factor2 = sqrt( fact_double(j1 + m1) * fact_double(j1 - m1) *
		     fact_double(j2 + m2) * fact_double(j2 - m2) *
		     fact_double(j3 + m3) * fact_double(j3 - m3) );

     /* calculate the minimum value of t */

     min = 0.;
     temp = - j3 + j2 - m1;
     if (temp > min)
	min = temp;
     temp = - j3 + j1 + m2;
     if (temp > min)
	min = temp;
     mint = ceil(min);

     /* calculate the maximum value of t */

     max = j1 + j2 - j3;
     temp = j1 - m1;
     if ( temp < max )
	max = temp;
     temp = j2 + m2;
     if (temp < max)
	max = temp;
     maxt = floor(max);

     for(t = mint; t <= maxt; t++)
     {
	sum += phase( (int)(j1 - j2 - m3 + t) ) /
	   ( fact_double(t) * fact_double(j1 + j2 - j3 - t) *
	     fact_double(j3 - j2 + m1 + t) * fact_double(j3 - j1 - m2 + t) *
	     fact_double(j1 - m1 - t) * fact_double(j2 + m2 - t) );

     }

     return (factor1 * factor2 * sum );
}


/* The function clebshgordon() returns the Clebsch-Gordon coefficient */

double clebschgordon( double j1, double m1, double j2, double m2,
		      double j3, double m3)
{
     return ( phase( (int)(j1 - j2 + m3) ) * sqrt(2. * j3 + 1.) *
	      threej( j1, j2, j3, m1, m2, - m3) );

}


/* The function sixj() returns the 6-j symbol */

double sixj( double j1, double j2, double j3,
	     double l1, double l2, double l3)
{
     double mj1, mj2, mj3, ml1, ml2, ml3;
     double sum = 0.;

     if ( (j3 < fabs(j1 - j2)) || (j3 > (j1 + j2)) ||
	  (l3 < fabs(j1 - l2)) || (l3 > (j1 + l2)) ||
	  (j3 < fabs(l1 - l2)) || (j3 > (l1 + l2)) ||
	  (l3 < fabs(l1 - j2)) || (l3 > (l1 + j2)) )
	  return 0.;

     for(mj1 = - j1; mj1 <= j1; mj1++)
     {
	for(mj2 = - j2; mj2 <= j2; mj2++)
	{
	   for(mj3 = - j3; mj3 <= j3; mj3++)
	   {
	      for(ml1 = - l1; ml1 <= l1; ml1++)
	      {
		 for(ml2 = - l2; ml2 <= l2; ml2++)
		 {
		    for(ml3 = - l3; ml3 <= l3; ml3++)
		    {
		    sum += phase((int)(j1 + j2 + j3 + l1 + l2 + l3 +
				       mj1 + mj2 + mj3 + ml1 + ml2 + ml3)) *
			   threej(j1, j2, j3, mj1, mj2, mj3) *
			   threej(j1, l2, l3, - mj1, ml2, - ml3) *
			   threej(l1, j2, l3, - ml1, - mj2, ml3) *
			   threej(l1, l2, j3, ml1, - ml2, - mj3);
		    }
		 }
	      }
	   }
	}
     }
     return sum;
}