Interpolation.m

%This program backs out the location of an impact based on the time delay
%in signal recieving between any two sensors.  Based on the paper: "Low
%Energy Impact Damage Monitoring of Composites Using Dynamic Strain Signals
%From FBG Sensors..."

testnum=9;  %number of reference tests

X=[-80,5,90,-80,5,90,-80,5,90]+81;
Y=[-50,-50,-50,0,0,0,50,50,50]+51;
t12=[144,-16,-149,102,-12,-110,57,-8,65];
t23=[29,37,73,24,4,7,-5,-27,-50];
t34=[-27,26,79,-61,41,153,-143,11,159];
t41=[-127,-60,-46,-25,-22,-13,95,53,17];

nx=max(X)-min(X); %grid dimensions
ny=max(Y)-min(Y);

%create surface of know time delays
plane12=ones(nx,ny);
plane23=ones(nx,ny);
plane34=ones(nx,ny);
plane41=ones(nx,ny);
for i=1:testnum
    plane12(X(i),Y(i))=t12(i);
    plane23(X(i),Y(i))=t23(i);
    plane34(X(i),Y(i))=t34(i);
    plane41(X(i),Y(i))=t41(i);
end

%interpolation calculations for time delays at arbitrary points
fxy1=zeros(mean(X),1);
fxy2=zeros(mean(X),1);
fxy=zeros(mean(X),mean(X));

%plane12
for x=min(X):mean(X)
    for y=min(X):mean(Y)
        if x<mean(X)+1 && x>min(X) && y<mean(Y)+1 && y>min(Y)
            if plane12(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane12(min(X),min(Y))+((x-min(X))/(mean(X)-min(X)))*plane12(mean(X),min(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane12(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane12(mean(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane12(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X)+1:max(X)
    for y=min(Y):mean(Y)
        if x<max(X) && x>mean(X) && y<mean(Y)+1 && y>min(Y)
            if plane12(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane12(mean(X),min(Y))+((x-mean(X))/(max(X)-mean(X)))*plane12(max(X),min(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane12(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane12(max(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane12(x,y)=fxy(x,y);
            end
        end
    end
end
for x=min(X):mean(X)
    for y=mean(Y):max(Y)
        if x<mean(X)+1 && x>min(X) && y<max(Y) && y>mean(Y)
            if plane12(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane12(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane12(mean(X),mean(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane12(min(X),max(Y))+((x-min(X))/(mean(X)-min(X)))*plane12(mean(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane12(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X):max(X)
    for y=mean(Y):max(Y)
        if x>mean(X) && x<max(X) && y<max(Y) && y>mean(Y)
            if plane12(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane12(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane12(max(X),mean(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane12(mean(X),max(Y))+((x-mean(X))/(max(X)-mean(X)))*plane12(max(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane12(x,y)=fxy(x,y);
            end
        end
    end
end

%plane23
for x=min(X):mean(X)
    for y=min(X):mean(Y)
        if x<mean(X)+1 && x>min(X) && y<mean(Y)+1 && y>min(Y)
            if plane23(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane23(min(X),min(Y))+((x-min(X))/(mean(X)-min(X)))*plane23(mean(X),min(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane23(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane23(mean(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane23(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X)+1:max(X)
    for y=min(Y):mean(Y)
        if x<max(X) && x>mean(X) && y<mean(Y)+1 && y>min(Y)
            if plane23(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane23(mean(X),min(Y))+((x-mean(X))/(max(X)-mean(X)))*plane23(max(X),min(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane23(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane23(max(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane23(x,y)=fxy(x,y);
            end
        end
    end
end
for x=min(X):mean(X)
    for y=mean(Y):max(Y)
        if x<mean(X)+1 && x>min(X) && y<max(Y) && y>mean(Y)
            if plane23(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane23(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane23(mean(X),mean(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane23(min(X),max(Y))+((x-min(X))/(mean(X)-min(X)))*plane23(mean(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane23(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X):max(X)
    for y=mean(Y):max(Y)
        if x>mean(X) && x<max(X) && y<max(Y) && y>mean(Y)
            if plane23(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane23(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane23(max(X),mean(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane23(mean(X),max(Y))+((x-mean(X))/(max(X)-mean(X)))*plane23(max(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane23(x,y)=fxy(x,y);
            end
        end
    end
end

%plane34
for x=min(X):mean(X)
    for y=min(X):mean(Y)
        if x<mean(X)+1 && x>min(X) && y<mean(Y)+1 && y>min(Y)
            if plane34(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane34(min(X),min(Y))+((x-min(X))/(mean(X)-min(X)))*plane34(mean(X),min(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane34(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane34(mean(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane34(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X)+1:max(X)
    for y=min(Y):mean(Y)
        if x<max(X) && x>mean(X) && y<mean(Y)+1 && y>min(Y)
            if plane34(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane34(mean(X),min(Y))+((x-mean(X))/(max(X)-mean(X)))*plane34(max(X),min(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane34(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane34(max(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane34(x,y)=fxy(x,y);
            end
        end
    end
end
for x=min(X):mean(X)
    for y=mean(Y):max(Y)
        if x<mean(X)+1 && x>min(X) && y<max(Y) && y>mean(Y)
            if plane34(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane34(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane34(mean(X),mean(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane34(min(X),max(Y))+((x-min(X))/(mean(X)-min(X)))*plane34(mean(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane34(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X):max(X)
    for y=mean(Y):max(Y)
        if x>mean(X) && x<max(X) && y<max(Y) && y>mean(Y)
            if plane34(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane34(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane34(max(X),mean(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane34(mean(X),max(Y))+((x-mean(X))/(max(X)-mean(X)))*plane34(max(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane34(x,y)=fxy(x,y);
            end
        end
    end
end

%plane41
for x=min(X):mean(X)
    for y=min(X):mean(Y)
        if x<mean(X)+1 && x>min(X) && y<mean(Y)+1 && y>min(Y)
            if plane41(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane41(min(X),min(Y))+((x-min(X))/(mean(X)-min(X)))*plane41(mean(X),min(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane41(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane41(mean(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane41(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X)+1:max(X)
    for y=min(Y):mean(Y)
        if x<max(X) && x>mean(X) && y<mean(Y)+1 && y>min(Y)
            if plane41(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane41(mean(X),min(Y))+((x-mean(X))/(max(X)-mean(X)))*plane41(max(X),min(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane41(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane41(max(X),mean(Y));
                fxy(x,y)=((mean(Y)-y)/(mean(Y)-min(Y)))*fxy1(x)+((y-min(Y))/(mean(Y)-min(Y)))*fxy2(x);
                plane41(x,y)=fxy(x,y);
            end
        end
    end
end
for x=min(X):mean(X)
    for y=mean(Y):max(Y)
        if x<mean(X)+1 && x>min(X) && y<max(Y) && y>mean(Y)
            if plane41(x,y)==1
                fxy1(x)=((mean(X)-x)/(mean(X)-min(X)))*plane41(min(X),mean(Y))+((x-min(X))/(mean(X)-min(X)))*plane41(mean(X),mean(Y));
                fxy2(x)=((mean(X)-x)/(mean(X)-min(X)))*plane41(min(X),max(Y))+((x-min(X))/(mean(X)-min(X)))*plane41(mean(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane41(x,y)=fxy(x,y);
            end
        end
    end
end
for x=mean(X):max(X)
    for y=mean(Y):max(Y)
        if x>mean(X) && x<max(X) && y<max(Y) && y>mean(Y)
            if plane41(x,y)==1
                fxy1(x)=((max(X)-x)/(max(X)-mean(X)))*plane41(mean(X),mean(Y))+((x-mean(X))/(max(X)-mean(X)))*plane41(max(X),mean(Y));
                fxy2(x)=((max(X)-x)/(max(X)-mean(X)))*plane41(mean(X),max(Y))+((x-mean(X))/(max(X)-mean(X)))*plane41(max(X),max(Y));
                fxy(x,y)=((max(Y)-y)/(max(Y)-mean(Y)))*fxy1(x)+((y-mean(Y))/(max(Y)-mean(Y)))*fxy2(x);
                plane41(x,y)=fxy(x,y);
            end
        end
    end
end


%find coordinates of constant time delays
[X1,Y1]=find(plane12<-77.5 & plane12>-78.5);
a=[X1,Y1]-90;
[X2,Y2]=find(plane23<-19.5 & plane23>-20.5);
b=[X2,Y2]-90;
[X3,Y3]=find(plane34<91.5 & plane34>90.5);
c=[X3,Y3]-90;
[X4,Y4]=find(plane41<7.5 & plane41>6.5);
d=[X4,Y4]-90;

%plot lines of constant time delays
figure(1)
scatter(a(:,1)',a(:,2)) %isoline in plane12
hold on
scatter(b(:,1)',b(:,2)) %isoline in plane23
hold on
scatter(c(:,1)',c(:,2)) %isoline in plane34
hold on
scatter(d(:,1)',d(:,2)) %isoline in plane41

test=ones(max(X)-min(X),max(Y)-min(Y));
for i=1:length(X1)
    test(X1(i),Y1(i))=200;
end
for i=1:length(X2)
    test(X2(i),Y2(i))=200;
end
for i=1:length(X3)
    test(X3(i),Y3(i))=200;
end
for i=1:length(X4)
    test(X4(i),Y4(i))=200;
end

%calculate intersections of the isolines
Isect12=[0,0];
for i=1:length(X1)
    for j=1:length(X2)
        if X1(i)==X2(j) && Y1(i)==Y2(j)
            Isect12=[X1(i),Y1(i)];
        end
    end
end

Isect13=[0,0];
for i=1:length(X1)
    for j=1:length(X3)
        if X1(i)==X3(j) && Y1(i)==Y3(j)
            Isect13=[X1(i),Y1(i)];
        end
    end
end

Isect14=[0,0];
for i=1:length(X1)
    for j=1:length(X4)
        if X1(i)==X4(j) && Y1(i)==Y4(j)
            Isect14=[X1(i),Y1(i)];
        end
    end
end

Isect23=[0,0];
for i=1:length(X2)
    for j=1:length(X3)
        if X2(i)==X3(j) && Y2(i)==Y3(j)
            Isect23=[X2(i),Y2(i)];
        end
    end
end

Isect24=[0,0];
for i=1:length(X2)
    for j=1:length(X4)
        if X2(i)==X4(j) && Y2(i)==Y4(j)
            Isect24=[X2(i),Y2(i)];
        end
    end
end

Isect34=[0,0];
for i=1:length(X3)
    for j=1:length(X4)
        if X3(i)==X4(j) && Y3(i)==Y4(j)
            Isect34=[X3(i),Y3(i)];
        end
    end
end

%count number of intersections
count=[Isect12(1),Isect13(1),Isect14(1),Isect23(1),Isect24(1),Isect34(1)];
count=count(count~=0);

%calculate predicted x,y coordinates
XX=(Isect12(1)+Isect13(1)+Isect14(1)+Isect23(1)+Isect24(1)+Isect34(1));
YY=(Isect12(2)+Isect13(2)+Isect14(2)+Isect23(2)+Isect24(2)+Isect34(2));

X_avg=XX/length(count)-90;
Y_avg=YY/length(count)-90;

test(round(X_avg)+90,round(Y_avg)+90)=300;  %predicted point of impact
test(50+90,-20+90)=400;                     %actual point of impact = (50,-20)

%graphs
figure(1)
contour(test)

figure(2)
surf(test)