Diff of /src/_op_geod2geoc.cc [000000] .. [d89be7] Maximize Restore

  Switch to side-by-side view

--- a
+++ b/src/_op_geod2geoc.cc
@@ -0,0 +1,141 @@
+/* -*- coding: utf-8 -*- */
+/* 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 this software; see the file COPYING.  If not, see
+ * <http://www.gnu.org/licenses/>.
+ */
+/******************************************************************************/
+/******************************************************************************/
+#define HELPTEXT "\
+-*- texinfo -*-\n\
+@deftypefn{Loadable Function}{[@var{X},@var{Y},@var{Z}] =}\
+_op_geod2geoc(@var{lon},@var{lat},@var{h},@var{a},@var{e2})\n\
+\n\
+@cindex Geodetic to geocentric coordinates.\n\
+\n\
+This function converts geodetic coordinates into cartesian tridimensional\n\
+geocentric coordinates using the PROJ.4 function pj_geodetic_to_geocentric().\n\
+\n\
+@var{lon} is a column vector containing the geodetic longitude, in radians.\n\
+@var{lat} is a column vector containing the geodetic latitude, in radians.\n\
+@var{h} is a column vector containing the ellipsoidal height.\n\
+@var{a} is a scalar containing the semi-major axis of the ellipsoid.\n\
+@var{e2} is a scalar containing the squared first eccentricity of the \
+ellipsoid.\n\
+\n\
+The coordinate vectors @var{lon}, @var{lat} and @var{h} must be all scalars\n\
+or all column vectors (of the same size).\n\
+The units of @var{h} and @var{a} must be the same.\n\
+\n\
+@var{X} is a column vector containing the X geocentric coordinate, in the same \
+units of @var{a}.\n\
+@var{Y} is a column vector containing the Y geocentric coordinate, in the same \
+units of @var{a}.\n\
+@var{Z} is a column vector containing the Z geocentric coordinate, in the same \
+units of @var{a}.\n\
+@seealso{_op_geoc2geod}\n\
+@end deftypefn"
+/******************************************************************************/
+/******************************************************************************/
+#include<octave/oct.h>
+#include<proj_api.h>
+#include<cstring>
+#include<cstdlib>
+/******************************************************************************/
+/******************************************************************************/
+#define ERRORTEXT 1000
+/******************************************************************************/
+/******************************************************************************/
+DEFUN_DLD(_op_geod2geoc,args,,HELPTEXT)
+{
+    //error message
+    char errorText[ERRORTEXT]="_op_geod2geoc:\n\t";
+    //output list
+    octave_value_list outputList;
+    ////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////
+    //checking input arguments
+    if(args.length()!=5)
+    {
+        //error text
+        sprintf(&errorText[strlen(errorText)],
+                "Incorrect number of input arguments\n\t"
+                "See help _op_geod2geoc");
+        //error message
+        error(errorText);
+    }
+    else
+    {
+        //error code
+        int* idError=NULL;
+        //loop index
+        size_t i=0;
+        //geodetic coordinates
+        ColumnVector lon=args(0).column_vector_value();
+        ColumnVector lat=args(1).column_vector_value();
+        ColumnVector h=args(2).column_vector_value();
+        //ellipsoidal parameters
+        double a=args(3).double_value();
+        double e2=args(4).double_value();
+        //number of elements
+        size_t nElem=static_cast<size_t>(lon.rows());
+        //geocentric coordinates
+        ColumnVector xOut(nElem);
+        ColumnVector yOut(nElem);
+        ColumnVector zOut(nElem);
+        //computation vectors
+        double* x=NULL;
+        double* y=NULL;
+        double* z=NULL;
+        ////////////////////////////////////////////////////////////////////////
+        ////////////////////////////////////////////////////////////////////////
+        //copy input data
+        for(i=0;i<nElem;i++)
+        {
+            xOut(i) = lon(i);
+            yOut(i) = lat(i);
+            zOut(i) = h(i);
+        }
+        //pointers to output data
+        x = const_cast<double*>(xOut.data());
+        y = const_cast<double*>(yOut.data());
+        z = const_cast<double*>(zOut.data());
+        ////////////////////////////////////////////////////////////////////////
+        ////////////////////////////////////////////////////////////////////////
+        //transformation
+        if(pj_geodetic_to_geocentric(a,e2,nElem,1,x,y,z))
+        {
+            //get error code
+            idError = pj_get_errno_ref();
+            //error text
+            sprintf(&errorText[strlen(errorText)],"%s",pj_strerrno(*idError));
+            //error message
+            error(errorText);
+            //exit
+            return outputList;
+        }
+        ////////////////////////////////////////////////////////////////////////
+        ////////////////////////////////////////////////////////////////////////
+        //output parameters list
+        outputList(0) = xOut;
+        outputList(1) = yOut;
+        outputList(2) = zOut;
+    }
+    ////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////
+    //exit
+    return outputList;
+}