## Copyright (C) 2018 Anonymous 
##
## 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
## @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 name 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
## Also can accept a unit argument
## @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{wgs84Ellipsoid}
## @end deftypefn

function Ell = referenceEllipsoid ( name , unit )
  
  if ! nargin 
    name = 'wgs84';
    unit = 'meter';
  endif
  
  if nargin < 2
    unit = 'meter';
  endif
 
  if ( isnumeric ( name ) && isreal ( name ) )
    name = num2str ( name );
  elseif ( ! ischar ( name )) 
    error ( "referenceEllipsoid: value must be a string or number" )
  endif
  
  switch 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 { 'grs80', '1980', '7019' } 
     Ell.Code = 7019;
     Ell.Name = 'Geodetic Reference System 1980';
     SemimajorAxis = 6378137;                             
     InverseFlattening = 298.257222101;
     
    case { 'int24' , '1924', '7022' }
     Ell.Code = 7022;
     Ell.Name = 'International 1924';
     SemimajorAxis = 6378388;                             
     InverseFlattening = 297;
     
    case  { 'wgs84' , '1984', '7030' }  
     Ell.Code = 7030;
     Ell.Name = 'World Geodetic System 1984';
     SemimajorAxis = 6378137;                             
     InverseFlattening = 298.2572235630;    
     
    otherwise
     error (" Ellipsoid %s has not been implemented", name )
  endswitch
   
 #Calculations
[ Ell.LengthUnit , Ell.SemimajorAxis , Ell.SemiminorAxis , ...
  Ell.InverseFlattening , Ell.Eccentricity , Ell.Flattening , ...
  Ell.ThirdFlattening , Ell.MeanRadius , Ell.SurfaceArea , ...
  Ell.Volume ]  =  param_calc( SemimajorAxis , InverseFlattening , unit );
 
 function [ Length, Major , Minor , Inverse ,  Ecc , ... 
            Flat , Third , Mean , Surface , Vol ] = param_calc(  SMA , InvF , unit )
            
   ratio = unitsratio (unit , 'Meters');
   Length = unit;
   Major = SemimajorAxis * ratio;
   Inverse = InvF;
   Ecc = flat2ecc( 1 / InvF );
   Minor = minaxis( Major , Ecc ); 
   Flat = 1 / InvF;
   Third = ( Major - Minor ) / ( Major + Minor );
   Mean = ( 2 * Major + Minor ) / 3;
   P = 1.6075; # From Knud Thomsen this results in a max error of 1.061%
   Surface = 4 * pi * ( ( Major ^ ( 2 * P) +  2 * ( Major * Minor ) ^ P ) / 3 ) ^ ( 1 / P );
   Vol = ( 4 * pi ) / 3 * Major ^ 2 * Minor;
  end
 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' )