/*************************    Program 5.5    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define NMAX 30
#define NTEL  5

void bssl (bj,by,n,x)
/* Subroutine to generate J_n(x) and Y_n(x) with given x and
   up to N=NMAX-NTEL.  Copyright (c) Tao Pang 1997. */
int n;
double x;
double bj[],by[];
{
int i;
double pi,sum,sum1,gamma;
double b1[NMAX+1];

if (n > (NMAX-NTEL))
  {
  printf("N is too large.\n");
  exit(1);
  }
pi = 4*atan(1);
gamma = 0.5772156649;

b1[NMAX] = 0;
b1[NMAX-1] = 1;

sum = 0;
for (i = NMAX-1; i > 0; i = i-1) 
  {
  b1[i-1] = 2*i*b1[i]/x-b1[i+1];
  if (i%2 == 0) sum = sum+2*b1[i];
  }
sum = sum+b1[0];

for (i = 0; i <= n; ++i) 
  {
  bj[i] = b1[i]/sum;
  }

/* Generating Y_n(x) starts here */

sum1 = 0;
for (i = 1; i <= NMAX/2; ++i) 
  {
  sum1 = sum1+pow(-1,i)*b1[2*i]/i;
  }
sum1 = -4*sum1/(pi*sum);
by[0] = 2*(log(x/2)+gamma)*bj[0]/pi+sum1;
by[1] = (bj[1]*by[0]-2/(pi*x))/bj[0];

for (i = 1; i <n-1; ++i) 
  {
  by[i+1] = 2*i*by[i]/x-by[i-1];
  }
}