function [a,r] = coval_yxt(xg,xd,lat_s,lon_s,t_s)
% 3-D lat,lon,time covariance function for mapping altimeter data
%
% a = coval_yxt(xg,xd,lat_s,lon_s,t_s)
%
% Input:
%      xg      m*3 model grid: columns are lat, long, time
%      xd      n*3 data  grid: columns are lat, long, time
%      lat_s   horizontal latitude scale  (degrees) 
%      lon_s   horizontal longitude scale (degrees)
%      t_s     time scale (days)
%
% Output:
%      a       m*n unscaled covariance function 
%
% John Wilkin 17 Jul 97
% modified from N.Bindoff/A.Walker versions

% separate the three coordinates into column vectors for grid and row
% vectors for the data 
xg_lat = xg(:,1);
xd_lat = xd(:,1)';
xg_lon = xg(:,2);
xd_lon = xd(:,2)';
xg_t   = xg(:,3);
xd_t   = xd(:,3)';

m = length(xg_lat); % number of grid points
n = length(xd_lat); % number of data points

% compute distance metrics in space (lat/lon) and time, normalizing by the
% input covariance scales 

% latitude/longitude term
rsq = ((xg_lat(:,ones([1 n]))-xd_lat(ones([m 1]),:))/lat_s).^2 + ...
    ((xg_lon(:,ones([1 n]))-xd_lon(ones([m 1]),:))/lon_s).^2;
r = sqrt(rsq);

% time term
t = abs(xg_t(:,ones([1 n]))-xd_t(ones([m 1]),:))/t_s;

% evaluate unscaled covariance function
% functional form is the altimeter covariance function used by De Mey
a = exp(-r).*(1+r+rsq/6-(r.*rsq)/6).*exp(-t.^2);