lats = [34.1989 34.0105 34.067894 34.1166 34.0151]; longs = [72.0231 71.9876 71.992783 72.0216 71.5249];
[x,y] = grn2eqa(lats,longs,[34.1166, 72.0216])
t1 = (sqrt((x(5)-x(4))^2+(y(5)-y(4))^2+(z_s-z(4))^2)-sqrt((x(5)-x(1))^2+(y(5)-y(1))^2+(z_s-z(1))^2))/c
t2 = (sqrt((x(5)-x(4))^2+(y(5)-y(4))^2+(z_s-z(4))^2)-sqrt((x(5)-x(2))^2+(y(5)-y(2))^2+(z_s-z(2))^2))/c
t3 = (sqrt((x(5)-x(4))^2+(y(5)-y(4))^2+(z_s-z(4))^2)-sqrt((x(5)-x(3))^2+(y(5)-y(3))^2+(z_s-z(3))^2))/c
eqn1 = sqrt((xs-x(4))^2+(ys-y(4))^2+(zs-z(4))^2)-sqrt((xs-x(1))^2+(ys-y(1))^2+(zs-z(1))^2)-(c*t1);
eqn2 = sqrt((xs-x(4))^2+(ys-y(4))^2+(zs-z(4))^2)-sqrt((xs-x(2))^2+(ys-y(2))^2+(zs-z(2))^2)-(c*t2);
eqn3 = sqrt((xs-x(4))^2+(ys-y(4))^2+(zs-z(4))^2)-sqrt((xs-x(3))^2+(ys-y(3))^2+(zs-z(3))^2)-(c*t3);
sol = solve([eqn1, eqn2, eqn3], [xs ys zs]);
fimplicit3(eqn1,[-0.008 0.001 -0.002 0.0015 -0.015 0.015],'EdgeColor','b','FaceAlpha',0)
fimplicit3(eqn2,[-0.008 0.001 -0.002 0.0015 -0.015 0.015],'EdgeColor','g','FaceAlpha',0)
fimplicit3(eqn3,[-0.008 0.001 -0.002 0.0015 -0.015 0.015],'EdgeColor','y','FaceAlpha',0)
possibleSol(1,m) = double(sol.xs(n));
possibleSol(2,m) = double(sol.ys(n));
possibleSol(3,m) = double(sol.zs(n))
idx = possibleSol(3,:) < 0 | any(imag(possibleSol) ~=0)
possibleSol(:, idx) = [];
[lat,long] = eqa2grn(possibleSol(1),possibleSol(2),[34.1166, 72.0216])
plot3(x(1:4),y(1:4),z, 'ro', 'LineWidth', 2, 'MarkerSize', 10);
plot3(possibleSol(1),possibleSol(2),possibleSol(3), 'b+', 'LineWidth', 4, 'MarkerSize', 10)
plot3(x(5),y(5),z_s,'g+', 'LineWidth', 2, 'MarkerSize', 10)
legend({'Receivers', 'Source','exactvalue'})