C语言实现matlab的interp2函数

 

 

/********一维差值函数*******/

double inter_linear(double x0, double x1, double y0, double y1, double x)
{
double a0, a1, y;
a0 = (x - x1) / (x0 - x1);
a1 = (x - x0) / (x1 - x0);
y = a0*y0 + a1*y1;

return (y);
}
/****************************/

/*****二维插值函数*****/
double interp2(double *x, double *y, double z[][500], int m, int n, double a, double b) //一维数组x、y。二维数组z，500可换。mn为向量大小。ab为插值点位置
{
int i, j, tmpi, tmpj;
double w, w1, w2;
tmpi = 0;
tmpj = 0;
w = 0.0;
for (i = 0; i<m; i++) //确定a在x轴的位置
{
if ((a <= x[i]) && (a>x[i + 1])) //
{
tmpi = i;
break;
}
if( (a>=x[i])&&(a<x[i+1]) ) //
{
tmpi = i;
break;
}
}
for (j = 0; j<n; j++) //确定b在y轴的位置
{
if ((b >= y[j]) && (b<y[j + 1]))
{
tmpj = j;
break;
}
}
/********x方向进行插值*****************/
if (x[tmpi] == a)
{
//
w1 = z[tmpi][tmpj];
w2 = z[tmpi][tmpj + 1];
/**********/
if (y[tmpj] == b)
{
//
w = w1;
}
else
{
//
w = inter_linear(y[tmpj], y[tmpj + 1], w1, w2, b);
}
}
else
{
//
w1 = inter_linear(x[tmpi], x[tmpi + 1], z[tmpi][tmpj], z[tmpi + 1][tmpj], a);
w2 = inter_linear(x[tmpi], x[tmpi + 1], z[tmpi][tmpj + 1], z[tmpi + 1][tmpj + 1], a);
if (y[tmpj] == b)
{
//
w = w1;
}
else
{
//
w = inter_linear(y[tmpj], y[tmpj + 1], w1, w2, b);
}
}
return (w);
}
/****************************/