function [ys,d1s,d2s] = ste_tan0( vs,grid )
% [ys,d1s,d2s] = ste_tan0( vs,grid )
%  The Stesjkal-Tanner (ST) equation for diffusion imaging
%   y_i = exp(-(v' x_i)^2)
%  Kernel function used for DWI mixture models
% 
% Inputs: 
%  - vs, 3xK grid of diffusion-weighted fiber directions
%  - grid, 3xn grid of measurements, x_1,...,x_n
%
% Outputs:
%  - ys, nxK matrix: the ith column is the ST kernel for the ith fiver in vs
%  - d1s, 1x3 cell:
%     d1s{1}, nxK matrix: derivative of ST kernel wrt x-coordinate of diffusion-weighted direction
%     d1s{2}, nxK matrix: derivative of ST kernel wrt y-coordinate of diffusion-weighted direction
%     d1s{3}, nxK matrix: derivative of ST kernel wrt z-coordinate of diffusion-weighted direction
%  - d2s, 3x3 cell:
%     d1s{1,1}, nxK matrix: d^2/dx^2 of ST kernel wrt diffusion-weighted direction
%     d1s{1,2}, nxK matrix: d^2/dxdy of ST kernel wrt diffusion-weighted direction
%     ...
%     d1s{3,3}, nxK matrix: d^2/dz^2 of ST kernel wrt diffusion-weighted direction

n=size(grid,2);
pp=size(vs,2);
raw = grid'*vs;
ys = exp(-raw.^2);
d1s = cell(1,3);
coords = cell(1,3);
for ii=1:3;
    coords{ii} = repmat(grid(ii,:)',1,pp);
end
for ii=1:3;
    d1s{ii} = -2.*coords{ii}.*raw.*ys;
end
d2s = cell(3);
for ii=1:3;
    for jj=1:3;
        d2s{ii,jj} = (4.*coords{ii}.*coords{jj}.*raw.^2 -2.*coords{ii}.*coords{jj}).*ys;
    end
end
end

