#include<stdio.h>
#include<math.h>

double b_laplace(double alpha, double s, double j);
double D_b_laplace();
double D2_b_laplace();
double D3_b_laplace();


/*************************************/
/* Laplace coefficients via summation 
   equation 6.68 of Murray and Dermott, page 237
   b_s^j (alpha) 
This program checked against expressions and values 
on page 262-263 in Murray and Dermott
*/
/*************************************/
#define ACC 1e-10  /* accuracy required for coefficient, size of last term */
double b_laplace(double alpha, double s, double j)
{
   double sum = 1.0;
   double fac = 1.0;
   double osum= 1.0;
   double diff = 100.0;
   double j_a;
   int i,k;

   j_a = fabs(j);

   if (alpha >1.0){
     printf("b_laplace: alpha > 1.0 error\n");
     return -1.0;
   }

/* outside the sum */
   for(k=0;k<j_a;k++) osum *= (s+k)/(k+1.0);
   osum *= pow(alpha,j_a);

   i = 0;
   while(diff>ACC){
     fac *= (s+i)*(s+j_a+i)/((i+1.0)*(j_a+i+1.0))*pow(alpha,2.0);
     sum += fac;
     diff = osum*fac;
     i++;
   }
   return sum*osum*2.0;
}




/*************************************/
/* Derivative of the laplace coeficient
   Db_s^j = d b_s^j/d alpha
   From equation 6.70 from Murray and Dermott page 237 */
/*************************************/

double D_b_laplace(double alpha, double s, double j)
{
  return s*(b_laplace(alpha, s+1.0, j-1.0)
 -2.0*alpha*b_laplace(alpha, s+1.0, j    )
       +    b_laplace(alpha, s+1.0, j+1.0));
}


/*************************************/
/* Second Derivative of the laplace coeficient
   D^2 b_s^j = d^2 b_s^j/d alpha^2
   From equation 6.71 from Murray and Dermott page 237  */
/*************************************/
double D2_b_laplace(double alpha, double s, double j)
{
  return s*(D_b_laplace(alpha, s+1.0, j-1.0)
 -2.0*alpha*D_b_laplace(alpha, s+1.0, j    )
       +    D_b_laplace(alpha, s+1.0, j+1.0)
         -2.0*b_laplace(alpha, s+1.0, j)
        );
}
  

/*************************************/
/* Third  Derivative of the laplace coeficient
   D^3 b_s^j = d^3 b_s^j/d alpha^3
   From equation 6.71 from Murray and Dermott page 237  */
/*************************************/
double D3_b_laplace(double alpha, double s, double j)
{
  return s*(D2_b_laplace(alpha, s+1.0, j-1.0)
 -2.0*alpha*D2_b_laplace(alpha, s+1.0, j    )
       +    D2_b_laplace(alpha, s+1.0, j+1.0)
        -2.0*2.0*D_b_laplace(alpha, s+1.0, j)
        );
}

double sign(double x)
{
  if (x>=0) return 1.0;
  else return -1.0;
}

/* all three routines checked with A5 on page 262 (M+D)*/
 
/* 
int main(int argc,char *argv[])
{
  double s,j,alpha,A5;
  j = 3.0;
  s = 0.5;
  while (1==1){
    printf("alpha=: ");
    scanf("%lf",&alpha);
    A5 = (21.0*b_laplace(alpha,s,j) 
           +10.0*alpha*D_b_laplace(alpha,s,j) 
           +alpha*alpha*D2_b_laplace(alpha,s,j) )/8.0;
    printf("%.5f\n",A5);
  }

}
*/