Verifying Spherical Coordinate Methods?
in [Math]
|
Prev: Does "iff" make sense in a definition?
Next: Al Gore, Jr. and Senior's key backer was Armand Hammer; nowadays, perhaps Beyond Petroleum
From: KBH on 9 Aug 2010 07:04 On Jul 26, 8:38 pm, KBH <emptyp...(a)hotmail.com> wrote: > Converting to local topocentric moves the pole to the user's location > and that might be related. > > Here is KBH code: > > {Latitude and longitude to ECEF xyz} > {KBH Code} > ra:= 180 / Pi; > a:= 6378137; > ee:= 0.006694379990; > Write(' Input latitude in degrees: '); > ReadLn(lat); > Write(' Input longitude in degrees: '); > ReadLn(lon); > Write(' Input ellipsoidal height in meters: '); > ReadLn(h); > lat:= lat / ra; > lon:= lon / ra; > v:= a / Sqrt(1 - (ee * Sqr(Sin(lat)))); > x:= (v + h) * Cos(lat) * Cos(lon); > y:= (v + h) * Cos(lat) * Sin(lon); > z:= ((1 - ee) * v + h) * Sin(lat); > > ---------------------------------------------------- > > {ECEF xyz to local topocentric} > {KBH Code} > nr:= (-x * Sin(lat) * Cos(lon)) + (-y * Sin(lat) * Sin(lon)) + (z * > Cos(lat)); > ea:= (-x * Sin(lon)) + (y * Cos(lon)); > ht:= (x * Cos(lat) * Cos(lon)) + (y * Cos(lat) * Sin(lon)) + (z * > Sin(lat)); > > {Note the first input is the lat and lon of the first xyz point > converted to > local topocentric and then that lat and lon is held in the > computations for > all the next xyz input points going to the same local topocentric > system} > > ----------------------------------------------------- > > {ECEF xyz to latitude and longitude} > {KBH Code} > ra:= 180 / Pi; > a:= 6378137; > f:= 0.003352811; > fm:= 1 - f; > ee:= 0.006694379990; > WriteLn; > Write(' Input x coordinate in meters: '); > ReadLn(x); > Write(' Input y coordinate in meters: '); > ReadLn(y); > Write(' Input z coordinate in meters: '); > ReadLn(z); > If (x = 0) Then x:= 0.0001; > {If (x > a) Or (y > a) Then Exit; > If (x < -a) Or (y < -a) Then Exit;} > p:= Sqrt(Sqr(x) + Sqr(y)); > r:= Sqrt(Sqr(p) + Sqr(z)); > u:= z / p * (fm + (ee * a / r)); > u:= ArcTan(u); > lon:= ArcTan(y / x); > la:= (z * fm) + (ee * a * Sin(u) * Sqr(Sin(u))); > lt:= fm * (p - (ee * a * Cos(u) * Sqr(Cos(u)))); > lat:= ArcTan(la / lt); > ha:= p * Cos(lat) + (z * Sin(lat)); > hb:= a * Sqrt(1 - (ee * Sqr(Sin(lat)))); > h:= ha - hb; > lon:= lon * ra; > lat:= lat * ra; > If (x < 0) Then lon:= 180 + lon; > If (lon > 180) Then lon:= lon - 360; > > ------------------------------------------------------ The code that I posted converts latitude and longitude to ECEF XYZ and then on to something one source calls a "Local Tangent Plane". It makes the ECEF y-coordinate of the first point as 0 and makes the z- coordinate near maximum. Using simple plane formulas there is very little direction error between two points in 135,000 meters while the distance error is about 1 meter in 13,500.
From: KBH on 9 Aug 2010 07:29 On Aug 9, 7:04 am, KBH <emptyp...(a)hotmail.com> wrote: > On Jul 26, 8:38 pm, KBH <emptyp...(a)hotmail.com> wrote: > > > > > > > Converting to local topocentric moves the pole to the user's location > > and that might be related. > > > Here is KBH code: > > > {Latitude and longitude to ECEF xyz} > > {KBH Code} > > ra:= 180 / Pi; > > a:= 6378137; > > ee:= 0.006694379990; > > Write(' Input latitude in degrees: '); > > ReadLn(lat); > > Write(' Input longitude in degrees: '); > > ReadLn(lon); > > Write(' Input ellipsoidal height in meters: '); > > ReadLn(h); > > lat:= lat / ra; > > lon:= lon / ra; > > v:= a / Sqrt(1 - (ee * Sqr(Sin(lat)))); > > x:= (v + h) * Cos(lat) * Cos(lon); > > y:= (v + h) * Cos(lat) * Sin(lon); > > z:= ((1 - ee) * v + h) * Sin(lat); > > > ---------------------------------------------------- > > > {ECEF xyz to local topocentric} > > {KBH Code} > > nr:= (-x * Sin(lat) * Cos(lon)) + (-y * Sin(lat) * Sin(lon)) + (z * > > Cos(lat)); > > ea:= (-x * Sin(lon)) + (y * Cos(lon)); > > ht:= (x * Cos(lat) * Cos(lon)) + (y * Cos(lat) * Sin(lon)) + (z * > > Sin(lat)); > > > {Note the first input is the lat and lon of the first xyz point > > converted to > > local topocentric and then that lat and lon is held in the > > computations for > > all the next xyz input points going to the same local topocentric > > system} > > > ----------------------------------------------------- > > > {ECEF xyz to latitude and longitude} > > {KBH Code} > > ra:= 180 / Pi; > > a:= 6378137; > > f:= 0.003352811; > > fm:= 1 - f; > > ee:= 0.006694379990; > > WriteLn; > > Write(' Input x coordinate in meters: '); > > ReadLn(x); > > Write(' Input y coordinate in meters: '); > > ReadLn(y); > > Write(' Input z coordinate in meters: '); > > ReadLn(z); > > If (x = 0) Then x:= 0.0001; > > {If (x > a) Or (y > a) Then Exit; > > If (x < -a) Or (y < -a) Then Exit;} > > p:= Sqrt(Sqr(x) + Sqr(y)); > > r:= Sqrt(Sqr(p) + Sqr(z)); > > u:= z / p * (fm + (ee * a / r)); > > u:= ArcTan(u); > > lon:= ArcTan(y / x); > > la:= (z * fm) + (ee * a * Sin(u) * Sqr(Sin(u))); > > lt:= fm * (p - (ee * a * Cos(u) * Sqr(Cos(u)))); > > lat:= ArcTan(la / lt); > > ha:= p * Cos(lat) + (z * Sin(lat)); > > hb:= a * Sqrt(1 - (ee * Sqr(Sin(lat)))); > > h:= ha - hb; > > lon:= lon * ra; > > lat:= lat * ra; > > If (x < 0) Then lon:= 180 + lon; > > If (lon > 180) Then lon:= lon - 360; > > > ------------------------------------------------------ > > The code that I posted converts latitude and longitude to ECEF XYZ and > then on to something one source calls a "Local Tangent Plane". It > makes the ECEF y-coordinate of the first point as 0 and makes the z- > coordinate near maximum. Using simple plane formulas there is very > little direction error between two points in 135,000 meters while the > distance error is about 1 meter in 13,500.- Hide quoted text - > > - Show quoted text - Here's a link: http://en.wikipedia.org/wiki/Geodetic_system Scroll down to "From ECEF to ENU" and they work with the difference in a pair of points at a time rather than transform one point at a time.
From: KBH on 10 Aug 2010 04:03 On Aug 9, 7:29 am, KBH <emptyp...(a)hotmail.com> wrote: > On Aug 9, 7:04 am, KBH <emptyp...(a)hotmail.com> wrote: > > > > > > > On Jul 26, 8:38 pm, KBH <emptyp...(a)hotmail.com> wrote: > > > > Converting to local topocentric moves the pole to the user's location > > > and that might be related. > > > > Here is KBH code: > > > > {Latitude and longitude to ECEF xyz} > > > {KBH Code} > > > ra:= 180 / Pi; > > > a:= 6378137; > > > ee:= 0.006694379990; > > > Write(' Input latitude in degrees: '); > > > ReadLn(lat); > > > Write(' Input longitude in degrees: '); > > > ReadLn(lon); > > > Write(' Input ellipsoidal height in meters: '); > > > ReadLn(h); > > > lat:= lat / ra; > > > lon:= lon / ra; > > > v:= a / Sqrt(1 - (ee * Sqr(Sin(lat)))); > > > x:= (v + h) * Cos(lat) * Cos(lon); > > > y:= (v + h) * Cos(lat) * Sin(lon); > > > z:= ((1 - ee) * v + h) * Sin(lat); > > > > ---------------------------------------------------- > > > > {ECEF xyz to local topocentric} > > > {KBH Code} > > > nr:= (-x * Sin(lat) * Cos(lon)) + (-y * Sin(lat) * Sin(lon)) + (z * > > > Cos(lat)); > > > ea:= (-x * Sin(lon)) + (y * Cos(lon)); > > > ht:= (x * Cos(lat) * Cos(lon)) + (y * Cos(lat) * Sin(lon)) + (z * > > > Sin(lat)); > > > > {Note the first input is the lat and lon of the first xyz point > > > converted to > > > local topocentric and then that lat and lon is held in the > > > computations for > > > all the next xyz input points going to the same local topocentric > > > system} > > > > ----------------------------------------------------- > > > > {ECEF xyz to latitude and longitude} > > > {KBH Code} > > > ra:= 180 / Pi; > > > a:= 6378137; > > > f:= 0.003352811; > > > fm:= 1 - f; > > > ee:= 0.006694379990; > > > WriteLn; > > > Write(' Input x coordinate in meters: '); > > > ReadLn(x); > > > Write(' Input y coordinate in meters: '); > > > ReadLn(y); > > > Write(' Input z coordinate in meters: '); > > > ReadLn(z); > > > If (x = 0) Then x:= 0.0001; > > > {If (x > a) Or (y > a) Then Exit; > > > If (x < -a) Or (y < -a) Then Exit;} > > > p:= Sqrt(Sqr(x) + Sqr(y)); > > > r:= Sqrt(Sqr(p) + Sqr(z)); > > > u:= z / p * (fm + (ee * a / r)); > > > u:= ArcTan(u); > > > lon:= ArcTan(y / x); > > > la:= (z * fm) + (ee * a * Sin(u) * Sqr(Sin(u))); > > > lt:= fm * (p - (ee * a * Cos(u) * Sqr(Cos(u)))); > > > lat:= ArcTan(la / lt); > > > ha:= p * Cos(lat) + (z * Sin(lat)); > > > hb:= a * Sqrt(1 - (ee * Sqr(Sin(lat)))); > > > h:= ha - hb; > > > lon:= lon * ra; > > > lat:= lat * ra; > > > If (x < 0) Then lon:= 180 + lon; > > > If (lon > 180) Then lon:= lon - 360; > > > > ------------------------------------------------------ > > > The code that I posted converts latitude and longitude to ECEF XYZ and > > then on to something one source calls a "Local Tangent Plane". It > > makes the ECEF y-coordinate of the first point as 0 and makes the z- > > coordinate near maximum. Using simple plane formulas there is very > > little direction error between two points in 135,000 meters while the > > distance error is about 1 meter in 13,500.- Hide quoted text - > > > - Show quoted text - > > Here's a link: > > http://en.wikipedia.org/wiki/Geodetic_system > > Scroll down to "From ECEF to ENU" and they work with the difference in > a pair of points at a time rather than transform one point at a time.- Hide quoted text - > > - Show quoted text - The Local Topocentric code that I posted works correctly with spherical ECEF coordinates. That's like v = 6366706.84 and z = 6366706.84 * Sin(lat) in the "Latitude and longitude to ECEF xyz". Then the Local Topocentric code makes the first local point the pole with an x-coordinate of 0 and a y-coordinate of 0. Then the z- coordinate is the maximum height. A forward x,y point can be produced on the Local Topocentric grid with simple plane grid formulas. An inverse distance of an existing point could be matched or a difference between a cord and arc on a 6366706.84 radius could be allowed for. (The grid directions can have spherical excess so the Local Topocentric should be as local as possible.) The z-coordinate of the forward point on the Local Topocentric grid would be the Square Root of [6366706.84^2 - x^2] . Then the inverse matrix could take the forward point to ECEF spherical xyz and from there to latitude and longitude.
From: KBH on 10 Aug 2010 04:37
On Aug 10, 4:03 am, KBH <emptyp...(a)hotmail.com> wrote: > On Aug 9, 7:29 am, KBH <emptyp...(a)hotmail.com> wrote: > > > > > > > On Aug 9, 7:04 am, KBH <emptyp...(a)hotmail.com> wrote: > > > > On Jul 26, 8:38 pm, KBH <emptyp...(a)hotmail.com> wrote: > > > > > Converting to local topocentric moves the pole to the user's location > > > > and that might be related. > > > > > Here is KBH code: > > > > > {Latitude and longitude to ECEF xyz} > > > > {KBH Code} > > > > ra:= 180 / Pi; > > > > a:= 6378137; > > > > ee:= 0.006694379990; > > > > Write(' Input latitude in degrees: '); > > > > ReadLn(lat); > > > > Write(' Input longitude in degrees: '); > > > > ReadLn(lon); > > > > Write(' Input ellipsoidal height in meters: '); > > > > ReadLn(h); > > > > lat:= lat / ra; > > > > lon:= lon / ra; > > > > v:= a / Sqrt(1 - (ee * Sqr(Sin(lat)))); > > > > x:= (v + h) * Cos(lat) * Cos(lon); > > > > y:= (v + h) * Cos(lat) * Sin(lon); > > > > z:= ((1 - ee) * v + h) * Sin(lat); > > > > > ---------------------------------------------------- > > > > > {ECEF xyz to local topocentric} > > > > {KBH Code} > > > > nr:= (-x * Sin(lat) * Cos(lon)) + (-y * Sin(lat) * Sin(lon)) + (z * > > > > Cos(lat)); > > > > ea:= (-x * Sin(lon)) + (y * Cos(lon)); > > > > ht:= (x * Cos(lat) * Cos(lon)) + (y * Cos(lat) * Sin(lon)) + (z * > > > > Sin(lat)); > > > > > {Note the first input is the lat and lon of the first xyz point > > > > converted to > > > > local topocentric and then that lat and lon is held in the > > > > computations for > > > > all the next xyz input points going to the same local topocentric > > > > system} > > > > > ----------------------------------------------------- > > > > > {ECEF xyz to latitude and longitude} > > > > {KBH Code} > > > > ra:= 180 / Pi; > > > > a:= 6378137; > > > > f:= 0.003352811; > > > > fm:= 1 - f; > > > > ee:= 0.006694379990; > > > > WriteLn; > > > > Write(' Input x coordinate in meters: '); > > > > ReadLn(x); > > > > Write(' Input y coordinate in meters: '); > > > > ReadLn(y); > > > > Write(' Input z coordinate in meters: '); > > > > ReadLn(z); > > > > If (x = 0) Then x:= 0.0001; > > > > {If (x > a) Or (y > a) Then Exit; > > > > If (x < -a) Or (y < -a) Then Exit;} > > > > p:= Sqrt(Sqr(x) + Sqr(y)); > > > > r:= Sqrt(Sqr(p) + Sqr(z)); > > > > u:= z / p * (fm + (ee * a / r)); > > > > u:= ArcTan(u); > > > > lon:= ArcTan(y / x); > > > > la:= (z * fm) + (ee * a * Sin(u) * Sqr(Sin(u))); > > > > lt:= fm * (p - (ee * a * Cos(u) * Sqr(Cos(u)))); > > > > lat:= ArcTan(la / lt); > > > > ha:= p * Cos(lat) + (z * Sin(lat)); > > > > hb:= a * Sqrt(1 - (ee * Sqr(Sin(lat)))); > > > > h:= ha - hb; > > > > lon:= lon * ra; > > > > lat:= lat * ra; > > > > If (x < 0) Then lon:= 180 + lon; > > > > If (lon > 180) Then lon:= lon - 360; > > > > > ------------------------------------------------------ > > > > The code that I posted converts latitude and longitude to ECEF XYZ and > > > then on to something one source calls a "Local Tangent Plane". It > > > makes the ECEF y-coordinate of the first point as 0 and makes the z- > > > coordinate near maximum. Using simple plane formulas there is very > > > little direction error between two points in 135,000 meters while the > > > distance error is about 1 meter in 13,500.- Hide quoted text - > > > > - Show quoted text - > > > Here's a link: > > >http://en.wikipedia.org/wiki/Geodetic_system > > > Scroll down to "From ECEF to ENU" and they work with the difference in > > a pair of points at a time rather than transform one point at a time.- Hide quoted text - > > > - Show quoted text - > > The Local Topocentric code that I posted works correctly with > spherical ECEF coordinates. That's like v = 6366706.84 and z = > 6366706.84 * Sin(lat) in the "Latitude and longitude to ECEF xyz". > Then the Local Topocentric code makes the first local point the pole > with an x-coordinate of 0 and a y-coordinate of 0. Then the z- > coordinate is the maximum height. > > A forward x,y point can be produced on the Local Topocentric grid with > simple plane grid formulas. An inverse distance of an existing point > could be matched or a difference between a cord and arc on a > 6366706.84 radius could be allowed for. (The grid directions can have > spherical excess so the Local Topocentric should be as local as > possible.) The z-coordinate of the forward point on the Local > Topocentric grid would be the Square Root of [6366706.84^2 - x^2] . > Then the inverse matrix could take the forward point to ECEF spherical > xyz and from there to latitude and longitude.- Hide quoted text - > > - Show quoted text - Oh, the conversion from spherical ECEF to Local Topocentric would be trivial if only one point were going to the Local Topocentric grid and then additional points produced on the grid itself. But if two or more points are being converted from spherical ECEF to Local Topocentric then the conversion matrix or formulas are needed. But convert each point at a time not a difference in two points. |