# include <math.h>
extern "C" {
# include "matrix.h"
# include "matrix2.h"
}
# include "s3glob.inc"
# include "s3mo.inc"
double  w_1(double x) 
{
return (x-2.0*x*x+x*x*x);
}
 
/* */
/* second spline:*/
double  w_2(double x) 
{
double dtmp;
double t1;
t1 = x*x;
dtmp = 3.0*t1-2.0*t1*x;
return (dtmp);
}
 
/*            */
/* the first spline:         */
double  coe_spline1(double x) 
  {
  double t1,t2,t4;
  double dtmp;
  t1 = pow(x,0.3333333333333333);
  t2 = t1*t1;
  t4 = x*x;
  dtmp = 0.1666666667/t2*(-16.0*x+33.0*t4-5.0)*sinc(-(-log(-x/(-1+x))+j*h)/h)*(
-1.0+x);
  return (dtmp);
  }
 
/*            */
/* the main sinc sum:         */
double  coe_sincs(double x) 
  {
  double t4,t7,t10,t11,t12,t17,t20,t23,t24,t26;
  double dtmp;
  t4 = log(-x/(-1+x));
  t7 = 1/h;
  t10 = pow(x,0.3333333333333333);
  t11 = t10*t10;
  t12 = t11*t11;
  t17 = (-t4+k*h)*t7;
  t20 = sinc_p(-t17)*h;
  t23 = h*h;
  t24 = sinc(-t17)*t23;
  t26 = x*x;
  dtmp = 0.1666666667*sinc(-(-t4+j*h)*t7)/t12*(-6.0*sinc_pp(-t17)+5.0*t20*x+t20
+7.0*t24-9.0*t24*x+2.0*t24*t26)/t23;
  return (dtmp);
  }
 
/*            */
/* the second spline:         */
double  coe_spline2(double x) 
  {
  double t12;
  double dtmp;
  t12 = pow(x,0.3333333333333333);
  dtmp = -(-4.0+11.0*x)*sinc(-(-log(-x/(-1+x))+j*h)/h)*t12*(-1.0+x);
  return (dtmp);
  }
 
/*            */
/*            */
double  psi(double x) 
  {
  double t1;
  double dtmp;
  t1 = exp(x);
  dtmp = t1/(t1+1.0);
  return (dtmp);
  }
 
/*            */
/*            */
double  Sphi(int j, double h, double x) 
  {
  
  double dtmp;
  dtmp = sinc((log(x/(1-x))-j*h)/h);
  return (dtmp);
  }
 
/*            */
/*            */
double  u(double x) 
{
  double dtmp;
{ 
double t2;
t2 = pow(-1.0+x,2.0);
dtmp = x*t2*u_vec->ve[0];
} 
{ 
double t1,t3;
t1 = x*x;
t3 = t1*x;
dtmp += a-3.0*a*t1+2.0*a*t3-b*t1+b*t3;
} 
{ 
double t1;
t1 = x*x;
dtmp += 3.0*t1*u_vec->ve[m+1]-2.0*t1*x*u_vec->ve[m+1];
} 
for(k=-M;k<=N;k++)
  {
  double t10,t11,t13,t14;
  t10 = u_vec->ve[k+M+1]*sinc(-(-log(-x/(-1+x))+k*h)/h);
  t11 = pow(x,0.3333333333333333);
  t13 = t11*t11;
  t14 = t13*t13;
  dtmp += t10*t11-t10*t14;
  }
return (dtmp);
}
 
double  u_p(double x) 
{
  double dtmp;
{ 
dtmp = (-1.0+3.0*x)*(-1.0+x)*(u_vec->ve[0]);
} 
{ 
double t2;
t2 = x*x;
dtmp += -6.0*a*x+6.0*a*t2-2.0*b*x+3.0*b*t2;
} 
{ 
double t2;
t2 = x*x;
dtmp += 6.0*x*(u_vec->ve[m+1])-6.0*t2*(u_vec->ve[m+1]);
} 
for(k=-M;k<=N;k++)
  {
  double t1,t3,t4,t5,t12,t17;
  t1 = 1/h;
  t3 = pow(x,0.3333333333333333);
  t4 = t3*t3;
  t5 = 1/t4;
  t12 = (-log(-x/(-1+x))+k*h)*t1;
  t17 = sinc(-t12);
  dtmp += 0.9999999999*(u_vec->ve[k+M+1])*t1*t5*sinc_p(-t12)+0.3333333333*
(u_vec->ve[k+M+1])*t5*t17-0.13333333332E1*(u_vec->ve[k+M+1])*t3*t17;
  }
return (dtmp);
}
 
double  coe_rhs(double x) 
  {
  double t1,t2,t3,t4,t5,t7,t8,t21,t49;
  double dtmp;
  t1 = x*x;
  t2 = t1*x;
  t3 = t1*t1;
  t4 = sqrt(x);
  t5 = t4*t1;
  t7 = 1.0-x;
  t8 = sqrt(t7);
  t21 = b-a;
  t49 = pow(x,0.3333333333333333);
  dtmp = (-0.8333333333E-1*(-14.0*t1+7.0*t2+16.0*t3+21.0*t5*a-10.0*t5*t8*b-12.0
*t4*x*a)*t4/t3/t8+a*(-8.0+18.0*x)+t21*(-2.0+6.0*x)+0.1666666667/x*(a*(-8.0*x+
9.0*t1)+t21*(-2.0*x+3.0*t1))-1/t1*(a*(1.0-4.0*t1+3.0*t2)+t21*(-t1+t2)))*sinc((
log(x/(1-x))-j*h)/h)*t49*t7;
  return (dtmp);
  }