| 1 |
dimitri |
1.1 |
function [range,A12,A21]=dist(lat,long,argu1,argu2); |
| 2 |
|
|
% DIST Computes distance and bearing between points on the earth |
| 3 |
|
|
% using various reference spheroids. |
| 4 |
|
|
% |
| 5 |
|
|
% [RANGE,AF,AR]=DIST(LAT,LONG) computes the ranges RANGE between |
| 6 |
|
|
% points specified in the LAT and LONG vectors (decimal degrees with |
| 7 |
|
|
% positive indicating north/east). Forward and reverse bearings |
| 8 |
|
|
% (degrees) are returned in AF, AR. |
| 9 |
|
|
% |
| 10 |
|
|
% [RANGE,GLAT,GLONG]=DIST(LAT,LONG,N) computes N-point geodesics |
| 11 |
|
|
% between successive points. Each successive geodesic occupies |
| 12 |
|
|
% it's own row (N>=2) |
| 13 |
|
|
% |
| 14 |
|
|
% [..]=DIST(...,'ellipsoid') uses the specified ellipsoid |
| 15 |
|
|
% to get distances and bearing. Available ellipsoids are: |
| 16 |
|
|
% |
| 17 |
|
|
% 'clarke66' Clarke 1866 |
| 18 |
|
|
% 'iau73' IAU 1973 |
| 19 |
|
|
% 'wgs84' WGS 1984 |
| 20 |
|
|
% 'sphere' Sphere of radius 6371.0 km |
| 21 |
|
|
% |
| 22 |
|
|
% The default is 'wgs84'. |
| 23 |
|
|
% |
| 24 |
|
|
% Ellipsoid formulas are recommended for distance d<2000 km, |
| 25 |
|
|
% but can be used for longer distances. |
| 26 |
|
|
|
| 27 |
|
|
%Notes: RP (WHOI) 3/Dec/91 |
| 28 |
|
|
% Mostly copied from BDC "dist.f" routine (copied from ....?), but |
| 29 |
|
|
% then wildly modified to bring it in line with Matlab vectorization. |
| 30 |
|
|
% |
| 31 |
|
|
% RP (WHOI) 6/Dec/91 |
| 32 |
|
|
% Feeping Creaturism! - added geodesic computations. This turned |
| 33 |
|
|
% out to be pretty hairy since there were a lot of branch problems |
| 34 |
|
|
% with asin, atan when computing geodesics subtending > 90 degrees |
| 35 |
|
|
% that were ignored in the original code! |
| 36 |
|
|
% RP (WHOI) 15/Jan/91 |
| 37 |
|
|
% Fixed some bothersome special cases, like when computing geodesics |
| 38 |
|
|
% and N=2, or LAT=0... |
| 39 |
|
|
|
| 40 |
|
|
%C GIVEN THE LATITUDES AND LONGITUDES (IN DEG.) IT ASSUMES THE IAU SPHERO |
| 41 |
|
|
%C DEFINED IN THE NOTES ON PAGE 523 OF THE EXPLANATORY SUPPLEMENT TO THE |
| 42 |
|
|
%C AMERICAN EPHEMERIS. |
| 43 |
|
|
%C |
| 44 |
|
|
%C THIS PROGRAM COMPUTES THE DISTANCE ALONG THE NORMAL |
| 45 |
|
|
%C SECTION (IN M.) OF A SPECIFIED REFERENCE SPHEROID GIVEN |
| 46 |
|
|
%C THE GEODETIC LATITUDES AND LONGITUDES OF THE END POINTS |
| 47 |
|
|
%C *** IN DECIMAL DEGREES *** |
| 48 |
|
|
%C |
| 49 |
|
|
%C IT USES ROBBIN'S FORMULA, AS GIVEN BY BOMFORD, GEODESY, |
| 50 |
|
|
%C FOURTH EDITION, P. 122. CORRECT TO ONE PART IN 10**8 |
| 51 |
|
|
%C AT 1600 KM. ERRORS OF 20 M AT 5000 KM. |
| 52 |
|
|
%C |
| 53 |
|
|
%C CHECK: SMITHSONIAN METEOROLOGICAL TABLES, PP. 483 AND 484, |
| 54 |
|
|
%C GIVES LENGTHS OF ONE DEGREE OF LATITUDE AND LONGITUDE |
| 55 |
|
|
%C AS A FUNCTION OF LATITUDE. (SO DOES THE EPHEMERIS ABOVE) |
| 56 |
|
|
%C |
| 57 |
|
|
%C PETER WORCESTER, AS TOLD TO BRUCE CORNUELLE...1983 MAY 27 |
| 58 |
|
|
%C |
| 59 |
|
|
|
| 60 |
|
|
spheroid='wgs84'; |
| 61 |
|
|
geodes=0; |
| 62 |
|
|
if (nargin >= 3), |
| 63 |
|
|
if (isstr(argu1)), |
| 64 |
|
|
spheroid=argu1; |
| 65 |
|
|
else |
| 66 |
|
|
geodes=1; |
| 67 |
|
|
Ngeodes=argu1; |
| 68 |
|
|
if (Ngeodes <2), error('Must have at least 2 points in a goedesic!');end; |
| 69 |
|
|
if (nargin==4), spheroid=argu2; end; |
| 70 |
|
|
end; |
| 71 |
|
|
end; |
| 72 |
|
|
|
| 73 |
|
|
if (spheroid(1:3)=='sph'), |
| 74 |
|
|
A = 6371000.0; |
| 75 |
|
|
B = A; |
| 76 |
|
|
E = sqrt(A*A-B*B)/A; |
| 77 |
|
|
EPS= E*E/(1-E*E); |
| 78 |
|
|
elseif (spheroid(1:3)=='cla'), |
| 79 |
|
|
A = 6378206.4E0; |
| 80 |
|
|
B = 6356583.8E0; |
| 81 |
|
|
E= sqrt(A*A-B*B)/A; |
| 82 |
|
|
EPS = E*E/(1.-E*E); |
| 83 |
|
|
elseif(spheroid(1:3)=='iau'), |
| 84 |
|
|
A = 6378160.e0; |
| 85 |
|
|
B = 6356774.516E0; |
| 86 |
|
|
E = sqrt(A*A-B*B)/A; |
| 87 |
|
|
EPS = E*E/(1.-E*E); |
| 88 |
|
|
elseif(spheroid(1:3)=='wgs'), |
| 89 |
|
|
|
| 90 |
|
|
%c on 9/11/88, Peter Worcester gave me the constants for the |
| 91 |
|
|
%c WGS84 spheroid, and he gave A (semi-major axis), F = (A-B)/A |
| 92 |
|
|
%c (flattening) (where B is the semi-minor axis), and E is the |
| 93 |
|
|
%c eccentricity, E = ( (A**2 - B**2)**.5 )/ A |
| 94 |
|
|
%c the numbers from peter are: A=6378137.; 1/F = 298.257223563 |
| 95 |
|
|
%c E = 0.081819191 |
| 96 |
|
|
A = 6378137.; |
| 97 |
|
|
E = 0.081819191; |
| 98 |
|
|
EPS= E*E/(1.-E*E); |
| 99 |
|
|
|
| 100 |
|
|
B = sqrt(A^2-(E*A)^2); % added by D Menemenlis, 4 nov 97 |
| 101 |
|
|
|
| 102 |
|
|
else |
| 103 |
|
|
error('dist: Unknown spheroid specified!'); |
| 104 |
|
|
end; |
| 105 |
|
|
|
| 106 |
|
|
|
| 107 |
|
|
NN=max(size(lat)); |
| 108 |
|
|
if (NN ~= max(size(long))), |
| 109 |
|
|
error('dist: Lat, Long vectors of different sizes!'); |
| 110 |
|
|
end |
| 111 |
|
|
|
| 112 |
|
|
if (NN==size(lat)), rowvec=0; % It is easier if things are column vectors, |
| 113 |
|
|
else rowvec=1; end; % but we have to fix things before returning! |
| 114 |
|
|
|
| 115 |
|
|
lat=lat(:)*pi/180; % convert to radians |
| 116 |
|
|
long=long(:)*pi/180; |
| 117 |
|
|
|
| 118 |
|
|
lat(lat==0)=eps*ones(sum(lat==0),1); % Fixes some nasty 0/0 cases in the |
| 119 |
|
|
% geodesics stuff |
| 120 |
|
|
|
| 121 |
|
|
PHI1=lat(1:NN-1); % endpoints of each segment |
| 122 |
|
|
XLAM1=long(1:NN-1); |
| 123 |
|
|
PHI2=lat(2:NN); |
| 124 |
|
|
XLAM2=long(2:NN); |
| 125 |
|
|
|
| 126 |
|
|
% wiggle lines of constant lat to prevent numerical probs. |
| 127 |
|
|
if (any(PHI1==PHI2)), |
| 128 |
|
|
for ii=1:NN-1, |
| 129 |
|
|
if (PHI1(ii)==PHI2(ii)), PHI2(ii)=PHI2(ii)+ 1e-14; end; |
| 130 |
|
|
end; |
| 131 |
|
|
end; |
| 132 |
|
|
% wiggle lines of constant long to prevent numerical probs. |
| 133 |
|
|
if (any(XLAM1==XLAM2)), |
| 134 |
|
|
for ii=1:NN-1, |
| 135 |
|
|
if (XLAM1(ii)==XLAM2(ii)), XLAM2(ii)=XLAM2(ii)+ 1e-14; end; |
| 136 |
|
|
end; |
| 137 |
|
|
end; |
| 138 |
|
|
|
| 139 |
|
|
|
| 140 |
|
|
|
| 141 |
|
|
%C COMPUTE THE RADIUS OF CURVATURE IN THE PRIME VERTICAL FOR |
| 142 |
|
|
%C EACH POINT |
| 143 |
|
|
|
| 144 |
|
|
xnu=A./sqrt(1.0-(E*sin(lat)).^2); |
| 145 |
|
|
xnu1=xnu(1:NN-1); |
| 146 |
|
|
xnu2=xnu(2:NN); |
| 147 |
|
|
|
| 148 |
|
|
%C*** COMPUTE THE AZIMUTHS. A12 (A21) IS THE AZIMUTH AT POINT 1 (2) |
| 149 |
|
|
%C OF THE NORMAL SECTION CONTAINING THE POINT 2 (1) |
| 150 |
|
|
|
| 151 |
|
|
TPSI2=(1.-E*E)*tan(PHI2) + E*E*xnu1.*sin(PHI1)./(xnu2.*cos(PHI2)); |
| 152 |
|
|
PSI2=atan(TPSI2); |
| 153 |
|
|
|
| 154 |
|
|
%C*** SOME FORM OF ANGLE DIFFERENCE COMPUTED HERE?? |
| 155 |
|
|
|
| 156 |
|
|
DPHI2=PHI2-PSI2; |
| 157 |
|
|
DLAM=XLAM2-XLAM1; |
| 158 |
|
|
CTA12=(cos(PHI1).*TPSI2 - sin(PHI1).*cos(DLAM))./sin(DLAM); |
| 159 |
|
|
A12=atan((1.)./CTA12); |
| 160 |
|
|
CTA21P=(sin(PSI2).*cos(DLAM) - cos(PSI2).*tan(PHI1))./sin(DLAM); |
| 161 |
|
|
A21P=atan((1.)./CTA21P); |
| 162 |
|
|
|
| 163 |
|
|
%C GET THE QUADRANT RIGHT |
| 164 |
|
|
DLAM2=(abs(DLAM)<pi).*DLAM + (DLAM>=pi).*(-2*pi+DLAM) + ... |
| 165 |
|
|
(DLAM<=-pi).*(2*pi+DLAM); |
| 166 |
|
|
A12=A12+(A12<-pi)*2*pi-(A12>=pi)*2*pi; |
| 167 |
|
|
A12=A12+pi*sign(-A12).*( sign(A12) ~= sign(DLAM2) ); |
| 168 |
|
|
A21P=A21P+(A21P<-pi)*2*pi-(A21P>=pi)*2*pi; |
| 169 |
|
|
A21P=A21P+pi*sign(-A21P).*( sign(A21P) ~= sign(-DLAM2) ); |
| 170 |
|
|
%%A12*180/pi |
| 171 |
|
|
%%A21P*180/pi |
| 172 |
|
|
|
| 173 |
|
|
|
| 174 |
|
|
SSIG=sin(DLAM).*cos(PSI2)./sin(A12); |
| 175 |
|
|
% At this point we are OK if the angle < 90...but otherwise |
| 176 |
|
|
% we get the wrong branch of asin! |
| 177 |
|
|
% This fudge will correct every case on a sphere, and *almost* |
| 178 |
|
|
% every case on an ellipsoid (wrong hnadling will be when |
| 179 |
|
|
% angle is almost exactly 90 degrees) |
| 180 |
|
|
dd2=[cos(long).*cos(lat) sin(long).*cos(lat) sin(lat)]; |
| 181 |
|
|
dd2=sum((diff(dd2).*diff(dd2))')'; |
| 182 |
|
|
if ( any(abs(dd2-2) < 2*((B-A)/A))^2 ), |
| 183 |
|
|
disp('dist: Warning...point(s) too close to 90 degrees apart'); |
| 184 |
|
|
end; |
| 185 |
|
|
bigbrnch=dd2>2; |
| 186 |
|
|
|
| 187 |
|
|
SIG=asin(SSIG).*(bigbrnch==0) + (pi-asin(SSIG)).*bigbrnch; |
| 188 |
|
|
|
| 189 |
|
|
SSIGC=-sin(DLAM).*cos(PHI1)./sin(A21P); |
| 190 |
|
|
SIGC=asin(SSIGC); |
| 191 |
|
|
A21 = A21P - DPHI2.*sin(A21P).*tan(SIG/2.0); |
| 192 |
|
|
|
| 193 |
|
|
%C COMPUTE RANGE |
| 194 |
|
|
|
| 195 |
|
|
G2=EPS*(sin(PHI1)).^2; |
| 196 |
|
|
G=sqrt(G2); |
| 197 |
|
|
H2=EPS*(cos(PHI1).*cos(A12)).^2; |
| 198 |
|
|
H=sqrt(H2); |
| 199 |
|
|
TERM1=-SIG.*SIG.*H2.*(1.0-H2)/6.0; |
| 200 |
|
|
TERM2=(SIG.^3).*G.*H.*(1.0-2.0*H2)/8.0; |
| 201 |
|
|
TERM3=(SIG.^4).*(H2.*(4.0-7.0*H2)-3.0*G2.*(1.0-7.0*H2))/120.0; |
| 202 |
|
|
TERM4=-(SIG.^5).*G.*H/48.0; |
| 203 |
|
|
|
| 204 |
|
|
range=xnu1.*SIG.*(1.0+TERM1+TERM2+TERM3+TERM4); |
| 205 |
|
|
|
| 206 |
|
|
|
| 207 |
|
|
if (geodes), |
| 208 |
|
|
|
| 209 |
|
|
%c now calculate the locations along the ray path. (for extra accuracy, could |
| 210 |
|
|
%c do it from start to halfway, then from end for the rest, switching from A12 |
| 211 |
|
|
%c to A21... |
| 212 |
|
|
%c started to use Rudoe's formula, page 117 in Bomford...(1980, fourth edition) |
| 213 |
|
|
%c but then went to Clarke's best formula (pg 118) |
| 214 |
|
|
|
| 215 |
|
|
%RP I am doing this twice because this formula doesn't work when we go |
| 216 |
|
|
%past 90 degrees! |
| 217 |
|
|
Ngd1=round(Ngeodes/2); |
| 218 |
|
|
|
| 219 |
|
|
% First time...away from point 1 |
| 220 |
|
|
if (Ngd1>1), |
| 221 |
|
|
wns=ones(1,Ngd1); |
| 222 |
|
|
CP1CA12 = (cos(PHI1).*cos(A12)).^2; |
| 223 |
|
|
R2PRM = -EPS.*CP1CA12; |
| 224 |
|
|
R3PRM = 3.0*EPS.*(1.0-R2PRM).*cos(PHI1).*sin(PHI1).*cos(A12); |
| 225 |
|
|
C1 = R2PRM.*(1.0+R2PRM)/6.0*wns; |
| 226 |
|
|
C2 = R3PRM.*(1.0+3.0*R2PRM)/24.0*wns; |
| 227 |
|
|
R2PRM=R2PRM*wns; |
| 228 |
|
|
R3PRM=R3PRM*wns; |
| 229 |
|
|
|
| 230 |
|
|
%c now have to loop over positions |
| 231 |
|
|
RLRAT = (range./xnu1)*([0:Ngd1-1]/(Ngeodes-1)); |
| 232 |
|
|
|
| 233 |
|
|
THETA = RLRAT.*(1 - (RLRAT.^2).*(C1 - C2.*RLRAT)); |
| 234 |
|
|
C3 = 1.0 - (R2PRM.*(THETA.^2))/2.0 - (R3PRM.*(THETA.^3))/6.0; |
| 235 |
|
|
DSINPSI =(sin(PHI1)*wns).*cos(THETA) + ... |
| 236 |
|
|
((cos(PHI1).*cos(A12))*wns).*sin(THETA); |
| 237 |
|
|
%try to identify the branch...got to other branch if range> 1/4 circle |
| 238 |
|
|
PSI = asin(DSINPSI); |
| 239 |
|
|
|
| 240 |
|
|
DCOSPSI = cos(PSI); |
| 241 |
|
|
DSINDLA = (sin(A12)*wns).*sin(THETA)./DCOSPSI; |
| 242 |
|
|
DTANPHI=(1.0+EPS)*(1.0 - (E^2)*C3.*(sin(PHI1)*wns)./DSINPSI).*tan(PSI); |
| 243 |
|
|
%C compute output latitude (phi) and long (xla) in radians |
| 244 |
|
|
%c I believe these are absolute, and don't need source coords added |
| 245 |
|
|
PHI = atan(DTANPHI); |
| 246 |
|
|
% fix branch cut stuff - |
| 247 |
|
|
otherbrcnh= sign(DLAM2*wns) ~= sign([sign(DLAM2) diff(DSINDLA')'] ); |
| 248 |
|
|
XLA = XLAM1*wns + asin(DSINDLA).*(otherbrcnh==0) + ... |
| 249 |
|
|
(pi-asin(DSINDLA)).*(otherbrcnh); |
| 250 |
|
|
else |
| 251 |
|
|
PHI=PHI1; |
| 252 |
|
|
XLA=XLAM1; |
| 253 |
|
|
end; |
| 254 |
|
|
|
| 255 |
|
|
% Now we do the same thing, but in the reverse direction from the receiver! |
| 256 |
|
|
if (Ngeodes-Ngd1>1), |
| 257 |
|
|
wns=ones(1,Ngeodes-Ngd1); |
| 258 |
|
|
CP2CA21 = (cos(PHI2).*cos(A21)).^2; |
| 259 |
|
|
R2PRM = -EPS.*CP2CA21; |
| 260 |
|
|
R3PRM = 3.0*EPS.*(1.0-R2PRM).*cos(PHI2).*sin(PHI2).*cos(A21); |
| 261 |
|
|
C1 = R2PRM.*(1.0+R2PRM)/6.0*wns; |
| 262 |
|
|
C2 = R3PRM.*(1.0+3.0*R2PRM)/24.0*wns; |
| 263 |
|
|
R2PRM=R2PRM*wns; |
| 264 |
|
|
R3PRM=R3PRM*wns; |
| 265 |
|
|
|
| 266 |
|
|
%c now have to loop over positions |
| 267 |
|
|
RLRAT = (range./xnu2)*([0:Ngeodes-Ngd1-1]/(Ngeodes-1)); |
| 268 |
|
|
|
| 269 |
|
|
THETA = RLRAT.*(1 - (RLRAT.^2).*(C1 - C2.*RLRAT)); |
| 270 |
|
|
C3 = 1.0 - (R2PRM.*(THETA.^2))/2.0 - (R3PRM.*(THETA.^3))/6.0; |
| 271 |
|
|
DSINPSI =(sin(PHI2)*wns).*cos(THETA) + ... |
| 272 |
|
|
((cos(PHI2).*cos(A21))*wns).*sin(THETA); |
| 273 |
|
|
%try to identify the branch...got to other branch if range> 1/4 circle |
| 274 |
|
|
PSI = asin(DSINPSI); |
| 275 |
|
|
|
| 276 |
|
|
DCOSPSI = cos(PSI); |
| 277 |
|
|
DSINDLA = (sin(A21)*wns).*sin(THETA)./DCOSPSI; |
| 278 |
|
|
DTANPHI=(1.0+EPS)*(1.0 - (E^2)*C3.*(sin(PHI2)*wns)./DSINPSI).*tan(PSI); |
| 279 |
|
|
%C compute output latitude (phi) and long (xla) in radians |
| 280 |
|
|
%c I believe these are absolute, and don't need source coords added |
| 281 |
|
|
PHI = [PHI fliplr(atan(DTANPHI))]; |
| 282 |
|
|
% fix branch cut stuff |
| 283 |
|
|
otherbrcnh= sign(-DLAM2*wns) ~= sign( [sign(-DLAM2) diff(DSINDLA')'] ); |
| 284 |
|
|
XLA = [XLA fliplr(XLAM2*wns + asin(DSINDLA).*(otherbrcnh==0) + ... |
| 285 |
|
|
(pi-asin(DSINDLA)).*(otherbrcnh))]; |
| 286 |
|
|
else |
| 287 |
|
|
PHI = [PHI PHI2]; |
| 288 |
|
|
XLA = [XLA XLAM2]; |
| 289 |
|
|
end; |
| 290 |
|
|
|
| 291 |
|
|
%c convert to degrees |
| 292 |
|
|
A12 = PHI*180/pi; |
| 293 |
|
|
A21 = XLA*180/pi; |
| 294 |
|
|
range=range*([0:Ngeodes-1]/(Ngeodes-1)); |
| 295 |
|
|
|
| 296 |
|
|
|
| 297 |
|
|
else |
| 298 |
|
|
|
| 299 |
|
|
%C*** CONVERT TO DECIMAL DEGREES |
| 300 |
|
|
A12=A12*180/pi; |
| 301 |
|
|
A21=A21*180/pi; |
| 302 |
|
|
if (rowvec), |
| 303 |
|
|
range=range'; |
| 304 |
|
|
A12=A12'; |
| 305 |
|
|
A21=A21'; |
| 306 |
|
|
end; |
| 307 |
|
|
end; |
| 308 |
|
|
|
| 309 |
|
|
|
| 310 |
|
|
|
| 311 |
|
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph