diff options
Diffstat (limited to 'inst/scxsc.m')
| -rw-r--r-- | inst/scxsc.m | 42 |
1 files changed, 21 insertions, 21 deletions
diff --git a/inst/scxsc.m b/inst/scxsc.m index cb29f9c..ad713be 100644 --- a/inst/scxsc.m +++ b/inst/scxsc.m @@ -1,4 +1,4 @@ -## Copyright (C) 2020 The Octave Project Developers +## 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 @@ -179,40 +179,40 @@ function [lat, lon, idl] = scxsc (varargin) ## Explore spherical distances between circle centers. ## Use haversine formula for approx.distances < 25..50 degrees - id = abs ((vect(:, 1) + 360) .- (vect(:, 4) + 360)) > 25; + id = abs ((vect(:, 1) + 360) - (vect(:, 4) + 360)) > 25; id = id | abs ((vect(:, 2) + 180) - (vect(:, 5) + 180)) > 25; ## Plain spherical cosine formula sphd = NaN (nv, 1); - sphd(id) = acosd (sind (vect(id, 1)) .* sind (vect(id, 4)) .+ ... + sphd(id) = acosd (sind (vect(id, 1)) .* sind (vect(id, 4)) + ... cosd (vect(id, 1)) .* cosd (vect(id, 4)) .* ... - cosd (vect(id, 2) .- vect(id, 5))); + cosd (vect(id, 2) - vect(id, 5))); ## Haversine formula sphd(! id) = 2 * asind (sqrt ( ... - (sind (abs (vect(! id, 1) .- vect(! id, 4)) / 2)) .^ 2 .+ ... + (sind (abs (vect(! id, 1) - vect(! id, 4)) / 2)) .^ 2 + ... cosd (vect(! id, 1)) .* cosd (vect(! id, 4)) .* ... - sind (abs (vect(! id, 2) .- vect(! id, 5)) / 2) .^2)); + sind (abs (vect(! id, 2) - vect(! id, 5)) / 2) .^2)); ## Init tracker for various reasons for issues with intersections idl = zeros (nv, 1); ## Find circle pairs that have no intersections. Case 1: sphd > sum of radii - idl(sphd .- (vect(:, 3) .+ vect(:, 6)) > 0) = 1; + idl(sphd - (vect(:, 3) + vect(:, 6)) > 0) = 1; ## Find circle pairs that have no intersections. Case 2: one circle ## completely lying in the other w/o tangent points - idl(vect(:, 3) .- sphd > vect(:, 6)) = 1; - idl(vect(:, 6) .- sphd > vect(:, 3)) = 1; + idl(vect(:, 3) - sphd > vect(:, 6)) = 1; + idl(vect(:, 6) - sphd > vect(:, 3)) = 1; ## Find circle pairs with coinciding (or antipodal) centers ==> ## Some may have no intersections (one enclosed in the other) [alat, alon] = antipode (vect(:, 4), vect(:, 5)); - idlc = find ((abs (vect(:, 1) .- vect (:, 4)) < 1e-13 & ... - abs (vect(:, 2) .- vect (:, 5)) < 1e-13) | ... - (abs (vect(:, 1) .- alat) < 1e-13 & ... - abs (vect(:, 2) .- alon) < 1e-13)); + idlc = find ((abs (vect(:, 1) - vect (:, 4)) < 1e-13 & ... + abs (vect(:, 2) - vect (:, 5)) < 1e-13) | ... + (abs (vect(:, 1) - alat) < 1e-13 & ... + abs (vect(:, 2) - alon) < 1e-13)); idl(idlc) = 1; ## Some of these may have coinciding circles - idlca = abs (vect(idlc, 3) .- vect(idlc, 6)) < 2 * eps | ... - abs (vect(idlc, 3) - 180 .+ vect(idlc, 6)) < 2 * eps; + idlca = abs (vect(idlc, 3) - vect(idlc, 6)) < 2 * eps | ... + abs (vect(idlc, 3) - 180 + vect(idlc, 6)) < 2 * eps; idl(idlc(idlca)) = 3; ## Set up pointer to valid pairs idv = idl < 1; @@ -255,17 +255,17 @@ function [lat, lon, idl] = scxsc (varargin) r1 = vect(idv, 3); r2 = vect(idv, 6); - a = (cosd (r1) .- q .* cosd (r2)) ./ (1 .- q2); - b = (cosd (r2) .- q .* cosd (r1)) ./ (1 .- q2); + a = (cosd (r1) - q .* cosd (r2)) ./ (1 - q2); + b = (cosd (r2) - q .* cosd (r1)) ./ (1 - q2); - x0 = a .* c1 .+ b .* c2; + x0 = a .* c1 + b .* c2; x02 = dot (x0, x0, 2); ## Dot product cannot be larger than 1 here x02 = sign (x02) .* min (1, abs (x02)); - t = sqrt ((1 .- x02) ./ n2); # <=== No chance n2 can be zero? (div by zero !!!) - p1 = x0 .+ (t .* n); - p2 = x0 .- (t .* n); + t = sqrt ((1 - x02) ./ n2); # <=== No chance n2 can be zero? (div by zero !!!) + p1 = x0 + (t .* n); + p2 = x0 - (t .* n); lat(idv, 1) = atan2d (p1(:, 3), hypot (p1(:, 1), p1(:, 2))); lon(idv, 1) = atan2d (p1(:, 2), p1(:, 1)); |