## Copyright (C) 2018-2020 Philip Nienhuis
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## This program is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSEll. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; see the file COPYING. If not, see
## .
## -*- texinfo -*-
## @deftypefn {Function File} {} referenceEllipsoid (@var{name}, @var{unit})
## Returns the parameters of an ellipsoid.
##
## @var{Name} can be the name (e.g., "wgs84") or (integer) EPSG code of a
## reference ellipoid. Case is not important. If the number or code 0
## (zero) is specified, referenceEllipsoid echoes a list of implemented
## ellipsoids to the screen.
##
## @var{unit} can be the name of any unit accepted by function
## validateLengthUnit.m. Also here case is not important.
##
## The output consists of a scalar struct with fields "Code" (EPSG code of the
## ellipsoid), "Name" (name of the ellipsoid), "LengthUnit", "SemimajorAxis",
## "SemiminorAxis", "InverseFlattening", "Eccentricity", "Flattening",
## "ThirdFlattening", "MeanRadius", "SurfaceArea", and "Volume".
##
## Examples:
##
## @example
## >> E = referenceEllipsoid ("wgs84")
## E =
##
## scalar structure containing the fields:
##
## Code = 7030
## Name = World Geodetic System 1984
## LengthUnit = meter
## SemimajorAxis = 6378137
## SemiminorAxis = 6.3568e+06
## InverseFlattening = 298.26
## Eccentricity = 0.081819
## Flattening = 0.0033528
## ThirdFlattening = 0.0016792
## MeanRadius = 6.3710e+06
## SurfaceArea = 5.1007e+14
## Volume = 1.0832e+21
## @end example
##
## The code number can be used:
##
## @example
## >> E = referenceEllipsoid (7019)
## E =
##
## scalar structure containing the fields:
##
## Code = 7019
## Name = Geodetic Reference System 1980
## LengthUnit = meter
## SemimajorAxis = 6.3781e+06
## SemiminorAxis = 6.3568e+06
## InverseFlattening = 298.26
## Eccentricity = 0.081819
## Flattening = 0.0033528
## ThirdFlattening = 0.0016792
## MeanRadius = 6.371e+06
## SurfaceArea = 5.1007e+14
## Volume = 1.0832e+21
## @end example
##
## @example
## >> E = referenceEllipsoid ("wgs84", "km")
## E =
##
## scalar structure containing the fields:
##
## Code = 7030
## Name = World Geodetic System 1984
## LengthUnit = km
## SemimajorAxis = 6378.1
## SemiminorAxis = 6356.8
## InverseFlattening = 298.26
## Eccentricity = 0.081819
## Flattening = 0.0033528
## ThirdFlattening = 0.0016792
## MeanRadius = 6371.0
## SurfaceArea = 5.1007e+08
## Volume = 1.0832e+12
## @end example
##
## @seealso{validateLengthUnit, wgs84Ellipsoid}
## @end deftypefn
## Function supplied by anonymous contributor, see:
## https://savannah.gnu.org/patch/index.php?9634
function Ell = referenceEllipsoid (name="wgs84", unit="meter")
## List of implemented codes. To be updated if codes are added.
persistent codes;
if (isempty (codes))
codes = {"7001", "Airy30", " "; ...
"7002", "Airy49", " "; ...
"7003", "Australia", " "; ...
"7004", "Bessel", " "; ...
"7008", "Clarke66", " "; ...
"7012", "Bessel 1880", " "; ...
"7015", "Everest 1830", " "; ...
"7019", "grs80", "1980"; ...
"7022", "int24", "1924"; ...
"7024", "Kras40", "Krasovsky"; ...
"7030", "Wgs84", "1984"; ...
"7043", "Wgs72", "1972"; ...
" ", "Sun", " "; ...
" ", "Mercury", " "; ...
" ", "Venus", " "; ...
" ", "Earth", " "; ...
" ", "Moon", " "; ...
" ", "Mars", " "; ...
" ", "Jupiter", " "; ...
" ", "Saturn", " "; ...
" ", "Uranus", " "; ...
" ", "Neptune", " "; ...
" ", "Pluto", " "};
endif
if (isnumeric (name) && isreal (name))
name = num2str (fix (name));
elseif (! ischar (name))
error ("referenceEllipsoid: value must be a string or real number");
elseif (strcmp (name, " "))
error ("referenceEllipsoid: name required");
endif
if (! ischar (unit))
error ("referenceEllipsoid: length name expected for input arg. #2");
endif
switch lower (name)
## Semimajor axis and Inverse flattening from
## USER's HANDBOOK ON DATUM TRANSFORMATIONS INVOLVING WGS 84
## 3rd Edition, July 2003, (Last correction August 2008),
## Special Publication No. 60
## Codenames are from https://epsg.io/
## Planet values from Report of the IAU Working Group on
## CartographicCoordinates and Rotational Elements: 2015
case lower (codes(1, :))
Code = 7001;
Name = "Airy 1830";
SemimajorAxis = 6377563.396;
InverseFlattening = 299.3249646;
case lower (codes (2, :))
Code = 7002;
Name = "Airy Modified 1849";
SemimajorAxis = 6377340.189;
InverseFlattening = 299.3249646;
case lower (codes (3, :))
Code = 7003;
Name = "Australian National";
SemimajorAxis = 6378160;
InverseFlattening = 298.25;
case lower (codes (4, :))
Code = 7004;
Name = "Bessel 1841";
SemimajorAxis = 6377397.155;
InverseFlattening = 299.1528128;
case lower (codes (5, :))
Code = 7008;
Name = "Clarke 1866";
SemimajorAxis = 6378206.4;
InverseFlattening = 294.9786982;
case lower (codes (6, :))
Code = 7012;
Name = "Bessel 1880";
SemimajorAxis = 6378249.145;
InverseFlattening = 293.465;
case lower (codes (7, :))
Code = 7015;
Name = "Everest 1830";
SemimajorAxis = 6377276.34518;
InverseFlattening = 300.8017;
case lower (codes (8, :))
Code = 7019;
Name = "Geodetic Reference System 1980";
SemimajorAxis = 6378137;
InverseFlattening = 298.257222101;
case lower (codes (9, :))
Code = 7022;
Name = "International 1924";
SemimajorAxis = 6378388;
InverseFlattening = 297;
case lower (codes (10, :))
Code = 7024;
Name = "Krasovsky 1940";
SemimajorAxis = 6378245;
InverseFlattening = 298.3;
case lower (codes (11, :))
Code = 7030;
Name = "World Geodetic System 1984";
SemimajorAxis = 6378137;
InverseFlattening = 298.2572235630;
case lower (codes (12, :))
Code = 7043;
Name = "World Geodetic System 1972";
SemimajorAxis = 6378135;
InverseFlattening = 298.26;
case lower (codes (13, :))
Code = " ";
Name = "Sun";
SemimajorAxis = 695700000;
InverseFlattening = 111111;
case lower (codes (14, :))
Code = " ";
Name = "Mercury";
SemimajorAxis = 2440530;
InverseFlattening = 1075;
case lower (codes (15, :))
Code = " ";
Name = "Venus";
SemimajorAxis = 6051800;
InverseFlattening = Inf;
case lower (codes (16, :))
Code = " ";
Name = "Earth";
SemimajorAxis = 6378137;
InverseFlattening = 298.2572235630;
case lower (codes (17, :))
Code = " ";
Name = "Moon";
SemimajorAxis = 1738100;
InverseFlattening = 833.33;
case lower (codes (18, :))
Code = " ";
Name = "Mars";
SemimajorAxis = 3396190;
InverseFlattening = 169.894;
case lower (codes (19, :))
Code = " ";
Name = "Jupiter";
SemimajorAxis = 71492000;
InverseFlattening = 15.4144;
case lower (codes (20, :))
Code = " ";
Name = "Saturn";
SemimajorAxis = 60268000;
InverseFlattening = 10.208;
case lower (codes (21, :))
Code = " ";
Name = "Uranus";
SemimajorAxis = 25559000;
InverseFlattening = 43.616;
case lower (codes (22, :))
Code = " ";
Name = "Neptune";
SemimajorAxis = 24764000;
InverseFlattening = 58.5437;
case lower (codes (23, :))
Code = " ";
Name = "Pluto";
SemimajorAxis = 1188300;
InverseFlattening = INF;
case "0"
## Show list of codes
printf ("\n referenceEllipsoid.m:\n list of implemented ellipsoids:\n");
printf (" Code Alias 1 Alias 2\n");
printf (" ==== ======= =======\n");
for ii=1:size (codes, 1)
printf ("%5s %15s %15s\n", codes (ii, :){:});
endfor
return
otherwise
error ("referenceEllipsoid: ellipsoid %s has not been implemented", name)
endswitch
## Calculations
Ell = param_calc (Code, Name, SemimajorAxis, InverseFlattening, unit);
endfunction
function ell = param_calc (Code, Name, SemimajorAxis, InvF, unit)
ell.Code = Code;
ell.Name = Name;
ratio = unitsratio (unit, "Meters");
ell.LengthUnit = unit;
Major = SemimajorAxis * ratio;
ell.SemimajorAxis = Major;
Inverse = InvF;
ell.InverseFlattening = InvF;
Ecc = flat2ecc (1 / InvF);
ell.Eccentricity = Ecc;
Minor = minaxis (Major, Ecc);
ell.SemiminorAxis = Minor;
ell.Flattening = 1 / InvF;
ell.ThirdFlattening = (Major - Minor) / (Major + Minor);
ell.MeanRadius = (2 * Major + Minor) / 3;
## From Knud Thomsen this results in a max error of 1.061 %:
P = 1.6075;
Surface = 4 * pi * ((Major^(2*P) + 2 * (Major * Minor)^P) / 3 )^(1/P);
ell.SurfaceArea = Surface;
ell.Volume = (4 * pi) / 3 * Major^2 * Minor;
endfunction
%!test
%!
%! E = referenceEllipsoid ("wgs84");
%! assert ( E.SemiminorAxis, 6356752.314245, 10e-7 )
%! assert ( E.Eccentricity, 0.081819221456, 10e-8)
%! assert ( E.Flattening, 1 / 298.2572235630, 10e-8 )
%!error referenceEllipsoid ( 7i )
%!error referenceEllipsoid ( "yy" )
%!error referenceEllipsoid ( " " )