I got 2 main functions.
*
Linear interpolation in 1D array:
Function is Y(X) and the result is Y(x).
*
Linear interpolation in 2D array:
Function is Y(X1, X2) and the result is Y(x, y).
And function to test it.
#include
#include
#include
double lin_int(int size, double* X, const double* Y, double x){
int i;
if (x >= *(X + size-1)){
return *(Y + size-1);
} else
if (x = x && *(X + i - 1) *(X1+m) && y > *(X2+n)) return *(Y+m*n-1); //по X1 и X2 больше оба значения x, y - возвращаем Y[m][n]
if (x < *(X1)) x = *(X1); //по X1 значение x меньше X1[0] - ставим x = X1[0]
if (y < *(X2)) y = *(X2); //по X2 значение y меньше X2[0] - ставим y = X2[0]
if (x > *(X1+m-1)) x = *(X1+m-1); //по X1 значение x больше X1[m] - ставим x = X1[m]
if (y > *(X2+n-1)) y = *(X2+n-1); //по X2 значение y больше X2[n] - ставим y = X2[n]
for (i = 0; i < n-1; i++){
if (x == *(X1 + i)){ //Проверка что мы в узле X1
in_X1 = true;
break;
} else
if (x > *(X1 + i) && x < *(X1 + i + 1)) break; //Если не в узле, фиксируем индекс предыдущего узла
}
for (j = 0; j < m-1; j++){
if (y == *(X2 + j)){ //Проверка что мы в узле X2
in_X2 = true;
break;
} else
if (y > *(X2 + j) && y < *(X2 + j + 1)) break; //Если не в узле, фиксируем индекс предыдущего узла
}
if (in_X1 && in_X2) { //Если попали в два узла сетки, возвращаем уже известное значение
return *((Y + i*n) + j);
} else
if (in_X1) { //Если в узле только по X1, возвращаем lin_int между узлами X2
return lin_int(n, X2, (Y+i*n), y);
} else
if (in_X2) { //Если в узле только по X2, возвращаем lin_int между узлами X1
x1 = *((Y+i*n)+j);
x2 = *((Y+(i+1)*n)+j);
if (*(X1+i+1) == *(X1+i)){
return x1 + (x2-x1) * (x-*(X1+i)) / (*(X1+i)*2); //два совпадающих значения в узлах
}
return x1 + (x2-x1) * (x-*(X1+i)) / (*(X1+i+1) - *(X1+i));
} else { //Если ни в какой узел не попали
if (*(X1+i+1) == *(X1+i)){
return x1 + (x2-x1) * (x-*(X1+i)) / (*(X1+i)*2); //два совпадающих значения в узлах
}
x1 = lin_int(n, X2, (Y+i*n), y); //Если не попали ни в какой узел, вычисляем x1, x2 и возвращаем lin_int между ними
x2 = lin_int(n, X2, (Y+(i+1)*n), y);
return x1 + (x2-x1) * (x-*(X1+i)) / (*(X1+i+1) - *(X1+i));
}
}
void test_lin_int_mass(void){
double X1[4] = { 0.3, 1.0, 4.0, 10.0 };
double X2[3] = { 0.0, 90.0, 180.0 };
double Y[4][3] =
{/* 0 90 180 */
/*0.3*/ {1.0, -0.107, -0.5},
/*1.0*/ {1.3, -0.3, -0.7},
/*4.0*/ {1.5, -0.2, -0.6},
/*10.*/ {1.53, -0.253, -0.5}
};
double i,j;
for (i=0; i
0 comments:
Post a Comment
Thanks