## Copyright (C) 2018 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
## <http://www.gnu.org/licenses/>.

## -*- 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.
##
## @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 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")

  if (isnumeric (name) && isreal (name))
    name = num2str (fix (name));
  elseif (! ischar (name))
    error ("referenceEllipsoid: value must be a string or real number");
  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/
    
    case {"airy30", "7001"}
      Code = 7001;
      Name = "Airy 1830";
      SemimajorAxis = 6377563.396;
      InverseFlattening = 299.3249646;
        
    case {"airy49", "7002"}
      Code = 7002;
      Name = "Airy Modified 1849";
      SemimajorAxis = 6377340.189;
      InverseFlattening = 299.3249646;

    case {"australia", "7003"}
      Code = 7003;
      Name = "Australian National";
      SemimajorAxis = 6378160;
      InverseFlattening = 298.25;

    case {"bessel", "7004"}
      Code = 7004;
      Name = "Bessel 1841";
      SemimajorAxis = 6377397.155;
      InverseFlattening = 299.1528128;
    
    case {"clarke66", "7008"}
      Code = 7008;
      Name = "Clarke 1866";
      SemimajorAxis = 6378206.4;
      InverseFlattening = 294.9786982;
    
    case {"clarke80", "7012"}
      Code = 7012;
      Name = "Bessel 1880";
      SemimajorAxis = 6378249.145;
      InverseFlattening = 293.465;
        
    case {"everest", "7015"}
      Code = 7015;
      Name = "Everest 1830";
      SemimajorAxis = 6377276.34518;
      InverseFlattening = 300.8017;
    
    case {"grs80", "1980", "7019"}
      Code = 7019;
      Name = "Geodetic Reference System 1980";
      SemimajorAxis = 6378137;
      InverseFlattening = 298.257222101;

    case {"int24" ,"1924", "7022"}
      Code = 7022;
      Name = "International 1924";
      SemimajorAxis = 6378388;
      InverseFlattening = 297;

    case {"kras40" ,"krasovsky", "7024"}
      Code = 7024;
      Name = "Krasovsky 1940";
      SemimajorAxis = 6378245;
      InverseFlattening = 298.3;

    case {"wgs84", "1984", "7030"}  
      Code = 7030;
      Name = "World Geodetic System 1984";
      SemimajorAxis = 6378137;
      InverseFlattening = 298.2572235630;

    case {"wgs72", "1972", "7043"}  
      Code = 7043;
      Name = "World Geodetic System 1972";
      SemimajorAxis = 6378135;
      InverseFlattening = 298.26;
    otherwise
      error (" 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 <value must be a string> referenceEllipsoid ( 7i )
%!error <not been implemented> referenceEllipsoid ( "yy" )
