00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <GeographicLib/GeodesicLineExact.hpp>
00030
00031 namespace GeographicLib {
00032
00033 using namespace std;
00034
00035 GeodesicLineExact::GeodesicLineExact(const GeodesicExact& g,
00036 real lat1, real lon1, real azi1,
00037 unsigned caps)
00038 : tiny_(g.tiny_)
00039 , _a(g._a)
00040 , _f(g._f)
00041 , _b(g._b)
00042 , _c2(g._c2)
00043 , _f1(g._f1)
00044 , _e2(g._e2)
00045 , _E(0, 0)
00046
00047 , _caps(caps | LATITUDE | AZIMUTH)
00048 {
00049 azi1 = GeodesicExact::AngRound(Math::AngNormalize(azi1));
00050 lon1 = Math::AngNormalize(lon1);
00051 _lat1 = lat1;
00052 _lon1 = lon1;
00053 _azi1 = azi1;
00054
00055 real alp1 = azi1 * Math::degree();
00056
00057
00058 _salp1 = azi1 == -180 ? 0 : sin(alp1);
00059 _calp1 = abs(azi1) == 90 ? 0 : cos(alp1);
00060 real cbet1, sbet1, phi;
00061 phi = lat1 * Math::degree();
00062
00063 sbet1 = _f1 * sin(phi);
00064 cbet1 = abs(lat1) == 90 ? tiny_ : cos(phi);
00065 GeodesicExact::SinCosNorm(sbet1, cbet1);
00066 _dn1 = (_f >= 0 ? sqrt(1 + g._ep2 * Math::sq(sbet1)) :
00067 sqrt(1 - _e2 * Math::sq(cbet1)) / _f1);
00068
00069
00070 _salp0 = _salp1 * cbet1;
00071
00072
00073 _calp0 = Math::hypot(_calp1, _salp1 * sbet1);
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083 _ssig1 = sbet1; _somg1 = _salp0 * sbet1;
00084 _csig1 = _comg1 = sbet1 != 0 || _calp1 != 0 ? cbet1 * _calp1 : 1;
00085
00086 _cchi1 = _f1 * _dn1 * _comg1;
00087 GeodesicExact::SinCosNorm(_ssig1, _csig1);
00088
00089
00090
00091 _k2 = Math::sq(_calp0) * g._ep2;
00092 _E.Reset(-_k2, -g._ep2, 1 + _k2, 1 + g._ep2);
00093
00094 if (_caps & CAP_E) {
00095 _E0 = _E.E() / (Math::pi() / 2);
00096 _E1 = _E.deltaE(_ssig1, _csig1, _dn1);
00097 real s = sin(_E1), c = cos(_E1);
00098
00099 _stau1 = _ssig1 * c + _csig1 * s;
00100 _ctau1 = _csig1 * c - _ssig1 * s;
00101
00102
00103 }
00104
00105 if (_caps & CAP_D) {
00106 _D0 = _E.D() / (Math::pi() / 2);
00107 _D1 = _E.deltaD(_ssig1, _csig1, _dn1);
00108 }
00109
00110 if (_caps & CAP_H) {
00111 _H0 = _E.H() / (Math::pi() / 2);
00112 _H1 = _E.deltaH(_ssig1, _csig1, _dn1);
00113 }
00114
00115 if (_caps & CAP_C4) {
00116 real eps = _k2 / (2 * (1 + sqrt(1 + _k2)) + _k2);
00117 g.C4f(eps, _C4a);
00118
00119 _A4 = Math::sq(_a) * _calp0 * _salp0 * _e2;
00120 _B41 = GeodesicExact::CosSeries(_ssig1, _csig1, _C4a, nC4_);
00121 }
00122 }
00123
00124 Math::real GeodesicLineExact::GenPosition(bool arcmode, real s12_a12,
00125 unsigned outmask,
00126 real& lat2, real& lon2, real& azi2,
00127 real& s12, real& m12,
00128 real& M12, real& M21,
00129 real& S12)
00130 const {
00131 outmask &= _caps & OUT_ALL;
00132 if (!( Init() && (arcmode || (_caps & DISTANCE_IN & OUT_ALL)) ))
00133
00134 return Math::NaN();
00135
00136
00137 real sig12, ssig12, csig12, E2 = 0, AB1 = 0;
00138 if (arcmode) {
00139
00140 sig12 = s12_a12 * Math::degree();
00141 real s12a = abs(s12_a12);
00142 s12a -= 180 * floor(s12a / 180);
00143 ssig12 = s12a == 0 ? 0 : sin(sig12);
00144 csig12 = s12a == 90 ? 0 : cos(sig12);
00145 } else {
00146
00147 real
00148 tau12 = s12_a12 / (_b * _E0),
00149 s = sin(tau12),
00150 c = cos(tau12);
00151
00152 E2 = - _E.deltaEinv(_stau1 * c + _ctau1 * s, _ctau1 * c - _stau1 * s);
00153 sig12 = tau12 - (E2 - _E1);
00154 ssig12 = sin(sig12);
00155 csig12 = cos(sig12);
00156 }
00157
00158 real lam12, lon12;
00159 real ssig2, csig2, sbet2, cbet2, salp2, calp2;
00160
00161 ssig2 = _ssig1 * csig12 + _csig1 * ssig12;
00162 csig2 = _csig1 * csig12 - _ssig1 * ssig12;
00163 real dn2 = _E.Delta(ssig2, csig2);
00164 if (outmask & (DISTANCE | REDUCEDLENGTH | GEODESICSCALE)) {
00165 if (arcmode) {
00166 E2 = _E.deltaE(ssig2, csig2, dn2);
00167 }
00168 AB1 = _E0 * (E2 - _E1);
00169 }
00170
00171 sbet2 = _calp0 * ssig2;
00172
00173 cbet2 = Math::hypot(_salp0, _calp0 * csig2);
00174 if (cbet2 == 0)
00175
00176 cbet2 = csig2 = tiny_;
00177
00178 salp2 = _salp0; calp2 = _calp0 * csig2;
00179
00180 if (outmask & DISTANCE)
00181 s12 = arcmode ? _b * (_E0 * sig12 + AB1) : s12_a12;
00182
00183 if (outmask & LONGITUDE) {
00184 real somg2 = _salp0 * ssig2, comg2 = csig2;
00185
00186 real cchi2 = _f1 * dn2 * comg2;
00187 lam12 = atan2(somg2 * _cchi1 - cchi2 * _somg1,
00188 cchi2 * _cchi1 + somg2 * _somg1) -
00189 _e2/_f1 * _salp0 * _H0 * (sig12 + _E.deltaH(ssig2, csig2, dn2) - _H1 );
00190 lon12 = lam12 / Math::degree();
00191
00192
00193 lon12 = Math::AngNormalize2(lon12);
00194 lon2 = Math::AngNormalize(_lon1 + lon12);
00195 }
00196
00197 if (outmask & LATITUDE)
00198 lat2 = atan2(sbet2, _f1 * cbet2) / Math::degree();
00199
00200 if (outmask & AZIMUTH)
00201
00202 azi2 = 0 - atan2(-salp2, calp2) / Math::degree();
00203
00204 if (outmask & (REDUCEDLENGTH | GEODESICSCALE)) {
00205 real J12 = _k2 * _D0 * (sig12 + _E.deltaD(ssig2, csig2, dn2) - _D1);
00206 if (outmask & REDUCEDLENGTH)
00207
00208
00209 m12 = _b * ((dn2 * (_csig1 * ssig2) - _dn1 * (_ssig1 * csig2))
00210 - _csig1 * csig2 * J12);
00211 if (outmask & GEODESICSCALE) {
00212 real t = _k2 * (ssig2 - _ssig1) * (ssig2 + _ssig1) / (_dn1 + dn2);
00213 M12 = csig12 + (t * ssig2 - csig2 * J12) * _ssig1 / _dn1;
00214 M21 = csig12 - (t * _ssig1 - _csig1 * J12) * ssig2 / dn2;
00215 }
00216 }
00217
00218 if (outmask & AREA) {
00219 real
00220 B42 = GeodesicExact::CosSeries(ssig2, csig2, _C4a, nC4_);
00221 real salp12, calp12;
00222 if (_calp0 == 0 || _salp0 == 0) {
00223
00224 salp12 = salp2 * _calp1 - calp2 * _salp1;
00225 calp12 = calp2 * _calp1 + salp2 * _salp1;
00226
00227
00228
00229 if (salp12 == 0 && calp12 < 0) {
00230 salp12 = tiny_ * _calp1;
00231 calp12 = -1;
00232 }
00233 } else {
00234
00235
00236
00237
00238
00239
00240
00241
00242 salp12 = _calp0 * _salp0 *
00243 (csig12 <= 0 ? _csig1 * (1 - csig12) + ssig12 * _ssig1 :
00244 ssig12 * (_csig1 * ssig12 / (1 + csig12) + _ssig1));
00245 calp12 = Math::sq(_salp0) + Math::sq(_calp0) * _csig1 * csig2;
00246 }
00247 S12 = _c2 * atan2(salp12, calp12) + _A4 * (B42 - _B41);
00248 }
00249
00250 return arcmode ? s12_a12 : sig12 / Math::degree();
00251 }
00252
00253 }