function [G,KT,KP,H,S]=linearise_model(T0,P0,p)
% This function implements the linear respons of the 0D fusion reactor to
% changes dT in temperature and dP in auxiliary power away from an
% equilibrium operating point T0, P0.

% What is obtained is the linear approximation
% ddT/dt [eV] = 1/(3en) ( dS/dT dT + dS/dP dP ) = KT dt + KP dP
% This is then converted into a transfer function
% dT/dP = G(s)

% Additionally, a transfer function  between neutron power deviation away
% from its equilibrium value and dT is derived
% dPn/dT = H


%% linearise sources

S=sources(T0,P0,p);

Salpha=S.S_alpha(end);
Sohm=S.S_ohm(end);
Saux=S.S_aux(end);
Srad=S.S_rad(end);
Scond=S.S_cond(end);
SH=Salpha+Sohm+Saux;    % total heating power


q_e = 1.60217657E-19; % electron charge [coulomb] or Boltzmann constant [J/eV]
a = [-21.38, -25.20, -7.101e-2 1.938e-4 4.925e-6 -3.984e-8];    % redefine some constants
alp = 0.2935;   % redefine some constants
Sexp = -alp*a(1)./((T0/1000).^alp) + a(3)*(T0/1000) + ...
         2*a(4)*(T0/1000).^2 + 3*a(5)*(T0/1000).^3 + 4*a(6).*(T0/1000).^4;    % recalculate quantity I need later
V = (2*pi*p.R_0)*(pi*p.kappa*p.a^2); % plasma volume [m^3]


% construct KT
KTohm=-3/2*Sohm/(3*p.n*q_e*T0);
KTrad=-Srad/2/(3*p.n*q_e*T0);
KTcond=(-1+3/2*0.69*Sohm/SH-0.69*Sexp*Salpha/SH)*Scond/(3*p.n*q_e*T0);
KTalpha=Sexp*Salpha/(3*p.n*q_e*T0);
KT=KTohm+KTrad+KTcond+KTalpha;

% contruct KP
KPaux=1/(3*p.n*q_e*V);
KPcond=-0.69*Scond/SH/(3*p.n*q_e*V);
KP=KPaux+KPcond;

% KT=1/(3*p.n*q_e*T0) ...
%     *( 3/2*(-1+0.69*Scond/SH)*Sohm-Srad/2-Scond+(1-0.69*Scond/SH)*Sexp*Salpha );
% KP=1/(3*p.n*q_e*V)*(1-0.69*Scond/SH);


%% construct transfer functions

s = tf('s');
G=KP/(s-KT);  % delta(Paux) to delta(T) transfer function

H=4*V*Sexp*Salpha/T0;    % delta(neutron power) to delta(T) transfer function
% H=4*V*KTalpha*(3*p.n*q_e);  % alternative calculation

end