%Program Awange-Grafarend Groebner Basis 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. clear all syms x1 x2 x3 x4 a0 a1 a2 a3 b0 b1 b2 b3 c0 c1 c2 c3 d0 d1 d2 d3 a03 a13 a23 b03 b13 b23 b3 c03 c13 c23 d30 d31 d32 syms e03 e13 e23 format long e % Symbolic computation of the stationary receiver coordinates {z=x3},{y=x2} and {x=x1} cx3=(-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23-a23*b13*d30*x4+a13*b23*d30*x4+... a23*b03*d31*x4-a03*b23*d31*x4-a13*b03*d32*x4+a03*b13*d32*x4)/(a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+... a13*b03*c23-a03*b13*c23); cx2=(-a23*e13+a13*e23-a23*c13*x3+a13*c23*x3-a23*d31*x4+a13*d32*x4)/(a23*b13-a13*b23); cx1=(e23+b23*x2+c23*x3+d32*x4)/(-a23); % Symbolic computation of the Coefficients of the quadratic equation cf2=((-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)^2/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)^2+... (-a23*c13*(-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+a13*c23*... (-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)-a23*d31+a13*d32)^2/... (a23*b13-a13*b23)^2+(b23*(-a23*c13*(-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+a13*c23*... (-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)-a23*d31+a13*d32)/... (a23*b13-a13*b23)+c23*(-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+d32)^2/a23^2-1);%*x4^2+ cf1=(2*((-a23*e13+a13*e23-a23*c13*(-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+a13*c23*... (-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23))/... (a23*b13-a13*b23)-b0)*(-a23*c13*(-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+a13*c23*... (-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)-a23*d31+a13*d32)/... (a23*b13-a13*b23)-2*(-(e23+b23*(-a23*e13+a13*e23-a23*c13*... (-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+a13*c23*... (-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23))/(a23*b13-a13*b23)+c23*... (-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23))/a23-a0)*... (b23*(-a23*c13*(-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+a13*c23*... (-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)-a23*d31+a13*d32)/... (a23*b13-a13*b23)+c23*(-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+d32)/a23+2*d0+2*... ((-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)-c0)*... (-a23*b13*d30+a13*b23*d30+a23*b03*d31-a03*b23*d31-a13*b03*d32+a03*b13*d32)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23));%*x4+ cf0=(-(e23+b23*(-a23*e13+a13*e23-a23*c13*(-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-... a13*b03*e23+a03*b13*e23)/(a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)+... a13*c23*(-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23))/(a23*b13-a13*b23)+c23*... (-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23))/a23-a0)^2-d0^2+... ((-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23)-c0)^2+... ((-a23*e13+a13*e23-a23*c13*(-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-... a13*b03*e23+a03*b13*e23)/(a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-... a03*b13*c23)+a13*c23*(-a23*b13*e03+a13*b23*e03+a23*b03*e13-a03*b23*e13-a13*b03*e23+a03*b13*e23)/... (a23*b13*c03-a13*b23*c03-a23*b03*c13+a03*b23*c13+a13*b03*c23-a03*b13*c23))/(a23*b13-a13*b23)-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))