Skip to content

Commit

Permalink
1.1: library
Browse files Browse the repository at this point in the history
  • Loading branch information
zvezdochiot committed Jul 10, 2022
1 parent 063187b commit 5a62c66
Show file tree
Hide file tree
Showing 5 changed files with 158 additions and 173 deletions.
14 changes: 11 additions & 3 deletions Makefile
Original file line number Diff line number Diff line change
@@ -1,11 +1,13 @@
PROJNAME = egm96
SRCS = src
PLIBNAME = lib$(PROJNAME).a
PROGNAME1 = $(PROJNAME)
PROGNAME2 = $(PROJNAME)_gen
PROGS = $(PROGNAME1) $(PROGNAME2)
CC = gcc
CFLAGS = -Wall
LDFLAGS = -lm -s
AR = ar
PREFIX = /usr/local
INSTALL = install
LN = ln -fs
Expand All @@ -17,12 +19,18 @@ RM = rm -f
all: $(PROGS)

clean:
$(RM) $(PROGS) man_*.pdf
$(RM) $(SRCS)/*.o $(PLIBNAME) $(PROGS) man_*.pdf

$(PROGNAME1): $(SRCS)/$(PROGNAME1).c
$(SRCS)/$(PROJNAME).o: $(SRCS)/$(PROJNAME).c
$(CC) $(CFLAGS) -c $^ -o $@

$(PLIBNAME): $(SRCS)/$(PROJNAME).o
$(AR) rcs $@ $^

$(PROGNAME1): $(SRCS)/$(PROGNAME1)_cli.c $(PLIBNAME)
$(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS)

$(PROGNAME2): $(SRCS)/$(PROGNAME2).c
$(PROGNAME2): $(SRCS)/$(PROGNAME2).c $(PLIBNAME)
$(CC) $(CFLAGS) $^ -o $@ $(LDFLAGS)

manual: man_$(PROGNAME1).pdf man_$(PROGNAME2).pdf
Expand Down
97 changes: 34 additions & 63 deletions src/egm96.c
Original file line number Diff line number Diff line change
Expand Up @@ -45,9 +45,12 @@
#include <math.h>
#include "egm96.h"

//! EGM96 correction and harmonic coefficients
static double egm96_data[_coeffs+1][4]={0};

/* ************************************************************************** */

double hundu(double p[_coeffs+1],
double egm96_hundu(double p[_coeffs+1],
double sinml[_361+1], double cosml[_361+1],
double gr, double re)
{
Expand Down Expand Up @@ -88,7 +91,7 @@ double hundu(double p[_coeffs+1],
return ((a * GM) / (gr * re)) + (ac / 100.0) - 0.53;
}

void dscml(double rlon, double sinml[_361+1], double cosml[_361+1])
void egm96_dscml(double rlon, double sinml[_361+1], double cosml[_361+1])
{
double a = sin(rlon);
double b = cos(rlon);
Expand Down Expand Up @@ -117,7 +120,7 @@ void dscml(double rlon, double sinml[_361+1], double cosml[_361+1])
* Original programmer: Oscar L. Colombo, Dept. of Geodetic Science the Ohio State University, August 1980.
* ineiev: I removed the derivatives, for they are never computed here.
*/
void legfdn(unsigned m, double theta, double rleg[_361+1])
void egm96_legfdn(unsigned m, double theta, double rleg[_361+1])
{
static double drts[1301], dirt[1301], cothet, sithet, rlnn[_361+1];
static int ir; // TODO 'ir' must be set to zero before the first call to this sub.
Expand Down Expand Up @@ -192,7 +195,7 @@ void legfdn(unsigned m, double theta, double rleg[_361+1])
* latitude, and an approximate value of normal gravity at the point based the
* constants of the WGS84(g873) system are used.
*/
void radgra(double lat, double lon, double *rlat, double *gr, double *re)
void egm96_radgra(double lat, double lon, double *rlat, double *gr, double *re)
{
const double a = 6378137.0;
const double e2 = 0.00669437999013;
Expand All @@ -216,36 +219,36 @@ void radgra(double lat, double lon, double *rlat, double *gr, double *re)
* \param lon: Longitude in radians.
* \return The geoid undulation / altitude offset (in meters).
*/
double undulation(double lat, double lon)
double egm96_undulation(double lat, double lon)
{
double p[_coeffs+1], sinml[_361+1], cosml[_361+1], rleg[_361+1];
double rlat, gr, re;
unsigned int j, m, i, nmax1 = _nmax + 1;

// compute the geocentric latitude, geocentric radius, normal gravity
radgra(lat, lon, &rlat, &gr, &re);
egm96_radgra(lat, lon, &rlat, &gr, &re);
rlat = (M_PI / 2) - rlat;

for (j = 1; j <= nmax1; j++)
{
m = j - 1;
legfdn(m, rlat, rleg);
egm96_legfdn(m, rlat, rleg);
for (i = j ; i <= nmax1; i++)
{
p[(((i - 1) * i) / 2) + m + 1] = rleg[i];
}
}
dscml(lon, sinml, cosml);
egm96_dscml(lon, sinml, cosml);

return hundu(p, sinml, cosml, gr, re);
return egm96_hundu(p, sinml, cosml, gr, re);
}

/* ************************************************************************** */

double egm96_compute_altitude_offset(double lat, double lon)
{
const double rad = (180.0 / M_PI);
return undulation(lat/rad, lon/rad);
return egm96_undulation(lat/rad, lon/rad);
}

/* ************************************************************************** */
Expand All @@ -257,7 +260,7 @@ double egm96_compute_altitude_offset(double lat, double lon)
* system of constants and are identical to those values used in the NIMA gridding procedure.
* Computed using subroutine grs written by N.K. PAVLIS.
*/
void dhcsin(FILE *f_12)
void egm96_dhcsin(FILE *f_12)
{
double c, s, ec, es;

Expand Down Expand Up @@ -289,7 +292,7 @@ void dhcsin(FILE *f_12)

/* ************************************************************************** */

void init_arrays(char* egmname, char* corname)
void egm96_init_arrays(char* egmname, char* corname)
{
FILE *f_1, *f_12;
int ig, n, m;
Expand All @@ -314,71 +317,39 @@ void init_arrays(char* egmname, char* corname)
// the correction coefficients are now read in
// the potential coefficients are now read in,
// and the reference even degree zonal harmonic coefficients removed to degree 6
dhcsin(f_12);
egm96_dhcsin(f_12);

fclose(f_1);
fclose(f_12);
}
}

/* ************************************************************************** */

/*!
* \brief Main function.
* \return 0 if success.
*
* The input files consist of:
* - correction coefficient set ("CORRCOEF") => unit = 1
* - potential coefficient set ("EGM96") => unit = 12
* - points at which to compute ("INPUT.dat") => unit = 14
* The output file is:
* - computed geoid heights ("OUTPUT.dat") => unit = 20
* - precomputed egm96_data.h (to use with the library)
* \brief Write precomputed EGM96 correction and harmonic coefficients to egm96_data.h
*/
int main(int argc, char *argv[])
void egm96_write_arrays(char* genname)
{
FILE *f_14, *f_20;
double flat, flon, u;
if (argc < 3)
{
printf("EGM96\n");
printf("Usage:\n");
printf(" %s EGM96-file COR-file [INPUT-file] [OUTPUT-file]\n", argv[0]);
}
else
FILE *precomp_out;
unsigned int i;

precomp_out = fopen(genname, "wb");
if (precomp_out)
{
init_arrays(argv[1], argv[2]);

if (argc > 3)
f_14 = fopen(argv[3], "rb");
else
f_14 = stdin;
if (argc > 4)
f_20 = fopen(argv[4], "wb");
else
f_20 = stdout;
if (f_14 && f_20)
fprintf(precomp_out, "#ifndef EGM96_DATA_H\n");
fprintf(precomp_out, "#define EGM96_DATA_H\n\n");
fprintf(precomp_out, "//! Precomputed EGM96 correction and harmonic coefficients\n");
fprintf(precomp_out, "static const double egm96_data[65342][4] = {\n");

for (i = 0; i <= _coeffs; i++)
{
// read geodetic latitude,longitude at point undulation is wanted
while (2 == fscanf(f_14, "%lg %lg", &flat, &flon))
{
// compute the geocentric latitude, geocentric radius, normal gravity
u = egm96_compute_altitude_offset(flat, flon);
fprintf(precomp_out, "{%g,%g,%g,%g},\n", egm96_data[i][0], egm96_data[i][1], egm96_data[i][2], egm96_data[i][3]);
}

// u is the geoid undulation from the EGM96 potential coefficient model
// including the height anomaly to geoid undulation correction term
// and a correction term to have the undulations refer to the
// WGS84 ellipsoid. the geoid undulation unit is meters.
fprintf(precomp_out, "};\n\n");
fprintf(precomp_out, "#endif // EGM96_DATA_H\n");

fprintf(f_20, "%14.7f %14.7f %16.7f\n", flat, flon, u);
}

fclose(f_14);
fclose(f_20);
}
fclose(precomp_out);
}

return 0;
}

/* ************************************************************************** */
5 changes: 2 additions & 3 deletions src/egm96.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,8 @@
#define _nmax (360) //!< Maximum degree and orders of harmonic coefficients.
#define _361 (361)

//! EGM96 correction and harmonic coefficients
static double egm96_data[_coeffs+1][4]={0};

void egm96_init_arrays(char* egmname, char* corname);
void egm96_write_arrays(char* genname);
/*!
* \brief Compute the geoid undulation from the EGM96 potential coefficient model, for a given latitude and longitude.
* \param latitude: Latitude in degrees.
Expand Down
108 changes: 108 additions & 0 deletions src/egm96_cli.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,108 @@
/*
* Copyright (c) 2006 D.Ineiev <[email protected]>
* Copyright (c) 2020 Emeric Grange <[email protected]>
*
* This software is provided 'as-is', without any express or implied warranty.
* In no event will the authors be held liable for any damages arising from
* the use of this software.
*
* Permission is granted to anyone to use this software for any purpose,
* including commercial applications, and to alter it and redistribute it
* freely, subject to the following restrictions:
*
* 1. The origin of this software must not be misrepresented; you must not
* claim that you wrote the original software. If you use this software
* in a product, an acknowledgment in the product documentation would be
* appreciated but is not required.
* 2. Altered source versions must be plainly marked as such, and must not be
* misrepresented as being the original software.
* 3. This notice may not be removed or altered from any source distribution.
*/

/*
* This program is designed for the calculation of a geoid undulation at a point
* whose latitude and longitude is specified.
*
* This program is designed to be used with the constants of EGM96 and those of
* the WGS84(g873) system. The undulation will refer to the WGS84 ellipsoid.
*
* It's designed to use the potential coefficient model EGM96 and a set of
* spherical harmonic coefficients of a correction term.
* The correction term is composed of several different components, the primary
* one being the conversion of a height anomaly to a geoid undulation.
* The principles of this procedure were initially described in the paper:
* - use of potential coefficient models for geoid undulation determination using
* a spherical harmonic representation of the height anomaly/geoid undulation
* difference by R.H. Rapp, Journal of Geodesy, 1996.
*
* This program is a modification of the program described in the following report:
* - a fortran program for the computation of gravimetric quantities from high
* degree spherical harmonic expansions, Richard H. Rapp, report 334, Department
* of Geodetic Science and Surveying, the Ohio State University, Columbus, 1982
*/

#include <stdio.h>
#include <math.h>
#include "egm96.h"

/* ************************************************************************** */

/*!
* \brief Main function.
* \return 0 if success.
*
* The input files consist of:
* - correction coefficient set ("CORRCOEF") => unit = 1
* - potential coefficient set ("EGM96") => unit = 12
* - points at which to compute ("INPUT.dat") => unit = 14
* The output file is:
* - computed geoid heights ("OUTPUT.dat") => unit = 20
* - precomputed egm96_data.h (to use with the library)
*/
int main(int argc, char *argv[])
{
FILE *f_14, *f_20;
double flat, flon, u;
if (argc < 3)
{
printf("EGM96\n");
printf("Usage:\n");
printf(" %s EGM96-file COR-file [INPUT-file] [OUTPUT-file]\n", argv[0]);
}
else
{
egm96_init_arrays(argv[1], argv[2]);

if (argc > 3)
f_14 = fopen(argv[3], "rb");
else
f_14 = stdin;
if (argc > 4)
f_20 = fopen(argv[4], "wb");
else
f_20 = stdout;
if (f_14 && f_20)
{
// read geodetic latitude,longitude at point undulation is wanted
while (2 == fscanf(f_14, "%lg %lg", &flat, &flon))
{
// compute the geocentric latitude, geocentric radius, normal gravity
u = egm96_compute_altitude_offset(flat, flon);

// u is the geoid undulation from the EGM96 potential coefficient model
// including the height anomaly to geoid undulation correction term
// and a correction term to have the undulations refer to the
// WGS84 ellipsoid. the geoid undulation unit is meters.

fprintf(f_20, "%14.7f %14.7f %16.7f\n", flat, flon, u);
}

fclose(f_14);
fclose(f_20);
}
}

return 0;
}

/* ************************************************************************** */
Loading

0 comments on commit 5a62c66

Please sign in to comment.