-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathgsw_distance.m
139 lines (122 loc) · 5 KB
/
gsw_distance.m
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
function distance = gsw_distance(long,lat,p)
% gsw_distance spherical earth distance between points
% in long, lat coordinates at given pressure
%==========================================================================
%
% USAGE:
% distance = gsw_distance(long,lat,{p})
%
% DESCRIPTION:
% Calculates the distance in metres between successive points in the
% vectors long and lat, computed using the Haversine formula on a
% spherical earth of radius 6,371 km, being the radius of a sphere having
% the same volume as Earth. For a sperical Earth of radius 6,371,000 m,
% one nautical mile is 1,853.2488 m, thus one degree of latitude is
% 111,194.93 m.
% Note. Distances are probably good to better than 1% of the "true"
% distance on the ellipsoidal earth.
%
% INPUT:
% long = longitude in decimal degress [ 0 ... +360 ]
% or [-180 ... +180 ]
% lat = latitude in decimal degress north [ -90 ... +90 ]
%
% OPTIONAL:
% p = sea pressure ( default is 0 ) [ dbar ]
% ( i.e. absolute pressure - 10.1325 dbar )
%
% lat and long need to have the same dimensions, Mx1 or 1xN or MxN.
% p, if provided, may have dimensions 1x1 or Mx1 or 1xN or MxN,
% where lat & long are Mx1 or 1xN or MxN.
%
% OUTPUT:
% distance = Distance between points on a spherical [ m ]
% Earth at pressure (p)
% Note. The output is in m not km.
%
% AUTHOR:
% 6th November, 2000 by Rich Pawlowicz [ [email protected] ]
% Note. This function was extracted from Rich Pawlowicz's m_map package,
% which is available from http://www.eos.ubc.ca/~rich/map.html
%
% MODIFIED:
% 4th April, 2011 by Paul Barker and Trevor McDougall.
%
% VERSION NUMBER: 3.05 (27th January 2015)
%
% REFERENCE:
% http://www.eos.ubc.ca/~rich/map.html
%
% The software is available from http://www.TEOS-10.org
%
%==========================================================================
%--------------------------------------------------------------------------
% Check variables
%--------------------------------------------------------------------------
if ~(nargin == 2 | nargin == 3)
error('gsw_distance: Requires either two or three inputs')
end %if
if nargin == 2
p = zeros(size(lat));
end %if
[mla,nla] = size(lat);
[mlo,nlo] = size(long);
[mp,np] = size(p);
if (mla ~= mlo | nla ~= nlo )
error('*** Input array dimensions in gsw_distance do not agree ***')
elseif mla == 1 & mlo == 1 & nla == 1 & nlo == 1
error('*** No, you need more than one point to find a distance! ***')
end
transposed = 0;
if (mla == 1) & (nla == 1) % lat is a scalar. Fill to size of p
lat = lat*ones(mp,np);
long = long*ones(mp,np);
elseif (nla == np) & (mla == 1) % lat is row vector,
lat = lat(ones(1,mp),:); % copy down each column.
long = long(ones(1,mp),:);
elseif (mla == mp) & (nla == 1) % lat is column vector.
lat = lat.';
long = long.';
transposed = 1;
lat = lat(ones(1,np),:);
long = long(ones(1,np),:);
elseif (mla == np) & (nla == 1) % lat is column vector and p is a column vector
lat = lat.';
long = long.';
transposed = 1;
elseif (mla == mp) & (nla == np)
% ok
end %if
[mla,nla] = size(lat);
if (mp == 1) & (np == 1) % p is a scalar. Fill to size of lat
p = p*ones(mla,nla);
elseif (np == nla) & (mp == 1) % p is row vector,
p = p(ones(1,mla),:); % copy down each column.
elseif (mp == nla) & (np == mla) % p is row vector,
p = p.';
elseif (np == nla) & (np == 1) % p is column vector.
p = p.'; % transpose.
p = p(:,ones(1,nla)); % copy along each row.
elseif (mp == mla) & (np == nla)
% ok
else
error('gsw_distance: p has wrong dimensions')
end %if
%--------------------------------------------------------------------------
% Start of the calculation
%--------------------------------------------------------------------------
pi180 = pi/180;
earth_radius = 6371000; % Earth's radius in metres.
dlong = pi180*(long(:,2:nla)-long(:,1:nla-1));
dlat = pi180*(lat(:,2:nla)-lat(:,1:nla-1));
a = (sin(0.5*dlat)).^2 + cos(lat(:,1:nla-1)*pi180).*cos(lat(:,2:nla)*pi180).*(sin(0.5*dlong)).^2;
angles = 2 * atan2(sqrt(a),sqrt(1-a));
p_mid = 0.5*(p(:,1:nla-1) + p(:,1:nla-1));
lat_mid = 0.5*(lat(:,1:nla-1) + lat(:,2:nla));
z = gsw_z_from_p(p_mid,lat_mid); % Note. z is height and is negative
% in the ocean.
distance = (earth_radius + z).*angles; % Note. The output is in m not km.
if transposed
distance = distance.';
end
end