function

Geodesy.latLonToLocalNorthEast

Great-circle inverse of localNorthEastToLatLon

Public declarations

Explicit source declarations; dimensions, modifications and defaults are shown as written, without evaluation.

Role	Name	Type	Dimensions / attributes / default	Description
input	`lat0_deg`	`Real`		Reference latitude [deg]
input	`lon0_deg`	`Real`		Reference longitude [deg]
input	`lat_deg`	`Real`		Target latitude [deg]
input	`lon_deg`	`Real`		Target longitude [deg]
input	`earth_radius_m`	`Real`	`= 6378137.0`	Spherical Earth radius [m]
output	`north_east_m`	`Real`	`[2]`	Local north and east displacement [m]

Modelica declarations and implementation

Show parameters, interfaces and source

Parameters, inputs, outputs, defaults and equations from the source file.

function latLonToLocalNorthEast
    "Great-circle inverse of localNorthEastToLatLon"
    input Real lat0_deg "Reference latitude [deg]";
    input Real lon0_deg "Reference longitude [deg]";
    input Real lat_deg "Target latitude [deg]";
    input Real lon_deg "Target longitude [deg]";
    input Real earth_radius_m = 6378137.0 "Spherical Earth radius [m]";
    output Real north_east_m[2] "Local north and east displacement [m]";
  protected
    constant Real pi = 3.141592653589793;
    Real lat0;
    Real lon0;
    Real lat;
    Real lon;
    Real dLon;
    Real haversine;
    Real centralAngle;
    Real distance_m;
    Real bearingY;
    Real bearingX;
    Real bearing;

  algorithm
    lat0 := lat0_deg * pi / 180.0;
    lon0 := lon0_deg * pi / 180.0;
    lat := lat_deg * pi / 180.0;
    lon := lon_deg * pi / 180.0;
    dLon := lon - lon0;

    // Central angle via the haversine, clamped to a valid arcsine domain.
    haversine := sin((lat - lat0) / 2.0) * sin((lat - lat0) / 2.0)
      + cos(lat0) * cos(lat) * sin(dLon / 2.0) * sin(dLon / 2.0);
    haversine := max(0.0, min(1.0, haversine));
    centralAngle := 2.0 * asin(sqrt(haversine));

    if centralAngle <= 1.0e-12 then
      north_east_m := {0.0, 0.0};
    else
      // Initial bearing from the origin to the target, then decompose the
      // arc length into north and east so that the forward projection above
      // recovers the target exactly.
      bearingY := sin(dLon) * cos(lat);
      bearingX := cos(lat0) * sin(lat) - sin(lat0) * cos(lat) * cos(dLon);
      bearing := atan2(bearingY, bearingX);
      distance_m := earth_radius_m * centralAngle;
      north_east_m := {distance_m * cos(bearing), distance_m * sin(bearing)};
    end if;
  end latLonToLocalNorthEast;