function   [XP YP]=chamberlin_p(latT,lonT,L0,lat,lon);
%
% function   [XP YP]=chamberlin_p(latT,lonT,L0,lat,lon);
%
% The Chamberlin Trimetric Projection - Analytic Form
%
% The formalism for the code here is due to John P. Snyder (or perhaps someone before him?).
% The code was extracted from Snyder's personal suite of projection routines that were
% coded in QBASIC.  This code was downloaded from this location:
%     http://www.galleryofmapprojections.com/snyder/snyder_full.zip
% Documentation associated with this code explicitly states that the routines are released
% to the public domain without restriction by Paul B. Anderson for John P. Snyder,
% dated October 7, 2007.  Paul collaborated with Snyder, who was the U.S. Geological
% Survey's mapping projection expert from 1978-1997.  Snyder passed away in 1997 at
% the age of 71.
%
% I know nothing of the theory behind this analytic form for this projection.  The code
% here was just transcoded from QBASIC to matlab, and then vectorized as much as possible.  
% The projection is for a sphere, rather than an ellipsoid.
%
% There is no known inverse formalism for this projection, that is, given an (X,Y) there is
% no analytic way to determine (lat,lon).  However, by defining grids (x,y) <-> (lat,lon),
% an inverse can easily be found by interpolation.
%
% Input Parameters:
% (latT,lonT) vectors of length 3 with the defining three points for this projection.
%   These points should be located near the boundaries of the map for an optimal map.
% L0 is the radius of the Earth
% (lat,lon) are the points to project to a plane.  
%   (lat,lon) can be either points, vectors, or arrays.
%   The routine is a little slow because it has to go point-by-point, I believe.
%
% Output Parameters
% (X,Y) the points on the plane associated with (lat,lon)
%
% B. Dushaw - January 2010
%

XP=zeros(size(lon));
YP=zeros(size(lon));

% Get set up - first the three defining points.
% The values K1 and K2 served two purposes in the QBASIC routine.  One was as dimension size for
% arrays, the other was as constants K1=1.0, K2=2.0.  I've set them to be 2, 3 because arrays
% in matlab do not start with 0, i.e., it is 1:K2 (1:3), rather than 0:K2 (0:2).
K1=2; 
K2=3;

% Define some constants.
PI=pi;
DG1=pi/180;
PI2=2*pi;
PIH=PI/2;
T=0.00015;
third=L0/3.0;

%%%%%%%%%%%%%%%  BEGIN SETUP - LOCATIONS OF TRIMETRIC POINTS AND WORKING VARIABLES  %%%%%%%%%%%%%%
CEN=zeros(2,3);
CEN(1,:)=latT;
CEN(2,:)=lonT;

CEN=CEN*DG1;
SCEN=sin(CEN);
CCEN=cos(CEN);

%%%   DISTANCES BETWEEN POINTS ON A UNIT SPHERE, I THINK.
E=zeros(3,1);
E(1)= SCEN(1,1) *SCEN(1,K1)+CCEN(1,1)* CCEN(1,K1)*(CCEN(K1,1)* CCEN(K1,K1)+SCEN(K1,1)* SCEN(K1,K1));
E(K1)=SCEN(1,K1)*SCEN(1,K2)+CCEN(1,K1)*CCEN(1,K2)*(CCEN(K1,K1)*CCEN(K1,K2)+SCEN(K1,K1)*SCEN(K1,K2));
E(K2)=SCEN(1,K2)*SCEN(1,1) +CCEN(1,K2)*CCEN(1,1)* (CCEN(K1,K2)*CCEN(K1,1)+ SCEN(K1,K2)*SCEN(K1,1));

E=PIH-atan(E./sqrt(abs(1.0-E.*E)));

A=zeros(2,1);
A(1)= 0.5*(E(1)*E(1)+E(K2)*E(K2)-E(K1)*E(K1))/(E(1)*E(K2));
A(K1)=0.5*(E(1)*E(1)+E(K1)*E(K1)-E(K2)*E(K2))/(E(1)*E(K1));

I=abs(A)==1;
A(I)=(PIH-PIH*sign(CEN(1,1)-CEN(1,K2)))*ones(size(A(I)));
A(~I)=PIH-atan(A(~I)./sqrt(abs(1.0-A(~I).*A(~I))))*sign(CEN(1,1)-CEN(1,K2));

XY(K1,1)=E(K2)*sin(A(1));
XY(K1,K1)=XY(K1,1);
XY(1,1)=-0.5*E(1);
XY(1,K1)=-XY(1,1);
XY(1,K2)=XY(1,1)+E(K2)*cos(A(1));

H=zeros(3,1);
H(1)= (SCEN(1,K2)-cos(E(K2))*SCEN(1,1))/ (sin(E(K2))*CCEN(1,1));
H(K1)=(SCEN(1,1)- cos(E(1))* SCEN(1,K1))/(sin(E(1))* CCEN(1,K1));
H(K2)=(SCEN(1,K1)-cos(E(K1))*SCEN(1,K2))/(sin(E(K1))*CCEN(1,K2));

I=abs(H)==1;
H(I)=zeros(size(H(I)));
H(~I)=PIH-atan(H(~I)./sqrt(abs(1.0-H(~I).*H(~I))));

%%%%%%%%%%%%%%%  END OF SETUP  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%  NOW PROJECT THE POINTS  %%%%%%%%%%%%%%%%%%%%

NNN=length(lat(:));
lat=lat*DG1;
lon=lon*DG1;
sP=sin(lat); cP=cos(lat);
sL=sin(lon); cL=cos(lon);

RR1= SCEN(1,1)* sP+CCEN(1,1)* cP.*(CCEN(K1,1)* cL+SCEN(K1,1)* sL);
RR2= SCEN(1,K1)*sP+CCEN(1,K1)*cP.*(CCEN(K1,K1)*cL+SCEN(K1,K1)*sL);
RR3= SCEN(1,K2)*sP+CCEN(1,K2)*cP.*(CCEN(K1,K2)*cL+SCEN(K1,K2)*sL);

R=zeros(3,1);

% There does not seem to be a way to vectorize this big loop.
% It must be point-by-point, I think.
for iki=1:NNN

P=lat(iki);
L=lon(iki);

if isnan(P) | isnan(L),
% nothing to do for NaN's
XP(iki)=NaN;
YP(iki)=NaN;
else

P1=sP(iki);
P2=cP(iki);
L1=sL(iki);
L2=cL(iki);

% Distances from the base points to the projection point.
R(1) = RR1(iki);
R(2) = RR2(iki);
R(3) = RR3(iki);

I=R==1;
R(I)=zeros(size(R(I)));
R(~I)=PIH-atan(R(~I)./sqrt(abs(1.0-R(~I).*R(~I))));

if R(K1)==0, 
   A(K2)=0;
   A2=0;
else
  A(K2)=(P1-cos(R(K1))*SCEN(1,K1))/(sin(R(K1))*CCEN(1,K1));
  if abs(A(K2))==1.0,
    A(K2)=0;
  else
    A(K2)=PIH-atan(A(K2)/sqrt(abs(1.0-A(K2)*A(K2))));
  end
  A2=0.5*(E(1)*E(1)+R(K1)*R(K1)-R(1)*R(1))/(E(1)*R(K1));
  if abs(A2)==1.0 ,
    A2=PIH*(1.0-sign(A2));
  else
    A2=PIH-atan(A2/sqrt(abs(1.0-A2*A2)));
  end
end

if abs(L-CEN(K1,K1))>PI, 
   L=L-PI2;
   if abs(L-CEN(K1,K1))>PI, 
     L=L+2*PI2;
   end
end

if (L<CEN(K1,K1)-T) & (A(K2)<H(K1)) | (L>CEN(K1,K1)+T) & (A(K2)<PI-H(K1)) | (abs(L-CEN(K1,K1))<2.0*T) & (P>CEN(1,K1)), 
  A2=-A2;
end

XA= XY(1,K1)-R(K1)*cos(A2);
YA=XY(K1,K1)-R(K1)*sin(A2);

%%%%%%%%%%%%%%%%%%%%  END FIRST INTERSECTION POINT  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if R(K2)==0, 
  A4=0;
  A3=0;
else
  A4=(P1-cos(R(K2))*SCEN(1,K2))/(sin(R(K2))*CCEN(1,K2));
  if abs(A4)==1.0,
     A4=PIH*(1.0-sign(A4));
  else
     A4=PIH-atan(A4/sqrt(abs(1.0-A4*A4)));
  end
  A3=0.5*(E(K1)*E(K1)+R(K2)*R(K2)-R(K1)*R(K1))/(E(K1)*R(K2));
  if abs(A3)==1.0, 
     A3=PIH*(1.0-sign(A3));
  else
     A3=PIH-atan(A3/sqrt(abs(1.0-A3*A3)));
  end
end

if (abs(L-CEN(K1,K2))>PI),
  L=L-PI2;
  if abs(L-CEN(K1,K2))>PI, 
    L=L+2.0*PI2;
  end
end
if (L<CEN(K1,K2)-T) & (A4>PI-H(K2)) | (L>CEN(K1,K2)+T) & (A4>H(K2)) | (abs(L-CEN(K1,K2))<2.0*T) & (P<CEN(1,K2)),
    A3=-A3;
end

XB=XY(1,K2)+R(K2)*cos(A3+A(K1));
YB=         R(K2)*sin(A3+A(K1));

%%%%%%%%%%%%%%%%%%%%  END SECOND INTERSECTION POINT  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if R(1)==0,
   A5=0;
   A1=0;
else
   A5=(P1-cos(R(1))*SCEN(1,1))/(sin(R(1))*CCEN(1,1));
   if abs(A5)==1.0,
      A5=PIH*(1.0-sign(A5));
   else
      A5=PIH-atan(A5/sqrt(abs(1.0-A5*A5)));
   end
   A1=0.5*(E(K2)*E(K2)+R(1)*R(1)-R(K2)*R(K2))/(E(K2)*R(1));
   if abs(A1)==1.0, 
     A1=PIH*(1.0-sign(A1));
   else
     A1=PIH-atan(A1/sqrt(abs(1.0-A1*A1)));
   end
end

if abs(L-CEN(K1,1))>PI, 
   L=L-PI2;
   if abs(L-CEN(K1,1))>PI, 
      L=L+2.0*PI2;
   end
end

if (L<CEN(K1,1)-T) & (A5>PI-H(1)) | (L>CEN(K1,1)+T) & (A5>H(1)) | (abs(L-CEN(K1,1))<2.0*T) & (P<CEN(1,1)),
     A1=-A1;
end

X=XY(1,1) +R(1)*cos(A(1)-A1);
Y=XY(K1,1)-R(1)*sin(A(1)-A1);

%%%%%%%%%%%%%%%%%%%%  END THIRD INTERSECTION POINT  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%  NOW AVERAGE THE THREE POINTS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X=third*(XA+XB+X);
Y=third*(YA+YB+Y-CEN(1,1)-CEN(1,K1)-CEN(1,K2));

XP(iki)=X;
YP(iki)=Y;
end   % end NaN if
end   % end iki
