From 5a62c6676f5ad13c9f1751adaef482e4d21b1231 Mon Sep 17 00:00:00 2001 From: zvezdochiot Date: Sun, 10 Jul 2022 03:52:16 +0300 Subject: [PATCH] 1.1: library --- Makefile | 14 +++++-- src/egm96.c | 97 +++++++++++++++---------------------------- src/egm96.h | 5 +-- src/egm96_cli.c | 108 ++++++++++++++++++++++++++++++++++++++++++++++++ src/egm96_gen.c | 107 ++--------------------------------------------- 5 files changed, 158 insertions(+), 173 deletions(-) create mode 100644 src/egm96_cli.c diff --git a/Makefile b/Makefile index 09f740b..a7e4bff 100644 --- a/Makefile +++ b/Makefile @@ -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 @@ -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 diff --git a/src/egm96.c b/src/egm96.c index fb52a22..c2836d4 100644 --- a/src/egm96.c +++ b/src/egm96.c @@ -45,9 +45,12 @@ #include #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) { @@ -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); @@ -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. @@ -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; @@ -216,28 +219,28 @@ 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); } /* ************************************************************************** */ @@ -245,7 +248,7 @@ double undulation(double lat, double lon) 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); } /* ************************************************************************** */ @@ -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; @@ -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; @@ -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; } /* ************************************************************************** */ diff --git a/src/egm96.h b/src/egm96.h index ccdca59..75eee09 100644 --- a/src/egm96.h +++ b/src/egm96.h @@ -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. diff --git a/src/egm96_cli.c b/src/egm96_cli.c new file mode 100644 index 0000000..e8252cb --- /dev/null +++ b/src/egm96_cli.c @@ -0,0 +1,108 @@ +/* + * Copyright (c) 2006 D.Ineiev + * Copyright (c) 2020 Emeric Grange + * + * 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 +#include +#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; +} + +/* ************************************************************************** */ diff --git a/src/egm96_gen.c b/src/egm96_gen.c index a2fc2bb..8892eea 100644 --- a/src/egm96_gen.c +++ b/src/egm96_gen.c @@ -25,107 +25,6 @@ /* ************************************************************************** */ -/*! - * \param f_12: EGM96 coefficients file. - * - * The even degree zonal coefficients given below were computed for the WGS84(g873) - * 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) -{ - double c, s, ec, es; - - const double j2 = 0.108262982131e-2; - const double j4 = -.237091120053e-05; - const double j6 = 0.608346498882e-8; - const double j8 = -0.142681087920e-10; - const double j10 = 0.121439275882e-13; - - unsigned n, m = (((_nmax+1) * (_nmax+2)) / 2); - for (n = 1; n <= m; n++) - { - egm96_data[n][2] = egm96_data[n][3] = 0; - } - - while (6 == fscanf(f_12, "%i %i %lf %lf %lf %lf", &n, &m, &c, &s, &ec, &es)) - { - if (n > _nmax) continue; - n = ((n * (n + 1)) / 2) + m + 1; - egm96_data[n][2] = c; - egm96_data[n][3] = s; - } - egm96_data[4][2] += j2 / sqrt(5); - egm96_data[11][2] += j4 / 3.0; - egm96_data[22][2] += j6 / sqrt(13); - egm96_data[37][2] += j8 / sqrt(17); - egm96_data[56][2] += j10 / sqrt(21); -} - -/* ************************************************************************** */ - -void init_arrays(char* egmname, char* corname) -{ - FILE *f_1, *f_12; - int ig, n, m; - double t1, t2; - unsigned int i; - - f_1 = fopen(corname, "rb"); // correction coefficient file: modified with 'sed -e"s/D/e/g"' to be read with fscanf - f_12 = fopen(egmname, "rb"); // potential coefficient file - - if (f_1 && f_12) - { - for (i = 1; i <= _coeffs; i++) - egm96_data[i][0] = egm96_data[i][1] = 0; - - while (4 == fscanf(f_1, "%i %i %lg %lg", &n, &m, &t1, &t2)) - { - ig = ((n * (n+1)) / 2) + m + 1; - egm96_data[ig][0] = t1; - egm96_data[ig][1] = t2; - } - - // 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); - - fclose(f_1); - fclose(f_12); - } -} - -/*! - * \brief Write precomputed EGM96 correction and harmonic coefficients to egm96_data.h - */ -void write_arrays(char* genname) -{ - FILE *precomp_out; - unsigned int i; - - precomp_out = fopen(genname, "wb"); - if (precomp_out) - { - 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++) - { - fprintf(precomp_out, "{%g,%g,%g,%g},\n", egm96_data[i][0], egm96_data[i][1], egm96_data[i][2], egm96_data[i][3]); - } - - fprintf(precomp_out, "};\n\n"); - fprintf(precomp_out, "#endif // EGM96_DATA_H\n"); - - fclose(precomp_out); - } -} - -/* ************************************************************************** */ - /*! * \brief Main function. * \return 0 if success. @@ -148,12 +47,12 @@ int main(int argc, char *argv[]) } else { - init_arrays(argv[1], argv[2]); + egm96_init_arrays(argv[1], argv[2]); if (argc > 3) - write_arrays(argv[3]); + egm96_write_arrays(argv[3]); else - write_arrays("EGM96_data.h"); + egm96_write_arrays("EGM96_data.h"); } return 0;