#
# Do all the initialisations
#
theta:=Array([161.2, 45.10, 309.0]):
sigma:= Array([0.8, 0.6, 1.3, 2.0]):
X:= Array([746, 629, 1571, 155]):
Y:= Array([1393, 375, 259, 987]):
d4:=864.3:
pi:=3.1415926540:
for i to 3 do
theta[i]:= 2*pi*theta[i]/360;
sigma[i]:= 2*pi*sigma[i]/360;
end do:
#
# Generate equations in appropriate form
# for solution
#
S:=add(((arctan(x-X[i],y-Y[i])-theta[i])/sigma[i])^2, i=1..3)
+((((x-X[4])^2+(y-Y[4])^2)^(1/2)-d4)/sigma[4])^2:
eqns:={diff(S,x)=0, diff(S,y)=0}:#
# Attempt solution
#
sol:=fsolve(eqns, {x, y});
whattype(sol);#
# Hmmm no solutions. Note the denominators: whenever y=Y[i]
# looks like one ends up with "0*Infinity" - pretty difficult
# to work with.
#
# Plot the surfaces just to see if there are regions where
# solutions look possible
#
plot3d( [ lhs(eqns[1]), lhs(eqns[2]) ],
x = -1000..1000,
y = -1000..1000
);
#
# Zoom in to a likely region - the big spikey bit
#
plot3d( [ lhs(eqns[1]), lhs(eqns[2]) ],
x = 100..400,
y = 200..500
);#
# Looks promising so test for a solution in this region
# (Probably many other solutions - unfortunately I suspect
# most of them are close to the areas where the functions
# become indeterminate
#
sol:=fsolve
( { diff(S,x) = 0, diff(S,y) = 0},
{ x = 100..300, y = 300..500}
);#
# *Looks* as if all solutions are for x>0, y>0, and probably(?)
# don't extend much beyond max(X), max(Y) so split this region
# into sub-regions and search within each subregion
#
t1:=time():
eqns:= { diff(S,x) = 0, diff(S,y) = 0 }:
for xr from 0 by 200 to 2000 do
for yr from 0 by 200 to 2000 do
sol:= fsolve
( eqns,
{ x = xr..xr+200, y = yr..yr+200}
);
if whattype(sol)=set
then printf(" sol=%a\134n", sol);
fi;
od;
od;
t2:=time()-t1;