function   [X Y LAT LON]=chamberlin(latA,lonA,latB,lonB,latC,lonC,DLt,DLn,minlat,minlon,maxlat,maxlon);
% function [X Y LAT LON]=chamberlin(latA,lonA,latB,lonB,latC,lonC,DLt,DLn,minlat,minlon,maxlat,maxlon);
%
% Chamberlin trimetric projection
%
% Brute-force projection using the m_map's/Pawlowicz's m_idist routine and 2-D linear interpolation.
% Through m_idist, this routine uses the WGS'84 ellipsoid by default.  Modify the calls to this
% routine to specify a different ellipsoid.
%
% We will call this particular projection the Chamberlin-Dushaw-Pawlowicz projection, or the
% CDP projection, for short.
%
% Only returns the corresponding grids of (lat,lon) <-> (x,y) for subsequent user interpolation.
% A strictly analytic solution does not seem possible with the design of this code.
%
% Input:
% (latA,lonA), (latB,lonB), (latC, lonC) of defining triad
% DLt,DLn is the desired lat,lon increment for the grid.  These should be ca. 1 degree.
% These values determine the interpolation accuracy for the projection - smaller values minimize 
% the interpolation error, but the calculation takes longer and you'll need more memory.
% (max,min)(lat,lon) define the rectangular grid area.  These must lie outside 
% the points of interest.  The calculated grid extends only to these values.
%
% The map is aligned so that points A and B are zonal with Point A at the origin.
% The lines between three points A, B, C are almost, but not quite, geodesics.
%
% Triangle          ---C---
% and               | / \ |
% Grid              |/   \|
% Domain            A_____B
%
% The procedure is to first calculate a grid projection, and then for the user
% to interpolate the (lat,lon) points using that grid.
%
% After calling this routine, 
%  x=interp2(LON,LAT,X,lon,lat);
%  y=interp2(LON,LAT,Y,lon,lat);
% will give you the projection (x,y) for (lat,lon).   
% Use similar calls for the inverse.
%
% => There are issues in this projection when the points are near the 
%    geodesics of the defining triad.  This interpolation routine is best, 
%    otherwise singleton points near these geodesics may be in error.
%
% Output:
% (X,Y) and (LAT,LON) are corresponding grids for projection and inverse.
%
% N.B. You can test for dislocations (errors) in your grid output by plotting 
% its 2nd derivative:
%   plot(diff(diff(Xp')),diff(diff(Yp')),'r'); 
%   hold on; 
%   plot(diff(diff(Xp'))',diff(diff(Yp'))','r')
% This figure should be smooth and well behaved - no glitches allowed!
%
% Derived from the description given in:
% Flattening the Earth: Two Thousand Years of Map Projections (Paperback)
% by John P. Snyder, University Of Chicago Press (December 5, 1997)
% 384 pp. ISBN-10: 0226767477  ISBN-13: 978-0226767475
%
% See also:  
% An album of map projections, by John P. Snyder and Philip M. Voxland
% USGS Professional Paper, Report 1453, 1989, 249 pp.
% URL: http://pubs.er.usgs.gov/usgspubs/pp/pp1453
%
% Brian Dushaw - December 2009-January 2010
% Applied Physics Laboratory, University of Washington
%

% First define the grid...
f=minlon:DLn:maxlon;
g=minlat:DLt:maxlat;
lon=ones(length(g),1)*f;
lat=g'*ones(1,length(f));
clear f g

[x1 y1 x2 y2 x3 y3 X1 X2 X3] = getstd(latA,lonA,latB,lonB,latC,lonC);
% Gets coordinates of the vertices and the lengths of the triangle sides

[x y xa ya] = getpts(lat,lon,x1,y1,x2,y2,x3,y3,X1,X2,X3,latA,lonA,latB,lonB,latC,lonC);
% (x,y) are the points
% (xa,ya) are possible alternate points; see below.

%%%%%%%%%%%%%%%%%%%%  FIX DISLOCATED POINTS ON THE GRID  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now fix the dislocated points along the triangle geodesics.
% These set of loops don't take very long to calculate - the main action is in the call to getpts.
% Takes about 1/3 the time of getpts.
% If you get dislocated points, try increasing the maximum size of GL - 8? 12? 16?
[n m]=size(x);
for GL=1:8,  
% Routine to identify which point (original or alternate) best fits the grid...
% Do the loops several times... First remove the gross errors, then do the fine tuning...
% Have to go point-by-point, but only need to look at those points for which (xa,ya) are not NaN's.
if GL==1 | GL==8 | GL==12 | GL==14, iq1=1:n;,    iq2=m:-1:1;, end
if GL==2 | GL==7 | GL==11 | GL==16, iq1=1:n;,    iq2=1:m;,    end
if GL==3 | GL==6 | GL==10 | GL==13, iq1=n:-1:1;, iq2=1:m;,    end
if GL==4 | GL==5 | GL==9  | GL==15, iq1=n:-1:1;, iq2=m:-1:1;, end
for indy=iq1
for indx=iq2
    hx=xa(indy,indx); 
    if ~isnan(hx),
      hy=ya(indy,indx);
      gx=x(indy,indx); 
      gy=y(indy,indx); 
      d1=0;
      d2=0;

      % d1 measures how well the recommended point fit the grid.
      % d2 measures how well the alternate point fits the grid.
      % Distances along X
      if indx>1 & indx < m, 
        d1=d1+(sqrt((x(indy,indx-1) - gx)^2 + (y(indy,indx-1)-gy)^2) - sqrt((gx-x(indy,indx+1))^2 + (gy-y(indy,indx+1))^2))^2;
        d2=d2+(sqrt((x(indy,indx-1) - hx)^2 + (y(indy,indx-1)-hy)^2) - sqrt((hx-x(indy,indx+1))^2 + (hy-y(indy,indx+1))^2))^2;
      elseif indx==1,
        % on the left edge
        d1=d1+(sqrt((x(indy,indx+1) - x(indy,indx+2))^2 + (y(indy,indx+1)-y(indy,indx+2))^2) - sqrt((gx-x(indy,indx+1))^2 + (gy-y(indy,indx+1))^2))^2;
        d2=d2+(sqrt((x(indy,indx+1) - x(indy,indx+2))^2 + (y(indy,indx+1)-y(indy,indx+2))^2) - sqrt((hx-x(indy,indx+1))^2 + (hy-y(indy,indx+1))^2))^2;
      elseif indx==m
        % on the right edge
        d1=d1+(sqrt((x(indy,indx-1) - gx)^2 + (y(indy,indx-1)-gy)^2) - sqrt((x(indy,indx-2)-x(indy,indx-1))^2 + (y(indy,indx-2)-y(indy,indx-1))^2))^2;
        d2=d2+(sqrt((x(indy,indx-1) - hx)^2 + (y(indy,indx-1)-hy)^2) - sqrt((x(indy,indx-2)-x(indy,indx-1))^2 + (y(indy,indx-2)-y(indy,indx-1))^2))^2;
      end
      
      % Distances along Y
      if indy>1 & indy < n, 
        d1=d1+(sqrt((x(indy-1,indx) - gx)^2 + (y(indy-1,indx)-gy)^2) - sqrt((gx-x(indy+1,indx))^2 + (gy-y(indy+1,indx))^2))^2;
        d2=d2+(sqrt((x(indy-1,indx) - hx)^2 + (y(indy-1,indx)-hy)^2) - sqrt((hx-x(indy+1,indx))^2 + (hy-y(indy+1,indx))^2))^2;
      elseif indy==1,
        % on the bottom edge
        d1=d1+(sqrt((x(indy+1,indx) - x(indy+2,indx))^2 + (y(indy+1,indx)-y(indy+2,indx))^2) - sqrt((gx-x(indy+1,indx))^2 + (gy-y(indy+1,indx))^2))^2;
        d2=d2+(sqrt((x(indy+1,indx) - x(indy+2,indx))^2 + (y(indy+1,indx)-y(indy+2,indx))^2) - sqrt((hx-x(indy+1,indx))^2 + (hy-y(indy+1,indx))^2))^2;
      elseif indy==n,
        % on the top edge
        d1=d1+(sqrt((x(indy-1,indx) - gx)^2 + (y(indy-1,indx)-gy)^2) - sqrt((x(indy-2,indx)-x(indy-1,indx))^2 + (y(indy-2,indx)-y(indy-1,indx))^2))^2;
        d2=d2+(sqrt((x(indy-1,indx) - hx)^2 + (y(indy-1,indx)-hy)^2) - sqrt((x(indy-2,indx)-x(indy-1,indx))^2 + (y(indy-2,indx)-y(indy-1,indx))^2))^2;
      end

      if d2 < d1,
% The alternate point fits the grid better, it seems.
        x(indy,indx)=hx;
        y(indy,indx)=hy;
% Save the original point in the event of erroneous corrections, which happens.
        xa(indy,indx)=gx;
        ya(indy,indx)=gy;
      end
    end  % end for "if"

end   % end for indx
end   % end for indy

end   % end for GL

%%%%%%%%%%%%%%%%%%%%  END OF FIX OF DISLOCATED POINTS ON THE GRID  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plot the grid - debugging for left over dislocated points
if 1==2
  plot(x,y,'r')
  hold on
  plot(x',y','r')
  plot(xa,ya,'b.')
end

X=x;
Y=y;
LAT=lat;
LON=lon;
%%%%%%%%%%%%%%%%%%%%  END OF MAIN ROUTINE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gets the projection - possible errors along geodesic triangle legs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x y xa ya] = getpts(lat,lon,x1,y1,x2,y2,x3,y3,X1,X2,X3,latA,lonA,latB,lonB,latC,lonC);

x=zeros(size(lat));
y=x;
xa=x;
ya=y;

% Calculate the three distances from vertices to (lat,lon) point in question.
% Most of the computation time for this routine is in these distances.
% Vectorized it is much faster!
[z1 a1 b1]=m_idist(lonA,latA,lon,lat);
[z2 a2 b2]=m_idist(lonB,latB,lon,lat);
[z3 a3 b3]=m_idist(lonC,latC,lon,lat);

for k=1:length(lat(:)),

  if lat(k)==latA & lon(k)==lonA,
     x(k)=x1; y(k)=y1;
  elseif lat(k)==latB & lon(k)==lonB,
     x(k)=x2; y(k)=y2;
  elseif lat(k)==latC & lon(k)==lonC,
     x(k)=x3; y(k)=y3;
  else
  
  z11=z1(k);
  z22=z2(k);
  z33=z3(k);

    if k==-11,
%    plots circle arcs for debugging
       a0=-100:100;
       a=-(a1-90)+a0; bx=x1+z1*cos(a*pi/180); by=y1+z1*sin(a*pi/180); plot(bx,by)
       hold on
       a=-(a2-90)+a0; bx=x2+z2*cos(a*pi/180); by=y2+z2*sin(a*pi/180); plot(bx,by)
       a=-(a3-90)+a0; bx=x3+z3*cos(a*pi/180); by=y3+z3*sin(a*pi/180); plot(bx,by)
    end

    % Solve for the point of intersection of the three arcs - its not a point of
    % intersection actually, so we take the average of the three points that lie
    % closest together.
    [x(k) y(k) xa(k) ya(k)] = solvetri(X1,X2,X3,x1,x2,x3,y1,y2,y3,z11,z22,z33);
   % x,y are the likely points; xa,ya are the alternate points.
  end

end

return

function [xs ys xa ya] = solvetri(X1,X2,X3,x1,x2,x3,y1,y2,y3,z1,z2,z3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve for intersection points of three circles.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% c.f., http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/
% Each pair of circles has two possible solutions, so need to 
% find those three solutions that are closest together.
% 
% When the points are near one of the geodesics between the defining triads,
% these two points of intersection are very close together.  This makes
% problems, since the correct point is not always found by the minimum.
% So we also return the alternate, possible points.
% Patched up in main calling routine, hopefully.
%
% (xs,ys) are the likely solutions
% (xa,ya) are the alternate solutions - (xs,ys) is often the wrong solution 
% along the geodesics of the triad.  The intersection of circles is not
% well-defined along these lines.

p3=zeros(6,2);
third=1/3;

% A,B
iX1=1/X1;
aa1=(z1^2 - z2^2 + X1^2)*0.5*iX1;
h1=iX1*sqrt(z1^2-aa1^2);
aa1=iX1*aa1;
p2x1=    aa1*(x2-x1); 
p2y1=    aa1*(y2-y1);
b1=h1*(y2-y1);
b2=h1*(x2-x1);
% intersection #1
p3(1,1)=p2x1 + b1;   p3(1,2)=p2y1 - b2;
% intersection #2
p3(2,1)=p2x1 - b1;   p3(2,2)=p2y1 + b2;

% A,C
iX2=1/X2;
aa2=(z1^2 - z3^2 + X2^2)*0.5*iX2;
h2=iX2*sqrt(z1^2-aa2^2);
aa2=iX2*aa2;
p2x2=    aa2*(x3-x1); 
p2y2=    aa2*(y3-y1);
b1=h2*(y3-y1);
b2=h2*(x3-x1);
% intersection #3
p3(3,1)=p2x2 + b1; p3(3,2)=p2y2 - b2; 
% intersection #4
p3(4,1)=p2x2 - b1; p3(4,2)=p2y2 + b2; 

% B,C
iX3=1/X3;
aa3=(z2^2 - z3^2 + X3^2)*0.5*iX3;
h3=iX3*sqrt(z2^2-aa3^2);
aa3=iX3*aa3;
p2x3=x2 + aa3*(x3-x2); 
p2y3=y2 + aa3*(y3-y2);
b1=h3*(y3-y2);
b2=h3*(x3-x2);
% intersection #5
p3(5,1)=p2x3 + b1; p3(5,2)=p2y3 - b2;
% intersection #6
p3(6,1)=p2x3 - b1; p3(6,2)=p2y3 + b2; 

% Now need to find and average the three points that are closest together

I1=[1 1 1 1 2 2 2 2 ];
I2=[3 3 4 4 3 3 4 4 ];
I3=[5 6 5 6 5 6 5 6 ]; 

% d is the length of the perimeter of the triangle, sort of - 
% the sum of the squares of the side lengths.
d=(p3(I1,1)-p3(I2,1)).^2 + (p3(I1,2)-p3(I2,2)).^2 + ...   
  (p3(I1,1)-p3(I3,1)).^2 + (p3(I1,2)-p3(I3,2)).^2 + ...
  (p3(I2,1)-p3(I3,1)).^2 + (p3(I2,2)-p3(I3,2)).^2;

[dmin id]=min(d);
xs1=p3(I1(id),1); ys1=p3(I1(id),2); 
xs2=p3(I2(id),1); ys2=p3(I2(id),2);
xs3=p3(I3(id),1); ys3=p3(I3(id),2);

xs=third*(xs1+xs2+xs3);
ys=third*(ys1+ys2+ys3);

% calculate the alternate positions, if any.
xa=NaN;
ya=NaN;
iq=d<(mean(d)/2);       % Possible bug in the making here...but only want to consider points
                        % that are close together, for which d is small.
                        % The alternate points are all near the base lines of the trimetric triangle.

if sum(iq) > 1,           
  d(id)=max(d);           % Neutralize the smallest value, to get at the next smallest value.
  [dmin id]=min(d);     % Get the next minimum value as the alternate
  xs1=p3(I1(id),1); ys1=p3(I1(id),2); 
  xs2=p3(I2(id),1); ys2=p3(I2(id),2);
  xs3=p3(I3(id),1); ys3=p3(I3(id),2);

  xa=third*(xs1+xs2+xs3);
  ya=third*(ys1+ys2+ys3);
end

return

function [x1 y1 x2 y2 x3 y3 X1 X2 X3] = getstd(latA,lonA,latB,lonB,latC,lonC);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DEFINITION OF STANDARD POINTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X1 = m_idist(lonA,latA,lonB,latB);
X2 = m_idist(lonA,latA,lonC,latC);
X3 = m_idist(lonB,latB,lonC,latC);

% (latA,lonA) is the origin.
x1=0;
y1=0;

% (latB,lonB) is to the right of (latA,lonA);
x2=X1;
y2=0;

% Use cosine rule to get the triangle angles.
cosA=(X3^2 - X1^2 - X2^2)/(-2*X1*X2);
cosB=(X2^2 - X1^2 - X3^2)/(-2*X1*X3);
cosC=(X1^2 - X2^2 - X3^2)/(-2*X2*X3);

% Get the location of third point.
x3=X2*cosA;
% Height of triangle using the formulas for Area
s=0.5*(X1 + X2 + X3);
Area=sqrt(s*(s-X1)*(s-X2)*(s-X3));
y3=2*Area/X1;     % A-B of length X1 is the horizontal base.

if latC > latA,
   y3=abs(y3);
else
   y3=-abs(y3);
end

return


