用C语言写的坡度算法

#include 
#include 
#include 
#include “hdf.h“
#include “mfhdf.h“
#include “HdfEosDef.h“

#define MAXSTRING			2048
#define MAXNAMEL			256			/* maximum length of a attribute name */
#define MAXATTR           16      /* maximum number of attributes */
#define MAXATTRL          65535    /* maximum length of attribute 64MB */
#define MAXSTRLEN    1024  /* maximum string length */
#define MAXDIMS      16   /* maximum number of dimensions */
#define MAXDIMSIZE   12150  /* maximum dimension size */
#define MAXATTRLEN   1048576  /* maximum attribute length 8MB */
#define NPROJPARAM   16     /* number of projection parameters */
#define XRES         -1.0  /* default grid column resolution */
#define YRES         -1.0  /* default grid row resolution */
#define FILLVALUE    0

#define NFIELD       2
#define NDATATYPE   43
#define NSPHEROID    20       /* number of spheroid types */

#define MAXGRID    128        /* number of maximum grids in a file */

#define RADIUS		6370997	/*radius 6370997m*/
/*#define PI			3.14159265358*/
#define NANVALUE	-32768.0

#define ELLIPSOID_A 6378137.0
#define ELLIPSOID_E2 0.006694385
#define Radians（x） （（x） * PI/180.0）

/*****************************/
/*    ____________________________
	  |c1      |c2      |c3      |
	  |        |        |        |
	  |        |  north |        |        
	  |        |        |        |
	  |________|________|________|          
	  |c4      |c5      |c6      |
	  |        |        |        |
	  |  east  | ns_med |  west  |
	  |        |        |        |
	  |________|________|________|
	  |c7      |c8      |c9      |
	  |        |        |        |
	  |        |  south |        |
	  |        |        |        |
	  |________|________|________|
*/

/* WARNING: this code is preliminary and may be changed
 * including calling sequences to any of the functions
 * defined here
 */

/* distance from point to point along a geodesic 
 * code from
 *   Paul D. Thomas
 *   “Spheroidal Geodesics Reference Systems and Local Geometry“
 *   U.S. Naval Oceanographic Office p. 162
 *   Engineering Library 526.3 T36s
 */

static float64 boa;
static float64 f;
static float64 ff64;
static float64 al;

/* must be called once to establish the ellipsoid */

/*!
 * \brief begin geodesic distance
 *
 * Initializes the distance calculations for the ellipsoid with
 * semi-major axis a （in meters） and ellipsoid eccentricity squared
 * e2. It is used only for the latitude-longitude projection.
 *
 *  \param a
 *  \param e2
 *  \return int
 */

int32 G_begin_geodesic_distance（float64 afloat64 e2）
{
    al = a;
    boa = sqrt （1 - e2）;
    f = 1 - boa;
    ff64 = f*f/64;

    return 0;
}

static float64 t1t2t3t4t1rt2r;

/* must be called first */

/*!
 * \brief set geodesic distance lat1
 *
 * Set the first latitude.
 *
 *  \param lat1
 *  \return int
 */

int32 G_set