--- a
+++ b/inst/op_geoc2geod.m
@@ -0,0 +1,183 @@
+## Copyright (C) 2009, José Luis García Pallero, <jgpallero@gmail.com>
+##
+## This file is part of OctPROJ.
+##
+## OctPROJ 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 PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File}{[@var{lon},@var{lat},@var{h}] =}op_geoc2geod(@var{X},@var{Y},@var{Z},@var{a},@var{f})
+##
+## This function converts cartesian tridimensional geocentric coordinates into
+## geodetic coordinates using the PROJ.4 function pj_geocentric_to_geodetic().
+##
+## @var{X} is a scalar or a vector containing the X geocentric coordinate.
+## @var{Y} is a scalar or a vector containing the Y geocentric coordinate.
+## @var{Z} is a scalar or a vector containing the Z geocentric coordinate.
+## @var{a} is a scalar containing the semi-major axis of the ellipsoid.
+## @var{f} is a scalar containing the flattening of the ellipsoid.
+##
+## If more than one of @var{X}, @var{Y} or @var{Z} are vectors, they must be of
+## the same size.
+## The units of @var{X}, @var{Y}, @var{Z} and @var{a} must be the same.
+##
+## @var{lon} is a scalar or a vector containing the geodetic longitude, in
+## radians.
+## @var{lat} is a scalar or a vector containing the geodetic latitude, in
+## radians.
+## @var{h} is a scalar or a vector containing the ellipsoidal height, in the
+## same units of @var{a}.
+##
+## @seealso{op_geod2geoc}
+## @end deftypefn
+
+
+
+
+function [lon,lat,h] = op_geoc2geod(X,Y,Z,a,f)
+
+try
+    functionName = 'op_geoc2geod';
+    argumentNumber = 5;
+
+%*******************************************************************************
+%NUMBER OF INPUT ARGUMENTS CHECKING
+%*******************************************************************************
+
+    %number of input arguments checking
+    if nargin~=argumentNumber
+        error(['Incorrect number of input arguments (%d)\n\t         ',...
+               'Correct number of input arguments = %d'],...
+              nargin,argumentNumber);
+    end
+
+%*******************************************************************************
+%INPUT ARGUMENTS CHECKING
+%*******************************************************************************
+
+    %X must be a column vector
+    if ismatrix(X)
+        %X dimensions
+        [row,col] = size(X);
+        %checking dimensions
+        if (row>1)&(col>1)
+            %error message
+            error('The first input argument must be scalar or column vector');
+        elseif col>1
+            %convert into column vector, if it was a row vector
+            X = X(:);
+        end
+    else
+        error('The first input argument is not numeric');
+    end
+    %Y must be a column vector
+    if ismatrix(Y)
+        %Y dimensions
+        [row,col] = size(Y);
+        %checking dimensions
+        if (row>1)&(col>1)
+            %error message
+            error('The second input argument must be scalar or column vector');
+        elseif col>1
+            %convert into column vector, if it was a row vector
+            Y = Y(:);
+        end
+    else
+        error('The second input argument is not numeric');
+    end
+    %Z must be a column vector
+    if ismatrix(Z)
+        %Z dimensions
+        [row,col] = size(Z);
+        %checking dimensions
+        if (row>1)&(col>1)
+            %error message
+            error('The third input argument must be scalar or column vector');
+        elseif col>1
+            %convert into column vector, if it was a row vector
+            Z = Z(:);
+        end
+    else
+        error('The third input argument is not numeric');
+    end
+    %checking number of rows
+    rowX = size(X,1);
+    rowY = size(Y,1);
+    rowZ = size(Z,1);
+    maxRow = max([rowX rowY rowZ]);
+    if ((rowX~=1)&(rowX~=maxRow))|((rowY~=1)&(rowY~=maxRow))|...
+       ((rowZ~=1)&(rowZ~=maxRow))
+        error('Incorrect dimensions in X, Y or Z');
+    end
+    %a must be a column vector
+    if ~isscalar(a)
+        error('The fourth input argument must be scalar');
+    end
+    %f must be a column vector
+    if ~isscalar(f)
+        error('The fifth input argument must be scalar');
+    end
+catch
+    %error message
+    error('\n\tIn function %s:\n\t -%s ',functionName,lasterr);
+end
+
+%*******************************************************************************
+%COMPUTATION
+%*******************************************************************************
+
+try
+    %first squared eccentricity of the ellipsoid
+    e2 = 2.0*f-f^2;
+    %calling oct function
+    [lon,lat,h] = _op_geoc2geod(X,Y,Z,a,e2);
+catch
+    %error message
+    error('\n\tIn function %s:\n\tIn function %s ',functionName,lasterr);
+end
+
+
+
+
+%*****END OF FUNCION*****
+
+
+
+
+%*****FUNCTION TESTS*****
+
+
+
+
+%!test
+%!  [lon,lat,h]=op_geoc2geod(2587045.819,1879598.809,5501461.606,6378388,1/297);
+%!  assert(lon,0.628318530616265,1e-11)
+%!  assert(lat,1.04719755124682,1e-11)
+%!  assert(h,999.999401183799,1e-5)
+%!error(op_geoc2geod)
+%!error(op_geoc2geod(1,2,3,4,5,6))
+%!error(op_geoc2geod('string',2,3,4,5))
+%!error(op_geoc2geod(1,'string',3,4,5))
+%!error(op_geoc2geod(1,2,'string',4,5))
+%!error(op_geoc2geod(1,2,3,[4 4],5))
+%!error(op_geoc2geod(1,2,3,4,[5 5]))
+%!error(op_geoc2geod([1 1;2 2],2,3,4,5))
+%!error(op_geoc2geod(1,[2 2;3 3],3,4,5))
+%!error(op_geoc2geod(1,2,[3 3;4 4],4,5))
+%!error(op_geoc2geod([1;1],[2 2 2],3,4,5))
+
+
+
+
+%*****END OF TESTS*****