function [grid,ergrid,X,Y]=objmap_latlon(depth,dist,data,xcor,ycor,err,phi,nx,ny);

% To plot a depth vs distance objective map
% 
% Input data and parameters:
% depth, vector with depths in m
% distance, vector in Km
lat=depth; % y-axis
lon=dist; % x-axis

% xcor and ycor; % decorrelation scales in km
% xcor=50; %input('Input x correlation length ');
% ycor=0.05; %input('Input y correlation length ');

% Error energy
% err=0.23; % input('Input error energy ')';

% phi=angle; % input('Input angle '); % "angle of the front relative to x-axis
% -45 for x-axis oriented /, 45 for x-axis oriented \.

% grid size
% nx=50; ny=100; 

%% %  Modified version of Dan Rudnick's code

% Matlab script to make an objective map using an objective function fit
% to approximate the mean.
%


%% Modified constants 
%Some constants.  Take care in changing these.
% nx=90; ny=100;
alon1=min(min(lon)); alon2=max(max(lon));
alat1=min(min(lat)); alat2=max(max(lat));
% tokm=111.324; tokmt=tokm*cos(pi/180*(alat1+alat2)/2);
xl=0.0; xr=(alon2-alon1); % in km already; % tokmt*(alon2-alon1);
yl=0.0; yr=(alat2-alat1)*10^(-3); % converting m to km; % tokm*(alat2-alat1);

xg=linspace(xl,xr,nx);
yg=linspace(yl,yr,ny);
nfunc=3;

%% 
% Ask for level, correlation lengths, angle, and error.
level=1; %input('Input level (1-50) ');
% xcor=input('Input x correlation length ');
% ycor=input('Input y correlation length ');
% phi=input('Input angle '); % "angle of the front relative to x-axis
% err=0.23; % input('Input error energy ')';
xcor2=xcor*xcor;
ycor2=ycor*ycor;
cosphi=cos(phi/180*pi);
sinphi=sin(phi/180*pi);


%Convert distances to km, and put data in vector d.
ii=find(~isnan(data(:,level)));
x=(lon(ii)-alon1); % already in Km % %tokmt*(lon(ii)-alon1);
y=(lat(ii)-alat1)*10^(-3); % from m to km; % %tokm*(lat(ii)-alat1); % depth*10^(-3); % 
d=data(ii,level);
nfact=length(d);
disp(sprintf('%5i data',nfact))

%Calculate the data-data covariance matrix E.
[X1,X2]=meshgrid(x);
[Y1,Y2]=meshgrid(y);
dx=X1-X2;
dy=Y1-Y2;
dxr=dx*cosphi-dy*sinphi;
dyr=dx*sinphi+dy*cosphi;
E=exp(-dxr.^2/xcor2-dyr.^2/ycor2);
E=E+err*eye(nfact);

%Calculate the data-grid points covariance matrix C.
C=zeros(nfact,nx*ny);
[Y,X]=meshgrid(yg,xg);
for n=1:nfact
   dx=X(:)'-x(n);
   dy=Y(:)'-y(n);
   dxr=dx*cosphi-dy*sinphi;
   dyr=dx*sinphi+dy*cosphi;
   C(n,:)=exp(-dxr.^2/xcor2-dyr.^2/ycor2);
end

%Calculate the matrix A to turn data into a grid.
m=1:nfunc;
F=fxy(m,x,y);
EF=E\F;
W=EF/(F'*EF);
EC=E\C;
f=fxy(m,X(:),Y(:))';
A=W*f+(eye(nfact)-W*F')*EC;

%Now calculate the grid.
grid=zeros(nx,ny);
grid(:)=A'*d;

%Calculate error grid.
ergrid=zeros(nx,ny);
EC=E*A;
for n=1:nx*ny   %this way saves memory
   ergrid(n)=A(:,n)'*(EC(:,n)-2*C(:,n))+1;
end

% to output grid
long=linspace(alon1,alon2,nx);
latg=linspace(alat1,alat2,ny);
X=long;
Y=latg;
