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
|
## Copyright (C) 2022 The Octave Project Developers
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {} {@var{lat}, @var{lon} =} gcxsc (@var{lat1}, @var{lon1}, @var{az}, @var{lat2}, @var{lon2}, @var{r1})
## @deftypefnx {} {@var{lat}, @var{lon} =} gcxsc (@var{lat1}, @var{lon1}, @var{az}, @var{lat2}, @var{lon2}, @var{r1}, @var{angleUnit})
## Determines the intersection points between a great circle and a small circle.
##
## Input:
## @itemize
## @item
## @var{lat1}, @var{lon1}, @var{az}: latitude, longitude, and azimuth of the
## great circle. These must be scalar values or vectors of equal length.
##
## @item
## @var{lat2}, @var{lon2}, @var{r2}: latitude, longitude, and range of the
## small circle. These must be scalar values or vectors of equal length.
##
## @item
## @var{angleUnit}: string for angular units ('degrees' or 'radians',
## case-insensitive, just the first character will do). Default is 'degrees'.
## @var{angleUnit} applies to all inputs and outputs.
## @end itemize
##
## Outputs:
##
## @itemize
## @item
## @var{lat} and @var{lon} are both Nx2 vectors of latitude(s) and longitude(s)
## of the intersection point(s). Circle pair(s) that have no intersection
## points or happen to lie on the same axis (identical or antipodal centers)
## NaN values are returned. @*
## If only one output vlues was requested
##
## @item
## Optional third output @var{istn}, if present, turns off warnings for
## coinciding circles or no intersections. It is an Nx1 vector indicating
## the intersection situation of each input pair of circles, with for each
## circle pair the values:
##
## @table @asis
## @item 0
## The pair of circles has two distinct intersection points.
##
## @item 1
## The circles have identical axis, so are either coinciding or don't have
## any intersection points.
##
## @item 2
## The pair of circles have just one common intersection point (tangent).
##
## @item 3
## The pair of circles are disjoint, have no intersection points.
## @end table
## @end itemize
##
## Example
## @example
## [newlat, newlon] = gcxsc (60, 25, 20, 55, 25, 2.5)
## newlat =
## 53.806 57.286
## newlon =
## 21.226 23.182
## @end example
## @seealso{gc2sc, gcxgc, scxsc}
## @end deftypefn
function [lat, lon, st] = gcxsc (varargin)
if (nargin < 6)
print_usage ();
elseif (nargin == 6)
angleUnit = "degrees";
else
angleUnit = varargin{7};
endif
if (! (all (cellfun ("isnumeric", varargin(1:6)) && ...
all (cellfun ("isreal", varargin(1:6))))))
error ("gcxsc: numeric values expected for first six inputs");
endif
isv = ! cellfun ("isscalar", varargin(1:6));
if (any (isv))
## At least one of the location inputs is a vector. Check sizes
numval = cellfun ("numel", varargin(isv));
if (any (diff (numval)))
error ("gcxgc: all vector inputs must have same lengths");
endif
nv = max (numval);
## Make sure all inputs are column vectors of same length
for ii=1:6
if (isv(ii))
varargin(ii) = {varargin{ii}(:)};
else
varargin(ii) = {(repmat (varargin{ii}, nv, 1))};
endif
endfor
else
nv = 1;
endif
vect = [varargin{1:6}];
if (! ischar (angleUnit))
error ("gcxsc: character value expected for 'angleUnit'");
elseif (strncmpi (angleUnit, "degrees", min (length (angleUnit), 7)))
## Latitude must be within (-90, 90)
if (any (abs (vect(:, [1 4])) >= 90))
error("gcxsc: latitude value(s) out of acceptable range (-90, 90)")
endif
vect = deg2rad (vect);
elseif (strncmpi (angleUnit, "radians", min (length (angleUnit), 7)))
## Latitude must be within (-pi/2, pi/2) as azimuth isn't defined there
if (any (abs (vect(:, [1 4])) >= pi / 2))
error("gcxsc: latitude value(s) out of acceptable range (-pi/2, pi/2)")
endif
else
error ("gcxsc: illegal input for 'angleUnit'");
endif
[slat, slong, srg] = gc2sc (vect(1), vect(2), vect(3), "r");
[lat, lon, st] = scxsc (slat, slong, srg, vect(4), vect(5), vect(6), "r");
if (nargout < 3)
if (any (st == 3))
warning ("Octave:coinciding-small-circles", ...
"scxsc: (some) circle pair(s) coincide.\n");
endif
if (any (st == 1))
warning ("Octave:no-intersecting-circles", ...
"scxsc: one or more circle pair(s) do not intersect.\n");
endif
endif
if (strncmpi (angleUnit, "degrees", length (angleUnit)))
lat = rad2deg (lat);
lon = rad2deg (lon);
endif
endfunction
%!test
%! [lat, lon] = gcxsc (60, 25, 20, 55, 25, 2.5);
%! assert (lat(1), 53.80564992, 1e-6);
%! assert (lon(1), 21.22598692, 1e-6);
%!error <numeric> gcxsc ("s", 0, 100, 10, 30, 0)
%!error <numeric> gcxsc (3i, 0, 100, 10, 30, 0)
%!error <numeric> gcxsc (50, "s", 100, 10, 30, 0)
%!error <numeric> gcxsc (50, 2i, 10, 10, 30, 0)
%!error <numeric> gcxsc (50, 0, "s", 10, 30, 0)
%!error <numeric> gcxsc (50, 0, 100i, 10, 30, 0)
%!error <numeric> gcxsc (50, 0, 100, "s", 30, 0)
%!error <numeric> gcxsc (50, 0, 100, 10i, 30, 0)
%!error <numeric> gcxsc (50, 0, 100, 10, "s", 0)
%!error <numeric> gcxsc (50, 0, 100, 10, 30i, 0)
%!error <numeric> gcxsc (50, 0, 100, 10, 30, "s")
%!error <numeric> gcxsc (50, 0, 100, 10, 30, 2i)
%!error <illegal> gcxsc (50, 0, 100, 10, 30, 0, "f")
%!error <illegal> gcxsc (50, 0, 100, 10, 30, 0, "degreef")
%!error <latitude value> gcxsc (190, 0, 90, -90.000001, 180, 80);
%!error <latitude value> gcxsc (190, 0, -90.001, -90.000001, 180, 80);
%!error <latitude value> gcxsc (pi/1.999, 0, pi/2, pi/2.0001, 2, 2*pi/3, "r");
|
