u[3]=2.0*u[2]-u[1];
            u[4]=2.0*u[3]-u[2];
          }
      }
    else
      { if (kk<=1)
          { u[3]=(y[kk+2]-y[kk+1])/h;
            if (kk==1)
              { u[1]=(y[1]-y[0])/h;
                u[0]=2.0*u[1]-u[2];
                if (n==4) u[4]=2.0*u[3]-u[2];
                else u[4]=(y[4]-y[3])/h;
              }
            else
              { u[1]=2.0*u[2]-u[3];
                u[0]=2.0*u[1]-u[2];
                u[4]=(y[3]-y[2])/h;
              }
          }
        else if (kk>=(n-3))
          { u[1]=(y[kk]-y[kk-1])/h;
            if (kk==(n-3))
              { u[3]=(y[n-1]-y[n-2])/h;
                u[4]=2.0*u[3]-u[2];
                if (n==4) u[0]=2.0*u[1]-u[2];
                else u[0]=(y[kk-1]-y[kk-2])/h;
              }
            else
              { u[3]=2.0*u[2]-u[1];
                u[4]=2.0*u[3]-u[2];
                u[0]=(y[kk-1]-y[kk-2])/h;
              }
          }
        else
          { u[1]=(y[kk]-y[kk-1])/h;
            u[0]=(y[kk-1]-y[kk-2])/h;
            u[3]=(y[kk+2]-y[kk+1])/h;
            u[4]=(y[kk+3]-y[kk+2])/h;
          }
      }
    s[0]=fabs(u[3]-u[2]);
    s[1]=fabs(u[0]-u[1]);
    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
         p=(u[1]+u[2])/2.0;
    else p=(s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);
    s[0]=fabs(u[3]-u[4]);
    s[1]=fabs(u[2]-u[1]);
    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
         q=(u[2]+u[3])/2.0;
    else q=(s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);
    s[0]=y[kk];
    s[1]=p;
    s[3]=h;
    s[2]=(3.0*u[2]-2.0*p-q)/s[3];
    s[3]=(q+p-2.0*u[2])/(s[3]*s[3]);
    if (k<0)
      { p=t-(x0+kk*h);
        s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;
      }
    return;
}
/////////////////////////////////////////////////////////////
void enspl(double x[],double y[],int n,int k,double t,double s[5])
{ 
	int kk,m,l;
    double u[5],p,q;
    s[4]=0.0; s[0]=0.0; s[1]=0.0; s[2]=0.0; s[3]=0.0;
    if (n<1) return;
    if (n==1) { s[0]=y[0]; s[4]=y[0]; return;}
    if (n==2)
      { s[0]=y[0]; s[1]=(y[1]-y[0])/(x[1]-x[0]);
        if (k<0)
          s[4]=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
        return;
      }
    if (k<0)
      { if (t<=x[1]) kk=0;
        else if (t>=x[n-1]) kk=n-2;
        else
          { kk=1; m=n;
            while (((kk-m)!=1)&&((kk-m)!=-1))
              { l=(kk+m)/2;
                if (t<x[l-1]) m=l;
                else kk=l;
              }
            kk=kk-1;
          }
      }
    else kk=k;
    if (kk>=n-1) kk=n-2;
    u[2]=(y[kk+1]-y[kk])/(x[kk+1]-x[kk]);
    if (n==3)
      { if (kk==0)
          { u[3]=(y[2]-y[1])/(x[2]-x[1]);
            u[4]=2.0*u[3]-u[2];
            u[1]=2.0*u[2]-u[3];
            u[0]=2.0*u[1]-u[2];
          }
        else
          { u[1]=(y[1]-y[0])/(x[1]-x[0]);
            u[0]=2.0*u[1]-u[2];
            u[3]=2.0*u[2]-u[1];
            u[4]=2.0*u[3]-u[2];
          }
      }
    else
      { if (kk<=1)
          { u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
            if (kk==1)
              { u[1]=(y[1]-y[0])/(x[1]-x[0]);
                u[0]=2.0*u[1]-u[2];
                if (n==4) u[4]=2.0*u[3]-u[2];
                else u[4]=(y[4]-y[3])/(x[4]-x[3]);
              }
            else
              { u[1]=2.0*u[2]-u[3];
                u[0]=2.0*u[1]-u[2];
                u[4]=(y[3]-y[2])/(x[3]-x[2]);
              }
          }
        else if (kk>=(n-3))
          { u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
            if (kk==(n-3))
              { u[3]=(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
                u[4]=2.0*u[3]-u[2];
                if (n==4) u[0]=2.0*u[1]-u[2];
                else u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
              }
            else
              { u[3]=2.0*u[2]-u[1];
                u[4]=2.0*u[3]-u[2];
                u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
              }
          }
        else
          { u[1]=(y[kk]-y[kk-1])/(x[kk]-x[kk-1]);
            u[0]=(y[kk-1]-y[kk-2])/(x[kk-1]-x[kk-2]);
            u[3]=(y[kk+2]-y[kk+1])/(x[kk+2]-x[kk+1]);
            u[4]=(y[kk+3]-y[kk+2])/(x[kk+3]-x[kk+2]);
          }
      }
    s[0]=fabs(u[3]-u[2]);
    s[1]=fabs(u[0]-u[1]);
    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
         p=(u[1]+u[2])/2.0;
    else p=(s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);
    s[0]=fabs(u[3]-u[4]);
    s[1]=fabs(u[2]-u[1]);
    if ((s[0]+1.0==1.0)&&(s[1]+1.0==1.0))
         q=(u[2]+u[3])/2.0;
    else q=(s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);
    s[0]=y[kk];
    s[1]=p;
    s[3]=x[kk+1]-x[kk];
    s[2]=(3.0*u[2]-2.0*p-q)/s[3];
    s[3]=(q+p-2.0*u[2])/(s[3]*s[3]);
    if (k<0)
      { p=t-x[kk];
        s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;
      }
    return;
}
/////////////////////////////////////////////////////////////
double espl1(double x[],double y[],int n,double dy[],double ddy[],double t[],int m,double z[],double dz[],double ddz[])
{ 
	int i,j;
    double h0,h1,alpha,beta,g,*s;
    s=(double*)malloc(n*sizeof(double));
    s[0]=dy[0]; dy[0]=0.0;
    h0=x[1]-x[0];
    for (j=1;j<=n-2;j++)
      { h1=x[j+1]-x[j];
        alpha=h0/(h0+h1);
        beta=(1.0-alpha)*(y[j]-y[j-1])/h0;
        beta=3.0*(beta+alpha*(y[j+1]-y[j])/h1);
        dy[j]=-alpha/(2.0+(1.0-alpha)*dy[j-1]);
        s[j]=(beta-(1.0-alpha)*s[j-1]);
        s[j]=s[j]/(2.0+(1.0-alpha)*dy[j-1]);
        h0=h1;
      }
    for (j=n-2;j>=0;j--)
      dy[j]=dy[j]*dy[j+1]+s[j];
    for (j=0;j<=n-2;j++) s[j]=x[j+1]-x[j];
    for (j=0;j<=n-2;j++)
      { h1=s[j]*s[j];
        ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
      }
    h1=s[n-2]*s[n-2];
    ddy[n-1]=6.*(y[n-2]-y[n-1])/h1+2.*(2.*dy[n-1]+dy[n-2])/s[n-2];
    g=0.0;
    for (i=0;i<=n-2;i++)
      { h1=0.5*s[i]*(y[i]+y[i+1]);
        h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
        g=g+h1;
      }
    for (j=0;j<=m-1;j++)
      { if (t[j]>=x[n-1]) i=n-2;
        else
          { i=0;
            while (t[j]>x[i+1]) i=i+1;
          }
        h1=(x[i+1]-t[j])/s[i];
        h0=h1*h1;
        z[j]=(3.0*h0-2.0*h0*h1)*y[i];
        z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
        dz[j]=6.0*(h0-h1)*y[i]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
        ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
        ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
        h1=(t[j]-x[i])/s[i];
        h0=h1*h1;
        z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
        z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
        dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
        ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s[i]*s[i]);
        ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
      }
    free(s);
    return(g);
}
/////////////////////////////////////////////////////////////
double espl2(double x[],double y[],int n,double dy[],double ddy[],double t[],int m,double z[],double dz[],double ddz[])
{ 
	int i,j;
    double h0,h1,alpha,beta,g,*s;
    s=(double*)malloc(n*sizeof(double));
    dy[0]=-0.5;
    h0=x[1]-x[0];
    s[0]=3.0*(y[1]-y[0])/(2.0*h0)-ddy[0]*h0/4.0;
    for (j=1;j<=n-2;j++)
      { h1=x[j+1]-x[j];
        alpha=h0/(h0+h1);
        beta=(1.0-alpha)*(y[j]-y[j-1])/h0;
        beta=3.0*(beta+alpha*(y[j+1]-y[j])/h1);
        dy[j]=-alpha/(2.0+(1.0-alpha)*dy[j-1]);
        s[j]=(beta-(1.0-alpha)*s[j-1]);
        s[j]=s[j]/(2.0+(1.0-alpha)*dy[j-1]);
        h0=h1;
      }
    dy[n-1]=(3.0*(y[n-1]-y[n-2])/h1+ddy[n-1]*h1/
            2.0-s[n-2])/(2.0+dy[n-2]);
    for (j=n-2;j>=0;j--)
      dy[j]=dy[j]*dy[j+1]+s[j];
    for (j=0;j<=n-2;j++) s[j]=x[j+1]-x[j];
    for (j=0;j<=n-2;j++)
      { h1=s[j]*s[j];
        ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
      }
    h1=s[n-2]*s[n-2];
    ddy[n-1]=6.*(y[n-2]-y[n-1])/h1+2.*(2.*dy[n-1]+dy[n-2])/s[n-2];
    g=0.0;
    for (i=0;i<=n-2;i++)
      { h1=0.5*s[i]*(y[i]+y[i+1]);
        h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
        g=g+h1;
      }
    for (j=0;j<=m-1;j++)
      { if (t[j]>=x[n-1]) i=n-2;
        else
          { i=0;
            while (t[j]>x[i+1]) i=i+1;
          }
        h1=(x[i+1]-t[j])/s[i];
        h0=h1*h1;
        z[j]=(3.0*h0-2.0*h0*h1)*y[i];
        z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
        dz[j]=6.0*(h0-h1)*y[i]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
        ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
        ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
        h1=(t[j]-x[i])/s[i];
        h0=h1*h1;
        z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
        z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
        dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
        ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s[i]*s[i]);
        ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
      }
    free(s);
    return(g);
}
/////////////////////////////////////////////////////////////
double espl3(double x[],double y[],int n,double dy[],double ddy[],double t[],int m,double z[],double dz[],double ddz[])
{ 
	int i,j;
    double h0,y0,h1,y1,alpha,beta,u,g,*s;
    s=(double*)malloc(n*sizeof(double));
    h0=x[n-1]-x[n-2];
    y0=y[n-1]-y[n-2];
    dy[0]=0.0; ddy[0]=0.0; ddy[n-1]=0.0;
    s[0]=1.0; s[n-1]=1.0;
    for (j=1;j<=n-1;j++)
      { h1=h0; y1=y0;
        h0=x[j]-x[j-1];
        y0=y[j]-y[j-1];
        alpha=h1/(h1+h0);
        beta=3.0*((1.0-alpha)*y1/h1+alpha*y0/h0);
        if (j<n-1)
          { u=2.0+(1.0-alpha)*dy[j-1];
            dy[j]=-alpha/u;
            s[j]=(alpha-1.0)*s[j-1]/u;
            ddy[j]=(beta-(1.0-alpha)*ddy[j-1])/u;
          }
      }
    for (j=n-2;j>=1;j--)
      { s[j]=dy[j]*s[j+1]+s[j];
        ddy[j]=dy[j]*ddy[j+1]+ddy[j];
      }
    dy[n-2]=(beta-alpha*ddy[1]-(1.0-alpha)*ddy[n-2])/
            (alpha*s[1]+(1.0-alpha)*s[n-2]+2.0);
    for (j=2;j<=n-1;j++)
        dy[j-2]=s[j-1]*dy[n-2]+ddy[j-1];
    dy[n-1]=dy[0];
    for (j=0;j<=n-2;j++) s[j]=x[j+1]-x[j];
    for (j=0;j<=n-2;j++)
      { h1=s[j]*s[j];
        ddy[j]=6.0*(y[j+1]-y[j])/h1-2.0*(2.0*dy[j]+dy[j+1])/s[j];
      }
    h1=s[n-2]*s[n-2];
    ddy[n-1]=6.*(y[n-2]-y[n-1])/h1+2.*(2.*dy[n-1]+dy[n-2])/s[n-2];
    g=0.0;
    for (i=0;i<=n-2;i++)
      { h1=0.5*s[i]*(y[i]+y[i+1]);
        h1=h1-s[i]*s[i]*s[i]*(ddy[i]+ddy[i+1])/24.0;
        g=g+h1;
      }
    for (j=0;j<=m-1;j++)
      { h0=t[j];
        while (h0>=x[n-1]) h0=h0-(x[n-1]-x[0]);
        while (h0<x[0]) h0=h0+(x[n-1]-x[0]);
        i=0;
        while (h0>x[i+1]) i=i+1;
        u=h0;
        h1=(x[i+1]-u)/s[i];
        h0=h1*h1;
        z[j]=(3.0*h0-2.0*h0*h1)*y[i];
        z[j]=z[j]+s[i]*(h0-h0*h1)*dy[i];
        dz[j]=6.0*(h0-h1)*y[i]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i];
        ddz[j]=(6.0-12.0*h1)*y[i]/(s[i]*s[i]);
        ddz[j]=ddz[j]+(2.0-6.0*h1)*dy[i]/s[i];
        h1=(u-x[i])/s[i];
        h0=h1*h1;
        z[j]=z[j]+(3.0*h0-2.0*h0*h1)*y[i+1];
        z[j]=z[j]-s[i]*(h0-h0*h1)*dy[i+1];
        dz[j]=dz[j]-6.0*(h0-h1)*y[i+1]/s[i];
        dz[j]=dz[j]+(3.0*h0-2.0*h1)*dy[i+1];
        ddz[j]=ddz[j]+(6.0-12.0*h1)*y[i+1]/(s[i]*s[i]);
        ddz[j]=ddz[j]-(2.0-6.0*h1)*dy[i+1]/s[i];
      }
    free(s);
    return(g);
}
/////////////////////////////////////////////////////////////
double eslq3(double x[],double y[],double z[],int n,int m,double u,double v)
{ 
	int nn,mm,ip,iq,i,j,k,l;
    double b[3],h,w;
    nn=3;
    if (n<=3) { ip=0;  nn=n;}
    else if (u<=x[1]) ip=0;
    else if (u>=x[n-2]) ip=n-3;
    else
      { i=1; j=n;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (u<x[l-1]) j=l;
            else i=l;
          }
        if (fabs(u-x[i-1])<fabs(u-x[j-1])) ip=i-2;
        else ip=i-1;
      }
    mm=3;
    if (m<=3) { iq=0; mm=m;}
    else if (v<=y[1]) iq=0;
    else if (v>=y[m-2]) iq=m-3;
    else
      { i=1; j=m;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (v<y[l-1]) j=l;
            else i=l;
          }
        if (fabs(v-y[i-1])<fabs(v-y[j-1])) iq=i-2;
        else iq=i-1;
      }
    for (i=0;i<=nn-1;i++)
      { b[i]=0.0;
        for (j=0;j<=mm-1;j++)
          { k=m*(ip+i)+(iq+j);
            h=z[k];
            for (k=0;k<=mm-1;k++)
              if (k!=j)
                h=h*(v-y[iq+k])/(y[iq+j]-y[iq+k]);
            b[i]=b[i]+h;
          }
      }
    w=0.0;
    for (i=0;i<=nn-1;i++)
      { h=b[i];
        for (j=0;j<=nn-1;j++)
          if (j!=i)
            h=h*(u-x[ip+j])/(x[ip+i]-x[ip+j]);
        w=w+h;
      }
    return(w);
}
/////////////////////////////////////////////////////////////
double eslgq(double x[],double y[],double z[],int n,int m,double u,double v)
{ 
	int ip,ipp,i,j,l,iq,iqq,k;
    double h,w,b[10];
    if (u<=x[0]) { ip=1; ipp=4;}
    else if (u>=x[n-1]) { ip=n-3; ipp=n;}
    else
      { i=1; j=n;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (u<x[l-1]) j=l;
            else i=l;
          }
        ip=i-3; ipp=i+4;
      }
    if (ip<1) ip=1;
    if (ipp>n) ipp=n;
    if (v<=y[0]) { iq=1; iqq=4;}
    else if (v>=y[m-1]) { iq=m-3; iqq=m;}
    else
      { i=1; j=m;
        while (((i-j)!=1)&&((i-j)!=-1))
          { l=(i+j)/2;
            if (v<y[l-1]) j=l;
            else i=l;
          }
        iq=i-3; iqq=i+4;
      }
    if (iq<1) iq=1;
    if (iqq>m) iqq=m;
    for (i=ip-1;i<=ipp-1;i++)
      { b[i-ip+1]=0.0;
        for (j=iq-1;j<=iqq-1;j++)
          { h=z[m*i+j];
            for (k=iq-1;k<=iqq-1;k++)
              if (k!=j) h=h*(v-y[k])/(y[j]-y[k]);
            b[i-ip+1]=b[i-ip+1]+h;
          }
      }
    w=0.0;
    for (i=ip-1;i<=ipp-1;i++)
      { h=b[i-ip+1];
        for (j=ip-1;j<=ipp-1;j++)
          if (j!=i) h=h*(u-x[j])/(x[i]-x[j]);
        w=w+h;
      }
    return(w);
}
/////////////////////////////////////////////////////////////
#endif