Triangulation.m

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))]