function LY = reprocess_Ohm_model(L,LY, varargin)
% Re-process the Ohm's law CDE type from input L, LY from liuqe or fge
% Useful to compute total plasma resistance Rp and total plasma self-inductance Lp 
%
% Original equal: dFydt = Myy dIydt + Mye dIedt = -Ry Iy  
% Multiplying scalarly by Iy/Ip one obtains
% 1) Lp dIp/dt + 1/2 Ip dLp/dt = -Rp Ip  - Mpe dIedt
% with 
% Rp = Iy'*(Ry.*Iy) /Ip^2 total plasma resistence
% Lp = (Iy(:)'*Myy*Iy(:))/Ip^2 = (Iy(:)'*Fy(:))/Ip^2 ttotal plasma self-inductance
% Mpe = (Iy(:)'*Mye )/Ip 
%
% It returns Lp(t) from analytic formula
% It computes Rp(t)  Rp =  (Iy*dFydt)/Ip^2. The estimation of Rp is more reliable during stationary phase where dIydt =0;
% Given Lp(t), Rp, Mpe solve the ode 1) and compare with dIp/dt inputs
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

%% Default parameters
P.dtwmvavg = 1e-3; 
P.doplot = true; % Make plots
P.idLpdt = true; % Include dLpdt in the model

%% Overwrite defaults with inputs
for k = 1:2:numel(varargin)
 P.(varargin{k}) = varargin{k+1};
end

%% Compute vacuum (Fy0) and plasma (Fyp) induced poloidal flux
nt = numel(LY.t);
LY.Fx0  = L.G.Mxa*LY.Ia + L.G.Mxu*LY.Iu;
LY.Fxp =  reshape(LY.Fx, [],nt) -  LY.Fx0;

%% Compute plasma self inductange Lp = (Iy' (Myy Iy))/Ip^2
LY.Lp = dot(LY.Fxp(L.lxy(:),:),reshape(LY.Iy,L.ny , nt))./LY.Ip.^2;

%%  Compute derivatives first order FD
fields = {'Ip', 'Ia', 'Iu', 'Fx0', 'Fxp', 'Lp'};
dFDo1 = @(x) [diff(x,1,2)./diff(LY.t), (x(:,end)-x(:,end-1))./(LY.t(end)-LY.t(end-1))];
for ii = 1:numel(fields)
 name_out = ['d' fields{ii} 'dt'];
 LY.(name_out) = dFDo1(LY.(fields{ii}));
end

%% Estimate Rp from (Iy*dFy/dt)/Ip = - Rp Ip. This is better estimated during stationary state where dFydt ~ dFy0dt
LY.Rp = -dot(LY.dFx0dt(L.lxy(:),:) + LY.dFxpdt(L.lxy(:),:) ,reshape(LY.Iy,L.ny , nt))./LY.Ip.^2;

%% Solve the ode for Ip
if P.idLpdt
dLpdt = LY.dLpdt;
else
dLpdt = 0;
end

Rp = LY.Rp; % Plasma resistence used
Lp = LY.Lp; % Internal inductance used

% Inputs
u = -dot(LY.dFx0dt(L.lxy(:),:),reshape(LY.Iy,L.ny , nt))./LY.Ip;

% Solve ODE
dIpdt = @(t, Ip)  (-interp1(LY.t, Rp + 1/2*dLpdt, t)*Ip + interp1(LY.t,u, t))/interp1(LY.t,Lp, t);
[t_out,Ip_out] = ode45(@(t,Ip) dIpdt(t, Ip), LY.t, LY.Ip(1)); % Solve ODE

%% Make plots
if P.doplot
 m = 3;
 n = 1;
 
 ind = 0;
 
 figure('Position', [100 100 400 800])
 ind = ind +1;
 subplot(m,n,ind)
 hold on
 plot(LY.t, LY.Ip)
 plot(t_out,Ip_out)
 legend({'input data', 'solved model'}, 'location', 'best')
 ylabel('Ip')
 xlabel('t [s]')
 title('Lp dIpdt = -(Rp + 1/2 dLpdt) Ip - Mpe dIedt')
 
 
 ind = ind +1;
 subplot(m,n,ind)
 hold on 
 plot(LY.t,Rp)
 plot(LY.t, 1/2*dLpdt)
 legend({'Rp  = -(Iy*dFy/dt)/Ip^2 ','1/2*dLpdt'}, 'location', 'best');
 xlabel('t [s]')
 title('Lp dIpdt = -(Rp + 1/2 dLpdt) Ip - Mpe dIedt')
 
 
 ind = ind +1;
 subplot(m,n,ind)
 hold on 
 plot(LY.t,(Rp + 1/2*dLpdt).*LY.Ip)
 plot(LY.t, u)
 plot(LY.t, Lp.*LY.dIpdt)
 legend({'Rp + 1/2*dLpdt', 'Mpe dIedt','Lp dIpdt'}, 'location','best');
 xlabel('t [s]')
 title('Lp dIpdt = -(Rp + 1/2 dLpdt) Ip - Mpe dIedt')
 
end

end



