%Program Awange-Grafarend Multipolynomial Resultant Algorithm % Department of Geodesy and GeoInformatics, University of Stuttgart, Germany-29th of August 2001 % Program computes the coefficients cs2 cs1 cs0 of the quadratic equation obtained from the Groebner basis approach % All that is required from the user is to replace the four satelite coordinates and their % respective pseudo-range measurements as explained in the input section of the program % The example considered is adopted from A. Kleusberg (1994), Zeitschrift für Vermessungswesen 119(1994)188-192 % and Grafarend and Shan (1996), ARTIFICIAL SATELLITES, Planetary Geodesy No 28 Vol 31 pp.133-147. %Multipolynomial solution of pseudorange equations based on the B. Sturmfel 1998 approach (Proceedings of Symposium in % applied Mathematics 53(1998) 25-29. clear all syms x1 x2 x3 x4 a03 b03 c03 d30 e03 a13 b13 syms c13 d31 e13 a23 b23 c23 d32 e23 a0 a1 a2 a3 b0 b1 b2 b3 c0 c1 c2 c3 d0 d1 d2 d3 % Hidding the variables by assigning them degree zero %Hidding x3 jc13=c03*x3+d30*x4+e03; jc23=c13*x3+d31*x4+e13; jc33=c23*x3+d32*x4+e23; % Hidding x2 jb12=b03*x2+d30*x4+e03; jb22=b13*x2+d31*x4+e13; jb32=b23*x2+d32*x4+e23; %Hidding x1 ja11=a03*x1+d30*x4+e03; ja21=a13*x1+d31*x4+e13; ja31=a23*x1+d32*x4+e23; % Computing the Jacobean determinants %determinant jc,jb,ja jc=det([a03 b03 jc13;a13 b13 jc23;a23 b23 jc33]); jb=det([a03 jb12 c03;a13 jb22 c13;a23 jb32 c23]); ja=det([ja11 b03 c03;ja21 b13 c13;ja31 b23 c23]); cjc=collect(jc,x3); cjb=collect(jb,x2); cja=collect(ja,x1); % Symbolic computation of the stationary receiver coordinates {z=x3},{y=x2} and {x=x1} cx3=-(a23*b03*d31*x4+a03*b13*d32*x4+a03*b13*e23-a23*b13*d30*x4-a03*b23*d31*x4-a03*b23*e13+... a13*b23*d30*x4-a13*b03*d32*x4-a13*b03*e23-a23*b13*e03+a23*b03*e13+a13*b23*e03)/... (a23*b03*c13+a13*b23*c03-a13*b03*c23-a23*b13*c03-a03*b23*c13+a03*b13*c23); cx2=-(a23*c13*d30*x4+a03*c23*d31*x4+a03*c23*e13-a23*c03*d31*x4-a03*c13*d32*x4-... a03*c13*e23+a13*c03*d32*x4-a13*c23*d30*x4-a13*c23*e03-a23*c03*e13+a23*c13*e03+a13*c03*e23)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13); cx1=-(e03*b13*c23+d32*x4*b03*c13+d30*x4*b13*c23-d30*x4*c13*b23-d31*x4*b03*c23-e03*c13*b23-... e13*b03*c23+e13*c03*b23+e23*b03*c13+d31*x4*c03*b23-d32*x4*c03*b13-e23*c03*b13)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13); % Symbolic computation of the Coefficients of the quadratic equation cf2=((-a23*b03*d31+a23*b13*d30+a03*b23*d31-a03*b13*d32-a13*b23*d30+a13*b03*d32)^2/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)^2+... (-a23*c13*d30+a23*c03*d31+a03*c13*d32-a03*c23*d31-a13*c03*d32+a13*c23*d30)^2/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)^2+... (-d32*b03*c13-d30*b13*c23+d31*b03*c23-d31*c03*b23+d32*c03*b13+d30*c13*b23)^2/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)^2-1);%*x4^2 cf1=(2*((a13*c23*e03+a03*c13*e23-a03*c23*e13-a13*c03*e23+a23*c03*e13-a23*c13*e03)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)-b0)*... (-a23*c13*d30+a23*c03*d31+a03*c13*d32-a03*c23*d31-a13*c03*d32+a13*c23*d30)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)+2*... ((-e03*b13*c23+e03*c13*b23+e13*b03*c23-e13*c03*b23-e23*b03*c13+e23*c03*b13)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)-a0)*... (-d32*b03*c13-d30*b13*c23+d31*b03*c23-d31*c03*b23+d32*c03*b13+d30*c13*b23)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)+2*d0+2*... ((a13*b03*e23+a03*b23*e13-a03*b13*e23-a13*b23*e03+a23*b13*e03-a23*b03*e13)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)-c0)*... (-a23*b03*d31+a23*b13*d30+a03*b23*d31-a03*b13*d32-a13*b23*d30+a13*b03*d32)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13));%*x4 cf0=((-e03*b13*c23+e03*c13*b23+e13*b03*c23-e13*c03*b23-e23*b03*c13+e23*c03*b13)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)-a0)^2-... d0^2+((a13*b03*e23+a03*b23*e13-a03*b13*e23-a13*b23*e03+a23*b13*e03-a23*b03*e13)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)-c0)^2+... ((a13*c23*e03+a03*c13*e23-a03*c23*e13-a13*c03*e23+a23*c03*e13-a23*c13*e03)/... (a23*c13*b03+a13*b23*c03-a13*c23*b03-a23*b13*c03-a03*c13*b23+a03*c23*b13)-b0)^2; % Input data as follows: % {a0,b0,c0}={X1,Y1,Z1} are the coordinates of the first satellite. d0 is the pseudo-range measured to this satellie % {a1,b1,c1}={X2,Y2,Z2} are the coordinates of the second satellite. d1 is the pseudo-range measured to this satellie % {a2,b2,c2}={X3,Y3,Z3} are the coordinates of the third satellite. d2 is the pseudo-range measured to this satellie % {a3,b3,c3}={X4,Y4,Z4} are the coordinates of the fourth satellite. d3 is the pseudo-range measured to this satellie a0=1.4832308660e+007;a1=-1.5799854050e+007;a2=1.984818910e+006;a3=-1.2480273190e+007; b0=-2.0466715890e+007;b1=-1.3301129170e+007;b2=-1.1867672960e+007;b3=-2.3382560530e+007; c0=-7.428634750e+006;c1=1.7133838240e+007;c2=2.3716920130e+007;c3=3.278472680e+006; d0=2.4310764064e+007;d1=2.2914600784e+007;d2=2.0628809405e+007;d3=2.3422377972e+007; f00=-a0^2-b0^2-c0^2+d0^2;f11=-a1^2-b1^2-c1^2+d1^2;f22=-a2^2-b2^2-c2^2+d2^2;f33=-a3^2-b3^2-c3^2+d3^2; a03=2*(a0-a3);b03=2*(b0-b3);c03=2*(c0-c3);d30=2*(d3-d0);e03=f00-f33; a13=2*(a1-a3);b13=2*(b1-b3);c13=2*(c1-c3);d31=2*(d3-d1);e13=f11-f33; a23=2*(a2-a3);b23=2*(b2-b3);c23=2*(c2-c3);d32=2*(d3-d2);e23=f22-f33; cs2=eval(cf2),cs1=eval(cf1),cs0=eval(cf0) % Computation of the stationary receiver range bias x4 sol1=vpa(solve(cs2*x4^2+cs1*x4+cs0),13);x4=sol1 % Computation of the stationary receiver coordinates {x1=x},{x2=y} and {x3=z} x3=eval(cx3) x2=eval(cx2) x1=eval(cx1) % Check by subtracting the measured from computed distances. check1=sqrt((a0-x1(2,1))^2+(b0-x2(2,1))^2+(c0-x3(2,1))^2)-(d0-x4(2,1)) check2=sqrt((a1-x1(2,1))^2+(b1-x2(2,1))^2+(c1-x3(2,1))^2)-(d1-x4(2,1)) check3=sqrt((a2-x1(2,1))^2+(b2-x2(2,1))^2+(c2-x3(2,1))^2)-(d2-x4(2,1)) check4=sqrt((a3-x1(2,1))^2+(b3-x2(2,1))^2+(c3-x3(2,1))^2)-(d3-x4(2,1))