From 95e17fce025e8fe688923563b331e641288c113d Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Wed, 5 May 2021 15:31:06 +0800 Subject: [PATCH] initial upload --- .gitignore | 4 + CMakeLists.txt | 25 + LICENSE.txt | 33 + README.md | 122 ++ lib/CMakeLists.txt | 49 + lib/constants.h | 60 + lib/geometry.cpp | 48 + lib/geometry.h | 37 + lib/glq.cpp | 265 +++++ lib/glq.h | 182 +++ lib/grav_tess.cpp | 628 ++++++++++ lib/grav_tess.h | 93 ++ lib/linalg.cpp | 142 +++ lib/linalg.h | 11 + lib/logger.cpp | 112 ++ lib/logger.h | 166 +++ lib/parsers.cpp | 678 +++++++++++ lib/parsers.h | 75 ++ lib/tessb_main.cpp | 385 ++++++ lib/tessb_main.h | 14 + lib/version.cpp | 17 + lib/version.h | 18 + toolkits/CMakeLists.txt | 47 + toolkits/tessbx.cpp | 11 + toolkits/tessby.cpp | 11 + toolkits/tessbz.cpp | 11 + toolkits/tessutil_combine_grids.cpp | 89 ++ toolkits/tessutil_gradient_calculator.cpp | 343 ++++++ toolkits/tessutil_magnetize_model.c | 1313 +++++++++++++++++++++ 29 files changed, 4989 insertions(+) create mode 100644 CMakeLists.txt create mode 100755 LICENSE.txt create mode 100755 README.md create mode 100644 lib/CMakeLists.txt create mode 100755 lib/constants.h create mode 100755 lib/geometry.cpp create mode 100755 lib/geometry.h create mode 100755 lib/glq.cpp create mode 100755 lib/glq.h create mode 100755 lib/grav_tess.cpp create mode 100755 lib/grav_tess.h create mode 100755 lib/linalg.cpp create mode 100755 lib/linalg.h create mode 100755 lib/logger.cpp create mode 100755 lib/logger.h create mode 100755 lib/parsers.cpp create mode 100755 lib/parsers.h create mode 100755 lib/tessb_main.cpp create mode 100755 lib/tessb_main.h create mode 100755 lib/version.cpp create mode 100755 lib/version.h create mode 100644 toolkits/CMakeLists.txt create mode 100755 toolkits/tessbx.cpp create mode 100755 toolkits/tessby.cpp create mode 100755 toolkits/tessbz.cpp create mode 100755 toolkits/tessutil_combine_grids.cpp create mode 100755 toolkits/tessutil_gradient_calculator.cpp create mode 100755 toolkits/tessutil_magnetize_model.c diff --git a/.gitignore b/.gitignore index 259148f..8242f61 100644 --- a/.gitignore +++ b/.gitignore @@ -30,3 +30,7 @@ *.exe *.out *.app + +build/ +.vs_code/ +.DS_Store \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..ec00c85 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,25 @@ +cmake_minimum_required(VERSION 3.15.2) +# 设置项目名称与语言 +project(LIBMAGTESS VERSION 1.0) +# 设置安装地址 +if(${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Linux") + message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME}) + set(CMAKE_INSTALL_PREFIX /usr/local) +elseif (${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Darwin") + message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME}) + set(CMAKE_INSTALL_PREFIX /usr/local) +elseif (${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Windows") + message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME}) + # 使用MinGW GCC编译时需取消注释 + #set(CMAKE_C_COMPILER gcc) + #set(CMAKE_CXX_COMPILER g++) + set(CMAKE_INSTALL_PREFIX D:/Library) +else() + message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME}) + set(CMAKE_INSTALL_PREFIX /usr/local) +endif() +message(STATUS "Install prefix: " ${CMAKE_INSTALL_PREFIX}) + +# 添加库源文件地址 +add_subdirectory(lib) +add_subdirectory(toolkits) \ No newline at end of file diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100755 index 0000000..fa9e309 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,33 @@ +BSD 2-Clause License + +Copyright (c) 2012-2015, Leonardo Uieda +Copyright (c) 2017-2020, Eldar Baykiev +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +3. Neither the name of Leonardo Uieda nor the names of any contributors* + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +* Contributors to the original tesseroids source and binary forms, which are +under BSD 3-Clause License diff --git a/README.md b/README.md new file mode 100755 index 0000000..60bd43a --- /dev/null +++ b/README.md @@ -0,0 +1,122 @@ +# Forward modeling of magnetic field in spherical coordinates + +Magnetic tesseroids is a collection of command-line tools for modelling of the magnetic field with spherical prisms (tesseroids) used as magnetic sources. + +To cite _magnetic tesseroids_ in publications, please use our paper published in Computers & Geosciences: + +>**Eldar Baykiev**, **Jörg Ebbing**, **Marco Brönner**, **Karl Fabian**, Forward modeling magnetic fields of induced and remanent magnetization in the lithosphere using tesseroids, _Computers & Geosciences_, Volume 96, November 2016, Pages 124-135, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2016.08.004. + +Article can also be found here http://goo.gl/x9g7gi (researchgate). + +## Usage of _magnetic tesseroids_ +_Magnetic tesseroids_ are based on the existing program of Leonardo Uieda called tesseroids (Uieda, 2013) of version 1.1 (https://github.com/leouieda/tesseroids/releases/tag/v1.1). It inherits the interface of tesseroids-1.1 but with several changes. Present appendix describes constants and units used be the program, as well as input and output format. +### Constants +1. Geocentric mean Earth's radius _R_E_ = 6378.137 km. +1. Magnetic permeability of a free space _µ_0_ = 4π × 10^-7 H·m^-1 + +### List of programs +The tessbx, tessby, tessbz are programs that calculate the corresponding components (x - north, y - east, **z - up**) of the magnetic field of the tesseroid model on the computational grid. + +### Input: tesseroid model +The input model file should be a text file where each line describe one tesseroid in such space separated format: +> `W E S N HEIGHT_OF_TOP HEIGHT_OF_BOTTOM DENSITY SUSCEPTIBILITY BX BY BZ` + +`W`, `E`, `S`, `N` correspond to the western, eastern, southern and northern edges of a tesseroid (_λ_1_, _λ_2_, _ϕ_1_, _ϕ_2_ respectively) expressed in decimal degrees [°]. +`HEIGHT_OF_TOP` and `HEIGHT_OF_BOTTOM` define the top and the bottom of tesseroid (_r_2_ and _r_1_ respectively). Both are counted from geocentric mean Earth's radius in meters [m]. If a tesseroid should be placed beneath the mean surface, than the values of these parameters should be negative. Note that `HEIGHT_OF_TOP` > `HEIGHT_OF_BOTTOM`. +`DENSITY` is the density _ρ_ of tesseroid in kilogram per cubic meter [kg/m^3] +`SUSCEPTIBILITY` is the susceptibility _χ_ of tesseroid in SI units. +`BX`, `BY` and `BZ` are the components of the magnetizing field in the local North-East-Up Cartesian coordinate system of a tesseroids' geometric center. They can be taken from any core field's model. Values are given in nanotesla [nT]. +In case of remanent magnetic field modeling, susceptibility must be set 1 SI and `BX`, `BY` and `BZ` values than would define the direction of remanent magnetization vector. +This example shows a model made of 3 neighboring tesseroids near the North Pole: +> `-74 -73 89 90 -1000.000000 -11650.000000 1.000000 1.000000 334.9504973176 -1969.9308033594 -56572.6324041700` + +> `-73 -72 89 90 -1000.000000 -11650.000000 1.000000 1.000000 370.1879538142 -1968.1093976826 -56571.2826313492` + +> `-72 -71 89 90 -1000.000000 -11650.000000 1.000000 1.000000 405.4388222633 -1965.6409379187 -56569.9502088641` + +### Input: computation grid +Computation grid can be regular or irregular and should be also a text file where each line describe the position of one computation point in such space separated format: +>`LON LAT ALT` + +`LON` and `LAT` correspond to the longitude and latitude of the point in decimal degrees [°]. +`ALT` corresponds to the altitude of the point above the mean surface in meters [m]. +Note that the program tessgrd from original tesseroids-1.1 can be used to create a regular computation grid (see Uieda, 2013). +This example shows a grid made of 6 points with the same latitude and the altitude of 400 km: +> `-6 51 400000 ` + +> `-5.8 51 400000 ` + +> `-5.6 51 400000 ` + +> `-5.4 51 400000 ` + +> `-5.2 51 400000 ` + +> `-5 51 400000` + +### Performing calculations +Example: to calculate the vertical component of the magnetic field of a model in file modelfile.txt on a grid from file gridpoints.txt one can simply use a console command: +``` +tessbz modelfile.txt < gridpoints.txt > gz_output.txt +``` + +The result would be written in the file gz_output.txt. +### Output format +The programs' output is a modified grid file where in the end of each line the calculated value of a corresponding magnetic field component would be written. Values are given in nanotesla [nT] in the local North-East-Up coordinate system of a computational point. +### Additional features +Magnetic tesseroids support features like piping and integration accuracy adjustment from tesseroids-1.1. Please, check sections in the tesseroids-1.1 manual (Uieda, 2013) relative to the gravity calculation programs to get more information. + +## Utilities +### tessutil_magnetize_model +This program is made to 'magnetize' any existing tesseroid model by any given main field spherical harmonic model. +Usage: +``` +tessutil_magnetize_model [SH coeff file] [input tesseroid model file] [day] [month] [year] [output tesseroid model file] +``` + +### tessutil_gradient_calculator +Gradient calculator (Baykiev et al., in press). +Usage: +``` +tessutil_gradient_calculator -bx[Bx grid file] -by[By grid file] -bz[Bx grid file] -o[output component] -c2 >> output_file.dat +``` + +All grid files should be in tessgrd format. With option `-c1` program reads input grid bz as its direction is upward, with option `-c2` - downward, just as in magnetic tesseroids output. Output of gradient calculator is always in North-East-Down coordinate system. + +Known issue: rounding error when processing grids with spacing equal or less than 0.2 degrees. + +### tessutil_combine_grids +Sums calculated grids. +Usage: +``` +tessutil_combine_grids [grid file1] [factor1] ... [grid fileN] [factorN] >> output_file.dat +``` + +Each grid is multiplied by factor (susceptibility) and then the sum of all grids is calculated. + +## Installation (version 1.1) +1. Download source code from [GitHub](https://github.com/eldarbaykiev/magnetic-tesseroids): + +``` +git clone https://github.com/eldarbaykiev/magnetic-tesseroids.git +``` + +2. On **Linux**, install [OpenBLAS](https://www.openblas.net/) library: + +``` +sudo apt-get install libopenblas-base libopenblas-dev +``` + +On **macOS**, make sure that [Xcode](https://developer.apple.com/xcode/) is installed and [Accelerate framework](https://developer.apple.com/documentation/accelerate) is available. + +3. Run **make** + +``` +make +``` + +To compile all utilities, run + +``` +make tools +``` diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt new file mode 100644 index 0000000..73398ab --- /dev/null +++ b/lib/CMakeLists.txt @@ -0,0 +1,49 @@ +if(${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Windows") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2") +endif() + +if(${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Linux") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -lopenblas -lm -O2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -lopenblas -lm -O2") +endif() + +if(${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Darwin") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -framework Accelerate -O2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -framework Accelerate -O2") +endif() + +# 设置库文件的输出地址 +set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) + +# 设定库源文件 +aux_source_directory(. LIBMAGTESS_SRC) + +# 以下部分为库的编译 +# 注意目标名必须唯一 所以不能直接生成相同名称的动态库与静态库 +# 注意此处不必为目标名称添加lib前缀和相应后缀，cmake会自行添加 +add_library(magtess SHARED ${LIBMAGTESS_SRC}) +# 首先添加静态库的生成命令 +add_library(magtess_static STATIC ${LIBMAGTESS_SRC}) +# 设置静态库的输出名称从而获得与动态库名称相同的静态库 +set_target_properties(magtess_static PROPERTIES OUTPUT_NAME "magtess") +# 设置输出目标属性以同时输出动态库与静态库 +set_target_properties(magtess PROPERTIES CLEAN_DIRECT_OUTPUT 1) +set_target_properties(magtess_static PROPERTIES CLEAN_DIRECT_OUTPUT 1) +# 设置动态库的版本号 +set_target_properties(magtess PROPERTIES VERSION 1.0 SOVERSION 1.0) + +# 库的安装命令 +if(WIN32) + install(TARGETS magtess DESTINATION lib) + install(TARGETS magtess_static DESTINATION lib) +else() + install(TARGETS magtess magtess_static + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib) +endif() + +# 头文件安装命令 +file(GLOB LIBMEGTESS_HEAD *.h) + +install(FILES ${LIBMEGTESS_HEAD} DESTINATION include/magtess) \ No newline at end of file diff --git a/lib/constants.h b/lib/constants.h new file mode 100755 index 0000000..cf04504 --- /dev/null +++ b/lib/constants.h @@ -0,0 +1,60 @@ +/* +Define constants used, like the gravitational constant and unit conversions. + +Values are assigned in file constants.c + +All values are in SI units! +*/ + +#ifndef _TESSEROIDS_CONSTANTS_H_ +#define _TESSEROIDS_CONSTANTS_H_ + +/* Mean Earth radius [\f$ m \f$] */ +const double MEAN_EARTH_RADIUS = 6378137.0; +const double EARTH_RADIUS_IGRF_KM = 6371.2; + +/* The gravitational constant [\f$ m^3*kg^{-1}*s^{-1} \f$] */ +const double G = 0.00000000006673; + + + +/* Conversion factor from SI units to Eotvos +[\f$ \frac{1}{s^2} = 10^9\ Eotvos \f$] */ +const double SI2EOTVOS = 1000000000.0; +const double EOTVOS2SI = 0.000000001; + +/* Conversion factor from SI units to mGal +[\f$ 1 \frac{m}{s^2} = 10^5\ mGal \f$] */ +const double SI2MGAL = 100000.0; + +/* Pi */ +#ifdef __cplusplus +const double PI = 3.1415926535897932384626433832795; +#else +#define PI 3.1415926535897932384626433832795 +#endif + +/* minimum distance-to-size ratio for potential computations to be accurate */ +const double TESSEROID_POT_SIZE_RATIO = 1.5; +/* Minimum distance-to-size ratio for gravity computations to be accurate */ +const double TESSEROID_GX_SIZE_RATIO = 3; +const double TESSEROID_GY_SIZE_RATIO = 3; +const double TESSEROID_GZ_SIZE_RATIO = 2; +/* Minimum distance-to-size ratio for gravity gradient computations to be +accurate */ +const double TESSEROID_GXX_SIZE_RATIO = 3; +const double TESSEROID_GXY_SIZE_RATIO = 4.5; +const double TESSEROID_GXZ_SIZE_RATIO = 4; +const double TESSEROID_GYY_SIZE_RATIO = 3; +const double TESSEROID_GYZ_SIZE_RATIO = 4; +const double TESSEROID_GZZ_SIZE_RATIO = 3; + +const double M_0 = 4 * (PI) * 0.0000001; + +const double DEG2RAD = (PI)/180.0; + + +#define FALSE 0 +#define TRUE 1 + +#endif diff --git a/lib/geometry.cpp b/lib/geometry.cpp new file mode 100755 index 0000000..97c1a5f --- /dev/null +++ b/lib/geometry.cpp @@ -0,0 +1,48 @@ +/* +Data structures for geometric elements and functions that operate on them. +Defines the TESSEROID, SPHERE, and PRISM structures. +*/ + + +#include +#include +#include +#include +#include "constants.h" +#include "logger.h" +#include "geometry.h" + + +/* Split a tesseroid into 8. */ +void split_tess(TESSEROID tess, TESSEROID *split) +{ + double dlon = 0.5*(tess.e - tess.w), + dlat = 0.5*(tess.n - tess.s), + dr = 0.5*(tess.r2 - tess.r1), + ws[2], ss[2], r1s[2]; + int i, j, k, t = 0; + + ws[0] = tess.w; + ws[1] = tess.w + dlon; + ss[0] = tess.s; + ss[1] = tess.s + dlat; + r1s[0] = tess.r1; + r1s[1] = tess.r1 + dr; + for(k = 0; k < 2; k++) + { + for(j = 0; j < 2; j++) + { + for(i = 0; i < 2; i++) + { + split[t].w = ws[i]; + split[t].e = ws[i] + dlon; + split[t].s = ss[j]; + split[t].n = ss[j] + dlat; + split[t].r1 = r1s[k]; + split[t].r2 = r1s[k] + dr; + split[t].density = tess.density; + t++; + } + } + } +} diff --git a/lib/geometry.h b/lib/geometry.h new file mode 100755 index 0000000..516bb21 --- /dev/null +++ b/lib/geometry.h @@ -0,0 +1,37 @@ +/* +Data structures for geometric elements and functions that operate on them. +Defines the TESSEROID, SPHERE, and PRISM structures. +*/ + +#ifndef _TESSEROIDS_GEOMETRY_H_ +#define _TESSEROIDS_GEOMETRY_H_ + + +/* Store information on a tesseroid */ +typedef struct tess_struct { + /* s, n, w, e in degrees. r1 and r2 are the smaller and larger radius */ + double density; /* in SI units */ + double w; /* western longitude border in degrees */ + double e; /* eastern longitude border in degrees */ + double s; /* southern latitude border in degrees */ + double n; /* northern latitude border in degrees */ + double r1; /* smallest radius border in SI units */ + double r2; /* largest radius border in SI units */ + double suscept; /* magnetic susceptibility */ + double Bx; /* x-component of ambient magnetic field */ + double By; /* y-component of ambient magnetic field */ + double Bz; /* z-component of ambient magnetic field */ + + double cos_a1; + double sin_a1; + double cos_b1; + double sin_b1; + //double Rx; + //double Ry; + //double Rz; +} TESSEROID; + + +void split_tess(TESSEROID tess, TESSEROID *split); + +#endif diff --git a/lib/glq.cpp b/lib/glq.cpp new file mode 100755 index 0000000..d95f3f7 --- /dev/null +++ b/lib/glq.cpp @@ -0,0 +1,265 @@ +/* +Functions for implementing a Gauss-Legendre Quadrature numerical integration. +*/ + + +#include +#include +#include "constants.h" +#include "logger.h" +#include "glq.h" + + +/* Make a new GLQ structure and set all the parameters needed */ +GLQ * glq_new(int order, double lower, double upper) +{ + GLQ *glq; + int rc; + + glq = (GLQ *)malloc(sizeof(GLQ)); + if(glq == NULL) + { + return NULL; + } + glq->order = order; + glq->nodes = (double *)malloc(sizeof(double)*order); + if(glq->nodes == NULL) + { + free(glq); + return NULL; + } + glq->nodes_unscaled = (double *)malloc(sizeof(double)*order); + if(glq->nodes_unscaled == NULL) + { + free(glq); + free(glq->nodes); + return NULL; + } + glq->weights = (double *)malloc(sizeof(double)*order); + if(glq->weights == NULL) + { + free(glq); + free(glq->nodes); + free(glq->nodes_unscaled); + return NULL; + } + rc = glq_nodes(order, glq->nodes_unscaled); + if(rc != 0 && rc != 3) + { + switch(rc) + { + case 1: + log_error("glq_nodes invalid GLQ order %d. Should be >= 2.", + order); + break; + case 2: + log_error("glq_nodes NULL pointer for nodes"); + break; + default: + log_error("glq_nodes unknown error code %g", rc); + break; + } + glq_free(glq); + return NULL; + } + else if(rc == 3) + { + log_warning("glq_nodes max iterations reached in root finder"); + log_warning("nodes might not have desired accuracy %g", GLQ_MAXERROR); + } + rc = glq_weights(order, glq->nodes_unscaled, glq->weights); + if(rc != 0) + { + switch(rc) + { + case 1: + log_error("glq_weights invalid GLQ order %d. Should be >= 2.", + order); + break; + case 2: + log_error("glq_weights NULL pointer for nodes"); + break; + case 3: + log_error("glq_weights NULL pointer for weights"); + break; + default: + log_error("glq_weights unknown error code %d\n", rc); + break; + } + glq_free(glq); + return NULL; + } + if(glq_set_limits(lower, upper, glq) != 0) + { + glq_free(glq); + return NULL; + } + return glq; +} + + +/* Free the memory allocated to make a GLQ structure */ +void glq_free(GLQ *glq) +{ + free(glq->nodes); + + free(glq->nodes_unscaled); + free(glq->weights); + free(glq); +} + + +/* Calculates the GLQ nodes using glq_next_root. */ +int glq_nodes(int order, double *nodes) +{ + register int i; + int rc = 0; + double initial; + + if(order < 2) + { + return 1; + } + if(nodes == NULL) + { + return 2; + } + for(i = 0; i < order; i++) + { + initial = cos(PI*(order - i - 0.25)/(order + 0.5)); + if(glq_next_root(initial, i, order, nodes) == 3) + { + rc = 3; + } + } + return rc; +} + + +/* Put the GLQ nodes to the integration limits IN PLACE. */ +int glq_set_limits(double lower, double upper, GLQ *glq) +{ + /* Only calculate once to optimize the code */ + double tmpplus = 0.5*(upper + lower), tmpminus = 0.5*(upper - lower); + register int i; + + if(glq->order < 2) + { + return 1; + } + if(glq->nodes == NULL) + { + return 2; + } + if(glq->nodes_unscaled == NULL) + { + return 2; + } + for(i = 0; i < glq->order; i++) + { + glq->nodes[i] = tmpminus*glq->nodes_unscaled[i] + tmpplus; + + } + return 0; +} + +/* Calculate the next Legendre polynomial root given the previous root found. */ +int glq_next_root(double initial, int root_index, int order, double *roots) +{ + double x1, x0, pn, pn_2, pn_1, pn_line, sum; + int it = 0; + register int n; + + if(order < 2) + { + return 1; + } + if(root_index < 0 || root_index >= order) + { + return 2; + } + x1 = initial; + do + { + x0 = x1; + + /* Calculate Pn(x0) */ + /* Starting from P0(x) and P1(x), */ + /* find the others using the recursive relation: */ + /* Pn(x)=(2n-1)xPn_1(x)/n - (n-1)Pn_2(x)/n */ + pn_1 = 1.; /* This is Po(x) */ + pn = x0; /* and this P1(x) */ + for(n = 2; n <= order; n++) + { + pn_2 = pn_1; + pn_1 = pn; + pn = ( ((2*n - 1)*x0*pn_1) - ((n - 1)*pn_2) )/n; + } + /* Now calculate Pn'(x0) using another recursive relation: */ + /* Pn'(x)=n(xPn(x)-Pn_1(x))/(x*x-1) */ + pn_line = order*(x0*pn - pn_1)/(x0*x0 - 1); + /* Sum the roots found so far */ + for(n = 0, sum = 0; n < root_index; n++) + { + sum += 1./(x0 - roots[n]); + } + /* Update the estimate for the root */ + x1 = x0 - (double)pn/(pn_line - pn*sum); + + /** Compute the absolute value of x */ + #define GLQ_ABS(x) ((x) < 0 ? -1*(x) : (x)) + } while(GLQ_ABS(x1 - x0) > GLQ_MAXERROR && ++it <= GLQ_MAXIT); + #undef GLQ_ABS + + roots[root_index] = x1; + + /* Tell the user if stagnation occurred */ + if(it > GLQ_MAXIT) + { + return 3; + } + return 0; +} + + +/* Calculates the weighting coefficients for the GLQ integration. */ +int glq_weights(int order, double *nodes, double *weights) +{ + register int i, n; + double xi, pn, pn_2, pn_1, pn_line; + + if(order < 2) + { + return 1; + } + if(nodes == NULL) + { + return 2; + } + if(weights == NULL) + { + return 3; + } + for(i = 0; i < order; i++){ + + xi = nodes[i]; + + /* Find Pn'(xi) with the recursive relation to find Pn and Pn-1: */ + /* Pn(x)=(2n-1)xPn_1(x)/n - (n-1)Pn_2(x)/n */ + /* Then use: Pn'(x)=n(xPn(x)-Pn_1(x))/(x*x-1) */ + + /* Find Pn and Pn-1 stating from P0 and P1 */ + pn_1 = 1; /* This is Po(x) */ + pn = xi; /* and this P1(x) */ + for(n = 2; n <= order; n++) + { + pn_2 = pn_1; + pn_1 = pn; + pn = ((2*n - 1)*xi*pn_1 - (n - 1)*pn_2)/n; + } + pn_line = order*(xi*pn - pn_1)/(xi*xi - 1.); + /* ith weight is: wi = 2/(1 - xi^2)(Pn'(xi)^2) */ + weights[i] = 2./((1 - xi*xi)*pn_line*pn_line); + } + return 0; +} diff --git a/lib/glq.h b/lib/glq.h new file mode 100755 index 0000000..84cb04a --- /dev/null +++ b/lib/glq.h @@ -0,0 +1,182 @@ +/* +Functions for implementing a Gauss-Legendre Quadrature numerical integration +(Hildebrand, 1987). + +Usage example +------------- + +To integrate the cossine function from 0 to 90 degrees: + + #include + #include + #include + #include "src/c/glq.h" + + int main(){ + // Create a new glq structure + GLQ *glq; + double result = 0, a = 0, b = 0.5*3.14; + int i; + + glq = glq_new(5, a, b); + + if(glq == NULL){ + printf("malloc error"); + return 1; + } + + // Calculate the integral + for(i = 0; i < glq->order; i++) + result += glq->weights[i]*cos(glq->nodes[i]); + + // Need to multiply by a scale factor of the integration limits + result *= 0.5*(b - a); + + printf("Integral of cossine from 0 to 90 degrees = %lf\n", result); + + // Free allocated memory + glq_free(glq); + + return 0; + } + +References +---------- + +* Hildebrand, F.B (1987): Introduction to numerical analysis. + Courier Dover Publications, 2. ed. +*/ + +#ifndef _TESSEROIDS_GLQ_H_ +#define _TESSEROIDS_GLQ_H_ + + +/** \var GLQ_MAXIT +Max iterations of the root-finder algorithm */ +const int GLQ_MAXIT = 1000; + + +/** \var GLQ_MAXERROR +Max error allowed for the root-finder algorithm */ +const double GLQ_MAXERROR = 0.000000000000001; + + +/** Store the nodes and weights needed for a GLQ integration */ +typedef struct glq_struct +{ + int order; /**< order of the quadrature, ie number of nodes */ + double *nodes; /**< abscissas or discretization points of the quadrature */ + double *weights; /**< weighting coefficients of the quadrature */ + double *nodes_unscaled; /**< nodes in [-1,1] interval */ +} GLQ; + + +/** Make a new GLQ structure and set all the parameters needed + +WARNING: Don't forget to free the memory malloced by this function using +glq_free()! + +Prints error and warning messages using the logging.h module. + +@param order order of the quadrature, ie number of nodes +@param lower lower integration limit +@param upper upper integration limit + +@return GLQ data structure with the nodes and weights calculated. NULL if there + was an error with allocation. +*/ +GLQ * glq_new(int order, double lower, double upper); + + +/** Free the memory allocated to make a GLQ structure + +@param glq pointer to the allocated memory +*/ +void glq_free(GLQ *glq); + + +/** Put the GLQ nodes to the integration limits IN PLACE. + +Will replace the values of glq.nodes with ones in the specified integration +limits. + +In case the GLQ structure was created with glq_new(), the integration limits can +be reset using this function. + +@param lower lower integration limit +@param upper upper integration limit +@param glq pointer to a GLQ structure created with glq_new() and with all + necessary memory allocated + +@return Return code: + - 0: if everything went OK + - 1: if invalid order + - 2: if NULL pointer for nodes or nodes_unscaled +*/ +int glq_set_limits(double lower, double upper, GLQ *glq); +/** Calculates the GLQ nodes using glq_next_root. + +Nodes will be in the [-1,1] interval. To convert them to the integration limits +use glq_scale_nodes + +@param order order of the quadrature, ie how many nodes. Must be >= 2. +@param nodes pre-allocated array to return the nodes. + +@return Return code: + - 0: if everything went OK + - 1: if invalid order + - 2: if NULL pointer for nodes + - 3: if number of maximum iterations was reached when calculating the root. + This usually means that the desired accuracy was not achieved. Default + desired accuracy is GLQ_MAXERROR. Default maximum iterations is + GLQ_MAXIT. +*/ +int glq_nodes(int order, double *nodes); + + +/** Calculate the next Legendre polynomial root given the previous root found. + +Uses the root-finder algorithm of: + + Barrera-Figueroa, V., Sosa-Pedroza, J. and López-Bonilla, J., 2006, + "Multiple root finder algorithm for Legendre and Chebyshev polynomials via + Newton's method", 2006, Annales mathematicae et Informaticae, 33, pp 3-13 + +@param initial initial estimate of the next root. I recommend the use of + \f$ \cos\left(\pi\frac{(N - i - 0.25)}{N + 0.5}\right) \f$, + where \f$ i \f$ is the index of the desired root +@param root_index index of the desired root, starting from 0 +@param order order of the Legendre polynomial, ie number of roots. +@param roots array with the roots found so far. Will return the next root in + roots[root_index], so make sure to malloc enough space. + +@return Return code: + - 0: if everything went OK + - 1: if order is not valid + - 2: if root_index is not valid (negative) + - 3: if number of maximum iterations was reached when calculating the root. + This usually means that the desired accuracy was not achieved. Default + desired accuracy is GLQ_MAXERROR. Default maximum iterations is + GLQ_MAXIT. +*/ +int glq_next_root(double initial, int root_index, int order, + double *roots); + + +/** Calculates the weighting coefficients for the GLQ integration. + +@param order order of the quadrature, ie number of nodes and weights. +@param nodes array containing the GLQ nodes calculated by glq_nodes. + IMPORTANT: needs the nodes in [-1,1] interval! Scaled nodes + will result in wrong weights! +@param weights pre-allocated array to return the weights + +@return Return code: + - 0: if everything went OK + - 1: if order is not valid + - 2: if nodes is a NULL pointer + - 3: if weights is a NULL pointer +*/ +int glq_weights(int order, double *nodes, double *weights); + +#endif diff --git a/lib/grav_tess.cpp b/lib/grav_tess.cpp new file mode 100755 index 0000000..0c93ba9 --- /dev/null +++ b/lib/grav_tess.cpp @@ -0,0 +1,628 @@ +/* +Functions that calculate the gravitational potential and its first and second +derivatives for the tesseroid. + +References +---------- + +* Grombein, T.; Seitz, K.; Heck, B. (2010): Untersuchungen zur effizienten +Berechnung topographischer Effekte auf den Gradiententensor am Fallbeispiel der +Satellitengradiometriemission GOCE. +KIT Scientific Reports 7547, ISBN 978-3-86644-510-9, KIT Scientific Publishing, +Karlsruhe, Germany. +*/ + + +#include +#include "logger.h" +#include "geometry.h" +#include "glq.h" +#include "constants.h" +#include "grav_tess.h" + + +/* Calculates the field of a tesseroid model at a given point. */ +double calc_tess_model(TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, double (*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ)) +{ + double res; + int tess; + + res = 0; + for(tess = 0; tess < size; tess++) + { + if(lonp >= model[tess].w && lonp <= model[tess].e && + latp >= model[tess].s && latp <= model[tess].n && + rp >= model[tess].r1 && rp <= model[tess].r2) + { + log_warning("Point (%g %g %g) is on tesseroid %d: %g %g %g %g %g %g %g. Can't guarantee accuracy.", + lonp, latp, rp - MEAN_EARTH_RADIUS, tess, + model[tess].w, model[tess].e, model[tess].s, + model[tess].n, model[tess].r2 - MEAN_EARTH_RADIUS, + model[tess].r1 - MEAN_EARTH_RADIUS, + model[tess].density); + } + glq_set_limits(model[tess].w, model[tess].e, glq_lon); + glq_set_limits(model[tess].s, model[tess].n, glq_lat); + glq_set_limits(model[tess].r1, model[tess].r2, glq_r); + res += field(model[tess], lonp, latp, rp, *glq_lon, *glq_lat, *glq_r); + } + return res; +} + +void calc_tess_model_triple(TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, void (*field_triple)(TESSEROID, double, double, double, GLQ, GLQ, GLQ, double*), double *res) +{ + double r1, r2, r3, ri[3]; + int tess; + + res[0] = 0; + res[1] = 0; + res[2] = 0; + for(tess = 0; tess < size; tess++) + { + ri[0] = 0; + ri[1] = 0; + ri[2] = 0; + + if(lonp >= model[tess].w && lonp <= model[tess].e && + latp >= model[tess].s && latp <= model[tess].n && + rp >= model[tess].r1 && rp <= model[tess].r2) + { + log_warning("Point (%g %g %g) is on tesseroid %d: %g %g %g %g %g %g %g. Can't guarantee accuracy.", + lonp, latp, rp - MEAN_EARTH_RADIUS, tess, + model[tess].w, model[tess].e, model[tess].s, + model[tess].n, model[tess].r2 - MEAN_EARTH_RADIUS, + model[tess].r1 - MEAN_EARTH_RADIUS, + model[tess].density); + } + glq_set_limits(model[tess].w, model[tess].e, glq_lon); + glq_set_limits(model[tess].s, model[tess].n, glq_lat); + glq_set_limits(model[tess].r1, model[tess].r2, glq_r); + field_triple(model[tess], lonp, latp, rp, *glq_lon, *glq_lat, *glq_r, ri); + + res[0] += ri[0]; + res[1] += ri[1]; + res[2] += ri[2]; + } + return; +} + + + +/* Adaptatively calculate the field of a tesseroid model at a given point */ +double calc_tess_model_adapt(TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, double (*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ), double ratio) +{ + double res, dist, lont, latt, rt, d2r = PI/180.; + int tess; + TESSEROID split[8]; + + res = 0; + for(tess = 0; tess < size; tess++) + { + rt = model[tess].r2; + lont = 0.5*(model[tess].w + model[tess].e); + latt = 0.5*(model[tess].s + model[tess].n); + dist = sqrt(rp*rp + rt*rt - 2*rp*rt*(sin(d2r*latp)*sin(d2r*latt) + + cos(d2r*latp)*cos(d2r*latt)*cos(d2r*(lonp - lont)))); + + /* Would get stuck in infinite loop if dist = 0 and get wrong results if + inside de tesseroid. Still do the calculation but warn user that it's + probably wrong. */ + if(lonp >= model[tess].w && lonp <= model[tess].e && + latp >= model[tess].s && latp <= model[tess].n && + rp >= model[tess].r1 && rp <= model[tess].r2) + { + log_warning("Point (%g %g %g) is on top of tesseroid %d: %g %g %g %g %g %g %g. Can't guarantee accuracy.", + lonp, latp, rp - MEAN_EARTH_RADIUS, tess, + model[tess].w, model[tess].e, model[tess].s, + model[tess].n, model[tess].r2 - MEAN_EARTH_RADIUS, + model[tess].r1 - MEAN_EARTH_RADIUS, + model[tess].density); + glq_set_limits(model[tess].w, model[tess].e, glq_lon); + glq_set_limits(model[tess].s, model[tess].n, glq_lat); + glq_set_limits(model[tess].r1, model[tess].r2, glq_r); + res += field(model[tess], lonp, latp, rp, *glq_lon, *glq_lat, + *glq_r); + } + /* Check if the computation point is at an acceptable distance. If not + split the tesseroid using the given ratio */ + else if( + dist < ratio*MEAN_EARTH_RADIUS*d2r*(model[tess].e - model[tess].w) || + dist < ratio*MEAN_EARTH_RADIUS*d2r*(model[tess].n - model[tess].s) || + dist < ratio*(model[tess].r2 - model[tess].r1)) + { + log_debug("Splitting tesseroid %d (%g %g %g %g %g %g %g) at point (%g %g %g) using ratio %g", + tess, model[tess].w, model[tess].e, model[tess].s, + model[tess].n, model[tess].r2 - MEAN_EARTH_RADIUS, + model[tess].r1 - MEAN_EARTH_RADIUS, model[tess].density, + lonp, latp, rp - MEAN_EARTH_RADIUS, ratio); + /* Do it recursively until ratio*size is smaller than distance */ + split_tess(model[tess], split); + res += calc_tess_model_adapt(split, 8, lonp, latp, rp, glq_lon, + glq_lat, glq_r, field, ratio); + } + else + { + glq_set_limits(model[tess].w, model[tess].e, glq_lon); + glq_set_limits(model[tess].s, model[tess].n, glq_lat); + glq_set_limits(model[tess].r1, model[tess].r2, glq_r); + res += field(model[tess], lonp, latp, rp, *glq_lon, *glq_lat, + *glq_r); + } + } + return res; +} + +/* Calculates gxx caused by a tesseroid. */ +double tess_gxx(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r) +{ + double d2r = PI/180., l_sqr, kphi, coslatp, coslatc, sinlatp, sinlatc, + coslon, rc, kappa, res; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + + l_sqr = rp*rp + rc*rc - 2*rp*rc*(sinlatp*sinlatc + + coslatp*coslatc*coslon); + + kphi = coslatp*sinlatc - sinlatp*coslatc*coslon; + + kappa = rc*rc*coslatc; + + res += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*rc*kphi*rc*kphi - l_sqr)/pow(l_sqr, 2.5); + } + } + } + + res *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + return res; +} + + +/* Calculates gxy caused by a tesseroid. */ +double tess_gxy(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r) +{ + double d2r = PI/180., l_sqr, kphi, coslatp, coslatc, sinlatp, sinlatc, + coslon, sinlon, rc, kappa, deltax, deltay, res; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + sinlon = sin(d2r*(glq_lon.nodes[k] - lonp)); + + l_sqr = rp*rp + rc*rc - 2*rp*rc*(sinlatp*sinlatc + + coslatp*coslatc*coslon); + + kphi = coslatp*sinlatc - sinlatp*coslatc*coslon; + + kappa = rc*rc*coslatc; + + deltax = rc*kphi; + + deltay = rc*coslatc*sinlon; + + res += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltax*deltay)/pow(l_sqr, 2.5); + } + } + } + + res *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + return res; +} + + +/* Calculates gxz caused by a tesseroid. */ +double tess_gxz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r) +{ + double d2r = PI/180., l_sqr, kphi, coslatp, coslatc, sinlatp, sinlatc, + coslon, cospsi, rc, kappa, deltax, deltaz, res; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + + cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon; + + l_sqr = rp*rp + rc*rc - 2*rp*rc*cospsi; + + kphi = coslatp*sinlatc - sinlatp*coslatc*coslon; + + kappa = rc*rc*coslatc; + + deltax = rc*kphi; + + deltaz = rc*cospsi - rp; + + res += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltax*deltaz)/pow(l_sqr, 2.5); + } + } + } + + res *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + return res; +} + + +/* Calculates gyy caused by a tesseroid. */ +double tess_gyy(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r) +{ + double d2r = PI/180., l_sqr, coslatp, coslatc, sinlatp, sinlatc, + coslon, sinlon, rc, kappa, deltay, res; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + sinlon = sin(d2r*(glq_lon.nodes[k] - lonp)); + + l_sqr = rp*rp + rc*rc - 2*rp*rc*(sinlatp*sinlatc + + coslatp*coslatc*coslon); + + kappa = rc*rc*coslatc; + + deltay = rc*coslatc*sinlon; + + res += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltay*deltay - l_sqr)/pow(l_sqr, 2.5); + } + } + } + + res *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + return res; +} + + +/* Calculates gyz caused by a tesseroid. */ +double tess_gyz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r) +{ + double d2r = PI/180., l_sqr, coslatp, coslatc, sinlatp, sinlatc, + coslon, sinlon, cospsi, rc, kappa, deltay, deltaz, res; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + sinlon = sin(d2r*(glq_lon.nodes[k] - lonp)); + + cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon; + + l_sqr = rp*rp + rc*rc - 2*rp*rc*cospsi; + + kappa = rc*rc*coslatc; + + deltay = rc*coslatc*sinlon; + + deltaz = rc*cospsi - rp; + + res += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltay*deltaz)/pow(l_sqr, 2.5); + } + } + } + + res *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + return res; +} + + +/* Calculates gzz caused by a tesseroid. */ +double tess_gzz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r) +{ + double d2r = PI/180., l_sqr, coslatp, coslatc, sinlatp, sinlatc, + coslon, cospsi, rc, kappa, deltaz, res; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + + cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon; + + l_sqr = rp*rp + rc*rc - 2*rp*rc*cospsi; + + kappa = rc*rc*coslatc; + + deltaz = rc*cospsi - rp; + + res += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltaz*deltaz - l_sqr)/pow(l_sqr, 2.5); + } + } + } + + res *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + return res; +} + +/*Calculate three gravity gradient components simultaneously*/ +void tess_gxz_gyz_gzz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r, double *res) +{ + double d2r = PI/180., l_sqr, coslatp, coslatc, sinlatp, sinlatc, sinlon, + coslon, cospsi, rc, kappa, deltaz, deltax, deltay, kphi, + res_gxz, res_gyz, res_gzz; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res_gxz = 0; + res_gyz = 0; + res_gzz = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + sinlon = sin(d2r*(glq_lon.nodes[k] - lonp)); + cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon; + + l_sqr = rp*rp + rc*rc - 2*rp*rc*cospsi; + kphi = coslatp*sinlatc - sinlatp*coslatc*coslon; + kappa = rc*rc*coslatc; + + deltax = rc*kphi; + deltay = rc*coslatc*sinlon; + deltaz = rc*cospsi - rp; + + res_gxz += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltax*deltaz)/pow(l_sqr, 2.5); + + res_gyz += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltay*deltaz)/pow(l_sqr, 2.5); + + res_gzz += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltaz*deltaz - l_sqr)/pow(l_sqr, 2.5); + } + } + } + + res_gxz *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res_gyz *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res_gzz *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res[0] = res_gxz; + res[1] = res_gyz; + res[2] = res_gzz; + + return; +} + +void tess_gxx_gxy_gxz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r, double *res) +{ + double d2r = PI/180., l_sqr, coslatp, coslatc, sinlatp, sinlatc, sinlon, + coslon, cospsi, rc, kappa, deltaz, deltax, deltay, kphi, + res_gxx, res_gxy, res_gxz; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res_gxx = 0; + res_gxy = 0; + res_gxz = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + sinlon = sin(d2r*(glq_lon.nodes[k] - lonp)); + cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon; + + l_sqr = rp*rp + rc*rc - 2*rp*rc*cospsi; + kphi = coslatp*sinlatc - sinlatp*coslatc*coslon; + kappa = rc*rc*coslatc; + + deltax = rc*kphi; + deltay = rc*coslatc*sinlon; + deltaz = rc*cospsi - rp; + + res_gxx += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*rc*kphi*rc*kphi - l_sqr)/pow(l_sqr, 2.5); + + res_gxy += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltax*deltay)/pow(l_sqr, 2.5); + + res_gxz += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltax*deltaz)/pow(l_sqr, 2.5); + } + } + } + + res_gxx *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res_gxy *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res_gxz *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res[0] = res_gxx; + res[1] = res_gxy; + res[2] = res_gxz; + + return; +} + + + +void tess_gxy_gyy_gyz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, + GLQ glq_lat, GLQ glq_r, double *res) +{ + double d2r = PI/180., l_sqr, coslatp, coslatc, sinlatp, sinlatc, sinlon, + coslon, cospsi, rc, kappa, deltaz, deltax, deltay, kphi, + res_gxy, res_gyy, res_gyz; + register int i, j, k; + + coslatp = cos(d2r*latp); + sinlatp = sin(d2r*latp); + + res_gxy = 0; + res_gyy = 0; + res_gyz = 0; + + for(k = 0; k < glq_lon.order; k++) + { + for(j = 0; j < glq_lat.order; j++) + { + for(i = 0; i < glq_r.order; i++) + { + rc = glq_r.nodes[i]; + sinlatc = sin(d2r*glq_lat.nodes[j]); + coslatc = cos(d2r*glq_lat.nodes[j]); + coslon = cos(d2r*(lonp - glq_lon.nodes[k])); + sinlon = sin(d2r*(glq_lon.nodes[k] - lonp)); + cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon; + + l_sqr = rp*rp + rc*rc - 2*rp*rc*cospsi; + kphi = coslatp*sinlatc - sinlatp*coslatc*coslon; + kappa = rc*rc*coslatc; + + deltax = rc*kphi; + deltay = rc*coslatc*sinlon; + deltaz = rc*cospsi - rp; + + res_gxy += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltax*deltay)/pow(l_sqr, 2.5); + + res_gyy += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltay*deltay - l_sqr)/pow(l_sqr, 2.5); + + res_gyz += glq_lon.weights[k]*glq_lat.weights[j]*glq_r.weights[i]* + kappa*(3*deltay*deltaz)/pow(l_sqr, 2.5); + + } + } + } + + res_gxy *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res_gyy *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res_gyz *= SI2EOTVOS*G*tess.density*d2r*(tess.e - tess.w)*d2r*(tess.n - tess.s)* + (tess.r2 - tess.r1)*0.125; + + res[0] = res_gxy; + res[1] = res_gyy; + res[2] = res_gyz; + + return; +} diff --git a/lib/grav_tess.h b/lib/grav_tess.h new file mode 100755 index 0000000..b0cf14f --- /dev/null +++ b/lib/grav_tess.h @@ -0,0 +1,93 @@ +/* +Functions that calculate the gravitational potential and its first and second +derivatives for the tesseroid. + +The gravity gradients can be calculated using the general formula of +Grombein et al. (2010). +The integrals are solved using the Gauss-Legendre Quadrature rule +(Asgharzadeh et al., 2007). + +The derivatives of the potential are made with respect to the local coordinate +system x->North, y->East, z->Up (away from center of the Earth). + +To maintain the standard convention, only for component gz the z axis is +inverted, so a positive density results in positive gz. + +Example +------- + +To calculate the gzz component due to a tesseroid on a regular grid: + + #include + #include "glq.h"r + #include "constants.h" + #include "grav_tess.h" + + int main() + { + TESSEROID tess = {1000, 44, 46, -1, 1, MEAN_EARTH_RADIUS - 100000, + MEAN_EARTH_RADIUS}; + GLQ *glqlon, *glqlat, *glqr; + double lon, lat, r = MEAN_EARTH_RADIUS + 1500000, res; + int order = 8; + + glqlon = glq_new(order, tess.w, tess.e); + glqlat = glq_new(order, tess.s, tess.n); + glqr = glq_new(order, tess.r1, tess.r2); + + for(lat = 20; lat <= 70; lat += 0.5) + { + for(lon = -25; lon <= 25; lon += 0.5) + { + res = tess_gzz(tess, lon, lat, r, *glqlon, *glqlat, *glqr); + printf("%g %g %g\n", lon, lat, res); + } + } + + glq_free(glqlon); + glq_free(glqlat); + glq_free(glqr); + + return 0; + } + +References +---------- + +Asgharzadeh, M.F., von Frese, R.R.B., Kim, H.R., Leftwich, T.E. & Kim, J.W. +(2007): Spherical prism gravity effects by Gauss-Legendre quadrature integration. +Geophysical Journal International, 169, 1-11. + +Grombein, T.; Seitz, K.; Heck, B. (2010): Untersuchungen zur effizienten +Berechnung topographischer Effekte auf den Gradiententensor am Fallbeispiel der +Satellitengradiometriemission GOCE. +KIT Scientific Reports 7547, ISBN 978-3-86644-510-9, KIT Scientific Publishing, +Karlsruhe, Germany. +*/ + +#ifndef _TESSEROIDS_GRAV_TESS_H_ +#define _TESSEROIDS_GRAV_TESS_H_ + + +/* Needed for definition of TESSEROID */ +#include "geometry.h" +/* Needed for definition of GLQ */ +#include "glq.h" + +double calc_tess_model(TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, double (*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ)); +void calc_tess_model_triple(TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, + void (*field_triple)(TESSEROID, double, double, double, GLQ, GLQ, GLQ, double*), double *res); +double calc_tess_model_adapt(TESSEROID *model, int size, double lonp, double latp, double rp, GLQ *glq_lon, GLQ *glq_lat, GLQ *glq_r, double (*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ), double ratio); + +double tess_gxx(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r); +double tess_gxy(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r); +double tess_gxz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r); +double tess_gyy(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r); +double tess_gyz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r); +double tess_gzz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r); + +void tess_gxz_gyz_gzz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r, double *res); +void tess_gxx_gxy_gxz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r, double *res); +void tess_gxy_gyy_gyz(TESSEROID tess, double lonp, double latp, double rp, GLQ glq_lon, GLQ glq_lat, GLQ glq_r, double *res); + +#endif diff --git a/lib/linalg.cpp b/lib/linalg.cpp new file mode 100755 index 0000000..82d028a --- /dev/null +++ b/lib/linalg.cpp @@ -0,0 +1,142 @@ +/* +Functions matrix and vector multiplications. + +*/ + +#include "linalg.h" +#include "constants.h" +#include + +//macOS only! +//#include + +#ifdef __linux__ // Debian, Ubuntu, Gentoo, Fedora, openSUSE, RedHat, Centos and other + #include +#elif defined(__APPLE__) && defined(__MACH__) + #include +#else + // +#endif + +/* Calculate magnetization vector in a coordinate system of a given point */ +void conv_vect_cblas(double *vect, double lon1, double lat1, double lon2, double lat2, double *res) +{ + double a1 = DEG2RAD*lat1; + double b1 = DEG2RAD*lon1; + double a2 = DEG2RAD*lat2; + double b2 = DEG2RAD*lon2; + + double cos_a1 = cos(PI/2.0-a1); + double sin_a1 = sin(PI/2.0-a1); + double cos_a2 = cos(PI/2.0-a2); + double sin_a2 = sin(PI/2.0-a2); + + double cos_b1 = cos(b1); + double sin_b1 = sin(b1); + double cos_b2 = cos(b2); + double sin_b2 = sin(b2); + + double Z1Y1[9] = {cos_a1*cos_b1, -sin_b1, cos_b1*sin_a1, cos_a1*sin_b1, cos_b1, sin_a1*sin_b1, -sin_a1, 0, cos_a1}; + double Z2Y2t[9] = {-cos_a2*cos_b2, -cos_a2*sin_b2, sin_a2, -sin_b2, cos_b2, 0, cos_b2*sin_a2, sin_a2*sin_b2, cos_a2}; + double R[9]; + + cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3, 1.0, Z1Y1, 3, Z2Y2t, 3, 0.0, R, 3); + + R[0] = -R[0]; + R[3] = -R[3]; + R[6] = -R[6]; + + cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1.0, R, 3, vect, 1, 0.0, res, 1); + + + return; +} + + +void conv_vect_cblas_precalc(double *vect, double cos_a1, double sin_a1, double cos_b1, double sin_b1, double cos_a2, double sin_a2, double cos_b2, double sin_b2, double *res) +{ +// double a1 = DEG2RAD*lat1; +// double b1 = DEG2RAD*lon1; +// double a2 = DEG2RAD*lat2; +// double b2 = DEG2RAD*lon2; + +// double cos_a1 = cos(PI/2.0-a1); +// double sin_a1 = sin(PI/2.0-a1); +// double cos_a2 = cos(PI/2.0-a2); +// double sin_a2 = sin(PI/2.0-a2); + +// double cos_b1 = cos(b1); +// double sin_b1 = sin(b1); +// double cos_b2 = cos(b2); +// double sin_b2 = sin(b2); + + double Z1Y1[9] = {cos_a1*cos_b1, -sin_b1, cos_b1*sin_a1, cos_a1*sin_b1, cos_b1, sin_a1*sin_b1, -sin_a1, 0, cos_a1}; + double Z2Y2t[9] = {-cos_a2*cos_b2, -cos_a2*sin_b2, sin_a2, -sin_b2, cos_b2, 0, cos_b2*sin_a2, sin_a2*sin_b2, cos_a2}; + double R[9]; + + cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3, 1.0, Z1Y1, 3, Z2Y2t, 3, 0.0, R, 3); + + R[0] = -R[0]; + R[3] = -R[3]; + R[6] = -R[6]; + + cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1.0, R, 3, vect, 1, 0.0, res, 1); + + + return; +} + +void from_loc_sphr_to_cart(double* columnvect_xyzloc, double colatitude, double longitude, double* columnvect_res) +{ +/*IMPORTANT: this subroutine is in the coordinate system NED*/ + + double phi = colatitude*DEG2RAD; + double lambda = longitude*DEG2RAD; + + double cos_phi = cos(phi); + double sin_phi = sin(phi); + double cos_lambda = cos(lambda); + double sin_lambda = sin(lambda); + + double columnvect_phi_unit[3] = {-sin_phi*cos_lambda, -sin_phi*sin_lambda, cos_phi}; + double columnvect_lambda_unit[3] = {-sin_lambda, cos_lambda, 0}; + double columnvect_r_unit[3] = {cos_phi*cos_lambda, cos_phi*sin_lambda, sin_phi}; + + columnvect_res[0] = columnvect_phi_unit[0]*columnvect_xyzloc[0]+columnvect_lambda_unit[0]*columnvect_xyzloc[1]+columnvect_r_unit[0]*columnvect_xyzloc[2]; + columnvect_res[1] = columnvect_phi_unit[1]*columnvect_xyzloc[0]+columnvect_lambda_unit[1]*columnvect_xyzloc[1]+columnvect_r_unit[1]*columnvect_xyzloc[2]; + columnvect_res[2] = columnvect_phi_unit[2]*columnvect_xyzloc[0]+columnvect_lambda_unit[2]*columnvect_xyzloc[1]+columnvect_r_unit[2]*columnvect_xyzloc[2]; + + return; +} + +void from_cart_to_loc_sphr(double* columnvect_xyzglob, double colatitude, double longitude, double* columnvect_res) +{ +/*IMPORTANT: this subroutine is in the coordinate system NED*/ + + double phi = colatitude*DEG2RAD; + double lambda = longitude*DEG2RAD; + + double cos_phi = cos(phi); + double sin_phi = sin(phi); + double cos_lambda = cos(lambda); + double sin_lambda = sin(lambda); + + double rowvect_phi_unit[3] = {-sin_phi*cos_lambda, -sin_phi*sin_lambda, cos_phi}; + double rowvect_lambda_unit[3] = {-sin_lambda, cos_lambda, 0}; + double rowvect_r_unit[3] = {cos_phi*cos_lambda, cos_phi*sin_lambda, sin_phi}; + + columnvect_res[0] = rowvect_phi_unit[0]*columnvect_xyzglob[0]+rowvect_phi_unit[1]*columnvect_xyzglob[1]+rowvect_phi_unit[2]*columnvect_xyzglob[2]; + columnvect_res[1] = rowvect_lambda_unit[0]*columnvect_xyzglob[0]+rowvect_lambda_unit[1]*columnvect_xyzglob[1]+rowvect_lambda_unit[2]*columnvect_xyzglob[2]; + columnvect_res[2] = rowvect_r_unit[0]*columnvect_xyzglob[0]+rowvect_r_unit[1]*columnvect_xyzglob[1]+rowvect_r_unit[2]*columnvect_xyzglob[2]; + + return; +} + +void from_loc_sphr_to_loc_sphr(double* columnvect_xyzloc, double colatitude1, double longitude1, double colatitude2, double longitude2, double* columnvect_res) +{ + double columnvect_xyzglob[3]; + from_loc_sphr_to_cart(columnvect_xyzloc, colatitude1, longitude1, columnvect_xyzglob); + from_cart_to_loc_sphr(columnvect_xyzglob, colatitude2, longitude2, columnvect_res); + return; + +} diff --git a/lib/linalg.h b/lib/linalg.h new file mode 100755 index 0000000..3e5ccf6 --- /dev/null +++ b/lib/linalg.h @@ -0,0 +1,11 @@ +#ifndef _LINALG_H_ +#define _LINALG_H_ + +void conv_vect_cblas(double *vect, double lon1, double lat1, double lon2, double lat2, double *res); +void conv_vect_cblas_precalc(double *vect, double cos_a1, double sin_a1, double cos_b1, double sin_b1, double cos_a2, double sin_a2, double cos_b2, double sin_b2, double *res); + +void from_loc_sphr_to_cart(double* columnvect_xyzloc, double colatitude, double longitude, double* columnvect_res); +void from_cart_to_loc_sphr(double* columnvect_xyzglob, double colatitude, double longitude, double* columnvect_res); +void from_loc_sphr_to_loc_sphr(double* columnvect_xyzloc, double colatitude1, double longitude1, double colatitude2, double longitude2, double* columnvect_res); + +#endif diff --git a/lib/logger.cpp b/lib/logger.cpp new file mode 100755 index 0000000..cb1b573 --- /dev/null +++ b/lib/logger.cpp @@ -0,0 +1,112 @@ +/* +Functions to set up logging. +*/ + + +#include +#include +#include +#include "logger.h" + +/* Initialize the logger so that it doesn't print by default */ +LOGGER logger = {100, 0, 100, NULL}; + + +/* Setup logging to stderr.*/ +void log_init(int level) +{ + logger.level = level; +} + + +/* Set logging to a file. */ +void log_tofile(FILE *logfile, int level) +{ + logger.filelogging = 1; + logger.logfile = logfile; + logger.file_level = level; +} + + +/* Log a message at debug level */ +void log_debug(const char *fmt, ...) +{ + char msg[10000]; + va_list args; + va_start(args, fmt); + vsprintf(msg, fmt, args); + va_end(args); + + if(logger.level <= LOG_DEBUG) + { + fprintf(stderr, "DEBUG: %s\n", msg); + } + + if(logger.filelogging && (logger.file_level <= LOG_DEBUG)) + { + fprintf(logger.logfile, "DEBUG: %s\n", msg); + } +} + + +/* Log a message at info level */ +void log_info(const char *fmt, ...) +{ + char msg[10000]; + va_list args; + va_start(args, fmt); + vsprintf(msg, fmt, args); + va_end(args); + + if(logger.level <= LOG_INFO) + { + fprintf(stderr, "%s\n", msg); + } + + if(logger.filelogging && logger.file_level <= LOG_INFO) + { + fprintf(logger.logfile, "%s\n", msg); + } +} + + +/* Log a message at warning level */ +void log_warning(const char *fmt, ...) +{ + char msg[10000]; + va_list args; + va_start(args, fmt); + vsprintf(msg, fmt, args); + va_end(args); + + if(logger.level <= LOG_WARNING) + { + fprintf(stderr, "WARNING: %s\n", msg); + } + + if(logger.filelogging && logger.file_level <= LOG_WARNING) + { + fprintf(logger.logfile, "WARNING: %s\n", msg); + } +} + + +/* Log a message at error level */ +void log_error(const char *fmt, ...) +{ + char msg[10000]; + va_list args; + va_start(args, fmt); + vsprintf(msg, fmt, args); + va_end(args); + + if(logger.level <= LOG_ERROR) + { + fprintf(stderr, "ERROR: %s\n", msg); + } + + if(logger.filelogging && logger.file_level <= LOG_ERROR) + { + fprintf(logger.logfile, "ERROR: %s\n", msg); + } +} diff --git a/lib/logger.h b/lib/logger.h new file mode 100755 index 0000000..8472fee --- /dev/null +++ b/lib/logger.h @@ -0,0 +1,166 @@ +/* +Functions to set up logging. + +Examples +-------- + +Logging to stderr: + + #include "logger.h" + + void my_func(){ + log_info("From my_func!\n"); + } + + int main(){ + // Enable logging to stderr in debug level + // will only print messages of level DEBUG or higher + log_init(LOG_DEBUG); + log_debug("debug line. The code is %d", LOG_DEBUG); + log_info("info line. The code is %d", LOG_INFO); + log_warning("warning line. The code is %d", LOG_WARNING); + log_error("error line. The code is %d", LOG_ERROR); + my_func(); + return 0; + } + +will print: + + DEBUG: debug line. The code is 0 + info line. The code is 1 + WARNING: warning line. The code is 2 + ERROR: error line. The code is 3 + From my_func! + +If function log_init() is not called than logging to stderr is disabled and no +messages will be printed. + +Logging to a file: + + #include + #include "logger.h" + + void my_func(){ + log_info("From my_func!\n"); + log_debug("Should not appear in log file\n"); + } + + int main(){ + // Enable logging to file "log.txt" in info level + // will not print DEBUG level messages + // since log_init was not called, there is no logging to stderr + FILE *logfile = fopen("log.txt", "w"); + log_tofile(logfile, LOG_INFO); + log_debug("debug line. The code is %d", LOG_DEBUG); + log_info("info line. The code is %d", LOG_INFO); + log_warning("warning line. The code is %d", LOG_WARNING); + log_error("error line. The code is %d", LOG_ERROR); + my_func(); + return 0; + } + +File log.txt will look like: + + info line. The code is 1 + WARNING: warning line. The code is 2 + ERROR: error line. The code is 3 + From my_func! + +Note that you can combine loggin to stderr and to a file with different +levels in the same program. +*/ + +#ifndef _TESSEROIDS_LOGGER_H_ +#define _TESSEROIDS_LOGGER_H_ + +/* Needed for definition of FILE */ +#include + + +/** Logging level for debug messages */ +#define LOG_DEBUG 1 +/** Logging level for general information */ +#define LOG_INFO 2 +/** Logging level for warning messages */ +#define LOG_WARNING 3 +/** Logging level for error messages */ +#define LOG_ERROR 4 + + +/** Keep the information on the global logger */ +typedef struct logger_struct +{ + int level; /**< level of logging */ + int filelogging; /**< flag to know wether loggint to a file is enabled */ + int file_level; /**< logging level for the file */ + FILE *logfile; /**< file to log to */ + +} LOGGER; + + +/** Global logger struct. Only declare in the main program! */ +//LOGGER logger; + + +/** Setup logging to stderr. + +@param level level of logging to be made. Can be one of: + - LOG_DEBUG + - LOG_INFO + - LOG_WARNING + - LOG_ERROR +*/ +void log_init(int level); + + +/** Set logging to a file. + +@param logfile FILE pointer to the already open file to log to. +@param level level of logging to be made to the file. Can be one of: + - LOG_DEBUG + - LOG_INFO + - LOG_WARNING + - LOG_ERROR +*/ +void log_tofile(FILE *logfile, int level); + + +/** Log a message at debug level. + +Pass parameters in the same format as printf() + +Prints a newline at the end. +*/ +void log_debug(const char *fmt, ...); + + +/** Log a message at info level. + +Pass parameters in the same format as printf() + +Does not print "INFO: " in front of the message when logging + +Prints a newline at the end. +*/ +void log_info(const char *fmt, ...); + + +/** Log a message at warning level. + +Pass parameters in the same format as printf() + +Prints a newline at the end. +*/ +void log_warning(const char *fmt, ...); + + +/** Log a message at error level. + +Pass parameters in the same format as printf() + +Prints a newline at the end. +*/ +void log_error(const char *fmt, ...); + + +#endif diff --git a/lib/parsers.cpp b/lib/parsers.cpp new file mode 100755 index 0000000..38a3a48 --- /dev/null +++ b/lib/parsers.cpp @@ -0,0 +1,678 @@ +/* +Input and output parsing tools. +*/ + + +#include +#include +#include +#include +#include "logger.h" +#include "version.h" +#include "parsers.h" +#include "constants.h" +#include "geometry.h" + +#include + + +/* Parse basic command line arguments for programs */ +int parse_basic_args(int argc, char **argv, const char *progname, + BASIC_ARGS *args, void (*print_help)(void)) +{ + int bad_args = 0, parsed_args = 0, total_args = 1, i; + char *params; + + /* Default values for options */ + args->verbose = 0; + args->logtofile = 0; + /* Parse arguments */ + for(i = 1; i < argc; i++) + { + if(argv[i][0] == '-') + { + switch(argv[i][1]) + { + case 'h': + if(argv[i][2] != '\0') + { + log_error("invalid argument '%s'", argv[i]); + bad_args++; + break; + } + print_help(); + return 2; + case 'v': + if(argv[i][2] != '\0') + { + log_error("invalid argument '%s'", argv[i]); + bad_args++; + break; + } + if(args->verbose) + { + log_error("repeated option -v"); + bad_args++; + break; + } + args->verbose = 1; + break; + case 'l': + { + if(args->logtofile) + { + log_error("repeated option -l"); + bad_args++; + break; + } + params = &argv[i][2]; + if(strlen(params) == 0) + { + log_error("bad input argument -l. Missing filename."); + bad_args++; + } + else + { + args->logtofile = 1; + args->logfname = params; + } + break; + } + case '-': + { + params = &argv[i][2]; + if(strcmp(params, "version")) + { + log_error("invalid argument '%s'", argv[i]); + bad_args++; + } + else + { + print_version(progname); + return 2; + } + break; + } + default: + log_error("invalid argument '%s'", argv[i]); + bad_args++; + break; + } + } + else + { + if(parsed_args == 0) + { + args->inputfname = argv[i]; + parsed_args++; + } + else + { + log_error("invalid argument '%s'. Already given model file %s", + argv[i], args->inputfname); + bad_args++; + } + } + } + /* Check if parsing went well */ + if(parsed_args > total_args) + { + log_error("%s: too many input arguments. given %d, max %d.", + progname, parsed_args, total_args); + } + if(bad_args > 0) + { + log_error("%d bad input argument(s)", bad_args); + return 1; + } + if(parsed_args < total_args) + { + return 3; + } + return 0; +} + + + +/* Parse command line arguments for tessh* programs */ +int parse_tessb_args(int argc, char **argv, const char *progname, + TESSB_ARGS *args, void (*print_help)(const char *)) +{ + int bad_args = 0, parsed_args = 0, total_args = 1, parsed_order = 0, + parsed_ratio1 = 0, parsed_ratio2 = 0, parsed_ratio3 = 0, i, nchar, nread; + char *params; + + /* Default values for options */ + args->verbose = 0; + args->logtofile = 0; + args->lon_order = 2; + args->lat_order = 2; + args->r_order = 2; + args->adaptative = 1; + args->ratio1 = 0; /* zero means use the default for the program */ + args->ratio2 = 0; + args->ratio3 = 0; + /* Parse arguments */ + for(i = 1; i < argc; i++) + { + if(argv[i][0] == '-') + { + switch(argv[i][1]) + { + case 'h': + if(argv[i][2] != '\0') + { + log_error("invalid argument '%s'", argv[i]); + bad_args++; + break; + } + print_help(progname); + return 2; + case 'v': + if(argv[i][2] != '\0') + { + log_error("invalid argument '%s'", argv[i]); + bad_args++; + break; + } + if(args->verbose) + { + log_error("repeated option -v"); + bad_args++; + break; + } + args->verbose = 1; + break; + case 'l': + { + if(args->logtofile) + { + log_error("repeated option -l"); + bad_args++; + break; + } + params = &argv[i][2]; + if(strlen(params) == 0) + { + log_error("bad input argument -l. Missing filename."); + bad_args++; + } + else + { + args->logtofile = 1; + args->logfname = params; + } + break; + } + case '-': + { + params = &argv[i][2]; + if(strcmp(params, "version")) + { + log_error("invalid argument '%s'", argv[i]); + bad_args++; + } + else + { + print_version(progname); + return 2; + } + break; + } + case 'a': + if(argv[i][2] != '\0') + { + log_error("invalid argument '%s'", argv[i]); + bad_args++; + break; + } + if(!args->adaptative) + { + log_error("repeated option -a"); + bad_args++; + break; + } + args->adaptative = 0; + break; + case 'o': + { + if(parsed_order) + { + log_error("repeated option -o"); + bad_args++; + break; + } + params = &argv[i][2]; + nchar = 0; + nread = sscanf(params, "%d/%d/%d%n", &(args->lon_order), + &(args->lat_order), &(args->r_order), &nchar); + if(nread != 3 || *(params + nchar) != '\0') + { + log_error("bad input argument '%s'", argv[i]); + bad_args++; + } + parsed_order = 1; + break; + } + case 't': + { + //ELDAR BAYKIEV/////////////////////////////////////////////////////////////////// + switch(argv[i][2]) + { + case '1': + { + if(parsed_ratio1) + { + log_error("repeated option -t1"); + bad_args++; + break; + } + params = &argv[i][3]; + nchar = 0; + nread = sscanf(params, "%lf%n", &(args->ratio1), &nchar); + if(nread != 1 || *(params + nchar) != '\0') + { + log_error("bad input argument '%s'", argv[i]); + bad_args++; + } + parsed_ratio1 = 1; + break; + } + + case '2': + { + if(parsed_ratio2) + { + log_error("repeated option -t2"); + bad_args++; + break; + } + params = &argv[i][3]; + nchar = 0; + nread = sscanf(params, "%lf%n", &(args->ratio2), &nchar); + if(nread != 1 || *(params + nchar) != '\0') + { + log_error("bad input argument '%s'", argv[i]); + bad_args++; + } + parsed_ratio2 = 1; + break; + } + + case '3': + { + if(parsed_ratio3) + { + log_error("repeated option -t3"); + bad_args++; + break; + } + params = &argv[i][3]; + nchar = 0; + nread = sscanf(params, "%lf%n", &(args->ratio3), &nchar); + if(nread != 1 || *(params + nchar) != '\0') + { + log_error("bad input argument '%s'", argv[i]); + bad_args++; + } + parsed_ratio3 = 1; + break; + } + } + //ELDAR BAYKIEV/////////////////////////////////////////////////////////////////// + } + default: + log_error("invalid argument '%s'", argv[i]); + bad_args++; + break; + } + } + else + { + if(parsed_args == 0) + { + args->modelfname = argv[i]; + parsed_args++; + } + else + { + log_error("invalid argument '%s'. Already given model file %s", + argv[i], args->modelfname); + bad_args++; + } + } + } + /* Check if parsing went well */ + if(bad_args > 0 || parsed_args != total_args) + { + if(parsed_args < total_args) + { + log_error("%s: missing input file.", + progname, parsed_args, total_args); + } + if(parsed_args > total_args) + { + log_error("%s: too many input arguments. given %d, max %d.", + progname, parsed_args, total_args); + } + if(bad_args > 0) + { + log_error("%d bad input argument(s)", bad_args); + } + return 1; + } + return 0; +} + +//parse arguments for gradient calculator +int parse_gradcalc_args(int argc, char **argv, const char *progname, GRADCALC_ARGS *args, void (*print_help)(const char *)) +{ + int bad_args = 0, parsed_args = 0, total_args = 5, i; + char *params; + + /* Default values for options */ + args->verbose = 0; + args->logtofile = 0; + + args->gridbx_set = FALSE; + args->gridby_set = FALSE; + args->gridbz_set = FALSE; + args->out_set = -1; + + args->bz_NEU_NED = -1; + args->bz_NEU_NED_set = FALSE; + + + /* Parse arguments */ + for(i = 1; i < argc; i++) + { + if(argv[i][0] == '-') + { + switch(argv[i][1]) + { + case 'h': + if(argv[i][2] != '\0') + { + //log_error("invalid argument '%s'", argv[i]); + printf("invalid argument '%s'\n", argv[i]); + bad_args++; + break; + } + print_help(progname); + return 2; + case '-': + { + params = &argv[i][2]; + if(strcmp(params, "version")) + { + printf("invalid argument '%s'\n", argv[i]); + bad_args++; + } + else + { + print_version(progname); + return 2; + } + break; + } + case 'b': + params = &argv[i][2]; + if(strlen(params) <= 1) + { + printf("bad input argument -b. Missing component and filename\n"); + bad_args++; + break; + } + else + { + switch(argv[i][2]) + { + case 'x': + if(args->gridbx_set) + { + printf("invalid argument '%s', gridfile bx already set\n", argv[i]); + bad_args++; + break; + } + else + { + args->gridbx_set = 1; + args->gridbx_fn = &argv[i][3]; + } + break; + case 'y': + if(args->gridby_set) + { + printf("invalid argument '%s', gridfile by already set\n", argv[i]); + bad_args++; + break; + } + else + { + args->gridby_set = 1; + args->gridby_fn = &argv[i][3]; + } + break; + case 'z': + if(args->gridbz_set) + { + printf("invalid argument '%s', gridfile by already set\n", argv[i]); + bad_args++; + break; + } + else + { + args->gridbz_set = 1; + args->gridbz_fn = &argv[i][3]; + } + break; + default: + printf("invalid argument '%s', wrong component\n", argv[i]); + bad_args++; + break; + } + } + break; + case 'c': + params = &argv[i][2]; + + if(args->bz_NEU_NED_set) + { + printf("invalid argument '%s', coordinate system is already set\n", argv[i]); + bad_args++; + break; + } + if(strlen(params) > 1) + { + printf("invalid argument '%s', specify coordinate system in the input grids\n", argv[i]); + bad_args++; + break; + } + if(argv[i][2] == '1') + { + args->bz_NEU_NED_set = 1; + args->bz_NEU_NED = 1; + break; + } + else if(argv[i][2] == '2') + { + args->bz_NEU_NED_set = 1; + args->bz_NEU_NED = -1; + break; + } + else + { + printf("invalid argument '%s', there are only NED (1) and NEU (2, default) coordinate systems\n", argv[i]); + bad_args++; + break; + } + break; + case 'o': + params = &argv[i][2]; + + if(args->out_set>=0) + { + printf("invalid argument '%s', output format is already set\n", argv[i]); + bad_args++; + break; + } + + if(strlen(params) != 1) + { + printf("invalid argument '%s', specify output format\n", argv[i]); + bad_args++; + break; + } + //TODO Add check if it is integer + args->out_set = atoi(params); + break; + default: + printf("invalid argument '%s'\n", argv[i]); + bad_args++; + break; + } + } + } + + if(parsed_args > total_args) + { + //log_error("%s: too many input arguments. given %d, max %d.", progname, parsed_args, total_args); + } + if(bad_args > 0) + { + //log_error("%d bad input argument(s)", bad_args); + return 1; + } + if(parsed_args < total_args) + { + return 3; + } + return 0; +} + +/* Strip trailing spaces and newlines from the end of a string */ +void strstrip(char *str) +{ + int i; + for(i = strlen(str) - 1; i >= 0; i--) + { + if(str[i] != ' ' && str[i] != '\n' && str[i] != '\r' && str[i] != '\0') + break; + } + str[i + 1] = '\0'; +} + + +/* Read a single tesseroid from a string */ +int gets_mag_tess(const char *str, TESSEROID *tess) +{ + double w, e, s, n, top, bot, dens, suscept, Bx, By, Bz, Rx, Ry, Rz; + int nread, nchars; + + nread = sscanf(str, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf%n", &w, &e, &s, + &n, &top, &bot, &dens, &suscept, &Bx, &By, &Bz, &nchars); + if(nread != 11 || str[nchars] != '\0') + { + return 1; + } + tess->w = w; + tess->e = e; + tess->s = s; + tess->n = n; + tess->r1 = MEAN_EARTH_RADIUS + bot; + tess->r2 = MEAN_EARTH_RADIUS + top; + tess->density = dens; + tess->suscept = suscept; + tess->Bx = Bx; + tess->By = By; + tess->Bz = Bz; + + tess->cos_a1 = cos(PI/2.0-DEG2RAD*(w+e)*0.5); + tess->sin_a1 = sin(PI/2.0-DEG2RAD*(w+e)*0.5); + tess->cos_b1 = cos(DEG2RAD*(s+n)*0.5); + tess->sin_b1 = sin(DEG2RAD*(s+n)*0.5); + return 0; +} + +//ELDAR BAYKIEV//////////////////////////////// +TESSEROID * read_mag_tess_model(FILE *modelfile, int *size) +{ + TESSEROID *model, *tmp; + int buffsize = 300, line, badinput = 0, error_exit = 0; + char sbuff[10000]; + + /* Start with a single buffer allocation and expand later if necessary */ + model = (TESSEROID *)malloc(buffsize*sizeof(TESSEROID)); + if(model == NULL) + { + log_error("problem allocating initial memory to load tesseroid model."); + return NULL; + } + *size = 0; + for(line = 1; !feof(modelfile); line++) + { + if(fgets(sbuff, 10000, modelfile) == NULL) + { + if(ferror(modelfile)) + { + log_error("problem encountered reading line %d.", line); + error_exit = 1; + break; + } + } + else + { + /* Check for comments and blank lines */ + if(sbuff[0] == '#' || sbuff[0] == '\r' || sbuff[0] == '\n') + { + continue; + } + if(*size == buffsize) + { + buffsize += buffsize; + tmp = (TESSEROID *)realloc(model, buffsize*sizeof(TESSEROID)); + if(tmp == NULL) + { + /* Need to free because realloc leaves unchanged in case of + error */ + free(model); + log_error("problem expanding memory for tesseroid model.\nModel is too big."); + return NULL; + } + model = tmp; + } + /* Remove any trailing spaces or newlines */ + strstrip(sbuff); + if(gets_mag_tess(sbuff, &model[*size])) + { + log_warning("bad/invalid tesseroid at line %d.", line); + badinput = 1; + continue; + } + (*size)++; + } + } + if(badinput || error_exit) + { + free(model); + return NULL; + } + /* Adjust the size of the model */ + if(*size != 0) + { + tmp = (TESSEROID *)realloc(model, (*size)*sizeof(TESSEROID)); + if(tmp == NULL) + { + /* Need to free because realloc leaves unchanged in case of + error */ + free(model); + log_error("problem freeing excess memory for tesseroid model."); + return NULL; + } + model = tmp; + } + return model; +} + + +/* Read a single rectangular prism from a string */ diff --git a/lib/parsers.h b/lib/parsers.h new file mode 100755 index 0000000..f351091 --- /dev/null +++ b/lib/parsers.h @@ -0,0 +1,75 @@ +/* +Input and output parsing tools. +*/ + + +#ifndef _TESSEROIDS_PARSERS_H_ +#define _TESSEROIDS_PARSERS_H_ + +/* Needed for definition of TESSEROID and PRISM */ +#include "geometry.h" +/* Need for the definition of FILE */ +#include + +/** Store basic input arguments and option flags */ +typedef struct basic_args +{ + char *inputfname; /**< name of the input file */ + int verbose; /**< flag to indicate if verbose printing is enabled */ + int logtofile; /**< flag to indicate if logging to a file is enabled */ + char *logfname; /**< name of the log file */ +} BASIC_ARGS; + +typedef struct tessh_args +{ + int lon_order; /**< glq order in longitude integration */ + int lat_order; /**< glq order in latitude integration */ + int r_order; /**< glq order in radial integration */ + char *modelfname; /**< name of the file with the tesseroid model */ + int verbose; /**< flag to indicate if verbose printing is enabled */ + int logtofile; /**< flag to indicate if logging to a file is enabled */ + char *logfname; /**< name of the log file */ + int adaptative; /**< flat to indicate wether to use the adaptative size + of tesseroid algorithm */ + double ratio1; /**< distance-size ratio used for recusive division */ + double ratio2; /**< distance-size ratio used for recusive division */ + double ratio3; /**< distance-size ratio used for recusive division */ +} TESSB_ARGS; + + +typedef struct gradcalc_args +{ + int gridbx_set; + int gridby_set; + int gridbz_set; + + char* gridbx_fn; + char* gridby_fn; + char* gridbz_fn; + + int out_set; + + + int bz_NEU_NED; + int bz_NEU_NED_set; + + int verbose; /**< flag to indicate if verbose printing is enabled */ + int logtofile; /**< flag to indicate if logging to a file is enabled */ + + + + +} GRADCALC_ARGS; + +int parse_basic_args(int argc, char **argv, const char *progname, BASIC_ARGS *args, void (*print_help)(void)); +int parse_tessb_args(int argc, char **argv, const char *progname, TESSB_ARGS *args, void (*print_help)(const char *)); + +int parse_gradcalc_args(int argc, char **argv, const char *progname, GRADCALC_ARGS *args, void (*print_help)(const char *)); + +void strstrip(char *str); + +int gets_mag_tess(const char *str, TESSEROID *tess); +TESSEROID * read_mag_tess_model(FILE *modelfile, int *size); + + +#endif diff --git a/lib/tessb_main.cpp b/lib/tessb_main.cpp new file mode 100755 index 0000000..4f93b9f --- /dev/null +++ b/lib/tessb_main.cpp @@ -0,0 +1,385 @@ +/* +Generic main function for the tessb* programs. +*/ + + + + +#include +#include +#include +#include +#include "logger.h" +#include "version.h" +#include "grav_tess.h" +#include "glq.h" +#include "constants.h" +#include "geometry.h" +#include "parsers.h" +#include "tessb_main.h" +#include "linalg.h" + +#include + +/* Print the help message for tessh* programs */ +void print_tessb_help(const char *progname) +{ + printf("MAGNETIC TESSEROIDS\n"); + printf("Usage: %s MODELFILE [OPTIONS]\n\n", progname); + if(strcmp(progname + 4, "pot") == 0) + { + printf("Calculate the potential due to a tesseroid model on\n"); + } + else + { + printf("Calculate the %s component due to a tesseroid model on\n", + progname + 4); + } + printf("TODO\n"); +} + + +/* Run the main for a generic tessh* program */ +int run_tessb_main(int argc, char **argv, const char *progname, + double (*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ), + double ratio1, double ratio2, double ratio3) +{ + TESSB_ARGS args; + GLQ *glq_lon, *glq_lat, *glq_r; + TESSEROID *model; + + int modelsize, rc, line, points = 0, error_exit = 0, bad_input = 0; + char buff[10000]; + + double lon, lat, height, res; + + /*variables for precalculation of trigonometrical functions for tesseroid centers*/ + double cos_a2, sin_a2, cos_b2, sin_b2; + + /*variables for precalculation of trigonometrical functions for grid points*/ + double lon_prev, lat_prev; + double cos_a2_prev, sin_a2_prev, cos_b2_prev, sin_b2_prev; + + FILE *logfile = NULL, *modelfile = NULL; + time_t rawtime; + clock_t tstart; + struct tm * timeinfo; + + double (*field1)(TESSEROID, double, double, double, GLQ, GLQ, GLQ); + double (*field2)(TESSEROID, double, double, double, GLQ, GLQ, GLQ); + double (*field3)(TESSEROID, double, double, double, GLQ, GLQ, GLQ); + void (*field_triple)(TESSEROID, double, double, double, GLQ, GLQ, GLQ, double*); + double ggt_1, ggt_2, ggt_3; + double gtt_v[3]; + int n_tesseroid; + + + log_init(LOG_INFO); + + rc = parse_tessb_args(argc, argv, progname, &args, &print_tessb_help); + if(rc == 2) + { + return 0; + } + if(rc == 1) + { + log_warning("Terminating due to bad input"); + log_warning("Try '%s -h' for instructions", progname); + + return 1; + } + + /* Set the appropriate logging level and log to file if necessary */ + if(!args.verbose) + { + log_init(LOG_WARNING); + } + if(args.logtofile) + { + logfile = fopen(args.logfname, "w"); + if(logfile == NULL) + { + log_error("unable to create log file %s", args.logfname); + log_warning("Terminating due to bad input"); + log_warning("Try '%s -h' for instructions", progname); + return 1; + } + log_tofile(logfile, LOG_DEBUG); + } + + /* Check if a custom distance-size ratio is given */ + if(args.ratio1 != 0) + { + ratio1 = args.ratio1; + } + if(args.ratio2 != 0) + { + ratio2 = args.ratio2; + } + if(args.ratio3 != 0) + { + ratio3 = args.ratio3; + } + + /* Print standard verbose */ + log_info("%s (Tesseroids project) %s", progname, tesseroids_version); + time(&rawtime); + timeinfo = localtime(&rawtime); + log_info("(local time) %s", asctime(timeinfo)); + log_info("Use recursive division of tesseroids: %s", + args.adaptative ? "True" : "False"); + log_info("Distance-size ratio1 for recusive division: %g", ratio1); + log_info("Distance-size ratio2 for recusive division: %g", ratio2); + log_info("Distance-size ratio3 for recusive division: %g", ratio3); + + /* Make the necessary GLQ structures */ + log_info("Using GLQ orders: %d lon / %d lat / %d r", args.lon_order, + args.lat_order, args.r_order); + glq_lon = glq_new(args.lon_order, -1, 1); + glq_lat = glq_new(args.lat_order, -1, 1); + glq_r = glq_new(args.r_order, -1, 1); + if(glq_lon == NULL || glq_lat == NULL || glq_r == NULL) + { + log_error("failed to create required GLQ structures"); + log_warning("Terminating due to bad input"); + log_warning("Try '%s -h' for instructions", progname); + if(args.logtofile) + fclose(logfile); + return 1; + } + + /* Read the tesseroid model file */ + log_info("Reading magnetic tesseroid model from file %s", args.modelfname); + modelfile = fopen(args.modelfname, "r"); + if(modelfile == NULL) + { + log_error("failed to open model file %s", args.modelfname); + log_warning("Terminating due to bad input"); + log_warning("Try '%s -h' for instructions", progname); + if(args.logtofile) + fclose(logfile); + return 1; + } + model = read_mag_tess_model(modelfile, &modelsize); + fclose(modelfile); + if(modelsize == 0) + { + log_error("tesseroid file %s is empty", args.modelfname); + log_warning("Terminating due to bad input"); + log_warning("Try '%s -h' for instructions", progname); + if(args.logtofile) + fclose(logfile); + return 1; + } + if(model == NULL) + { + log_error("failed to read model from file %s", args.modelfname); + log_warning("Terminating due to bad input"); + log_warning("Try '%s -h' for instructions", progname); + if(args.logtofile) + fclose(logfile); + return 1; + } + log_info("Total of %d tesseroid(s) read", modelsize); + + /* Print a header on the output with provenance information */ + if(strcmp(progname + 4, "pot") == 0) + { + printf("# Potential calculated with %s %s:\n", progname, + tesseroids_version); + } + else + { + printf("# %s component calculated with %s %s:\n", progname+4, progname, + tesseroids_version); + } + printf("# local time: %s", asctime(timeinfo)); + printf("# model file: %s (%d tesseroids)\n", args.modelfname, modelsize); + printf("# GLQ order: %d lon / %d lat / %d r\n", args.lon_order, + args.lat_order, args.r_order); + printf("# Use recursive division of tesseroids: %s\n", + args.adaptative ? "True" : "False"); + printf("# Distance-size ratio1 for recusive division: %g\n", ratio1); + printf("# Distance-size ratio2 for recusive division: %g\n", ratio2); + printf("# Distance-size ratio3 for recusive division: %g\n", ratio3); + + /////////////ELDAR BAYKIEV////////////// + /* Assign pointers to functions that calculate gravity gradient tensor components */ + if (!strcmp("tessbx", progname)) + { + field1 = &tess_gxx; + field2 = &tess_gxy; + field3 = &tess_gxz; + + field_triple = &tess_gxx_gxy_gxz; + } + + if (!strcmp("tessby", progname)) + { + field1 = &tess_gxy; + field2 = &tess_gyy; + field3 = &tess_gyz; + + field_triple = &tess_gxy_gyy_gyz; + } + + if (!strcmp("tessbz", progname)) + { + field1 = &tess_gxz; + field2 = &tess_gyz; + field3 = &tess_gzz; + + field_triple = &tess_gxz_gyz_gzz; + } + /////////////ELDAR BAYKIEV////////////// + + /* Read each computation point from stdin and calculate */ + log_info("Calculating (this may take a while)..."); + tstart = clock(); + + /////////////// + lon_prev = 0; + lat_prev = 0; + + for(line = 1; !feof(stdin); line++) + { + if(fgets(buff, 10000, stdin) == NULL) + { + if(ferror(stdin)) + { + log_error("problem encountered reading line %d", line); + error_exit = 1; + break; + } + } + else + { + /* Check for comments and blank lines */ + if(buff[0] == '#' || buff[0] == '\r' || buff[0] == '\n') + { + printf("%s", buff); + continue; + } + if(sscanf(buff, "%lf %lf %lf", &lon, &lat, &height) != 3) + { + log_warning("bad/invalid computation point at line %d", line); + log_warning("skipping this line and continuing"); + bad_input++; + continue; + } + /* Need to remove \n and \r from end of buff first to print the + result in the end */ + strstrip(buff); + + + + /////////////ELDAR BAYKIEV////////////// + res = 0; + if(args.adaptative) + { + for(n_tesseroid = 0; n_tesseroid < modelsize; n_tesseroid++) + { + gtt_v[0] = 0; + gtt_v[1] = 0; + gtt_v[2] = 0; + double B_to_H = model[n_tesseroid].suscept/(M_0);//IMPORTANT + double M_vect[3] = {model[n_tesseroid].Bx * B_to_H, model[n_tesseroid].By * B_to_H, model[n_tesseroid].Bz * B_to_H}; + double M_vect_p[3] = {0, 0, 0}; + + conv_vect_cblas(M_vect, (model[n_tesseroid].w + model[n_tesseroid].e)*0.5, (model[n_tesseroid].s + model[n_tesseroid].n)*0.5, lon, lat, M_vect_p); + + ggt_1 = calc_tess_model_adapt(&model[n_tesseroid], 1, lon, lat, height + MEAN_EARTH_RADIUS, glq_lon, glq_lat, glq_r, field1, ratio1); + ggt_2 = calc_tess_model_adapt(&model[n_tesseroid], 1, lon, lat, height + MEAN_EARTH_RADIUS, glq_lon, glq_lat, glq_r, field2, ratio2); + ggt_3 = calc_tess_model_adapt(&model[n_tesseroid], 1, lon, lat, height + MEAN_EARTH_RADIUS, glq_lon, glq_lat, glq_r, field3, ratio3); + + res = res + M_0*EOTVOS2SI*(ggt_1 * M_vect_p[0] + ggt_2 * M_vect_p[1] + ggt_3 * M_vect_p[2]) /(G*model[n_tesseroid].density*4*PI); + + //printf("res %g\n", res); + } + + ///////////////////////////////////////////////////////////////////////////////////////////////////////// + } + else + { + //precalculate trigonometrical functions + + if(lon == lon_prev) + { + cos_b2 = cos_b2_prev; + sin_b2 = sin_b2_prev; + } + else + { + cos_b2 = cos(DEG2RAD*lon); + sin_b2 = sin(DEG2RAD*lon); + } + + if(lat == lat_prev) + { + cos_a2 = cos_a2_prev; + sin_a2 = sin_a2_prev; + } + else + { + cos_a2 = cos(PI/2.0-DEG2RAD*lat); + sin_a2 = sin(PI/2.0-DEG2RAD*lat); + } + ///////////////////////////////////////////////////////////////////////////////////////////////////////// + for(n_tesseroid = 0; n_tesseroid < modelsize; n_tesseroid++) + { + gtt_v[0] = 0; + gtt_v[1] = 0; + gtt_v[2] = 0; + //ELDAR: TODO: PRECALCULATE SIC COSINE TABLES + double B_to_H = model[n_tesseroid].suscept/(M_0);//IMPORTANT + double M_vect[3] = {model[n_tesseroid].Bx * B_to_H, model[n_tesseroid].By * B_to_H, model[n_tesseroid].Bz * B_to_H}; + double M_vect_p[3] = {0, 0, 0}; + + //conv_vect_cblas(M_vect, (model[n_tesseroid].w + model[n_tesseroid].e)*0.5, (model[n_tesseroid].s + model[n_tesseroid].n)*0.5, lon, lat, M_vect_p); + conv_vect_cblas_precalc(M_vect, model[n_tesseroid].cos_a1, model[n_tesseroid].sin_a1, model[n_tesseroid].cos_b1, model[n_tesseroid].sin_b1, cos_a2, sin_a2, cos_b2, sin_b2, M_vect_p); + + calc_tess_model_triple(&model[n_tesseroid], 1, lon, lat, height + MEAN_EARTH_RADIUS, glq_lon, glq_lat, glq_r, field_triple, gtt_v); + + res = res + M_0*EOTVOS2SI*(gtt_v[0] * M_vect_p[0] + gtt_v[1] * M_vect_p[1] + gtt_v[2] * M_vect_p[2]) /(G*model[n_tesseroid].density*4*PI); + + } + + lon_prev = lon; + lat_prev = lat; + + cos_a2_prev = cos_a2; + sin_a2_prev = sin_a2; + + cos_b2_prev = cos_b2; + sin_b2_prev = sin_b2; + ///////////////////////////////////////////////////////////////////////////////////////////////////////// + } + printf("%s %.15g\n", buff, res); + points++; + } + } + if(bad_input) + { + log_warning("Encountered %d bad computation points which were skipped", + bad_input); + } + if(error_exit) + { + log_warning("Terminating due to error in input"); + log_warning("Try '%s -h' for instructions", progname); + } + else + { + log_info("Calculated on %d points in %.5g seconds", points, + (double)(clock() - tstart)/CLOCKS_PER_SEC); + } + /* Clean up */ + free(model); + glq_free(glq_lon); + glq_free(glq_lat); + glq_free(glq_r); + log_info("Done"); + if(args.logtofile) + fclose(logfile); + return 0; +} diff --git a/lib/tessb_main.h b/lib/tessb_main.h new file mode 100755 index 0000000..7a0ac83 --- /dev/null +++ b/lib/tessb_main.h @@ -0,0 +1,14 @@ +/* +Generic main function for the tessb* programs. +*/ + +#ifndef _TESSEROIDS_TESSH_MAIN_H_ +#define _TESSEROIDS_TESSH_MAIN_H_ + +#include "glq.h" +#include "geometry.h" + +void print_tessb_help(const char *progname); +int run_tessb_main(int argc, char **argv, const char *progname, double (*field)(TESSEROID, double, double, double, GLQ, GLQ, GLQ), double ratio1, double ratio2, double ratio3); + +#endif diff --git a/lib/version.cpp b/lib/version.cpp new file mode 100755 index 0000000..0f383b0 --- /dev/null +++ b/lib/version.cpp @@ -0,0 +1,17 @@ +/* +Hold the version number of the project. +*/ + +#include +#include "version.h" + + +/* Print version and license information */ +void print_version(const char *progname) +{ + printf("%s (Magnetic Tesseroids) %s\n\n", progname, tesseroids_version); + printf("Copyright (C) 2015, Eldar Baykiev.\n"); + printf("This program is based on the code of Leonardo Uieda.\n"); + printf("\n"); + printf("Developed by Eldar Baykiev.\n"); +} diff --git a/lib/version.h b/lib/version.h new file mode 100755 index 0000000..6ab6396 --- /dev/null +++ b/lib/version.h @@ -0,0 +1,18 @@ +/* +Hold the version number of the project. +*/ + +#ifndef _TESSEROIDS_VERSION_H_ +#define _TESSEROIDS_VERSION_H_ + +/** Current project version number */ +const char tesseroids_version[] = "1.1"; + + +/** Print version and license information + +@param progname name of the program +*/ +void print_version(const char *progname); + +#endif diff --git a/toolkits/CMakeLists.txt b/toolkits/CMakeLists.txt new file mode 100644 index 0000000..73b9fdd --- /dev/null +++ b/toolkits/CMakeLists.txt @@ -0,0 +1,47 @@ +# 设置编译选项 +if(${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Windows") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2") +endif() + +if(${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Linux") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -lopenblas -lm -O2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -lopenblas -lm -O2") +endif() + +if(${CMAKE_HOST_SYSTEM_NAME} STREQUAL "Darwin") + set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -framework Accelerate -O2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -framework Accelerate -O2") +endif() + +# 设置可执行文件的输出地址 +set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin/toolkits) + +# 下面设置tools的编译命令 首先设置一个宏 +macro(add_tools name) + # 添加可执行程序名称 + add_executable(${name} ${name}.cpp) + # 设置安装后的动态库调用地址 + set_target_properties(${name} PROPERTIES INSTALL_RPATH ${CMAKE_INSTALL_PREFIX}/lib) + # 链接动态库 + target_link_libraries(${name} PUBLIC magtess) + # 将可执行程序安装到/usr/local/sbin + install(TARGETS ${name} RUNTIME DESTINATION sbin/magtess) +endmacro() + +# 添加tools +add_tools(tessbx) +add_tools(tessby) +add_tools(tessbz) +add_tools(tessutil_combine_grids) +add_tools(tessutil_gradient_calculator) + +#add_tools(tessutil_magnetize_model) +# 添加可执行程序名称 +add_executable(tessutil_magnetize_model tessutil_magnetize_model.c) +# 设置安装后的动态库调用地址 +set_target_properties(tessutil_magnetize_model PROPERTIES INSTALL_RPATH ${CMAKE_INSTALL_PREFIX}/lib) +# 链接动态库 +target_link_libraries(tessutil_magnetize_model PUBLIC magtess) +# 将可执行程序安装到/usr/local/sbin +install(TARGETS tessutil_magnetize_model RUNTIME DESTINATION sbin/magtess) \ No newline at end of file diff --git a/toolkits/tessbx.cpp b/toolkits/tessbx.cpp new file mode 100755 index 0000000..8e6602a --- /dev/null +++ b/toolkits/tessbx.cpp @@ -0,0 +1,11 @@ +#include "../lib/constants.h" +#include "../lib/grav_tess.h" +#include "../lib/tessb_main.h" + + +/** Main tessbx*/ +int main(int argc, char **argv) +{ + return run_tessb_main(argc, argv, "tessbx", 0, TESSEROID_GXX_SIZE_RATIO, TESSEROID_GXY_SIZE_RATIO, TESSEROID_GXZ_SIZE_RATIO); + +} diff --git a/toolkits/tessby.cpp b/toolkits/tessby.cpp new file mode 100755 index 0000000..7d1115e --- /dev/null +++ b/toolkits/tessby.cpp @@ -0,0 +1,11 @@ +#include "../lib/constants.h" +#include "../lib/grav_tess.h" +#include "../lib/tessb_main.h" + + +/** Main tessby*/ +int main(int argc, char **argv) +{ + return run_tessb_main(argc, argv, "tessby", 0, TESSEROID_GXX_SIZE_RATIO, TESSEROID_GXY_SIZE_RATIO, TESSEROID_GXZ_SIZE_RATIO); + +} diff --git a/toolkits/tessbz.cpp b/toolkits/tessbz.cpp new file mode 100755 index 0000000..b398b26 --- /dev/null +++ b/toolkits/tessbz.cpp @@ -0,0 +1,11 @@ +#include "../lib/constants.h" +#include "../lib/grav_tess.h" +#include "../lib/tessb_main.h" + + +/** Main tessbz*/ +int main(int argc, char **argv) +{ + return run_tessb_main(argc, argv, "tessbz", 0, TESSEROID_GXX_SIZE_RATIO, TESSEROID_GXY_SIZE_RATIO, TESSEROID_GXZ_SIZE_RATIO); + +} diff --git a/toolkits/tessutil_combine_grids.cpp b/toolkits/tessutil_combine_grids.cpp new file mode 100755 index 0000000..53ebff9 --- /dev/null +++ b/toolkits/tessutil_combine_grids.cpp @@ -0,0 +1,89 @@ +#include +#include +#include +#include + +#define MAX_GRID_POINTS 16000 + +#define GRID_FORMAT "%lf %lf %f %lf" + +void printresult_withalt(double* longitudes, double* latitudes, float* altitudes, double* values, int n_values) +{ + + + for (int h = 0; h < n_values; h++) + printf( "%lf %lf %f %lf\n", longitudes[h], latitudes[h], altitudes[h], values[h]); + + return; +} + +void printresult(double* longitudes, double* latitudes, double* values, int n_values) +{ + + + for (int h = 0; h < n_values; h++) + printf( "%lf %lf %lf\n", longitudes[h], latitudes[h], values[h]); + + return; +} + + +int main(int argc, char**argv) +{ + int n_files = (argc-1)/2; + + double lons[MAX_GRID_POINTS]; + double lats[MAX_GRID_POINTS]; + float alts[MAX_GRID_POINTS]; + + double vals[MAX_GRID_POINTS]; + + double factor = 1; + + char * line = NULL; + size_t len = 0; + ssize_t read; + + int n_lines = 0; + + for (int i = 0; i < n_files; i++) + { + FILE * fp = fopen(argv[1+2*i], "r"); + sscanf(argv[1+2*i+1], "%lf", &factor); + + if (fp == NULL) + { + printf("ERROR: Can not open file with grid values.\n"); + exit(EXIT_FAILURE); + } + + n_lines = 0; + while ((read = getline(&line, &len, fp )) != -1) + { + + if ((line[0] != '#') && (strlen(line) > 2)) + { + n_lines++; + if (n_lines>MAX_GRID_POINTS) + { + printf("ERROR: Too many grid points (> %d) in the input. Recompile program with a bigger value of MAX_GRID_POINTS.\n", n_lines); + exit(EXIT_FAILURE); + } + double value; + + sscanf(line, GRID_FORMAT, &lons[n_lines-1], &lats[n_lines-1], &alts[n_lines-1], &value); + vals[n_lines-1] = vals[n_lines-1] + value*factor; + + } + } + fclose(fp); + } + + int no_alt = 0; + + if (no_alt) + printresult(lons, lats, vals, n_lines); + else + printresult_withalt(lons, lats, alts, vals, n_lines); + +} diff --git a/toolkits/tessutil_gradient_calculator.cpp b/toolkits/tessutil_gradient_calculator.cpp new file mode 100755 index 0000000..0f2b87b --- /dev/null +++ b/toolkits/tessutil_gradient_calculator.cpp @@ -0,0 +1,343 @@ +#include +#include +#include +#include + +#include "../lib/constants.h" +#include "../lib/parsers.h" +#include "../lib/linalg.h" + +#define MAX_GRID_POINTS 50000 + +#define GRID_FORMAT "%lf %lf %f %lf" + + +// TODO conversion of input/output units nT/km pT/km nT/m pT/m + + +void printcomp(double* longitudes, double* latitudes, double* values, int n_values) +{ + + + for (int h = 0; h < n_values; h++) + printf("%lf %lf %lf\n", longitudes[h], latitudes[h], values[h]); + + return; +} + +void printall(double* longitudes, double* latitudes, int n_values, double* values1, double* values2, double* values3, double* values4, double* values5, double* values6, double* values7) +{ + + + for (int h = 0; h < n_values; h++) + printf("%lf %lf %lf %lf %lf %lf %lf %lf %lf\n", longitudes[h], latitudes[h], values1[h], values2[h], values3[h], values4[h], values5[h], values6[h], values7[h]); + + return; +} + + +void print_gradcalc_help(const char *progname) +{ + printf("MAGNETIC TESSEROIDS: Gradient Calculator\n"); + printf("Usage: %s [OPTIONS]\n\n", progname); + printf("Options:\n"); + printf("\t-h\t\t\t Help\n\n"); + + printf("\t-bx[GRID FILENAME]\t Grid filename with Bx component\n"); + printf("\t-by[GRID FILENAME]\t Grid filename with By component\n"); + printf("\t-bz[GRID FILENAME]\t Grid filename with Bz component\n"); + printf("\tNOTE:\tall grids must be in format LON LAT ALT B*,\n\t\tstart from West-South corner and longitudes must increment first.\n\t\tLON, LAT in [deg], B* in [nT] and ALT in [m].\n\n"); + + printf("\t-c[COORD SYSTEM]\t Coordinate system in input grids. 1 - North-East-Down, 2 - North-East-Up\n"); + printf("\t-o[COMPONENT]\t\t If 0, then output format is LON LAT BXX BYX BZX BXY BYY BZY BZZ, if 1-7, then \n"); + printf("\tonly corresponding component would be printed with format LON LAT B**.\n"); + printf("\tNOTE: output is always in North-East-Down coordinate system.\n\t\tLON, LAT in [deg], B** in [nT/km].\n"); + + + + + //printf("TODO\n"); +} + + + + +int main(int argc, char**argv) +{ + const char * progname = "tessutil_gradient_calculator"; + GRADCALC_ARGS args; + int rc; + + + //log_init(LOG_INFO); + + rc = parse_gradcalc_args(argc, argv, progname, &args, &print_gradcalc_help); + if(rc == 2) //help + { + return 0; + } + if(rc == 1) + { + printf("Terminating due to bad input\n"); + printf("Try '%s -h' for instructions\n", progname); + + return 1; + } + + if ((args.gridbx_set == 0) || (args.gridbx_set == 0) || (args.gridbx_set == 0)) + { + printf("no input grids\n"); + exit(EXIT_FAILURE); + } + + if (args.bz_NEU_NED == 1) + printf("#Coordinate system in input grids: North-East-Down\n"); + else + printf("#Coordinate system in input grids: North-East-Up\n"); + printf("#Coordinate system in output grid: North-East-Down\n"); + + + double lons[MAX_GRID_POINTS]; + double lats[MAX_GRID_POINTS]; + float alts[MAX_GRID_POINTS]; + + double bx[MAX_GRID_POINTS]; + double by[MAX_GRID_POINTS]; + double bz[MAX_GRID_POINTS]; + //read file with bx + + char * line = NULL; + size_t len = 0; + ssize_t read; + + + + FILE * bxfp = fopen(args.gridbx_fn, "r"); + if (bxfp == NULL) + { + printf("ERROR: Can not open file with Bx values.\n"); + exit(EXIT_FAILURE); + } + + int n_lines = 0; + + + + while ((read = getline(&line, &len, bxfp )) != -1) + { + + if ((line[0] != '#') && (strlen(line) > 2)) + { + n_lines++; + if (n_lines>MAX_GRID_POINTS) + { + printf("ERROR: Too many grid points (> %d) in the input. Recompile program with a bigger value of MAX_GRID_POINTS.\n", n_lines); + exit(EXIT_FAILURE); + } + + sscanf(line, GRID_FORMAT, &lons[n_lines-1], &lats[n_lines-1], &alts[n_lines-1], &bx[n_lines-1]); + + } + } + fclose(bxfp); + + + /*number of grid points*/ + printf("#Number of grid points: %d\n", n_lines); + + /*grid spacing*/ + + double lon_min = lons[0]; + double lon_max = lons[n_lines-1]; + + double lon_step = lons[1]-lons[0]; + if ((lon_step == 0) || (lon_step <= 0) || (lon_max < lon_min)) + { + printf("ERROR: Wrong grid format. Longitudes must increment first. Use the format of tessgrd.\n"); + exit(EXIT_FAILURE); + } + + int lon_n = 1; + while ((lats[lon_n] == lats[0]) && (lon_n < n_lines)) + { + lon_n++; + } + double lat_step = lats[lon_n]-lats[0]; + + int lat_n = n_lines / lon_n; + + double lat_min = lats[0]; + double lat_max = lats[n_lines-1]; + + + + if ((lat_step <= 0) || (lat_step <= 0) || (lat_max < lat_min)) + { + printf("ERROR: Wrong grid format. Latitudes must increment. Use the format of tessgrd.\n"); + exit(EXIT_FAILURE); + } + + printf("#Longitudinal step: %lf, latitudinal step: %lf \n", lon_step, lat_step); + printf("#Longitudinal points: %d, latitudinal points: %d \n", lon_n, lat_n); + printf("#Edges: W %lf, E %lf, S %lf, N %lf \n", lon_min, lon_max, lat_min, lat_max); + + /* read other grids */ + // By + FILE * byfp = fopen(args.gridby_fn, "r"); + if (byfp == NULL) + { + printf("ERROR: Can not open file with Bx values.\n"); + exit(EXIT_FAILURE); + } + + int n_lines2 = 0; + while ((read = getline(&line, &len, byfp )) != -1) + { + + if ((line[0] != '#') && (strlen(line) > 2)) + { + n_lines2++; + //printf("%s", line); + double dummy1, dummy2; + float dummy3; + sscanf(line, GRID_FORMAT , &dummy1, &dummy2, &dummy3, &by[n_lines2-1]); + } + } + fclose(byfp); + + if (n_lines2 != n_lines) + { + printf("ERROR: Grid points of Bx and By do not coincide.\n"); + exit(EXIT_FAILURE); + } + + // Bz + FILE * bzfp = fopen(args.gridbz_fn, "r"); + if (bzfp == NULL) + { + printf("ERROR: Can not open file with Bx values.\n"); + exit(EXIT_FAILURE); + } + + n_lines2 = 0; + while ((read = getline(&line, &len, bzfp )) != -1) + { + if ((line[0] != '#') && (strlen(line) > 2)) + { + n_lines2++; + //printf("%s", line); + double dummy1, dummy2; + float dummy3; + double bz_curr; + sscanf(line, GRID_FORMAT, &dummy1, &dummy2, &dummy3,&bz_curr); + + bz[n_lines2-1] = args.bz_NEU_NED* bz_curr; //COORDINATE SYSTEM NEU or NED + } + } + fclose(byfp); + + if (n_lines2 != n_lines) + { + printf("ERROR: Grid points of Bx and Bz do not coincide.\n"); + exit(EXIT_FAILURE); + } + + + /* calculate gradients */ + + + double bxx[MAX_GRID_POINTS]; + double byx[MAX_GRID_POINTS]; + double bzx[MAX_GRID_POINTS]; + + double bxy[MAX_GRID_POINTS]; + double byy[MAX_GRID_POINTS]; + double bzy[MAX_GRID_POINTS]; + + double bzz[MAX_GRID_POINTS]; + + int cent_ind, west_ind, east_ind, south_ind, north_ind; + + for (int i = 1; i +#include +#include +#include +#include + +#include "../lib/constants.h" + +int my_isnan(double d) +{ + return (d != d); /* IEEE: only NaN is not equal to itself */ +} + +#define NaN log(-1.0) +#define FT2KM (1.0/0.0003048) +#define RAD2DEG (180.0/PI) + +#ifndef SEEK_SET +#define SEEK_SET 0 +#define SEEK_CUR 1 +#define SEEK_END 2 +#endif + +#define IEXT 0 +#define FALSE 0 +#define TRUE 1 /* constants */ +#define RECL 81 + +#define MAXINBUFF RECL+14 +/** Max size of in buffer **/ + +#define MAXREAD MAXINBUFF-2 +/** Max to read 2 less than total size (just to be safe) **/ + +#define MAXMOD 30 +/** Max number of models in a file **/ + +#define PATH MAXREAD +/** Max path and filename length **/ + +#define EXT_COEFF1 (double)0 +#define EXT_COEFF2 (double)0 +#define EXT_COEFF3 (double)0 + +#define MAXDEG 13 +#define MAXCOEFF (MAXDEG*(MAXDEG+2)+1) /* index starts with 1!, (from old Fortran?) */ +double gh1[MAXCOEFF]; +double gh2[MAXCOEFF]; +double gha[MAXCOEFF]; /* Geomag global variables */ +double ghb[MAXCOEFF]; +double d=0,f=0,h=0,i=0; +double dtemp,ftemp,htemp,itemp; +double x=0,y=0,z=0; +double xtemp,ytemp,ztemp; + +FILE *stream = NULL; /* Pointer to specified model data file */ + + +int main(int argc, char**argv) +{ + /* Variable declaration */ + + /* Variables related to tesseroids */ + FILE * tessfp; + FILE * tessoutfp; + char * line = NULL; + size_t len = 0; + ssize_t read; + + float W, E, S, N; + float HOT, HOB; + float DENSITY; + float SUSCEPT, OBX, OBY, OBZ; + float lon_c, lat_c, alt_c; + + int i, count; + + char tessfilename[PATH]; + char tessoutfilename[PATH]; + + + /* Control variables */ + int again = 1; + int decyears = 3; + int units = 4; + int decdeg = 3; + int range = -1; + int counter = 0; + int warn_H, warn_H_strong, warn_P; + + int modelI; /* Which model (Index) */ + int nmodel; /* Number of models in file */ + int max1[MAXMOD]; + int max2[MAXMOD]; + int max3[MAXMOD]; + int nmax; + int igdgc=3; + int isyear=-1; + int ismonth=-1; + int isday=-1; + int ieyear=-1; + int iemonth=-1; + int ieday=-1; + int ilat_deg=200; + int ilat_min=200; + int ilat_sec=200; + int ilon_deg=200; + int ilon_min=200; + int ilon_sec=200; + int fileline; + long irec_pos[MAXMOD]; + + int coords_from_file = 0; + int arg_err = 0; + int need_to_read_model = 1; + + char mdfile[PATH]; + char inbuff[MAXINBUFF]; + char model[MAXMOD][9]; + char *begin; + char *rest; + char args[7][MAXREAD]; + int iarg; + + char coord_fname[PATH]; + char out_fname[PATH]; + FILE *coordfile,*outfile; + int iline=0; + int read_flag; + + double epoch[MAXMOD]; + double yrmin[MAXMOD]; + double yrmax[MAXMOD]; + double minyr; + double maxyr; + double altmin[MAXMOD]; + double altmax[MAXMOD]; + double minalt; + double maxalt; + double alt=-999999; + double sdate=-1; + double step=-1; + double syr; + double edate=-1; + double latitude=200; + double longitude=200; + double ddot; + double fdot; + double hdot; + double idot; + double xdot; + double ydot; + double zdot; + double warn_H_val, warn_H_strong_val; + + /* Subroutines used */ + + + double degrees_to_decimal(); + double julday(); + int interpsh(); + int extrapsh(); + int shval3(); + int dihf(); + int safegets(char *buffer,int n); + int getshc(); + + + /* Initializations. */ + + inbuff[MAXREAD+1]='\0'; /* Just to protect mem. */ + inbuff[MAXINBUFF-1]='\0'; /* Just to protect mem. */ + + for (iarg=0; iarg MAXMOD) /* If too many headers */ + { + printf("ERROR: Too many models in file %s on line %d.", mdfile, fileline); + fclose(stream); + exit(6); + } + + irec_pos[modelI]=ftell(stream); + /* Get fields from buffer into individual vars. */ + sscanf(inbuff, "%s%lg%d%d%d%lg%lg%lg%lg", model[modelI], &epoch[modelI], &max1[modelI], &max2[modelI], &max3[modelI], &yrmin[modelI], &yrmax[modelI], &altmin[modelI], &altmax[modelI]); + + /* Compute date range for all models */ + if (modelI == 0) /*If first model */ + { + minyr=yrmin[0]; + maxyr=yrmax[0]; + } + else + { + if (yrmin[modelI]maxyr){ + maxyr=yrmax[modelI]; + } + } /* if modelI != 0 */ + + } /* If 1st 3 chars are spaces */ + + } /* While not end of model file */ + + nmodel = modelI + 1; + fclose(stream); + + decdeg=1; + decyears=2; + range=1; + igdgc=1; + units = 1; + + + sdate = julday(ismonth,isday,isyear); + + + /* Pick model */ + for (modelI=0; modelI 2 ? leap_year : 0)); + + return ((double)year + (day_in_year / (365.0 + leap_year))); +} + + +/****************************************************************************/ +/* */ +/* Subroutine getshc */ +/* */ +/****************************************************************************/ +/* */ +/* Reads spherical harmonic coefficients from the specified */ +/* model into an array. */ +/* */ +/* Input: */ +/* stream - Logical unit number */ +/* iflag - Flag for SV equal to ) or not equal to 0 */ +/* for designated read statements */ +/* strec - Starting record number to read from model */ +/* nmax_of_gh - Maximum degree and order of model */ +/* */ +/* Output: */ +/* gh1 or 2 - Schmidt quasi-normal internal spherical */ +/* harmonic coefficients */ +/* */ +/* FORTRAN */ +/* Bill Flanagan */ +/* NOAA CORPS, DESDIS, NGDC, 325 Broadway, Boulder CO. 80301 */ +/* */ +/* C */ +/* C. H. Shaffer */ +/* Lockheed Missiles and Space Company, Sunnyvale CA */ +/* August 15, 1988 */ +/* */ +/****************************************************************************/ + + +int getshc(file, iflag, strec, nmax_of_gh, gh) +char file[PATH]; +int iflag; +long int strec; +int nmax_of_gh; +int gh; +{ + char inbuff[MAXINBUFF]; + char irat[9]; + int ii,m,n,mm,nn; + int ios; + int line_num; + double g,hh; + double trash; + + stream = fopen(file, "rt"); + if (stream == NULL) + { + printf("\nError on opening file %s", file); + } + else + { + ii = 0; + ios = 0; + fseek(stream,strec,SEEK_SET); + for ( nn = 1; nn <= nmax_of_gh; ++nn) + { + for (mm = 0; mm <= nn; ++mm) + { + if (iflag == 1) + { + fgets(inbuff, MAXREAD, stream); + sscanf(inbuff, "%d%d%lg%lg%lg%lg%s%d", + &n, &m, &g, &hh, &trash, &trash, irat, &line_num); + } + else + { + fgets(inbuff, MAXREAD, stream); + sscanf(inbuff, "%d%d%lg%lg%lg%lg%s%d", + &n, &m, &trash, &trash, &g, &hh, irat, &line_num); + } + if ((nn != n) || (mm != m)) + { + ios = -2; + fclose(stream); + return(ios); + } + ii = ii + 1; + switch(gh) + { + case 1: gh1[ii] = g; + break; + case 2: gh2[ii] = g; + break; + default: printf("\nError in subroutine getshc"); + break; + } + if (m != 0) + { + ii = ii+ 1; + switch(gh) + { + case 1: gh1[ii] = hh; + break; + case 2: gh2[ii] = hh; + break; + default: printf("\nError in subroutine getshc"); + break; + } + } + } + } + } + fclose(stream); + return(ios); +} + + +/****************************************************************************/ +/* */ +/* Subroutine extrapsh */ +/* */ +/****************************************************************************/ +/* */ +/* Extrapolates linearly a spherical harmonic model with a */ +/* rate-of-change model. */ +/* */ +/* Input: */ +/* date - date of resulting model (in decimal year) */ +/* dte1 - date of base model */ +/* nmax1 - maximum degree and order of base model */ +/* gh1 - Schmidt quasi-normal internal spherical */ +/* harmonic coefficients of base model */ +/* nmax2 - maximum degree and order of rate-of-change model */ +/* gh2 - Schmidt quasi-normal internal spherical */ +/* harmonic coefficients of rate-of-change model */ +/* */ +/* Output: */ +/* gha or b - Schmidt quasi-normal internal spherical */ +/* harmonic coefficients */ +/* nmax - maximum degree and order of resulting model */ +/* */ +/* FORTRAN */ +/* A. Zunde */ +/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */ +/* */ +/* C */ +/* C. H. Shaffer */ +/* Lockheed Missiles and Space Company, Sunnyvale CA */ +/* August 16, 1988 */ +/* */ +/****************************************************************************/ + + +int extrapsh(date, dte1, nmax1, nmax2, gh) +double date; +double dte1; +int nmax1; +int nmax2; +int gh; +{ + int nmax; + int k, l; + int ii; + double factor; + + factor = date - dte1; + if (nmax1 == nmax2) + { + k = nmax1 * (nmax1 + 2); + nmax = nmax1; + } + else + { + if (nmax1 > nmax2) + { + k = nmax2 * (nmax2 + 2); + l = nmax1 * (nmax1 + 2); + switch(gh) + { + case 3: for ( ii = k + 1; ii <= l; ++ii) + { + gha[ii] = gh1[ii]; + } + break; + case 4: for ( ii = k + 1; ii <= l; ++ii) + { + ghb[ii] = gh1[ii]; + } + break; + default: printf("\nError in subroutine extrapsh"); + break; + } + nmax = nmax1; + } + else + { + k = nmax1 * (nmax1 + 2); + l = nmax2 * (nmax2 + 2); + switch(gh) + { + case 3: for ( ii = k + 1; ii <= l; ++ii) + { + gha[ii] = factor * gh2[ii]; + } + break; + case 4: for ( ii = k + 1; ii <= l; ++ii) + { + ghb[ii] = factor * gh2[ii]; + } + break; + default: printf("\nError in subroutine extrapsh"); + break; + } + nmax = nmax2; + } + } + switch(gh) + { + case 3: for ( ii = 1; ii <= k; ++ii) + { + gha[ii] = gh1[ii] + factor * gh2[ii]; + } + break; + case 4: for ( ii = 1; ii <= k; ++ii) + { + ghb[ii] = gh1[ii] + factor * gh2[ii]; + } + break; + default: printf("\nError in subroutine extrapsh"); + break; + } + return(nmax); +} + +/****************************************************************************/ +/* */ +/* Subroutine interpsh */ +/* */ +/****************************************************************************/ +/* */ +/* Interpolates linearly, in time, between two spherical harmonic */ +/* models. */ +/* */ +/* Input: */ +/* date - date of resulting model (in decimal year) */ +/* dte1 - date of earlier model */ +/* nmax1 - maximum degree and order of earlier model */ +/* gh1 - Schmidt quasi-normal internal spherical */ +/* harmonic coefficients of earlier model */ +/* dte2 - date of later model */ +/* nmax2 - maximum degree and order of later model */ +/* gh2 - Schmidt quasi-normal internal spherical */ +/* harmonic coefficients of internal model */ +/* */ +/* Output: */ +/* gha or b - coefficients of resulting model */ +/* nmax - maximum degree and order of resulting model */ +/* */ +/* FORTRAN */ +/* A. Zunde */ +/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */ +/* */ +/* C */ +/* C. H. Shaffer */ +/* Lockheed Missiles and Space Company, Sunnyvale CA */ +/* August 17, 1988 */ +/* */ +/****************************************************************************/ + + +int interpsh(date, dte1, nmax1, dte2, nmax2, gh) + double date; + double dte1; + int nmax1; + double dte2; + int nmax2; + int gh; +{ + int nmax; + int k, l; + int ii; + double factor; + + factor = (date - dte1) / (dte2 - dte1); + if (nmax1 == nmax2) + { + k = nmax1 * (nmax1 + 2); + nmax = nmax1; + } + else + { + if (nmax1 > nmax2) + { + k = nmax2 * (nmax2 + 2); + l = nmax1 * (nmax1 + 2); + switch(gh) + { + case 3: for ( ii = k + 1; ii <= l; ++ii) + { + gha[ii] = gh1[ii] + factor * (-gh1[ii]); + } + break; + case 4: for ( ii = k + 1; ii <= l; ++ii) + { + ghb[ii] = gh1[ii] + factor * (-gh1[ii]); + } + break; + default: printf("\nError in subroutine extrapsh"); + break; + } + nmax = nmax1; + } + else + { + k = nmax1 * (nmax1 + 2); + l = nmax2 * (nmax2 + 2); + switch(gh) + { + case 3: for ( ii = k + 1; ii <= l; ++ii) + { + gha[ii] = factor * gh2[ii]; + } + break; + case 4: for ( ii = k + 1; ii <= l; ++ii) + { + ghb[ii] = factor * gh2[ii]; + } + break; + default: printf("\nError in subroutine extrapsh"); + break; + } + nmax = nmax2; + } + } + switch(gh) + { + case 3: for ( ii = 1; ii <= k; ++ii) + { + gha[ii] = gh1[ii] + factor * (gh2[ii] - gh1[ii]); + } + break; + case 4: for ( ii = 1; ii <= k; ++ii) + { + ghb[ii] = gh1[ii] + factor * (gh2[ii] - gh1[ii]); + } + break; + default: printf("\nError in subroutine extrapsh"); + break; + } + return(nmax); +} + + + + + +/****************************************************************************/ +/* */ +/* Subroutine shval3 */ +/* */ +/****************************************************************************/ +/* */ +/* Calculates field components from spherical harmonic (sh) */ +/* models. */ +/* */ +/* Input: */ +/* igdgc - indicates coordinate system used; set equal */ +/* to 1 if geodetic, 2 if geocentric */ +/* latitude - north latitude, in degrees */ +/* longitude - east longitude, in degrees */ +/* elev - WGS84 altitude above ellipsoid (igdgc=1), or */ +/* radial distance from earth's center (igdgc=2) */ +/* a2,b2 - squares of semi-major and semi-minor axes of */ +/* the reference spheroid used for transforming */ +/* between geodetic and geocentric coordinates */ +/* or components */ +/* nmax - maximum degree and order of coefficients */ +/* iext - external coefficients flag (=0 if none) */ +/* ext1,2,3 - the three 1st-degree external coefficients */ +/* (not used if iext = 0) */ +/* */ +/* Output: */ +/* x - northward component */ +/* y - eastward component */ +/* z - vertically-downward component */ +/* */ +/* based on subroutine 'igrf' by D. R. Barraclough and S. R. C. Malin, */ +/* report no. 71/1, institute of geological sciences, U.K. */ +/* */ +/* FORTRAN */ +/* Norman W. Peddie */ +/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */ +/* */ +/* C */ +/* C. H. Shaffer */ +/* Lockheed Missiles and Space Company, Sunnyvale CA */ +/* August 17, 1988 */ +/* */ +/****************************************************************************/ + + +int shval3(igdgc, flat, flon, elev, nmax, gh, iext, ext1, ext2, ext3) + int igdgc; + double flat; + double flon; + double elev; + int nmax; + int gh; + int iext; + double ext1; + double ext2; + double ext3; +{ + double earths_radius = 6371.2; + double dtr = 0.01745329; + double slat; + double clat; + double ratio; + double aa, bb, cc, dd; + double sd; + double cd; + double r; + double a2; + double b2; + double rr; + double fm,fn; + double sl[14]; + double cl[14]; + double p[119]; + double q[119]; + int ii,j,k,l,m,n; + int npq; + int ios; + double argument; + double power; + a2 = 40680631.59; /* WGS84 */ + b2 = 40408299.98; /* WGS84 */ + ios = 0; + r = elev; + argument = flat * dtr; + slat = sin( argument ); + if ((90.0 - flat) < 0.001) + { + aa = 89.999; /* 300 ft. from North pole */ + } + else + { + if ((90.0 + flat) < 0.001) + { + aa = -89.999; /* 300 ft. from South pole */ + } + else + { + aa = flat; + } + } + argument = aa * dtr; + clat = cos( argument ); + argument = flon * dtr; + sl[1] = sin( argument ); + cl[1] = cos( argument ); + switch(gh) + { + case 3: x = 0; + y = 0; + z = 0; + break; + case 4: xtemp = 0; + ytemp = 0; + ztemp = 0; + break; + default: printf("\nError in subroutine shval3"); + break; + } + sd = 0.0; + cd = 1.0; + l = 1; + n = 0; + m = 1; + npq = (nmax * (nmax + 3)) / 2; + if (igdgc == 1) + { + aa = a2 * clat * clat; + bb = b2 * slat * slat; + cc = aa + bb; + argument = cc; + dd = sqrt( argument ); + argument = elev * (elev + 2.0 * dd) + (a2 * aa + b2 * bb) / cc; + r = sqrt( argument ); + cd = (elev + dd) / r; + sd = (a2 - b2) / dd * slat * clat / r; + aa = slat; + slat = slat * cd - clat * sd; + clat = clat * cd + aa * sd; + } + ratio = earths_radius / r; + argument = 3.0; + aa = sqrt( argument ); + p[1] = 2.0 * slat; + p[2] = 2.0 * clat; + p[3] = 4.5 * slat * slat - 1.5; + p[4] = 3.0 * aa * clat * slat; + q[1] = -clat; + q[2] = slat; + q[3] = -3.0 * clat * slat; + q[4] = aa * (slat * slat - clat * clat); + for ( k = 1; k <= npq; ++k) + { + if (n < m) + { + m = 0; + n = n + 1; + argument = ratio; + power = n + 2; + rr = pow(argument,power); + fn = n; + } + fm = m; + if (k >= 5) + { + if (m == n) + { + argument = (1.0 - 0.5/fm); + aa = sqrt( argument ); + j = k - n - 1; + p[k] = (1.0 + 1.0/fm) * aa * clat * p[j]; + q[k] = aa * (clat * q[j] + slat/fm * p[j]); + sl[m] = sl[m-1] * cl[1] + cl[m-1] * sl[1]; + cl[m] = cl[m-1] * cl[1] - sl[m-1] * sl[1]; + } + else + { + argument = fn*fn - fm*fm; + aa = sqrt( argument ); + argument = ((fn - 1.0)*(fn-1.0)) - (fm * fm); + bb = sqrt( argument )/aa; + cc = (2.0 * fn - 1.0)/aa; + ii = k - n; + j = k - 2 * n + 1; + p[k] = (fn + 1.0) * (cc * slat/fn * p[ii] - bb/(fn - 1.0) * p[j]); + q[k] = cc * (slat * q[ii] - clat/fn * p[ii]) - bb * q[j]; + } + } + switch(gh) + { + case 3: aa = rr * gha[l]; + break; + case 4: aa = rr * ghb[l]; + break; + default: printf("\nError in subroutine shval3"); + break; + } + if (m == 0) + { + switch(gh) + { + case 3: x = x + aa * q[k]; + z = z - aa * p[k]; + break; + case 4: xtemp = xtemp + aa * q[k]; + ztemp = ztemp - aa * p[k]; + break; + default: printf("\nError in subroutine shval3"); + break; + } + l = l + 1; + } + else + { + switch(gh) + { + case 3: bb = rr * gha[l+1]; + cc = aa * cl[m] + bb * sl[m]; + x = x + cc * q[k]; + z = z - cc * p[k]; + if (clat > 0) + { + y = y + (aa * sl[m] - bb * cl[m]) * + fm * p[k]/((fn + 1.0) * clat); + } + else + { + y = y + (aa * sl[m] - bb * cl[m]) * q[k] * slat; + } + l = l + 2; + break; + case 4: bb = rr * ghb[l+1]; + cc = aa * cl[m] + bb * sl[m]; + xtemp = xtemp + cc * q[k]; + ztemp = ztemp - cc * p[k]; + if (clat > 0) + { + ytemp = ytemp + (aa * sl[m] - bb * cl[m]) * + fm * p[k]/((fn + 1.0) * clat); + } + else + { + ytemp = ytemp + (aa * sl[m] - bb * cl[m]) * + q[k] * slat; + } + l = l + 2; + break; + default: printf("\nError in subroutine shval3"); + break; + } + } + m = m + 1; + } + if (iext != 0) + { + aa = ext2 * cl[1] + ext3 * sl[1]; + switch(gh) + { + case 3: x = x - ext1 * clat + aa * slat; + y = y + ext2 * sl[1] - ext3 * cl[1]; + z = z + ext1 * slat + aa * clat; + break; + case 4: xtemp = xtemp - ext1 * clat + aa * slat; + ytemp = ytemp + ext2 * sl[1] - ext3 * cl[1]; + ztemp = ztemp + ext1 * slat + aa * clat; + break; + default: printf("\nError in subroutine shval3"); + break; + } + } + switch(gh) + { + case 3: aa = x; + x = x * cd + z * sd; + z = z * cd - aa * sd; + break; + case 4: aa = xtemp; + xtemp = xtemp * cd + ztemp * sd; + ztemp = ztemp * cd - aa * sd; + break; + default: printf("\nError in subroutine shval3"); + break; + } + return(ios); +} + + +/****************************************************************************/ +/* */ +/* Subroutine dihf */ +/* */ +/****************************************************************************/ +/* */ +/* Computes the geomagnetic d, i, h, and f from x, y, and z. */ +/* */ +/* Input: */ +/* x - northward component */ +/* y - eastward component */ +/* z - vertically-downward component */ +/* */ +/* Output: */ +/* d - declination */ +/* i - inclination */ +/* h - horizontal intensity */ +/* f - total intensity */ +/* */ +/* FORTRAN */ +/* A. Zunde */ +/* USGS, MS 964, box 25046 Federal Center, Denver, CO. 80225 */ +/* */ +/* C */ +/* C. H. Shaffer */ +/* Lockheed Missiles and Space Company, Sunnyvale CA */ +/* August 22, 1988 */ +/* */ +/****************************************************************************/ + +int dihf (gh) + int gh; +{ + int ios; + int j; + double sn; + double h2; + double hpx; + double argument, argument2; + + ios = gh; + sn = 0.0001; + + switch(gh) + { + case 3: for (j = 1; j <= 1; ++j) + { + h2 = x*x + y*y; + argument = h2; + h = sqrt(argument); /* calculate horizontal intensity */ + argument = h2 + z*z; + f = sqrt(argument); /* calculate total intensity */ + if (f < sn) + { + d = NaN; /* If d and i cannot be determined, */ + i = NaN; /* set equal to NaN */ + } + else + { + argument = z; + argument2 = h; + i = atan2(argument,argument2); + if (h < sn) + { + d = NaN; + } + else + { + hpx = h + x; + if (hpx < sn) + { + d = PI; + } + else + { + argument = y; + argument2 = hpx; + d = 2.0 * atan2(argument,argument2); + } + } + } + } + break; + case 4: for (j = 1; j <= 1; ++j) + { + h2 = xtemp*xtemp + ytemp*ytemp; + argument = h2; + htemp = sqrt(argument); + argument = h2 + ztemp*ztemp; + ftemp = sqrt(argument); + if (ftemp < sn) + { + dtemp = NaN; /* If d and i cannot be determined, */ + itemp = NaN; /* set equal to 999.0 */ + } + else + { + argument = ztemp; + argument2 = htemp; + itemp = atan2(argument,argument2); + if (htemp < sn) + { + dtemp = NaN; + } + else + { + hpx = htemp + xtemp; + if (hpx < sn) + { + dtemp = PI; + } + else + { + argument = ytemp; + argument2 = hpx; + dtemp = 2.0 * atan2(argument,argument2); + } + } + } + } + break; + default: printf("\nError in subroutine dihf"); + break; + } + return(ios); +}