@Carl Love Thanks, I understand that. My difficulty was the sensitivity of the numerical result

to relatively minor changes to the input values.

## bearing.mpl, solve the target motion problem with bearings only.

##

## Consider a sensor platform moving through points (x,y) at times

## t[1..4] with the target bearings, Brg[1..4] taken at times t[1..4]

## with the target proceeding along a constant course and speed.

##

## time t, bearing line slope m, sensor position (x,y) are known

## values.

##

## Since this is a generated problem the target position at time t is

## provided to compare with the results.

##

#########################################################################

##

restart;

##

genKnownValues := proc()

description "set the known values",

"t - relative time",

"x - sensor x location at time t[i]",

"y - sensor y location at time t[i]",

"m - slope of the bearing lines at time t[i]",

"tgtPosit - target position at time t[i]";

global t, mBrg, mCalc, x, y, tgtPosit;

local dt, Cse, Spd, Brg, A, B, C, R, X;

local tgtX, tgtY, tgtRange, tgtCse, tgtSpd;

## relative and delta time

t := [0, 1+1/2, 3, 3+1/2];

dt := [0, seq(t[idx]-t[idx-1],idx=2..4)];

## sensor motion

Cse := [90, 90, 90, 50] *~ Pi/180; ## true heading

Spd := [15, 15, 15, 22]; ## knots

## bearings to the target at time t

Brg := [10, 358, 340, 330] *~ (Pi/180);

## calculate the sensor position vs time

x := ListTools[PartialSums](dt *~ Spd *~ map(cos, Cse));

y := ListTools[PartialSums](dt *~ Spd *~ map(sin, Cse));

## target values start the target at a known (x,y) position at a

## constant course and speed

tgtRange := 95+25/32; ## miles at t1, match octave value...

tgtCse := 170 * Pi/180; ## course

tgtSpd := 10; ## knots

## tgtX := tgtRange*cos(Brg[1]);

tgtX := 94 + 1/(3 + 1/(3 + 1/(2 + 1/3)));

tgtX := tgtX +~ ListTools[PartialSums](dt *~ tgtSpd *~ cos(tgtCse));

## tgtY := tgtRange*sin(Brg[1]);

tgtY := 17 + 1/(-4 + 1/(-3 + 1/10));

tgtY := tgtY +~ ListTools[PartialSums](dt *~ tgtSpd *~ sin(tgtCse));

## return target position vs time as a matrix

tgtPosit:=Matrix(4,2,[seq([tgtX[idx],tgtY[idx]][],idx=1..4)]);

## slope of the bearing lines

mBrg := map(tan,Brg);

mCalc := [seq((tgtY[idx]-y[idx])/(tgtX[idx] - x[idx]),idx=1..4)];

end proc:

##

#########################################################################

## t[], m[], x[], and y[] are known values

##

## equation of the bearing lines

eq1 := tgtY[1] - y[1] = m[1]*(tgtX[1]-x[1]):

eq2 := tgtY[2] - y[2] = m[2]*(tgtX[2]-x[2]):

eq3 := tgtY[3] - y[3] = m[3]*(tgtX[3]-x[3]):

eq4 := tgtY[4] - y[4] = m[4]*(tgtX[4]-x[4]):

## target X motion along the target line

eq5 := tgtX[2] - tgtX[1] = tgtVx*(t[2]-t[1]):

eq6 := tgtX[3] - tgtX[2] = tgtVx*(t[3]-t[2]):

eq7 := tgtX[4] - tgtX[3] = tgtVx*(t[4]-t[3]):

## target Y motion along the target line

eq8 := tgtY[2] - tgtY[1] = tgtVy*(t[2]-t[1]):

eq9 := tgtY[3] - tgtY[2] = tgtVy*(t[3]-t[2]):

eq10:= tgtY[4] - tgtY[3] = tgtVy*(t[4]-t[3]):

##

#########################################################################

##

## solve the equations

eqs := {eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10}:

Sol:= solve(eqs, {tgtVx, tgtVy, seq([tgtX[k], tgtY[k]][], k= 1..4)}):

##

SolSet := [tgtX[4]=eval(tgtX[4], Sol),

tgtY[4]=eval(tgtY[4], Sol),

tgtVx=eval(tgtVx, Sol),

tgtVy=eval(tgtVy, Sol)]:

genKnownValues():

with(geometry):

m:= mBrg:

point(P1,[eval(tgtX[4],SolSet),eval(tgtY[4],SolSet)]):

m := mCalc:

point(P2,[eval(tgtX[4],SolSet),eval(tgtY[4],SolSet)]):

evalf(distance(P1,P2));