/usr/share/doc/geographiclib/html/GeodesicLine_8cpp_source.html is in geographiclib-tools 1.8-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 | <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<title>GeographicLib: GeodesicLine.cpp Source File</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<link href="doxygen.css" rel="stylesheet" type="text/css"/>
</head>
<body>
<!-- Generated by Doxygen 1.7.1 -->
<div class="navigation" id="top">
<div class="tabs">
<ul class="tablist">
<li><a href="index.html"><span>Main Page</span></a></li>
<li><a href="pages.html"><span>Related Pages</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="annotated.html"><span>Classes</span></a></li>
<li class="current"><a href="files.html"><span>Files</span></a></li>
<li><a href="dirs.html"><span>Directories</span></a></li>
</ul>
</div>
<div class="tabs2">
<ul class="tablist">
<li><a href="files.html"><span>File List</span></a></li>
<li><a href="globals.html"><span>File Members</span></a></li>
</ul>
</div>
<div class="navpath">
<ul>
<li><a class="el" href="dir_b6dff8fdfd5d70bd61ba562840823df9.html">src</a> </li>
</ul>
</div>
</div>
<div class="header">
<div class="headertitle">
<h1>GeodesicLine.cpp</h1> </div>
</div>
<div class="contents">
<a href="GeodesicLine_8cpp.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">/**</span>
<a name="l00002"></a>00002 <span class="comment"> * \file GeodesicLine.cpp</span>
<a name="l00003"></a>00003 <span class="comment"> * \brief Implementation for GeographicLib::GeodesicLine class</span>
<a name="l00004"></a>00004 <span class="comment"> *</span>
<a name="l00005"></a>00005 <span class="comment"> * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com></span>
<a name="l00006"></a>00006 <span class="comment"> * and licensed under the LGPL. For more information, see</span>
<a name="l00007"></a>00007 <span class="comment"> * http://geographiclib.sourceforge.net/</span>
<a name="l00008"></a>00008 <span class="comment"> *</span>
<a name="l00009"></a>00009 <span class="comment"> * This is a reformulation of the geodesic problem. The notation is as</span>
<a name="l00010"></a>00010 <span class="comment"> * follows:</span>
<a name="l00011"></a>00011 <span class="comment"> * - at a general point (no suffix or 1 or 2 as suffix)</span>
<a name="l00012"></a>00012 <span class="comment"> * - phi = latitude</span>
<a name="l00013"></a>00013 <span class="comment"> * - beta = latitude on auxiliary sphere</span>
<a name="l00014"></a>00014 <span class="comment"> * - omega = longitude on auxiliary sphere</span>
<a name="l00015"></a>00015 <span class="comment"> * - lambda = longitude</span>
<a name="l00016"></a>00016 <span class="comment"> * - alpha = azimuth of great circle</span>
<a name="l00017"></a>00017 <span class="comment"> * - sigma = arc length along greate circle</span>
<a name="l00018"></a>00018 <span class="comment"> * - s = distance</span>
<a name="l00019"></a>00019 <span class="comment"> * - tau = scaled distance (= sigma at multiples of pi/2)</span>
<a name="l00020"></a>00020 <span class="comment"> * - at northwards equator crossing</span>
<a name="l00021"></a>00021 <span class="comment"> * - beta = phi = 0</span>
<a name="l00022"></a>00022 <span class="comment"> * - omega = lambda = 0</span>
<a name="l00023"></a>00023 <span class="comment"> * - alpha = alpha0</span>
<a name="l00024"></a>00024 <span class="comment"> * - sigma = s = 0</span>
<a name="l00025"></a>00025 <span class="comment"> * - a 12 suffix means a difference, e.g., s12 = s2 - s1.</span>
<a name="l00026"></a>00026 <span class="comment"> * - s and c prefixes mean sin and cos</span>
<a name="l00027"></a>00027 <span class="comment"> **********************************************************************/</span>
<a name="l00028"></a>00028
<a name="l00029"></a>00029 <span class="preprocessor">#include "<a class="code" href="GeodesicLine_8hpp.html" title="Header for GeographicLib::GeodesicLine class.">GeographicLib/GeodesicLine.hpp</a>"</span>
<a name="l00030"></a>00030
<a name="l00031"></a><a class="code" href="GeodesicLine_8cpp.html#a8c325bf3ebfa8126888b62b16495d64e">00031</a> <span class="preprocessor">#define GEOGRAPHICLIB_GEODESICLINE_CPP "$Id: GeodesicLine.cpp 6921 2010-12-31 14:34:50Z karney $"</span>
<a name="l00032"></a>00032 <span class="preprocessor"></span>
<a name="l00033"></a>00033 <a class="code" href="Constants_8hpp.html#af90fa899707a2ac513d5e4c76853bbf5">RCSID_DECL</a>(<a class="code" href="GeodesicLine_8cpp.html#a8c325bf3ebfa8126888b62b16495d64e">GEOGRAPHICLIB_GEODESICLINE_CPP</a>)
<a name="l00034"></a>00034 <a class="code" href="Constants_8hpp.html#af90fa899707a2ac513d5e4c76853bbf5">RCSID_DECL</a>(<a class="code" href="GeodesicLine_8hpp.html#a45305d595b59bfb061a5a0c6e1044485">GEOGRAPHICLIB_GEODESICLINE_HPP</a>)
<a name="l00035"></a>00035
<a name="l00036"></a>00036 namespace GeographicLib {
<a name="l00037"></a>00037
<a name="l00038"></a>00038 <span class="keyword">using namespace </span>std;
<a name="l00039"></a>00039
<a name="l00040"></a><a class="code" href="classGeographicLib_1_1GeodesicLine.html#a5669be85b4a574258f4a136f12d3f1bb">00040</a> GeodesicLine::GeodesicLine(<span class="keyword">const</span> <a class="code" href="classGeographicLib_1_1Geodesic.html" title="Geodesic calculations">Geodesic</a>& g,
<a name="l00041"></a>00041 real lat1, real lon1, real azi1,
<a name="l00042"></a>00042 <span class="keywordtype">unsigned</span> caps) <span class="keywordflow">throw</span>()
<a name="l00043"></a>00043 : _a(g._a)
<a name="l00044"></a>00044 , _r(g._r)
<a name="l00045"></a>00045 , _b(g._b)
<a name="l00046"></a>00046 , _c2(g._c2)
<a name="l00047"></a>00047 , _f1(g._f1)
<a name="l00048"></a>00048 <span class="comment">// Always allow latitude and azimuth</span>
<a name="l00049"></a>00049 , _caps(caps | LATITUDE | AZIMUTH)
<a name="l00050"></a>00050 {
<a name="l00051"></a>00051 azi1 = Geodesic::AngNormalize(azi1);
<a name="l00052"></a>00052 <span class="comment">// Guard against underflow in salp0</span>
<a name="l00053"></a>00053 azi1 = Geodesic::AngRound(azi1);
<a name="l00054"></a>00054 lon1 = Geodesic::AngNormalize(lon1);
<a name="l00055"></a>00055 _lat1 = lat1;
<a name="l00056"></a>00056 _lon1 = lon1;
<a name="l00057"></a>00057 _azi1 = azi1;
<a name="l00058"></a>00058 <span class="comment">// alp1 is in [0, pi]</span>
<a name="l00059"></a>00059 real alp1 = azi1 * Math::degree<real>();
<a name="l00060"></a>00060 <span class="comment">// Enforce sin(pi) == 0 and cos(pi/2) == 0. Better to face the ensuing</span>
<a name="l00061"></a>00061 <span class="comment">// problems directly than to skirt them.</span>
<a name="l00062"></a>00062 _salp1 = azi1 == -180 ? 0 : sin(alp1);
<a name="l00063"></a>00063 _calp1 = abs(azi1) == 90 ? 0 : cos(alp1);
<a name="l00064"></a>00064 real cbet1, sbet1, phi;
<a name="l00065"></a>00065 phi = lat1 * Math::degree<real>();
<a name="l00066"></a>00066 <span class="comment">// Ensure cbet1 = +epsilon at poles</span>
<a name="l00067"></a>00067 sbet1 = _f1 * sin(phi);
<a name="l00068"></a>00068 cbet1 = abs(lat1) == 90 ? Geodesic::eps2 : cos(phi);
<a name="l00069"></a>00069 Geodesic::SinCosNorm(sbet1, cbet1);
<a name="l00070"></a>00070
<a name="l00071"></a>00071 <span class="comment">// Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),</span>
<a name="l00072"></a>00072 _salp0 = _salp1 * cbet1; <span class="comment">// alp0 in [0, pi/2 - |bet1|]</span>
<a name="l00073"></a>00073 <span class="comment">// Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following</span>
<a name="l00074"></a>00074 <span class="comment">// is slightly better (consider the case salp1 = 0).</span>
<a name="l00075"></a>00075 _calp0 = Math::hypot(_calp1, _salp1 * sbet1);
<a name="l00076"></a>00076 <span class="comment">// Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).</span>
<a name="l00077"></a>00077 <span class="comment">// sig = 0 is nearest northward crossing of equator.</span>
<a name="l00078"></a>00078 <span class="comment">// With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).</span>
<a name="l00079"></a>00079 <span class="comment">// With bet1 = pi/2, alp1 = -pi, sig1 = pi/2</span>
<a name="l00080"></a>00080 <span class="comment">// With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2</span>
<a name="l00081"></a>00081 <span class="comment">// Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).</span>
<a name="l00082"></a>00082 <span class="comment">// With alp0 in (0, pi/2], quadrants for sig and omg coincide.</span>
<a name="l00083"></a>00083 <span class="comment">// No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.</span>
<a name="l00084"></a>00084 <span class="comment">// With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.</span>
<a name="l00085"></a>00085 _ssig1 = sbet1; _somg1 = _salp0 * sbet1;
<a name="l00086"></a>00086 _csig1 = _comg1 = sbet1 != 0 || _calp1 != 0 ? cbet1 * _calp1 : 1;
<a name="l00087"></a>00087 Geodesic::SinCosNorm(_ssig1, _csig1); <span class="comment">// sig1 in (-pi, pi]</span>
<a name="l00088"></a>00088 Geodesic::SinCosNorm(_somg1, _comg1);
<a name="l00089"></a>00089
<a name="l00090"></a>00090 _k2 = sq(_calp0) * g._ep2;
<a name="l00091"></a>00091 real eps = _k2 / (2 * (1 + sqrt(1 + _k2)) + _k2);
<a name="l00092"></a>00092
<a name="l00093"></a>00093 <span class="keywordflow">if</span> (_caps & CAP_C1) {
<a name="l00094"></a>00094 _A1m1 = Geodesic::A1m1f(eps);
<a name="l00095"></a>00095 Geodesic::C1f(eps, _C1a);
<a name="l00096"></a>00096 _B11 = Geodesic::SinCosSeries(<span class="keyword">true</span>, _ssig1, _csig1, _C1a, nC1);
<a name="l00097"></a>00097 real s = sin(_B11), c = cos(_B11);
<a name="l00098"></a>00098 <span class="comment">// tau1 = sig1 + B11</span>
<a name="l00099"></a>00099 _stau1 = _ssig1 * c + _csig1 * s;
<a name="l00100"></a>00100 _ctau1 = _csig1 * c - _ssig1 * s;
<a name="l00101"></a>00101 <span class="comment">// Not necessary because C1pa reverts C1a</span>
<a name="l00102"></a>00102 <span class="comment">// _B11 = -SinCosSeries(true, _stau1, _ctau1, _C1pa, nC1p);</span>
<a name="l00103"></a>00103 }
<a name="l00104"></a>00104
<a name="l00105"></a>00105 <span class="keywordflow">if</span> (_caps & CAP_C1p)
<a name="l00106"></a>00106 Geodesic::C1pf(eps, _C1pa);
<a name="l00107"></a>00107
<a name="l00108"></a>00108 <span class="keywordflow">if</span> (_caps & CAP_C2) {
<a name="l00109"></a>00109 _A2m1 = Geodesic::A2m1f(eps);
<a name="l00110"></a>00110 Geodesic::C2f(eps, _C2a);
<a name="l00111"></a>00111 _B21 = Geodesic::SinCosSeries(<span class="keyword">true</span>, _ssig1, _csig1, _C2a, nC2);
<a name="l00112"></a>00112 }
<a name="l00113"></a>00113
<a name="l00114"></a>00114 <span class="keywordflow">if</span> (_caps & CAP_C3) {
<a name="l00115"></a>00115 g.C3f(eps, _C3a);
<a name="l00116"></a>00116 _A3c = -g._f * _salp0 * g.A3f(eps);
<a name="l00117"></a>00117 _B31 = Geodesic::SinCosSeries(<span class="keyword">true</span>, _ssig1, _csig1, _C3a, nC3-1);
<a name="l00118"></a>00118 }
<a name="l00119"></a>00119
<a name="l00120"></a>00120 <span class="keywordflow">if</span> (_caps & CAP_C4) {
<a name="l00121"></a>00121 g.C4f(_k2, _C4a);
<a name="l00122"></a>00122 <span class="comment">// Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)</span>
<a name="l00123"></a>00123 _A4 = sq(g._a) * _calp0 * _salp0 * g._e2;
<a name="l00124"></a>00124 _B41 = Geodesic::SinCosSeries(<span class="keyword">false</span>, _ssig1, _csig1, _C4a, nC4);
<a name="l00125"></a>00125 }
<a name="l00126"></a>00126 }
<a name="l00127"></a>00127
<a name="l00128"></a><a class="code" href="classGeographicLib_1_1GeodesicLine.html#ad9522110fe6a5a0c7f912a6b1a23094d">00128</a> <a class="code" href="classGeographicLib_1_1Math.html#aeee4778d7cf2f9fb9648efe4911da59d">Math::real</a> GeodesicLine::GenPosition(<span class="keywordtype">bool</span> arcmode, real s12_a12,
<a name="l00129"></a>00129 <span class="keywordtype">unsigned</span> outmask,
<a name="l00130"></a>00130 real& lat2, real& lon2, real& azi2,
<a name="l00131"></a>00131 real& s12, real& m12,
<a name="l00132"></a>00132 real& M12, real& M21,
<a name="l00133"></a>00133 real& S12)
<a name="l00134"></a>00134 <span class="keyword">const</span> <span class="keywordflow">throw</span>() {
<a name="l00135"></a>00135 outmask &= _caps & OUT_ALL;
<a name="l00136"></a>00136 <span class="keywordflow">if</span> (!( Init() && (arcmode || (_caps & DISTANCE_IN & OUT_ALL)) ))
<a name="l00137"></a>00137 <span class="comment">// Uninitialized or impossible distance calculation requested</span>
<a name="l00138"></a>00138 <span class="keywordflow">return</span> Math::NaN();
<a name="l00139"></a>00139
<a name="l00140"></a>00140 <span class="comment">// Avoid warning about uninitialized B12.</span>
<a name="l00141"></a>00141 real sig12, ssig12, csig12, B12 = 0, AB1 = 0;
<a name="l00142"></a>00142 <span class="keywordflow">if</span> (arcmode) {
<a name="l00143"></a>00143 <span class="comment">// Interpret s12_a12 as spherical arc length</span>
<a name="l00144"></a>00144 sig12 = s12_a12 * Math::degree<real>();
<a name="l00145"></a>00145 real s12a = abs(s12_a12);
<a name="l00146"></a>00146 s12a -= 180 * floor(s12a / 180);
<a name="l00147"></a>00147 ssig12 = s12a == 0 ? 0 : sin(sig12);
<a name="l00148"></a>00148 csig12 = s12a == 90 ? 0 : cos(sig12);
<a name="l00149"></a>00149 } <span class="keywordflow">else</span> {
<a name="l00150"></a>00150 <span class="comment">// Interpret s12_a12 as distance</span>
<a name="l00151"></a>00151 real
<a name="l00152"></a>00152 tau12 = s12_a12 / (_b * (1 + _A1m1)),
<a name="l00153"></a>00153 s = sin(tau12),
<a name="l00154"></a>00154 c = cos(tau12);
<a name="l00155"></a>00155 <span class="comment">// tau2 = tau1 + tau12</span>
<a name="l00156"></a>00156 B12 = - Geodesic::SinCosSeries(<span class="keyword">true</span>, _stau1 * c + _ctau1 * s,
<a name="l00157"></a>00157 _ctau1 * c - _stau1 * s,
<a name="l00158"></a>00158 _C1pa, nC1p);
<a name="l00159"></a>00159 sig12 = tau12 - (B12 - _B11);
<a name="l00160"></a>00160 ssig12 = sin(sig12);
<a name="l00161"></a>00161 csig12 = cos(sig12);
<a name="l00162"></a>00162 }
<a name="l00163"></a>00163
<a name="l00164"></a>00164 real omg12, lam12, lon12;
<a name="l00165"></a>00165 real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2;
<a name="l00166"></a>00166 <span class="comment">// sig2 = sig1 + sig12</span>
<a name="l00167"></a>00167 ssig2 = _ssig1 * csig12 + _csig1 * ssig12;
<a name="l00168"></a>00168 csig2 = _csig1 * csig12 - _ssig1 * ssig12;
<a name="l00169"></a>00169 <span class="keywordflow">if</span> (outmask & (DISTANCE | REDUCEDLENGTH | GEODESICSCALE)) {
<a name="l00170"></a>00170 <span class="keywordflow">if</span> (arcmode)
<a name="l00171"></a>00171 B12 = Geodesic::SinCosSeries(<span class="keyword">true</span>, ssig2, csig2, _C1a, nC1);
<a name="l00172"></a>00172 AB1 = (1 + _A1m1) * (B12 - _B11);
<a name="l00173"></a>00173 }
<a name="l00174"></a>00174 <span class="comment">// sin(bet2) = cos(alp0) * sin(sig2)</span>
<a name="l00175"></a>00175 sbet2 = _calp0 * ssig2;
<a name="l00176"></a>00176 <span class="comment">// Alt: cbet2 = hypot(csig2, salp0 * ssig2);</span>
<a name="l00177"></a>00177 cbet2 = Math::hypot(_salp0, _calp0 * csig2);
<a name="l00178"></a>00178 <span class="keywordflow">if</span> (cbet2 == 0)
<a name="l00179"></a>00179 <span class="comment">// I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case</span>
<a name="l00180"></a>00180 cbet2 = csig2 = Geodesic::eps2;
<a name="l00181"></a>00181 <span class="comment">// tan(omg2) = sin(alp0) * tan(sig2)</span>
<a name="l00182"></a>00182 somg2 = _salp0 * ssig2; comg2 = csig2; <span class="comment">// No need to normalize</span>
<a name="l00183"></a>00183 <span class="comment">// tan(alp0) = cos(sig2)*tan(alp2)</span>
<a name="l00184"></a>00184 salp2 = _salp0; calp2 = _calp0 * csig2; <span class="comment">// No need to normalize</span>
<a name="l00185"></a>00185 <span class="comment">// omg12 = omg2 - omg1</span>
<a name="l00186"></a>00186 omg12 = atan2(somg2 * _comg1 - comg2 * _somg1,
<a name="l00187"></a>00187 comg2 * _comg1 + somg2 * _somg1);
<a name="l00188"></a>00188
<a name="l00189"></a>00189 <span class="keywordflow">if</span> (outmask & DISTANCE)
<a name="l00190"></a>00190 s12 = arcmode ? _b * ((1 + _A1m1) * sig12 + AB1) : s12_a12;
<a name="l00191"></a>00191
<a name="l00192"></a>00192 <span class="keywordflow">if</span> (outmask & LONGITUDE) {
<a name="l00193"></a>00193 lam12 = omg12 + _A3c *
<a name="l00194"></a>00194 ( sig12 + (Geodesic::SinCosSeries(<span class="keyword">true</span>, ssig2, csig2, _C3a, nC3-1)
<a name="l00195"></a>00195 - _B31));
<a name="l00196"></a>00196 lon12 = lam12 / Math::degree<real>();
<a name="l00197"></a>00197 <span class="comment">// Can't use AngNormalize because longitude might have wrapped multiple</span>
<a name="l00198"></a>00198 <span class="comment">// times.</span>
<a name="l00199"></a>00199 lon12 = lon12 - 360 * floor(lon12/360 + <a class="code" href="Geod_8cpp.html#a5caf95d46b184d9ca1d3764b3781b3c9">real</a>(0.5));
<a name="l00200"></a>00200 lon2 = Geodesic::AngNormalize(_lon1 + lon12);
<a name="l00201"></a>00201 }
<a name="l00202"></a>00202
<a name="l00203"></a>00203 <span class="keywordflow">if</span> (outmask & LATITUDE)
<a name="l00204"></a>00204 lat2 = atan2(sbet2, _f1 * cbet2) / Math::degree<real>();
<a name="l00205"></a>00205
<a name="l00206"></a>00206 <span class="keywordflow">if</span> (outmask & AZIMUTH)
<a name="l00207"></a>00207 <span class="comment">// minus signs give range [-180, 180). 0- converts -0 to +0.</span>
<a name="l00208"></a>00208 azi2 = 0 - atan2(-salp2, calp2) / Math::degree<real>();
<a name="l00209"></a>00209
<a name="l00210"></a>00210 <span class="keywordflow">if</span> (outmask & (REDUCEDLENGTH | GEODESICSCALE)) {
<a name="l00211"></a>00211 real
<a name="l00212"></a>00212 ssig1sq = sq(_ssig1),
<a name="l00213"></a>00213 ssig2sq = sq( ssig2),
<a name="l00214"></a>00214 w1 = sqrt(1 + _k2 * ssig1sq),
<a name="l00215"></a>00215 w2 = sqrt(1 + _k2 * ssig2sq),
<a name="l00216"></a>00216 B22 = Geodesic::SinCosSeries(<span class="keyword">true</span>, ssig2, csig2, _C2a, nC2),
<a name="l00217"></a>00217 AB2 = (1 + _A2m1) * (B22 - _B21),
<a name="l00218"></a>00218 J12 = (_A1m1 - _A2m1) * sig12 + (AB1 - AB2);
<a name="l00219"></a>00219 <span class="keywordflow">if</span> (outmask & REDUCEDLENGTH)
<a name="l00220"></a>00220 <span class="comment">// Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure</span>
<a name="l00221"></a>00221 <span class="comment">// accurate cancellation in the case of coincident points.</span>
<a name="l00222"></a>00222 m12 = _b * ((w2 * (_csig1 * ssig2) - w1 * (_ssig1 * csig2))
<a name="l00223"></a>00223 - _csig1 * csig2 * J12);
<a name="l00224"></a>00224 <span class="keywordflow">if</span> (outmask & GEODESICSCALE) {
<a name="l00225"></a>00225 M12 = csig12 + (_k2 * (ssig2sq - ssig1sq) * ssig2 / (w1 + w2)
<a name="l00226"></a>00226 - csig2 * J12) * _ssig1 / w1;
<a name="l00227"></a>00227 M21 = csig12 - (_k2 * (ssig2sq - ssig1sq) * _ssig1 / (w1 + w2)
<a name="l00228"></a>00228 - _csig1 * J12) * ssig2 / w2;
<a name="l00229"></a>00229 }
<a name="l00230"></a>00230 }
<a name="l00231"></a>00231
<a name="l00232"></a>00232 <span class="keywordflow">if</span> (outmask & AREA) {
<a name="l00233"></a>00233 real
<a name="l00234"></a>00234 B42 = Geodesic::SinCosSeries(<span class="keyword">false</span>, ssig2, csig2, _C4a, nC4);
<a name="l00235"></a>00235 real salp12, calp12;
<a name="l00236"></a>00236 <span class="keywordflow">if</span> (_calp0 == 0 || _salp0 == 0) {
<a name="l00237"></a>00237 <span class="comment">// alp12 = alp2 - alp1, used in atan2 so no need to normalized</span>
<a name="l00238"></a>00238 salp12 = salp2 * _calp1 - calp2 * _salp1;
<a name="l00239"></a>00239 calp12 = calp2 * _calp1 + salp2 * _salp1;
<a name="l00240"></a>00240 <span class="comment">// The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz</span>
<a name="l00241"></a>00241 <span class="comment">// salp12 = -0 and alp12 = -180. However this depends on the sign being</span>
<a name="l00242"></a>00242 <span class="comment">// attached to 0 correctly. The following ensures the correct behavior.</span>
<a name="l00243"></a>00243 <span class="keywordflow">if</span> (salp12 == 0 && calp12 < 0) {
<a name="l00244"></a>00244 salp12 = Geodesic::eps2 * _calp1;
<a name="l00245"></a>00245 calp12 = -1;
<a name="l00246"></a>00246 }
<a name="l00247"></a>00247 } <span class="keywordflow">else</span> {
<a name="l00248"></a>00248 <span class="comment">// tan(alp) = tan(alp0) * sec(sig)</span>
<a name="l00249"></a>00249 <span class="comment">// tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)</span>
<a name="l00250"></a>00250 <span class="comment">// = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)</span>
<a name="l00251"></a>00251 <span class="comment">// If csig12 > 0, write</span>
<a name="l00252"></a>00252 <span class="comment">// csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)</span>
<a name="l00253"></a>00253 <span class="comment">// else</span>
<a name="l00254"></a>00254 <span class="comment">// csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1</span>
<a name="l00255"></a>00255 <span class="comment">// No need to normalize</span>
<a name="l00256"></a>00256 salp12 = _calp0 * _salp0 *
<a name="l00257"></a>00257 (csig12 <= 0 ? _csig1 * (1 - csig12) + ssig12 * _ssig1 :
<a name="l00258"></a>00258 ssig12 * (_csig1 * ssig12 / (1 + csig12) + _ssig1));
<a name="l00259"></a>00259 calp12 = sq(_salp0) + sq(_calp0) * _csig1 * csig2;
<a name="l00260"></a>00260 }
<a name="l00261"></a>00261 S12 = _c2 * atan2(salp12, calp12) + _A4 * (B42 - _B41);
<a name="l00262"></a>00262 }
<a name="l00263"></a>00263
<a name="l00264"></a>00264 <span class="keywordflow">return</span> arcmode ? s12_a12 : sig12 / Math::degree<real>();
<a name="l00265"></a>00265 }
<a name="l00266"></a>00266 } <span class="comment">// namespace GeographicLib</span>
<a name="l00267"></a>00267
</pre></div></div>
</div>
<hr class="footer"/><address class="footer"><small>Generated on Tue Feb 22 2011 for GeographicLib by
<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.1 </small></address>
</body>
</html>
|