function NDEL = Ndiff(fname,argi,vcell)
% computes the numerical gradient of the tensor function 'fname', with respect to argument #argi
% the "point" at which it's evaluated is specified by the tensors in vcell
% syntax:
% NDEL = Ndiff('fname',argi,{x1,x2,...,x_n});
% output:
% if Y = f(xi1,xi2,...,xi_NX), then
% NDEL(yi1,yi2,...,yi_orderY,x_argi_1,x_argi_2,...,x_argi_orderX) = 
% dY(yi1,yi2,...,yi_orderY)/dX_argi(x_argi_1,x_argi_2,...,x_argi_orderX)
% for example:
% NLTA = Ndiff('myLL',1,{ThA,ThB,U1})
% NLTB = Ndiff('myLL',2,{ThA,ThB,U1})

dx = 1e-8;
%dx = 1e-12;

vstr = cell_ind_str('vcell',length(vcell));
% vstr = 'vcell{1},vcell{2},...,vcell{narg}'

y0 = eval([fname,'(',vstr,')']);
X0 = vcell{argi};
dimx = mydim(X0);
dimy = mydim(y0);
for ix = 1:prod(dimx)
  cx = vec2str(ind2coord(ix,dimx));
  % cx = 'i1,i2,...,i_n', which are the coordinates of cx
  X = X0;
  eval(['X(',cx,') = X(',cx,') + dx;']);
  vcell{argi} = X;
  y1 = eval([fname,'(',vstr,')']);
  Dy = y1 - y0;
  for iy = 1:prod(dimy)
    cy = vec2str(ind2coord(iy,dimy));
    dy = eval(['Dy(',cy,')']);
    del = dy/dx;
    cdel = concatc(cy,cx);
    eval(['NDEL(',cdel,') = del;']);
  end
end

return


function vstr = cell_ind_str(cname,dim)
vstr = '';
for j=1:dim
  vstr = [vstr,',',cname,'{',num2str(j),'}'];
end
vstr = vstr(2:end);
return


function vstr = vec2str(vec)
vstr = '';
for j=1:length(vec)
  vstr = [vstr,',',num2str(vec(j))];
end
vstr = vstr(2:end);
if isempty(vstr)
  vstr = '1';
end
return

function c = ind2coord(ind,dim)
ind = ind-1;
c = zeros(size(dim));
for ij = 1:length(dim)
  c(ij) = mod(ind,dim(ij));
  ind = (ind - c(ij))/dim(ij);
end
c = c + 1;
return


function ind = coord2ind(c,dim)
ind = 0;
c = c - 1;
f = 1;
for ij = 1:length(dim)
  ind = ind + c(ij)*f;
  f = f*dim(ij);
end
ind = ind+1;
return


function d = mydim(x)
d = [];
for j=1:length(size(x))
  if (size(x,j)>1)
    d(end+1) = size(x,j);
  end
end
return

function c3 = concatc(c1,c2)
if isempty(c1)
  c3=c2;
else
  if isempty(c2)
    c3=c1;
  else
    c3 = [c1,',',c2];
  end
end

return


function y = myrat(K,phi)
%yF = fnorm(K);

yphi = phi*K;

yPsi = psinorm(K);

y = yphi/(yPsi+realmin);

return