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=1
warning('Fill in linearised model for S_{Ohm} here')
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
warning('Fill in linearised model for P_{aux} here')
KPcond=-0.69*Scond/SH/(3*p.n*q_e*V);
KP=KPaux+KPcond;


%% construct transfer functions

s = tf('s');
G=tf(1)   % delta(T) to delta(Paux) transfer function
warning('Fill in transfer function between dP and dT here')

H=tf(1)    % delta(T) to delta(neutron power) transfer function
warning('Fill in transfer function between dPn and dT here (only needed in queston 3d)')


end