restart;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;for i to 3 do
theta[i] := evalf( 2*Pi*theta[i] / 360 );
if theta[i] > evalf(Pi) then
theta[i] := theta[i] - evalf(2*Pi) fi;
sigma[i] := evalf( 2*Pi*sigma[i] / 360 );
od:
i := 'i':
print( theta ); print( sigma );S := sum( ( (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;Sx := diff(S, x):
Sy := diff(S, y):Rxy := x = 900..1000, y = 600..800:Gr := grid=[40, 40]:PX:=plots[implicitplot](Sx, Rxy, Gr, thickness=1):
PY:=plots[implicitplot](Sy, Rxy, Gr, thickness=4):plots[display]({PX ,PY});sol := fsolve( { diff(S,x)=0, diff(S,y)=0 }, {x,y}, {Rxy} );S0 := evalf( subs( sol, S ));S1 := evalf( subs( sol, linalg[grad]( S, [x,y] )));S2 := evalf( subs( sol, linalg[hessian]( S, [x,y] ) ));
pmp0 := [ x-subs(sol,x), y-subs(sol,y) ];
Sapprox := expand( S0 + evalm( transpose(pmp0) &* S2 &* pmp0 ));stats[ statevalf, icdf, chisquare[4] ] ( 0.95 );ellips := { seq( stats[ statevalf, icdf, chisquare[4]](c) =
Sapprox, c = [0.5, 0.95, 0.999] ) }:
plots[implicitplot]( ellips, x = 950..1000, y = 700..750,
grid=[50,50], view=[950..1000, 700..750] );ev := linalg[eigenvects]( S2 );