WaBis

walter.bislins.ch

Method to calculate the Radius of the Earth from 3 Points

This page describes how the measured radius of the earth R3pt in the GNSS Data Viewer is calculated. The two calculation variants on this page differ by the calculated plane in which the circle given by 3 points is computed. The Geo Data App currently only uses the method Plane defined by the Perpendiculars at 2 Points.

Notation: On this page I notate vectors by bold letters like , unit vectors as .

Plane through Center of the Earth

Given three 3D points , and in cartesian coordinates, find the center and radius of the circle through this points so that the center of the earth and the points and lie in the same plane.

Plane defined by the Perpendiculars at 2 Points

ZoomConstruction of Arc Plane

Given three 3D points , and in cartesian coordinates, find the center and radius of the circle through this points so the points and lie in the plane. The plane is defined by the normalized vectors and perpendicular to the ellipsoid at the points and respectively so that the deviation angle of this normals from the plane is minimal.

In the figure the green vector is the mean vector between and and lies in the plane. The blue vector is the planes normal and is perpendicular to the green vector and the vector from to .

Note: The center of the earth lies normally not on this plane.

Problem with 3D Points

To find the center of a circle defined by 3 points, all points must lay on the same plane. The plane is calculated so that the points and lie on the plane. The middle point does not have to be exactly on the plane. It is projected perpendicular to the plane and becomes the 2D point .

Most of the calculations are the same for both plane variants. The only difference is how the local coordinate system of the plane is calculated. For the origin of the plane in both variants the point is chosen.

a) Plane through the Center of the Earth

First we define the plane by a local coordinate system with the origin at , X-axes from the center of the earth in the direction to point , Z-axes perpendicular to the vectors and and the Y-axes perpendicular to X and Z in a right handed system.

(1)

The hat above the variables like denotes unit vectors. gives the unit vector in the direction of and is the vector cross product.

b) Plane defined by the Perpendiculars at 2 Points

First we calculate the perpendicular vector for the plane as the mean vector of a perpendicular plane vector at point and at point . This z-vectors are the cross product of the vector from to and the perpendiculars and respectively.

The X-axis direction is then perpendicular to and and the Y-axes perpendicular to and in a right handed system.

Lets calculate some intermediate vectors:

(2)

Now we can calculate the plane coordinate system:

(3)

Common calculations

Construction Arc through 3 Points

Now we transform , and into this plane 2D coordinate system so that is the origin. Note that by this process the vector gets automatically projected onto the plane because the Z-component is ignored. The 2D plane coordinates of the points are denoted as , and :

(4)

The center of the circle lies on the two perpendiculars through the middle points between A'-B' and B'-C'. The middle points can be calculated by adding the vectors and scale it by a factor 1/2:

(5)

The direction of the perpendiculars through and can be calculted from the 90° rotated unit vectors from to and to :

(6)

Rotating by 90° counter clockwise using the function gives:

(7)

Now we can describe the equations for the two middle perpendicular lines. The center of the circle is the intersection of the two lines:

(8)

Because each vector in the equation (8) has 2 components, we have 2 equations for the unknown factors a and b.

(9)

We only need one of this factors to calculate the center of the circle on the plane. We can solve both equations (9) for a and subtract the results. This gives us one equation for the unknown b. Solving for b gives:

(10)

We can now insert b into (8) and get the center of the circle in plane coordinates:

(11)

The radius of the circle is the length of the vector from to , or :

(12)

The center of the circle in the original 3D coordinate system can be calculated as follows:

(13)

JavaScript Code

  CircleFrom3PointsAndNormals: function( A_3d, B_3d, C_3d, NA_3d, NC_3d ) {
    // returns { R: radius, Center: [ x, y z ] } coords in ECEF
    //
    // if NA_3d and NC_3d are defined, then
    //   the circle plane is such that A and B lie on the plane
    //   and the plane normal is calculated from the ellispoid perpendicular vectors NA_3d, NC_3d
    //   such that the plane lies between this vectors.
    // else
    //   the plane is defined by the points A_3d, B_3d, and the center of the earth

    // calculate plane local coord sys
    var x, y, z;
    if (V3.Ok(NA_3d) && V3.Ok(NC_3d)) {

      // the circle plane is such that A and B lie on the plane
      // and the plane normal is calculated from the ellispoid perpendicular vectors NA_3d, NC_3d
      // such that the plane lies between this vectors.
      // This is due to the fact that this vectors may be angled perpendicular to the connection of A and C
      var D_3d = V3.Norm( V3.Sub( C_3d, A_3d ) );
      var na = V3.Norm( V3.Mult( NA_3d, D_3d ) );
      var nc = V3.Norm( V3.Mult( NC_3d, D_3d ) );
      z = V3.Norm( V3.Add( na, nc ) );
      x = V3.Norm( V3.Mult( D_3d, z ) );
      y = V3.Mult( z, x );

    } else {

      // circle plane contains the center of the earth
      z = V3.Norm( V3.Mult( A_3d, C_3d ) );
      x = V3.Norm( A_3d );
      y = V3.Mult( z, x );

    }

    // transform points into plane coord system
    var B0_3d = V3.Sub( B_3d, A_3d );
    var C0_3d = V3.Sub( C_3d, A_3d );
    var A = [ 0, 0 ];
    var B = [ V3.ScalarProd( B0_3d, x ), V3.ScalarProd( B0_3d, y ) ];
    var C = [ V3.ScalarProd( C0_3d, x ), V3.ScalarProd( C0_3d, y ) ];

    // calculate mid points
    var Ma = V2.Scale( V2.Add( A, B ), 0.5 );
    var Mb = V2.Scale( V2.Add( B, C ), 0.5 );

    // calculate middle perpendicular directions
    var ra = V2.Norm( V2.Sub( B, A ) );
    var rb = V2.Norm( V2.Sub( C, B ) );
    // rotate counter clockwise
    var sa = [ -ra[1], ra[0] ];
    var sb = [ -rb[1], rb[0] ];

    // equations of middle perpendiculars are:
    // Z = Ma + a * sa = Mb + b * sb
    // solve for b (or a)
    var b = (sa[0] * (Mb[1] - Ma[1]) - sa[1] * (Mb[0] - Ma[0])) / (sb[0] * sa[1] - sb[1] * sa[0]);

    // center of circle on plane
    var Z = V2.Add( Mb, V2.Scale( sb, b ) );

    // radius of circle
    var R = V2.Length( V2.Sub( A, Z ) );
    if (b < 0) R *= -1;

    // transform center to 3D
    var Z0_3d = V3.Add( V3.Scale(x, Z[0]), V3.Scale(y, Z[1]) );
    Z_3d = V3.Add( Z0_3d, A_3d );

    return { R: R, Center: Z_3d };
  }

Weitere Infos zur Seite
Erzeugt Wednesday, July 17, 2019
von wabis
Zum Seitenanfang
Geändert Wednesday, July 17, 2019
von wabis