function [diststruct,nd,nw] = build_dist(model)
% build disturbance models, add relevant terms to model.psi, model.te etc
% spline matrix for additive source disturbance model

scal = 1e6; % disturbance enters through current/power density, typical MW/m^3

% repeat for state equations
iinit = 0; % init for indices
diststruct.psi = build_distmodel(model.psi.dist,model.rgrid,scal,iinit);

iinit = iinit + numel(diststruct.psi.ind);
diststruct.te  = build_distmodel(model.te.dist,model.rgrid,scal,iinit);

nd = iinit + numel(diststruct.te.ind); % total number of disturbances.

nw = size(diststruct.psi.Lamgaussw,2) + ...
  size(diststruct.te.Lamgaussw ,2);

return

function distmodel = build_distmodel(distparams,rgrid,scal,iinit)
% build spline matrices
ndist = distparams.ndist;
rhogauss = rgrid.rhogauss;
rho = rgrid.rho;

switch distparams.disttype
  case 'equidist'
    
    if ndist == 0
      distmodel.Lamgauss = zeros(numel(rhogauss),[]);
      distmodel.Lam      = zeros(numel(rho),[]);
    elseif ndist == 1
      distmodel.Lamgauss = scal*ones(numel(rhogauss),1);
      distmodel.Lam      = scal*ones(numel(rho),1);
    else
      xkts = linspace(0,1,ndist)';
      % choose order 1 splines for disturbance
      Lamgauss = full(splineval(rhogauss,xkts,1)');
      Lam      = full(splineval(rho,xkts,1)');
      % disturbance model
      % source = dist.Lamgauss*dk(dist.ind); % to be added in state equation
      distmodel.Lamgauss = scal*Lamgauss;
      distmodel.Lam      = scal*Lam;
      
    end
    
  case 'gaussian'
    if ndist~=1
      error('only single gaussian supported')
    end
    distmodel.Lamgauss = scal*exp(-4*((rhogauss-distparams.distgaussrpos).^2)./(distparams.distgausswidth^2));
    distmodel.Lam      = scal*exp(-4*((rho-distparams.distgaussrpos).^2)./(distparams.distgausswidth^2));
    
  otherwise
    error('Unknown disturbance type specified')
end
distmodel.ind = iinit+[1:ndist];

% matrix for noise model
switch distparams.noisemodel
  case 'default'
    distmodel.Lamgaussw = scal * full(splineval(rhogauss,rho,1)');
  case 'sameasdist'
    distmodel.Lamgaussw = distmodel.Lamgauss;
  otherwise
    error('Unknown noise model specified')
end
return