function SS = gencov(sm,cov_0,cov_e,cov_rhopow,cov_usegauss,cov_gauss_width,rgrid,T,varargin)
% function SS = gencov(sm,cov_0,cov_e,rgrid,T,varargin)
% generates covariance matrix representing smooth profiles

%%
smr = exp(-rgrid.^2/sm.^2); % smoothing vs rho
if cov_usegauss
  gauss_shape=exp(-(rgrid.^2)./(cov_gauss_width^2)).*(1-rgrid.^2); % gauss + force edge to zero withh power scal
  lco = sqrt(cov_e) + (sqrt(cov_0)-sqrt(cov_e))*gauss_shape;
else
  lco = sqrt(cov_e) + (sqrt(cov_0)-sqrt(cov_e))*(1-rgrid.^cov_rhopow); % lower covariance close to edge
end
STo = diag(lco)*toeplitz(smr)*diag(lco); % Covariance matrix after transformation

if ~isempty(varargin)
  % extra constraint for psi_e, special case.
  cov_psi_edge = varargin{1};
  ST = blkdiag(STo,cov_psi_edge);
else
  ST = STo;
end
%%
pT = pinv(T);
SS = pT*ST*pT';

return