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