00001 /** 00002 * \file Geoid.cpp 00003 * \brief Implementation for GeographicLib::Geoid class 00004 * 00005 * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> 00006 * and licensed under the LGPL. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #include "GeographicLib/Geoid.hpp" 00011 #include <sstream> 00012 #include <cstdlib> 00013 00014 #define GEOGRAPHICLIB_GEOID_CPP "$Id: Geoid.cpp 6785 2010-01-05 22:15:42Z karney $" 00015 00016 RCSID_DECL(GEOGRAPHICLIB_GEOID_CPP) 00017 RCSID_DECL(GEOGRAPHICLIB_GEOID_HPP) 00018 00019 #if !defined(GEOID_DEFAULT_PATH) 00020 #if defined(_MSC_VER) 00021 #define GEOID_DEFAULT_PATH "C:/cygwin/usr/local/share/GeographicLib/geoids" 00022 #else 00023 #define GEOID_DEFAULT_PATH "/usr/local/share/GeographicLib/geoids" 00024 #endif 00025 #endif 00026 00027 #if defined(_MSC_VER) 00028 // Squelch warnings about unsafe use of getenv 00029 #pragma warning (disable: 4996) 00030 #endif 00031 00032 namespace GeographicLib { 00033 00034 using namespace std; 00035 00036 // This is the transfer matrix for a 3rd order fit with a 12-point stencil 00037 // with weights 00038 // 00039 // \ x -1 0 1 2 00040 // y 00041 // -1 . 1 1 . 00042 // 0 1 2 2 1 00043 // 1 1 2 2 1 00044 // 2 . 1 1 . 00045 // 00046 // A algorithm for n-dimensional polynomial fits is described in 00047 // F. H. Lesh, 00048 // Multi-dimensional least-squares polynomial curve fitting, 00049 // CACM 2, 29-30 (1959). 00050 // 00051 // Here's the Maxima code to generate this matrix: 00052 // 00053 // /* The stencil and the weights */ 00054 // xarr:[ 00055 // 0, 1, 00056 // -1, 0, 1, 2, 00057 // -1, 0, 1, 2, 00058 // 0, 1]$ 00059 // yarr:[ 00060 // -1,-1, 00061 // 0, 0, 0, 0, 00062 // 1, 1, 1, 1, 00063 // 2, 2]$ 00064 // warr:[ 00065 // 1, 1, 00066 // 1, 2, 2, 1, 00067 // 1, 2, 2, 1, 00068 // 1, 1]$ 00069 // 00070 // /* [x exponent, y exponent] for cubic fit */ 00071 // pows:[ 00072 // [0,0], 00073 // [1,0],[0,1], 00074 // [2,0],[1,1],[0,2], 00075 // [3,0],[2,1],[1,2],[0,3]]$ 00076 // 00077 // basisvec(x,y,pows):=map(lambda([ex],(if ex[1]=0 then 1 else x^ex[1])* 00078 // (if ex[2]=0 then 1 else y^ex[2])),pows)$ 00079 // addterm(x,y,f,w,pows):=block([a,b,bb:basisvec(x,y,pows)], 00080 // a:w*(transpose(bb).bb), 00081 // b:(w*f) * bb, 00082 // [a,b])$ 00083 // 00084 // c3row(k):=block([a,b,c,pows:pows,n], 00085 // n:length(pows), 00086 // a:zeromatrix(n,n), 00087 // b:copylist(part(a,1)), 00088 // c:[a,b], 00089 // for i:1 thru length(xarr) do 00090 // c:c+addterm(xarr[i],yarr[i],if i=k then 1 else 0,warr[i],pows), 00091 // a:c[1],b:c[2], 00092 // part(transpose( a^^-1 . transpose(b)),1))$ 00093 // c3:[]$ 00094 // for k:1 thru length(warr) do c3:endcons(c3row(k),c3)$ 00095 // c3:apply(matrix,c3)$ 00096 // c0:part(ratsimp( 00097 // genmatrix(yc,1,length(warr)).abs(c3).genmatrix(yd,length(pows),1)),2)$ 00098 // c3:c0*c3$ 00099 00100 const Math::real Geoid::c0 = 240; // Common denominator 00101 const Math::real Geoid::c3[stencilsize * nterms] = { 00102 9, -18, -88, 0, 96, 90, 0, 0, -60, -20, 00103 -9, 18, 8, 0, -96, 30, 0, 0, 60, -20, 00104 9, -88, -18, 90, 96, 0, -20, -60, 0, 0, 00105 186, -42, -42, -150, -96, -150, 60, 60, 60, 60, 00106 54, 162, -78, 30, -24, -90, -60, 60, -60, 60, 00107 -9, -32, 18, 30, 24, 0, 20, -60, 0, 0, 00108 -9, 8, 18, 30, -96, 0, -20, 60, 0, 0, 00109 54, -78, 162, -90, -24, 30, 60, -60, 60, -60, 00110 -54, 78, 78, 90, 144, 90, -60, -60, -60, -60, 00111 9, -8, -18, -30, -24, 0, 20, 60, 0, 0, 00112 -9, 18, -32, 0, 24, 30, 0, 0, -60, 20, 00113 9, -18, -8, 0, -24, -30, 0, 0, 60, 20, 00114 }; 00115 00116 // Like c3, but with the coeffs of x, x^2, and x^3 constrained to be zero. 00117 // Use this at the N pole so that the height in independent of the longitude 00118 // there. 00119 // 00120 // Here's the Maxima code to generate this matrix (continued from above). 00121 // 00122 // /* figure which terms to exclude so that fit is indep of x at y=0 */ 00123 // mask:part(zeromatrix(1,length(pows)),1)+1$ 00124 // for i:1 thru length(pows) do 00125 // if pows[i][1]>0 and pows[i][2]=0 then mask[i]:0$ 00126 // 00127 // /* Same as c3row but with masked pows. */ 00128 // c3nrow(k):=block([a,b,c,powsa:[],n,d,e], 00129 // for i:1 thru length(mask) do if mask[i]>0 then 00130 // powsa:endcons(pows[i],powsa), 00131 // n:length(powsa), 00132 // a:zeromatrix(n,n), 00133 // b:copylist(part(a,1)), 00134 // c:[a,b], 00135 // for i:1 thru length(xarr) do 00136 // c:c+addterm(xarr[i],yarr[i],if i=k then 1 else 0,warr[i],powsa), 00137 // a:c[1],b:c[2], 00138 // d:part(transpose( a^^-1 . transpose(b)),1), 00139 // e:[], 00140 // for i:1 thru length(mask) do 00141 // if mask[i]>0 then (e:endcons(first(d),e),d:rest(d)) else e:endcons(0,e), 00142 // e)$ 00143 // c3n:[]$ 00144 // for k:1 thru length(warr) do c3n:endcons(c3nrow(k),c3n)$ 00145 // c3n:apply(matrix,c3n)$ 00146 // c0n:part(ratsimp( 00147 // genmatrix(yc,1,length(warr)).abs(c3n).genmatrix(yd,length(pows),1)),2)$ 00148 // c3n:c0n*c3n$ 00149 00150 const Math::real Geoid::c0n = 372; // Common denominator 00151 const Math::real Geoid::c3n[stencilsize * nterms] = { 00152 0, 0, -131, 0, 138, 144, 0, 0, -102, -31, 00153 0, 0, 7, 0, -138, 42, 0, 0, 102, -31, 00154 62, 0, -31, 0, 0, -62, 0, 0, 0, 31, 00155 124, 0, -62, 0, 0, -124, 0, 0, 0, 62, 00156 124, 0, -62, 0, 0, -124, 0, 0, 0, 62, 00157 62, 0, -31, 0, 0, -62, 0, 0, 0, 31, 00158 0, 0, 45, 0, -183, -9, 0, 93, 18, 0, 00159 0, 0, 216, 0, 33, 87, 0, -93, 12, -93, 00160 0, 0, 156, 0, 153, 99, 0, -93, -12, -93, 00161 0, 0, -45, 0, -3, 9, 0, 93, -18, 0, 00162 0, 0, -55, 0, 48, 42, 0, 0, -84, 31, 00163 0, 0, -7, 0, -48, -42, 0, 0, 84, 31, 00164 }; 00165 00166 // Like c3n, but y -> 1-y so that h is independent of x at y = 1. Use this 00167 // at the S pole so that the height in independent of the longitude there. 00168 // 00169 // Here's the Maxima code to generate this matrix (continued from above). 00170 // 00171 // /* Transform c3n to c3s by transforming y -> 1-y */ 00172 // vv:[ 00173 // v[11],v[12], 00174 // v[7],v[8],v[9],v[10], 00175 // v[3],v[4],v[5],v[6], 00176 // v[1],v[2]]$ 00177 // poly:expand(vv.(c3n/c0n).transpose(basisvec(x,1-y,pows)))$ 00178 // c3sf[i,j]:=coeff(coeff(coeff(poly,v[i]),x,pows[j][1]),y,pows[j][2])$ 00179 // c3s:genmatrix(c3sf,length(vv),length(pows))$ 00180 // c0s:part(ratsimp( 00181 // genmatrix(yc,1,length(warr)).abs(c3s).genmatrix(yd,length(pows),1)),2)$ 00182 // c3s:c0s*c3s$ 00183 00184 const Math::real Geoid::c0s = 372; // Common denominator 00185 const Math::real Geoid::c3s[stencilsize * nterms] = { 00186 18, -36, -122, 0, 120, 135, 0, 0, -84, -31, 00187 -18, 36, -2, 0, -120, 51, 0, 0, 84, -31, 00188 36, -165, -27, 93, 147, -9, 0, -93, 18, 0, 00189 210, 45, -111, -93, -57, -192, 0, 93, 12, 93,