function [x,y]=ll2xy(lat, lon, ref);
%ll2xy.m
%
% [x,y]=ll2xy(lat, lon, ref);
%
%Converts latitude and longitude coordinates to x,y in meters relative
%to some reference position.  If no ref=[lat,lon]
%is provided, the position of Boston is used (42.35 N, -71.0264).



if (nargin<3)
	ref=[42.3493, -71.0264];       %default ref is Boston
end;
n=length(lat);

R= 6378.155*1000;    %Radius at equator is 6378.155 km

R2=R*cos(ref(1)*pi/180);

x=(lon-ref(2))/360 * 2*pi*R2;  	%X
y=(lat-ref(1))/360 * 2*pi*R;		%Y
