1function [x, y, azi, rk] = cassini_fwd(lat0, lon0, lat, lon, ellipsoid)
2%CASSINI_FWD  Forward Cassini-Soldner projection
3%
4%   [x, y] = CASSINI_FWD(lat0, lon0, lat, lon)
5%   [x, y, azi, rk] = CASSINI_FWD(lat0, lon0, lat, lon, ellipsoid)
6%
7%   performs the forward Cassini-Soldner projection of points (lat,lon) to
8%   (x,y) using (lat0,lon0) as the center of projection.  These input
9%   arguments can be scalars or arrays of equal size.  The ellipsoid vector
10%   is of the form [a, e], where a is the equatorial radius in meters, e is
11%   the eccentricity.  If ellipsoid is omitted, the WGS84 ellipsoid (more
12%   precisely, the value returned by defaultellipsoid) is used.  projdoc
13%   defines the projection and gives the restrictions on the allowed ranges
14%   of the arguments.  The inverse projection is given by cassini_inv.
15%
16%   azi and rk give metric properties of the projection at (lat,lon); azi
17%   is the azimuth of the easting (x) direction and rk is the reciprocal of
18%   the northing (y) scale.  The scale in the easting direction is 1.
19%
20%   lat0, lon0, lat, lon, azi are in degrees.  The projected coordinates x,
21%   y are in meters (more precisely the units used for the equatorial
22%   radius).  rk is dimensionless.
23%
24%   See also PROJDOC, CASSINI_INV, GEODDISTANCE, DEFAULTELLIPSOID,
25%     FLAT2ECC.
26
27% Copyright (c) Charles Karney (2012-2015) <charles@karney.com>.
28
29  narginchk(4, 5)
30  if nargin < 5, ellipsoid = defaultellipsoid; end
31  try
32    Z = zeros(size(lat0 + lon0 + lat + lon));
33  catch
34    error('lat0, lon0, lat, lon have incompatible sizes')
35  end
36  if length(ellipsoid(:)) ~= 2
37    error('ellipsoid must be a vector of size 2')
38  end
39
40  degree = pi/180;
41  f = ecc2flat(ellipsoid(2));
42  dlon = AngDiff(lon0, lon) + Z;
43  [s12, azi1, azi2, ~, ~, ~, ~, sig12] = ...
44      geoddistance(lat, -abs(dlon), lat, abs(dlon), ellipsoid);
45  sig12 = 0.5 * sig12;
46  s12 = 0.5 * s12;
47  c = s12 == 0;
48  da = AngDiff(azi1, azi2)/2;
49  s = abs(dlon) <= 90;
50  azi1(c & s) = 90 - da(c & s);
51  azi2(c & s) = 90 + da(c & s);
52  s = ~s;
53  azi1(c & s) = -90 - da(c & s);
54  azi2(c & s) = -90 + da(c & s);
55  c = dlon < 0;
56  azi2(c) = azi1(c);
57  s12(c) = -s12(c);
58  sig12(c) = -sig12(c);
59  x = s12;
60  azi = AngNormalize(azi2);
61  [~, ~, ~, ~, ~, ~, rk] = ...
62      geodreckon(lat, dlon, -sig12, azi, ellipsoid, 1);
63  [sbet, cbet] = sincosdx(lat);
64  [sbet, cbet] = norm2((1-f) * sbet, cbet);
65  [sbet0, cbet0] = sincosdx(lat0);
66  [sbet0, cbet0] = norm2((1-f) * sbet0, cbet0);
67  [salp, calp] = sincosdx(azi);
68  salp0 = salp .* cbet;
69  calp0 = hypot(calp, salp .* sbet);
70  sbet1 = calp0;
71  c = lat + Z >= 0;
72  sbet1(~c) = -sbet1(~c);
73  cbet1 = abs(salp0);
74  c = abs(dlon) <= 90;
75  cbet1(~c) = -cbet1(~c);
76  sbet01 = sbet1 .* cbet0 - cbet1 .* sbet0;
77  cbet01 = cbet1 .* cbet0 + sbet1 .* sbet0;
78  sig01 = atan2(sbet01, cbet01) / degree;
79  [~, ~, ~, ~, ~, ~, ~, y] = geodreckon(lat0, 0, sig01, 0, ellipsoid, 1);
80end
81