function [Qx,Qd,Sx0,Sd0] = KFcovariances(KFm,KFp,model) 
%#codegen
% Generate covariance matrices from parameters

x_cov_j0   = KFp.jx_cov(1); % covariance on axis for current
x_cov_je   = KFp.jx_cov(2); % covariance at edge for current
x_cov_pe0  = KFp.px_cov(1);  % covariance at axis for power
x_cov_pee  = KFp.px_cov(2);  % covariance at edge for power

d_cov_j0   = KFp.jd_cov(1); % disturbance power covariance at axis
d_cov_je   = KFp.jd_cov(2); % disturbance power covariance at edge
d_cov_pe0  = KFp.pd_cov(1); % disturbance power covariance at axis
d_cov_pee  = KFp.pd_cov(2); % disturbance power covariance at edge

cov_rhopow = KFp.cov_rhopow; % convariance rho scaling: (1-rho^cov_rhopow)
cov_usegauss = KFp.cov_usegauss; % option to choose gauss shape
cov_gauss_width = KFp.cov_gauss_width; % gauss width

x_smooth_jp  = KFp.jx_sm; % smoothness
x_smooth_pe  = KFp.px_sm;
d_smooth_jp  = KFp.jd_sm;
d_smooth_pe  = KFp.pd_sm;

rgrid = model.rgrid.rhogauss;

% Specify covariance for x and d in terms of covariance of:
% for Psi: fp = [jpar;psiedge];
% for Te:  ft = [power density profile];

% transformation from source to state  noise representation
Tpx = blkdiag(model.dist.psi.Lamgaussw);
Ttx = blkdiag(model.dist.te.Lamgaussw);

Mps = gencov(x_smooth_jp ,x_cov_j0,x_cov_je,cov_rhopow,cov_usegauss,cov_gauss_width,rgrid,Tpx); %,x_cov_psie);
Mte = gencov(x_smooth_pe,x_cov_pe0,x_cov_pee,cov_rhopow,cov_usegauss,cov_gauss_width,rgrid,Ttx);
Qx = blkdiag(Mps,Mte);

% transformation from source to disturbance noise representation
Tpd = blkdiag(model.dist.psi.Lamgauss);
Ttd = blkdiag(model.dist.te.Lamgauss);

Mps = gencov(d_smooth_jp ,d_cov_j0,d_cov_je,cov_rhopow,cov_usegauss,cov_gauss_width,rgrid,Tpd);
Mte = gencov(d_smooth_pe,d_cov_pe0,d_cov_pee,cov_rhopow,cov_usegauss,cov_gauss_width,rgrid,Ttd);
Qd = blkdiag(Mps,Mte);

% initialize state error covariance matrices
Sx0 = KFp.s0*eye(KFm.nx);
Sd0 = zeros(KFm.nd);

return