## 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} {@var{r} =} rcurve (@var{type} , @var{spheroid} , @var{lat}, @var{angleUnit})
## This returns the length of the curve based on which type: meridian, parallel (default), or transverse
## With a spheroid (default: wgs84) and a latitude (default: degrees)
## 
##
## Examples:
##
## @example
## r = rcurve ( 'parallel', 'wgs84' , 45)
## => r =
## 4.5176e+06
## Note: this is in meters
## @end example 
## 
## @example
## r = rcurve ( '', '' , 45)
## => r =
## 4.5176e+06
## @end example
##
## @example
## r = rcurve ( '', '' , pi/4, 'radians')
## => r =
## 4.5176e+06
## @end example
##
## @example
## r = rcurve ( 'meridian', 'wgs84' , 45)
## => r =
##  6.3514e+06
## @end example
##
## @example
## r = rcurve ( 'transverse', 'wgs84' , 45)
## => r =
##  6.3888e+06
## @end example
## @end deftypefn
function r = rcurve ( type , spheroid , lat , angleUnit )
  
  if isempty ( type )
    type = 'parallel';
  endif
  
  if isempty ( spheroid )
    E = wgs84Ellipsoid;
  else
    E = referenceEllipsoid ( spheroid );
  endif 
  
  if (nargin < 4)
    angleUnit = "degrees";  
  elseif ( nargin < 3 )
    print_usage ();
  endif
  
  if ( strncmpi (lower (angleUnit), "r" , 1 ) == 0 && strncmpi (lower (angleUnit), "d", 1 ) == 0 )
    error ( "rcurve.m: angleUnit %s is not supported", angleUnit)
  endif
  
  if (! isnumeric ( lat ) || ! isreal ( lat ) )
    error ( "rcurve.m : numeric input expected");
  elseif ( strncmpi (lower (angleUnit), "r", 1) == 1 )
    c_l = cos ( lat );
    s_l = sin ( lat );
  else
    c_l = cosd ( lat );
    s_l = sind ( lat );
  endif
  
  #Insight From: Algorithms for Global Positioning pg 370-372
  
  e2 = E.Eccentricity ^ 2;
  R = E.SemimajorAxis;
  e_p = e2 / ( 1 - e2 );
  
  N = ( R * sqrt ( 1 + e_p ) / ( sqrt ( 1 + e_p * c_l ^ 2)));
  switch type
    
    case { 'meridian' }
    w = sqrt ( 1 - e2 * s_l ^ 2 );
    r = R * ( 1 - e2 ) / ( w ^ ( 3 / 2 ) );
    case { 'parallel' }
    r = N * c_l;
    case { 'transverse' }    
    r = N;
    otherwise
     error (" Type %s does not exist", type )
  endswitch
  
endfunction

%!test
%!
%! assert ( rcurve ( '' , '' , 45) ,  4517590.87885 , 10e-6)
%! assert ( rcurve ( 'transverse' , '' , 45) , 6388838.29012 , 10e-6)
%! assert ( rcurve ( 'meridian' , '' , 45) , 6351390.490795 , 10e-6)

%!error <angleUnit> rcurve ('', '' , 45 ,  "km" )
%!error <numeric input expected>  rcurve ( '', '', "A")
%!error <numeric input expected>  rcurve ( '', '', 45i)
%!error <does not exist> rcurve ( 'All' , '' , 45 )