From 908ab640896a93a8233e92dfeb5e0504878b4cb1 Mon Sep 17 00:00:00 2001 From: Eric J Alexandrescu Date: Mon, 6 Aug 2018 01:05:52 -0400 Subject: [PATCH] Add files via upload --- AcousticSourceErrorMin.m | 84 ++++++++++ Interpolation.m | 344 +++++++++++++++++++++++++++++++++++++++ Triangulation.m | 60 +++++++ 3 files changed, 488 insertions(+) create mode 100644 AcousticSourceErrorMin.m create mode 100644 Interpolation.m create mode 100644 Triangulation.m diff --git a/AcousticSourceErrorMin.m b/AcousticSourceErrorMin.m new file mode 100644 index 0000000..1966662 --- /dev/null +++ b/AcousticSourceErrorMin.m @@ -0,0 +1,84 @@ + +x1=0; y1=0; %m +x2=0; y2=0.20; +x3=0.10; y3=0; +X=[x1,x2,x3]; +Y=[y1,y2,y3]; +t12=2.7589e-5; %3.92322e-5; %s +t21=t12; +t13=-1.0425e-5; %3.37*40*10^-6; +t31=t13; +t23=-3.8014e-5; %-2.65811e-5; +t32=t23; +int=[0.1,0.1]; +c=1110; +c1=1110; +c2=1110; +c3=1110; + + +E_xy=@(C) (t23*(sqrt((x1-C(1)).^2+(y1-C(2)).^2)-sqrt((x2-C(1)).^2+(y2-C(2)).^2))-... + t12*(sqrt((x2-C(1)).^2+(y2-C(2)).^2)-sqrt((x3-C(1)).^2+(y3-C(2)).^2))).^2+... + (t31*(sqrt((x2-C(1)).^2+(y2-C(2)).^2)-sqrt((x3-C(1)).^2+(y3-C(2)).^2))-... + t23*(sqrt((x3-C(1)).^2+(y3-C(2)).^2)-sqrt((x1-C(1)).^2+(y1-C(2)).^2))).^2+... + (t12*(sqrt((x3-C(1)).^2+(y3-C(2)).^2)-sqrt((x1-C(1)).^2+(y1-C(2)).^2))-... + t31*(sqrt((x1-C(1)).^2+(y1-C(2)).^2)-sqrt((x2-C(1)).^2+(y2-C(2)).^2))).^2; + +E2_xy= @(C) (c*t12-(sqrt((x1-C(1)).^2+(y1-C(2)).^2)-sqrt((x2-C(1)).^2+(y2-C(2)).^2))).^2+... + (c*t13-(sqrt((x1-C(1)).^2+(y1-C(2)).^2)-sqrt((x3-C(1)).^2+(y3-C(2)).^2))).^2+... + (c*t23-(sqrt((x2-C(1)).^2+(y2-C(2)).^2)-sqrt((x3-C(1)).^2+(y3-C(2)).^2))).^2; + + coord4=fminsearch(E2_xy,int); + coord4=fliplr(coord4); + + +%{ +E_xy=@(C) (t23.*c3.*(c2.*sqrt((x1-C(1)).^2+(y1-C(2)).^2)-c1.*sqrt((x2-C(1)).^2+(y2-C(2)).^2))-... + t12.*c1.*(c3.*sqrt((x2-C(1)).^2+(y2-C(2)).^2)-c2.*sqrt((x3-C(1)).^2+(y3-C(2)).^2))).^2+... + (t31.*c1.*(c3.*sqrt((x2-C(1)).^2+(y2-C(2)).^2)-c2.*sqrt((x3-C(1)).^2+(y3-C(2)).^2))-... + t23.*c2.*(c1.*sqrt((x3-C(1)).^2+(y3-C(2)).^2)-c3.*sqrt((x1-C(1)).^2+(y1-C(2)).^2))).^2+... + (t12.*c2.*(c1.*sqrt((x3-C(1)).^2+(y3-C(2)).^2)-c3.*sqrt((x1-C(1)).^2+(y1-C(2)).^2))-... + t31.*c3.*(c2.*sqrt((x1-C(1)).^2+(y1-C(2)).^2)-c1.*sqrt((x2-C(1)).^2+(y2-C(2)).^2))).^2; + %} + +coord2=fminsearch(E_xy,int); + + +[x,y]=meshgrid(0:0.002:max(X),0:0.002:max(Y)); +%[x,y]=meshgrid(0.13:0.00005:0.132,0.022:0.00005:0.026); +%d1=sqrt((x1-x).^2+(y1-y).^2); +%d2=sqrt((x2-x).^2+(y2-y).^2); +%d3=sqrt((x3-x).^2+(y3-y).^2); +%z=(t23*(d1-d2)-t12*(d2-d1)).^2+(t31*(d2-d3)-t23*(d3-d1)).^2+(t12*(d3-d1)-t31*(d1-d2)).^2; +%surf(x,y,z) + +%{ +z2=(t23.*(sqrt((x1-x).^2+(y1-y).^2)-sqrt((x2-x).^2+(y2-y).^2))-... + t12.*(sqrt((x2-x).^2+(y2-y).^2)-sqrt((x3-x).^2+(y3-y).^2))).^2+... + (t31.*(sqrt((x2-x).^2+(y2-y).^2)-sqrt((x3-x).^2+(y3-y).^2))-... + t23.*(sqrt((x3-x).^2+(y3-y).^2)-sqrt((x1-x).^2+(y1-y).^2))).^2+... + (t12.*(sqrt((x3-x).^2+(y3-y).^2)-sqrt((x1-x).^2+(y1-y).^2))-... + t31.*(sqrt((x1-x).^2+(y1-y).^2)-sqrt((x2-x).^2+(y2-y).^2))).^2; +%} + +%An improved algorithm for detecting point of impact in anisotropic inhomogeneous plates +z3=(c*t12-(sqrt((x1-x).^2+(y1-y).^2)-sqrt((x2-x).^2+(y2-y).^2))).^2+... + (c*t13-(sqrt((x1-x).^2+(y1-y).^2)-sqrt((x3-x).^2+(y3-y).^2))).^2+... + (c*t23-(sqrt((x2-x).^2+(y2-y).^2)-sqrt((x3-x).^2+(y3-y).^2))).^2; + + +surf(x,y,log(z3)) +xlabel('x-axis (meters)') +ylabel('y-axis (meters)') +zlabel('log (Error Function)') + +Min=min(z3(:)); +[row,col]=find(z3==Min); +s=size(z3); +Ind=find(z3==Min); +[i,j]=ind2sub(s,Ind); + + +coord3=[x(row,col),y(row,col)] + + diff --git a/Interpolation.m b/Interpolation.m new file mode 100644 index 0000000..11ecc1f --- /dev/null +++ b/Interpolation.m @@ -0,0 +1,344 @@ +%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 xmin(X) && ymin(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 xmean(X) && ymin(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 xmin(X) && ymean(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) && xmean(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 xmin(X) && ymin(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 xmean(X) && ymin(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 xmin(X) && ymean(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) && xmean(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 xmin(X) && ymin(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 xmean(X) && ymin(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 xmin(X) && ymean(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) && xmean(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 xmin(X) && ymin(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 xmean(X) && ymin(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 xmin(X) && ymean(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) && xmean(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) diff --git a/Triangulation.m b/Triangulation.m new file mode 100644 index 0000000..df01fe3 --- /dev/null +++ b/Triangulation.m @@ -0,0 +1,60 @@ +x1=0; y1=0; +x2=0; y2=0.1; +x3=0.1; y3=0; + +c=1110; +t12=2.7589*10^-5; t21=t12; +t23=-3.8014*10^-5; t32=t23; +t13=-1.0425*10^-5; t31=t13; + +d12=abs(c*t12); d21=d12; +d23=abs(c*t23); d32=d23; +d13=abs(c*t13); d31=d13; + +dx=0.001; +eps=dx/10; +[x,y]=meshgrid(0:dx:x3:dx:y2); +z1=zeros(size(x)); +z2=zeros(size(x)); +z3=zeros(size(x)); + +f1=(x1-x).^2+(y1-y).^2; +f2=(x2-x).^2+(y2-y).^2; +f3=(x3-x).^2+(y3-y).^2; + +h=0; +dh=dx; +surface=z1+z2+z3; +iter=0; +while surface<3 & iter<500 + for i=1:length(x) + for j=1:length(y) + if f1(i,j)>(h+d12)^2-eps && f1(i,j)<(h+d12)^2+eps + z1(i,j)=1; + else + z1(i,j)=0; + end + if f2(i,j)>(h)^2-eps && f2(i,j)<(h)^2+eps + z2(i,j)=1; + else + z2(i,j)=0; + end + if f3(i,j)>(d23+h)^2-eps && f3(i,j)<(d23+h)^2+eps + z3(i,j)=1; + else + z3(i,j)=0; + end + end + end + surface=z1+z2+z3; + h=h+dh; + iter=iter+1; +end +contourf(linspace(0,x3,length(x)),linspace(0,y2,length(y)),surface) + +[row,col]=find(surface==3); +location=[col,row]; + +xcoord=mean(location(:,1)); +ycoord=mean(location(:,2)); +location=[x(round(xcoord),round(xcoord)),y(round(ycoord),round(ycoord))]