双三次样条插值算法代码

双三次样条插值算法代码 该代码是c语言完成的，已经调试过，可以直接嵌入程序中运行。

#include 
#include 
#include 
#define NR_END 1
#define FREE_ARG char*

void splie2（double x1a[]double x2a[]double** yaint mint ndouble** y2a）;
void spline（double x[]double y[]int ndouble yp1double ypndouble y2[]）;
double *vector（long nllong nh）;
void nrerror（char error_text[]）;
void free_vector（double* vlong nllong nh）;
void splin2（double x1a[]double x2a[]double** yadouble** y2aint mint ndouble x1double x2double*y）;
void splint（double xa[]double ya[]double y2a[]int ndouble xdouble* y）;


void splie2（double x1a[]double x2a[]double** yaint mint ndouble** y2a）
//给定列表函数ya[1..m][1..n]及独立变量x1a[1..m]和x2a[1..n]本程序构造ya的行的一维自然三次样条
//返回的二阶导数存于数组y2a[1..m][1..n]中。
{
	for（int j=1;j<=m;j++）
		spline（x2aya[j]n1.0e301.0e30y2a[j]）;//值1.0e30标志自然样条
}

void splin2（double x1a[]double x2a[]double** yadouble** y2aint mint ndouble x1double x2double*y）
//给定列表函数ya[1..m][1..n]及独立变量x1a[1..m]和x2a[1..n]并由该程序生成y2a，对给定所求的插值点x1x2，本程序
//用双三次样条插值求得插值函数值y。
{
	int j;
	double *ytmp*yytmp;

	ytmp=vector（1m）;
	yytmp=vector（1m）;
	for（j=1;j<=m;j++）
		splint（x2aya[j]y2a[j]nx2&yytmp[j]）;
	spline（x1ayytmpm1.0e301.0e30ytmp）;
	splint（x1ayytmpytmpmx1y）;
	free_vector（yytmp1m）;
	free_vector（ytmp1m）;
}

void spline（double x[]double y[]int ndouble yp1double ypndouble y2[]）
//给出数组x[1..n]和y]1..n]构成一个列表函数，即yi=f（xi）其中x1//第n个点处的一阶导致值yp1和ypn。相应地，本程序返回数组y2[1..n]，其值为插值函数在列表中xi点的二阶导数。
//如果yp1和（或）ypn大于或等于1.0e30则本程序将相应的边界条件设为自然样条，即边界的二阶导数为零。
{
	int ik;
    double pqnsigun*u;

	u=vector（1n-1）;
	if（yp1>0.99e30）
	    y2[1]=u[1]=0.0;   
	else
	{
		y2[1]=-0.5;
		u[1]=（3.0/（x[2]-x[1]））*（（y[2]-y[1]）/（x[2]-x[1]）-yp1）;	   
	}

	for（i=2;i<=n-1;i++）
	{
		sig=（x[i]-x[i-1]）/（x[i+1]-x[i-1]）;
		p=sig*y2[i-1]+2.0;
		y2[i]=（sig-1.0）/p;
		u[i]=（y[i+1]-y[i]）/（x[i+1]-x[i]）-（y[i]-y[i-1]）/（x[i]-x[i-1]）;
		u[i]=（6.0*u[i]/（x[i+1]-x[i-1]）-sig*u[i-1]）/p;
	}

	if（ypn>0.99e30）
		qn=un=0.0;
	else
	{
		qn=0.5;
		un=（3.0/（x[n]-x[n-1]））*（ypn-（y[n]-y[n-1]）/（x[n]-x[n-1]））;
	}

	y2[n]=（un-qn*u[n-1]）/（qn*y2[n-1]+1.0）;
	for（k=n-1;k>=1;k--）
		y2[k]=y2[k]*y2[k+1]+u[k];

	free_vector（u1n-1）;
}

double* vector（long nllong nh）
{
	double* v;
	v=（double*）malloc（（size_t）（（nh-nl+1+NR_END）*sizeof（double）））;
	if（!v）
		nrerror（“allocation failure in vector（）“）;

	return v-nl+NR_END;
}

void nrerror（char error_text[]）
{
	fprintf（stderr“Numerical Recipes run-time error...\n“）;
	fprintf（stderr“%s\n“error_text）;
	fprintf（stderr“...now exiting to system...\n“）;
	exit（1）;
}

void free_vector（double* vlong nllong nh）
{
	free（（FREE_ARG） （v+nl-NR_END））;
}

void splint（double xa[]double ya[]double y2a[]int ndouble xdouble* y）
{
	int klokhik;
	double hba;

	klo=1;
	khi=n;
	while（（khi-klo）>1）
	{
		k=（khi+klo）>>1;
		if（xa[k]>x）
			khi=k;
		else
			klo=k;
	}

	h=xa[khi]-xa[klo];
	if（h==0.0）
		nrerror

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       4304  2008-11-05 17:17  双三次样条插值\bi_cubic.dsp

     文件        539  2008-11-05 15:57  双三次样条插值\bi_cubic.dsw

     文件      50176  2009-01-06 18:27  双三次样条插值\bi_cubic.ncb

     文件      48640  2009-01-06 18:27  双三次样条插值\bi_cubic.opt

     文件       1272  2008-11-17 20:20  双三次样条插值\bi_cubic.plg

     文件      11260  2008-11-17 20:20  双三次样条插值\main.cpp

     目录          0  2009-01-06 18:27  双三次样条插值

----------- ---------  ---------- -----  ----

               116191                    7