function [mats, varargout] = build_matrices(model,config,geom,g,rho,rhogauss,varargin)
%function model = build_matrices(model,config,geom,g,rho,rhogauss)
% FE matrices and boundary conditions are prepared for state transport
% equations.
%
%% Build FE matrices

% to indicate projection on H mode spline basis
hmode_modeltype_imposed = false;
if nargout == 2
  hmode_modeltype_imposed = true;
end

% constants
mu0 = 4e-7*pi;
qelectron = 1.6e-19;

model.geom = geom; % temp placeholder to get right input for eval functions
model.rgrid.rho = rho;
model.rgrid.rhogauss = rhogauss;

% knots points for matrices
xkts = rho;
% spline order
sporder = config.numerics.sporder;
ngauss = config.numerics.ngauss;

%% PSI equation
% time-invariant vectors to go into FE matrices
psi = model.psi;

% psi equation
mpsi.fixed = 16*pi^2*mu0*rhogauss    * model.psi.scal;
apsi.fixed = 8*pi^2 *mu0*rhogauss.^2 * model.psi.scal;
dpsi.fixed = ones(size(rhogauss)) * model.psi.scal;
hpsi.fixed = ones(size(rhogauss)) * model.psi.scal;
spsi.fixed = -8*pi^2*mu0 * ones(size(rhogauss));

% Psi equation
di=0; dj=0; % zeroth derivative
[Mpsi.XI,Mpsi.XIindex,Mpsi.nsp] = FEconstruct(xkts,mpsi.fixed,sporder,ngauss,di,dj,[1 0 0]);

di=0; dj=1;
[Apsi.XI,Apsi.XIindex,Apsi.nsp] = FEconstruct(xkts,apsi.fixed,sporder,ngauss,di,dj,[1 0 0]);

di=1; dj=1; % first derivative
[Dpsi.XI,Dpsi.XIindex,Dpsi.nsp] = FEconstruct(xkts,dpsi.fixed,sporder,ngauss,di,dj,[1 0 0]);

di=0; dj=0;
[Hpsi.XI,Hpsi.XIindex,Hpsi.nsp] = FEconstruct(xkts,hpsi.fixed,sporder,ngauss,di,dj,[1 0 0]);

di=0; dj=[]; % for source term;
[Spsi.XI,Spsi.XIindex,Spsi.nsp] = FEconstruct(xkts,spsi.fixed,sporder,ngauss,di,dj,[1 0 0]);

%%%% Boundary conditions
Dpsi_BC = zeros(psi.nsp,psi.nsp,1);
Spsi_BC = zeros(psi.nsp,1,1);

switch psi.BC.type
  case 'Ip'
    Spsi_BC(end) = 16*pi^3*mu0;
    Dpsi_BC(end,end) = 0;
  case 'Vloop'
    error('not implemented yet');
  case 'PsiOH'
    error('yet to adapt this expression to new equations')
    %Dpsi_BC(end,end) = mu0/eq_config.Lext;
    %Spsi_BC(end)   = mu0/eq_config.Lext;
  otherwise
    error('unknown BC option');
end

% scaling:
Spsi_BC = Spsi_BC/model.psi.scal;

% Store in model structure
% Psi equation
mats.Mpsi = Mpsi;
mats.Apsi = Apsi;
mats.Dpsi = Dpsi;
mats.Hpsi = Hpsi;
mats.Spsi = Spsi;

mats.Dpsi_BC = Dpsi_BC;
mats.Spsi_BC = Spsi_BC;

%% Te equation
te=model.te;
% time-invariant vectors to go into FE matrices

% flag to indicate whether H mode spline basis needs to be generated
if hmode_modeltype_imposed
  nknots = numel(xkts);
  rhoped = config.hmode.params.rhoped;
  xkts_core = linspace(0, rhoped, nknots)';
end

% geometric profiles are now time-varying, so assign reduced fixed part
% only (rho and physical constants)

% Te equation
mTe.fixed = 3/2 * qelectron * model.te.scal;
aTe.fixed = 3/4 * qelectron * rhogauss * model.te.scal;
dTe.fixed = qelectron  * model.te.scal; % q electron is since T must be in Kelvin!
hTe.fixed = 1/2 * qelectron * model.te.scal;
sTe.fixed = ones(size(rhogauss));

di=0; dj=0; % zeroth derivative: MASS matrix
[MTe.XI,MTe.XIindex,MTe.nsp] = FEconstruct(xkts,mTe.fixed,sporder,ngauss,di,dj,[1 0 0]);
if hmode_modeltype_imposed
  [MTe_hmode.XI,MTe_hmode.XIindex,MTe_hmode.nsp] = FEconstruct(xkts,mTe.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
end

di=1; dj=0;
[ATe.XI,ATe.XIindex,ATe.nsp] = FEconstruct(xkts,aTe.fixed,sporder,ngauss,di,dj,[1 0 0]);
if hmode_modeltype_imposed
  [ATe_hmode.XI,ATe_hmode.XIindex,ATe_hmode.nsp] = FEconstruct(xkts,aTe.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
end

di=1; dj=1; % first derivative: STIFFNESS matrix
[DTe.XI,DTe.XIindex,DTe.nsp] = FEconstruct(xkts,dTe.fixed,sporder,ngauss,di,dj,[1 0 0]);
if hmode_modeltype_imposed
  [DTe_hmode.XI,DTe_hmode.XIindex,DTe_hmode.nsp] = FEconstruct(xkts,dTe.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
end

di=0; dj=0;
[HTe.XI,HTe.XIindex,HTe.nsp] = FEconstruct(xkts,hTe.fixed,sporder,ngauss,di,dj,[1 0 0]);
if hmode_modeltype_imposed
  [HTe_hmode.XI,HTe_hmode.XIindex,HTe_hmode.nsp] = FEconstruct(xkts,hTe.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
end

di=0; dj=[]; % for SOURCE term;
[STe.XI,STe.XIindex,STe.nsp] = FEconstruct(xkts,sTe.fixed,sporder,ngauss,di,dj,[1 0 0]);
if hmode_modeltype_imposed
  [STe_hmode.XI,STe_hmode.XIindex,STe_hmode.nsp] = FEconstruct(xkts,sTe.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
end

%%%% Boundary condition

% Te equation
% *Dirichlet condition: substitute last equation with Tehat = BC
% *Fixed normalized gradient: dTe/drho = c*Te
%   then boundary term becomes Lamb(end)*c*Te which is to be added to DTe

% DTe_BC = zeros(te.nsp,te.nsp,1);
% STe_BC = zeros(te.nsp,1,1);

% Te equation
mats.MTe = MTe;
mats.ATe = ATe;
mats.DTe = DTe;
mats.HTe = HTe;
mats.STe = STe;

% mats.DTe_BC = DTe_BC;
% mats.STe_BC = STe_BC;
if hmode_modeltype_imposed
  mats_hmode.MTe = MTe_hmode;
  mats_hmode.ATe = ATe_hmode;
  mats_hmode.DTe = DTe_hmode;
  mats_hmode.HTe = HTe_hmode;
  mats_hmode.STe = STe_hmode;

%   mats_hmode.DTe_BC = DTe_BC;
%   mats_hmode.STe_BC = STe_BC;
end

%% ne equation
if strcmp(model.ne.method,'state')
  ne=model.ne;
  % time-invariant vectors to go into FE matrices
  
  mne.fixed = model.ne.scal;
  ane.fixed = model.ne.scal;
  dne.fixed = model.ne.scal;
  hne.fixed = model.ne.scal;
    fne.fixed = 1;
  sne.fixed = ones(size(rhogauss));
  
  di=0; dj=0; % zeroth derivative: MASS matrix
  [Mne.XI,Mne.XIindex,Mne.nsp] = FEconstruct(xkts,mne.fixed,sporder,ngauss,di,dj,[1 0 0]);
  if hmode_modeltype_imposed
    [Mne_hmode.XI,Mne_hmode.XIindex,Mne_hmode.nsp] = FEconstruct(xkts,mne.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
  end
  
  di=1; dj=0;
  [Ane.XI,Ane.XIindex,Ane.nsp] = FEconstruct(xkts,ane.fixed,sporder,ngauss,di,dj,[1 0 0]);
  if hmode_modeltype_imposed
    [Ane_hmode.XI,Ane_hmode.XIindex,Ane_hmode.nsp] = FEconstruct(xkts,ane.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
  end
  
  di=1; dj=1; % first derivative: STIFFNESS matrix
  [Dne.XI,Dne.XIindex,Dne.nsp] = FEconstruct(xkts,dne.fixed,sporder,ngauss,di,dj,[1 0 0]);
  if hmode_modeltype_imposed
    [Dne_hmode.XI,Dne_hmode.XIindex,Dne_hmode.nsp] = FEconstruct(xkts,dne.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
  end
    
  di=0; dj=0;
  [Hne.XI,Hne.XIindex,Hne.nsp] = FEconstruct(xkts,hne.fixed,sporder,ngauss,di,dj,[1 0 0]);
  if hmode_modeltype_imposed
    [Hne_hmode.XI,Hne_hmode.XIindex,Hne_hmode.nsp] = FEconstruct(xkts,hne.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
  end
  
  di=0; dj=[]; % for SOURCE term;
  [Sne.XI,Sne.XIindex,Sne.nsp] = FEconstruct(xkts,sne.fixed,sporder,ngauss,di,dj,[1 0 0]);
  if hmode_modeltype_imposed
    [Sne_hmode.XI,Sne_hmode.XIindex,Sne_hmode.nsp] = FEconstruct(xkts,sne.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
  end
  di=1; dj=[]; % for FLUX term;
  [Fne.XI,Fne.XIindex,Fne.nsp] = FEconstruct(xkts,fne.fixed,sporder,ngauss,di,dj,[1 0 0]);
  if hmode_modeltype_imposed
    [Fne_hmode.XI,Fne_hmode.XIindex,Fne_hmode.nsp] = FEconstruct(xkts,fne.fixed,sporder,ngauss,di,dj,[1 0 0], xkts_core);
  end
  
  %%%% Boundary condition
  
  % ne equation
  % *Dirichlet condition: substitute last equation with nehat = BC
  
  Dne_BC = zeros(ne.nsp,ne.nsp,1);
  Sne_BC = zeros(ne.nsp,1,1);
  
  % ne equation
  mats.Mn = Mne;
  mats.An = Ane;
  mats.Dn = Dne;
  mats.Hn = Hne;
  mats.Sn = Sne;
  mats.Fn = Fne;
  
  mats.Dn_BC = Dne_BC;
  mats.Sn_BC = Sne_BC;
  if hmode_modeltype_imposed
    mats_hmode.Mn = Mne_hmode;
    mats_hmode.An = Ane_hmode;
    mats_hmode.Dn = Dne_hmode;
    mats_hmode.Hn = Hne_hmode;
    mats_hmode.Sn = Sne_hmode;
    mats_hmode.Fn = Fne_hmode;

    mats_hmode.Dn_BC = Dne_BC;
    mats_hmode.Sn_BC = Sne_BC;
  end
  
end
if hmode_modeltype_imposed
  varargout{1} = mats_hmode;
end
return

end
