clear all
clf
% This test case is the same as that from "An Album of Map Projections" by
% Snyder and Voxland, USGS.
%
% => EDIT AT BOTTOM TO INCLUDE YOUR OWN COASTLINE DATA 
%
% B. Dushaw - December 2009

Ln=(-200:10:25);
Lt=0:10:89;

LN=ones(length(Lt),1)*Ln;
LT=Lt'*ones(1,length(Ln));

% Define the three points for the "trimetric"
latA=55;
latB=55;
latC=10;

lonA=-150;
lonB=-35;
lonC=-92.5;

% Define the grid interval to return
DLT=5;
DLN=5;

minlat=0;
minlon=-200;
maxlat=89;
maxlon=30;

[Xp Yp LAT LON] = chamberlin(latA,lonA,latB,lonB,latC,lonC,DLT,DLN,minlat,minlon,maxlat,maxlon);
% Just returns the grid of the projection.  Interpolate to desired points manually.

% plot longitude lines
for k=1:length(Lt),
  ln=Ln(1):1:Ln(end);
  lt=Lt(k)*ones(size(ln));
  X=interp2(LON,LAT,Xp,ln,lt);
  Y=interp2(LON,LAT,Yp,ln,lt);
  plot(X,Y,'k')
  hold on
end
% plot latitude lines
for k=1:length(Ln),
  lt=Lt(1):1:Lt(end);
  ln=Ln(k)*ones(size(lt));
  X=interp2(LON,LAT,Xp,ln,lt);
  Y=interp2(LON,LAT,Yp,ln,lt);
  plot(X,Y,'k')
  hold on
end


if 1==2,
% plot coastlines
load high_res.mat
I=latC==999;
lt=latC/1000;
ln=lonC/1000;
lt(I)=NaN*lt(I);
ln(I)=NaN*ln(I);

minlat=min(Lt);
maxlat=max(Lt);
minlon=min(Ln);
maxlon=max(Ln);

%I=(-360 < ln) & (ln < -250);
I=(ln < -250);
ln(I)=ln(I)+360;

I1=(minlon < ln) & (ln < maxlon);
I2=(minlat < lt) & (lt < maxlat);
I=~(I1 & I2);

lt(I)=[];
ln(I)=[];
I=abs(diff(ln))>0.5;
ln(I)=1000*ln(I);
I=abs(diff(lt))>0.5;
lt(I)=1000*lt(I);


xc=interp2(LON,LAT,Xp,ln,lt);
yc=interp2(LON,LAT,Yp,ln,lt);

plot(xc,yc,'b','LineWidth',1)
end

axis([-2000000 8300000 -6500000 3700000])
axis('equal')

