/* to evaluate Eq (2.43) of Li's thesis (generalized to coupling) :
 *
 * xi[N][M]=(-1)^(M+N) ( l1  l2   s )  / (l1 l2 s)
 *                     ( -N  -M  M+N) /  ( 0  0 0)
 *
 *         =(-1)^(M+N) ( l1  s   l2 )  / (l1 s l2)
 *                     ( -N  M+N -M ) /  ( 0 0  0)
 *
 *
 *  for all M : -2<=M<=2
 *  and all N : -2<=N<=2
 */
 
float **eq2_43(l1,l2,s)
int l1,l2,s;
{
static float **xi;
float **farray2();
int dim1,dim2,M,N,d0;
float sign;
double *work,*darray1();
static first=1;
  if(first)
  { first=0;
    xi=farray2(-2,2,-2,2);
  }
  
  dim1=2*l1+1;
  dim2=2*l2+1;
  work=darray1(0,dim1*dim2-1);
  w3j_(&l1,&s,&l2,work,&dim1);
  d0=l1+dim1*l2;
  
  for(M=-2;M<=2;M++)
  for(N=-2;N<=2;N++) 
    xi[M][N]=0.0;
  
  sign=1.0;
  for(M=-2;M<=2;M++)
  { if(M>=-l2 && M<=l2)
    { for(N=-2;N<=2;N++)
      { if(N>=-l1&&N<=l1) 
         xi[N][M]=sign*work[N-dim1*M+d0]/work[d0];
      }
    }
    sign=-sign;
  }
  
  free_darray1(work,0);
   
  return xi;
}