#include <stdio.h>
#include <windows.h>
#include "float.h"
void Interp1(double * x, double * y,int n, double *xn,double *yn,int len)
{
double max,min;
//to find the max&min value of x-series
if(x[0]>x[n-1]){
max=x[0];
min=x[n-1];
}else{
min=x[0];
max=x[n-1];
}
for (int i = 0; i < len; i++) {
int j;
//in case the interoplation point alraedy exists in the array
for (j = 0; j < n; j++)
if(xn[i]==x[j])
break;
if(j!=n)
yn[i] = x[j];
else{
// we need linear interpolation/extrapolation.
if((xn[i]>=min)&&(xn[i]<=max)){
for (int idx = 0; idx < n; idx++)
{
if (xn[i] >= x[idx] && xn[i] <= x[idx+1])
{
//y = y0 + (y1-y0)*(x-x0)/(x1-x0);
double x0 = x[idx];
double x1 = x[idx + 1];
double y0 = y[idx];
double y1 = y[idx + 1];
yn[i] = y0 + (y1 - y0)*(xn[i] - x0) / (x1 - x0);
}
}
}else{//extrapolate
int idx;
//y = y0 + ((x - x0) / (x1 - x0)) * (y1 - y0)
if (xn[i] >= x[n-1])
{
idx = n-1;
}
else
{
idx = 0;
}
double y0 = y[idx];
double y1 = y[idx + 1];
double x0 = x[idx];
double x1 = x[idx + 1];
yn[i] = y0 + ((xn[i] - x0) / (x1 - x0)) * (y1 - y0);
}
}
}
}
double array[188<<2];
double yn_c[188];
int main()
{
FILE *fid;
fid = fopen("F:\\osphotonics\\matlab\\sync\\interp_data.bin","rb");
fread(array,8,188<<2,fid);
fclose(fid);
double *x;
double *y;
double *xn;
double *yn;
x = &array[188*0];
y = &array[188*1];
xn = &array[188*2];
yn = &array[188*3];
Interp1(x,y,188,xn,yn_c,188);
for (int i = 0; i < 188; ++i) {
printf("[%d]x=%.8llf y=%.8llf xn=%.8llf yn_mat=%.8llf yn_calc=%.8llf (diff)=%.8llf\n",i,x[i],y[i],xn[i],yn[i],yn_c[i],yn[i]-yn_c[i]);
}
return 0;
}