/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/***** #include "CH5.h"  插值 *****/
#ifndef CH5_H_
#define CH5_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"
//*******************************************************************
double enlgr(double x[],double y[],int n,double t);//一元全區間不等距插值
double eelgr(double x0,double h,int n,double y[],double t);//一元全區間等距插值
double enlg3(double x[],double y[],int n,double t);//一元三點不等距插值
double eelg3(double x0,double h,int n,double y[],double t);//一元三點等距插值
double enpqs(double x[],double y[],int n,double t);//連分式不等距插值
double eepqs(double x0,double h,int n,double y[],double t);//連分式等距插值
double enhmt(double x[],double y[],double dy[],int n,double t);//埃特金不等距插值
double eehmt(double x0,double h,int n,double y[],double dy[],double t);//埃特金等距插值
double enatk(double x[],double y[],int n,double t,double eps);//埃特金不等距逐步插值
double eeatk(double x0,double h,int n,double y[],double t,double eps);//埃特金等距逐步插值
void eespl(double x0,double h,int n,double y[],int k,double t,double s[5]);//光滑不等距插值
void enspl(double x[],double y[],int n,int k,double t,double s[5]);//光滑等距插值
double espl1(double x[],double y[],int n,double dy[],double ddy[],double t[],int m,double z[],double dz[],double ddz[]);//第一種邊界條件的三次樣條函數插值、微商與積分
double espl2(double x[],double y[],int n,double dy[],double ddy[],double t[],int m,double z[],double dz[],double ddz[]);//第二種邊界條件的三次樣條函數插值、微商與積分
double espl3(double x[],double y[],int n,double dy[],double ddy[],double t[],int m,double z[],double dz[],double ddz[]);//第三種邊界條件的三次樣條函數插值、微商與積分
double eslq3(double x[],double y[],double z[],int n,int m,double u,double v);//二元三點插值
double eslgq(double x[],double y[],double z[],int n,int m,double u,double v);//二元全區間插值

//*******************************************************************

double enlgr(double x[],double y[],int n,double t)
{ 
	int i,j,k,m;
    double z,s;
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0];return(z);}
    if (n==2)
      { z=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
        return(z);
      }
    i=0;
    while ((x[i]<t)&&(i<n)) i=i+1;
    k=i-4;
    if (k<0) k=0;
    m=i+3;
    if (m>n-1) m=n-1;
    for (i=k;i<=m;i++)
      { s=1.0;
        for (j=k;j<=m;j++)
          if (j!=i) s=s*(t-x[j])/(x[i]-x[j]);
        z=z+s*y[i];
      }
    return(z);
}
/////////////////////////////////////////////////////////////
double eelgr(double x0,double h,int n,double y[],double t)
{
	int i,j,k,m;
    double z,s,xi,xj;
    double p,q;//float p,q;
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0]; return(z);}
    if (n==2)
      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
        return(z);
      }
    if (t>x0)

      { p=(t-x0)/h; i=(int)p; q=(float)i;
        if (p>q) i=i+1;
      }
    else i=0;
    k=i-4;
    if (k<0) k=0;
    m=i+3;
    if (m>n-1) m=n-1;
    for (i=k;i<=m;i++)
      { s=1.0; xi=x0+i*h;
        for (j=k; j<=m; j++)
          if (j!=i)
            { xj=x0+j*h;
              s=s*(t-xj)/(xi-xj);
            }
        z=z+s*y[i];
      }
    return(z);
}
/////////////////////////////////////////////////////////////
double enlg3(double x[],double y[],int n,double t)
{ 
	int i,j,k,m;
    double z,s;
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0]; return(z);}
    if (n==2)
      { z=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
        return(z);
      }
    if (t<=x[1]) { k=0; m=2;}
    else if (t>=x[n-2]) { k=n-3; m=n-1;}
    else
      { k=1; m=n;
        while (m-k!=1)
          { i=(k+m)/2;
            if (t<x[i-1]) m=i;
            else k=i;
          }
        k=k-1; m=m-1;
        if (fabs(t-x[k])<fabs(t-x[m])) k=k-1;
        else m=m+1;
      }
    z=0.0;
    for (i=k;i<=m;i++)
      { s=1.0;
        for (j=k;j<=m;j++)
          if (j!=i) s=s*(t-x[j])/(x[i]-x[j]);
        z=z+s*y[i];
      }
    return(z);
}
/////////////////////////////////////////////////////////////
double eelg3(double x0,double h,int n,double y[],double t)
{ 
	int i,j,k,m;
    double z,s,xi,xj;
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0]; return(z);}
    if (n==2)
      { z=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
        return(z);
      }
    if (t<=x0+h) { k=0; m=2;}
    else if (t>=x0+(n-3)*h) { k=n-3; m=n-1;}
    else
      { i=(int)((t-x0)/h)+1;
        if (fabs(t-x0-i*h)>=fabs(t-x0-(i-1)*h))
          { k=i-2; m=i;}
        else {k=i-1; m=i+1;}
      }
    z=0.0;
    for (i=k;i<=m;i++)
      { s=1.0; xi=x0+i*h;
        for (j=k;j<=m;j++)
          if (j!=i)
            { xj=x0+j*h; s=s*(t-xj)/(xi-xj);}
        z=z+s*y[i];
      }
    return(z);
}
/////////////////////////////////////////////////////////////
double enpqs(double x[],double y[],int n,double t)
{ 
	int i,j,k,m,l;
    double z,h,b[8];
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0]; return(z);}
    if (n<=8) { k=0; m=n;}
    else if (t<x[4]) { k=0; m=8;}
    else if (t>x[n-5]) { k=n-8; m=8;}
    else
      { k=1; j=n;
        while (j-k!=1)
          { i=(k+j)/2;
            if (t<x[i-1]) j=i;
            else k=i;
          }
        k=k-4; m=8;
      }
    b[0]=y[k];
    for (i=2;i<=m;i++)
      { h=y[i+k-1]; l=0; j=1;
        while ((l==0)&&(j<=i-1))
            { if (fabs(h-b[j-1])+1.0==1.0) l=1;
              else h=(x[i+k-1]-x[j+k-1])/(h-b[j-1]);
              j=j+1;
            }
        b[i-1]=h;
        if (l!=0) b[i-1]=1.0e+35;
      }
    z=b[m-1];
    for (i=m-1;i>=1;i--) z=b[i-1]+(t-x[i+k-1])/z;
    return(z);
}
/////////////////////////////////////////////////////////////
double eepqs(double x0,double h,int n,double y[],double t)
{
	int i,j,k,m,l;
    double z,hh,xi,xj,b[8];
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0]; return(z);}
    if (n<=8) { k=0; m=n;}
    else if (t<(x0+4.0*h)) { k=0; m=8;}
    else if (t>(x0+(n-5)*h)) { k=n-8; m=8;}
    else { k=(int)((t-x0)/h)-3; m=8;}
    b[0]=y[k];
    for (i=2;i<=m;i++)
      { hh=y[i+k-1]; l=0; j=1;
        while ((l==0)&&(j<=i-1))
            { if (fabs(hh-b[j-1])+1.0==1.0) l=1;
              else
                { xi=x0+(i+k-1)*h;
                  xj=x0+(j+k-1)*h;
                  hh=(xi-xj)/(hh-b[j-1]);
                }
              j=j+1;
            }
        b[i-1]=hh;
        if (l!=0) b[i-1]=1.0e+35;
      }
    z=b[m-1];
    for (i=m-1;i>=1;i--)
      z=b[i-1]+(t-(x0+(i+k-1)*h))/z;
    return(z);
}
/////////////////////////////////////////////////////////////
double enhmt(double x[],double y[],double dy[],int n,double t)
{ 
	int i,j;
    double z,p,q,s;
    z=0.0;
    for (i=1;i<=n;i++)
      { s=1.0;
        for (j=1;j<=n;j++)
          if (j!=i) s=s*(t-x[j-1])/(x[i-1]-x[j-1]);
        s=s*s;
        p=0.0;
        for (j=1;j<=n;j++)
          if (j!=i) p=p+1.0/(x[i-1]-x[j-1]);
        q=y[i-1]+(t-x[i-1])*(dy[i-1]-2.0*y[i-1]*p);
        z=z+q*s;
      }
    return(z);
}
/////////////////////////////////////////////////////////////
double eehmt(double x0,double h,int n,double y[],double dy[],double t)
{ 
	int i,j;
    double z,s,p,q;
    z=0.0;
    for (i=1;i<=n;i++)
      { s=1.0; q=x0+(i-1)*h;
        for (j=1;j<=n;j++)
          { p=x0+(j-1)*h;
            if (j!=i) s=s*(t-p)/(q-p);
          }
        s=s*s;
        p=0.0;
        for (j=1;j<=n;j++)
          if (j!=i) p=p+1.0/(q-(x0+(j-1)*h));
        q=y[i-1]+(t-q)*(dy[i-1]-2.0*y[i-1]*p);
        z=z+q*s;
      }
    return(z);
}
/////////////////////////////////////////////////////////////
double enatk(double x[],double y[],int n,double t,double eps)
{ 
	int i,j,k,m,l;
    double z,xx[10],yy[10];
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0]; return(z);}
    m=10;
    if (m>n) m=n;
    if (t<=x[0]) k=1;
    else if (t>=x[n-1]) k=n;
    else
      { k=1; j=n;
        while ((k-j!=1)&&(k-j!=-1))
          { l=(k+j)/2;
            if (t<x[l-1]) j=l;
            else k=l;
          }
        if (fabs(t-x[l-1])>fabs(t-x[j-1])) k=j;
      }
    j=1; l=0;
    for (i=1;i<=m;i++)
      { k=k+j*l;
        if ((k<1)||(k>n))
          { l=l+1; j=-j; k=k+j*l;}
        xx[i-1]=x[k-1]; yy[i-1]=y[k-1];
        l=l+1; j=-j;
      }
    i=0;
    do
      { i=i+1; z=yy[i];
        for (j=0;j<=i-1;j++)
          z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx[i]);
        yy[i]=z;
      }
    while ((i!=m-1)&&(fabs(yy[i]-yy[i-1])>eps));
    return(z);
}
/////////////////////////////////////////////////////////////
double eeatk(double x0,double h,int n,double y[],double t,double eps)
{ 
	int i,j,k,m,l;
    double z,xx[10],yy[10];
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0]; return(z);}
    m=10;
    if (m>n) m=n;
    if (t<=x0) k=1;
    else if (t>=x0+(n-1)*h) k=n;
    else
      { k=1; j=n;
        while ((k-j!=1)&&(k-j!=-1))
          { l=(k+j)/2;
            if (t<x0+(l-1)*h) j=l;
            else k=l;
          }
        if (fabs(t-(x0+(l-1)*h))>fabs(t-(x0+(j-1)*h))) k=j;
      }
    j=1; l=0;
    for (i=1;i<=m;i++)
      { k=k+j*l;
        if ((k<1)||(k>n))
          { l=l+1; j=-j; k=k+j*l;}
        xx[i-1]=x0+(k-1)*h; yy[i-1]=y[k-1];
        l=l+1; j=-j;
      }
    i=0;
    do
      { i=i+1; z=yy[i];
        for (j=0;j<=i-1;j++)
          z=yy[j]+(t-xx[j])*(yy[j]-z)/(xx[j]-xx[i]);
        yy[i]=z;
      }
    while ((i!=m-1)&&(fabs(yy[i]-yy[i-1])>eps));
    return(z);
}
/////////////////////////////////////////////////////////////
void eespl(double x0,double h,int n,double y[],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])/h;
        if (k<0)
          s[4]=(y[1]*(t-x0)-y[0]*(t-x0-h))/h;
        return;
      }
    if (k<0)
      { if (t<=x0+h) kk=0;
        else if (t>=x0+(n-1)*h) kk=n-2;
        else
          { kk=1; m=n;
            while (((kk-m)!=1)&&((kk-m)!=-1))
              { l=(kk+m)/2;
                if (t<x0+(l-1)*h) 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])/h;
    if (n==3)
      { if (kk==0)
          { u[3]=(y[2]-y[1])/h;
            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])/h;
            u[0]=2.0*u[1]-u[2];