function [L,LX,LY] = LLXLY_IDS_corsica(L_LIU, a)
%% Load data drom IDS


% Get geomtrical information from the existing LIUQE
L = struct();
L.P = L_LIU.P;
fieldlist = {'rl';'zl'; 'rw';'zw'; 'Twa'; 'rv';'zv'; 'Tvs'; 'rf'; 'zf';'rm'; 'zm'; 'am'};
for  field = fieldlist.'
   L.G.(field{1})  = L_LIU.G.(field{1});
end

% Set new grid
L.G.rx = a.equilibrium.time_slice{1}.profiles_2d{1}.grid.dim1;
L.G.zx = a.equilibrium.time_slice{1}.profiles_2d{1}.grid.dim2;
L.P.pq = sqrt((a.equilibrium.time_slice{1}.profiles_1d.psi - a.equilibrium.time_slice{1}.global_quantities.psi_axis)/(a.equilibrium.time_slice{1}.global_quantities.psi_boundary - a.equilibrium.time_slice{1}.global_quantities.psi_axis));

% Recompute the geometrical quantities for the CORSICA grid
L.G = liug(L.G,L.P);
L = liuc(L.P, L.G); 

[rrx, zzx] = meshgrid(L.G.rx, L.G.zx);

%%  Get LX structure
% {Ia, Ip, rBt} available
% Ff, Bm, Ft recomputed

LX.t = a.pf_active.coil{1}.current.time;
LX.nt = numel(LX.t);

LX.Ip = zeros(LX.nt,1);
for ii=1:LX.nt
    LX.Ip(ii) = a.equilibrium.time_slice{ii}.global_quantities.ip;
end

% Ia is the current per number of turns
LX.Ia = zeros(L.G.na, LX.nt);
for ii= 1:L.G.na
   LX.Ia(ii,:) = a.pf_active.coil{ii}.current.data/sum(L.G.Twa(:,ii));
end

LX.rBt = zeros(LX.nt,1);
LX.rBt(:) = a.equilibrium.vacuum_toroidal_field.r0*a.equilibrium.vacuum_toroidal_field.b0;


%% Synthetic data that needs to be recomputed
LY.Ix = zeros( L.nzx, L.nrx ,LX.nt);
for tt=1:LX.nt
    LY.Ix(:,:,tt) = a.equilibrium.time_slice{tt}.profiles_2d{1}.j_tor./rrx*L.drx*L.dzx*1e-2; % Wrong definition in IDS Corsica
   % Ix_Corsica(:,:,tt) = a.equilibrium.time_slice{tt}.profiles_2d{1}.j_tor; % Wrong definition of CORSICA
   % Ix_Corsica(:,:,tt) = -Ix_Corsica(:,:,tt)*abs(LX.Ip(tt)/sum(sum(Ix_Corsica(:,:,tt))));
end


LX.Bm = zeros(L.G.nm, LX.nt);
for tt= 1:LX.nt
    LX.Bm(:,tt) = L.G.Bmx*reshape(LY.Ix(:,:,tt), L.nrx*L.nzx,1 )  + L.G.Bma*LX.Ia(:,tt);
end

LX.Ff = zeros(L.G.nf, LX.nt);
for tt= 1:LX.nt
    LX.Ff(:,tt) =  L.G.Mfx*reshape(LY.Ix(:,:,tt), L.nrx*L.nzx,1 )  + L.G.Mfa*LX.Ia(:,tt);
end

LX.Is = zeros(0,LX.nt);

LX.shot = 5000;

LY.Fx = zeros(L.nzx,L.nrx ,LX.nt);
for tt = 1:LX.nt
   LY.Fx(:,:,tt) =  a.equilibrium.time_slice{tt}.profiles_2d{1}.psi;
end

% Normalized flux of the R,Z grid
LY.Opy = zeros(L.nzy, L.nry, LX.nt);
LY.FNy = zeros(L.nzy, L.nry, LX.nt);
for tt=1:LX.nt
[LY.rA(tt),LY.zA(tt),LY.FA(tt),dr2FA,dz2FA,drzFA,rX,zX,~,~,~,~,rl,zl,...
    LY.rB(tt),LY.zB(tt),LY.FB(tt),LY.lB(tt),LY.lX(tt),LY.Opy(:,:,tt)] = meqpdom(LY.Fx(:,:,tt),LX.Ip(tt),L.P.isaddl,L);
LY.FNy(:,:,tt) =  (LY.Fx(2:end-1,2:end-1,tt) - LY.FA(tt))/(LY.FB(tt) - LY.FA(tt));
end



%% Recompute the plasma current as a check. There is a bit of a difference
psiN = (a.equilibrium.time_slice{tt}.profiles_1d.psi - a.equilibrium.time_slice{tt}.global_quantities.psi_axis)/(a.equilibrium.time_slice{tt}.global_quantities.psi_boundary - a.equilibrium.time_slice{tt}.global_quantities.psi_axis);
for tt=1:LX.nt
pprime = a.equilibrium.time_slice{tt}.profiles_1d.dpressure_dpsi;
p = spline(psiN,pprime);
pprimey = ppval(p, LY.FNy(:,:,tt));
pprimey = pprimey.*LY.Opy(:,:,tt);

ttprime = a.equilibrium.time_slice{tt}.profiles_1d.f_df_dpsi;
T2d2 = cumtrapz(a.equilibrium.time_slice{tt}.profiles_1d.psi,ttprime);
T2d2 = T2d2 - T2d2(end);

p = spline(psiN,ttprime);
ttprimey = ppval(p, LY.FNy(:,:,tt));
ttprimey = ttprimey.*LY.Opy(:,:,tt);

p = spline(psiN,T2d2);
Ty = ppval(p, LY.FNy(:,:,tt));
Ty = Ty.*LY.Opy(:,:,tt);

L.rymap = repmat(L.ry',numel(L.zy),1); % Matriz map ry
LX.Ft(tt) = sum(sum( Ty/LX.rBt(tt)*L.drx*L.dzx./L.rymap ));

Iy = 2*pi*(L.rymap.*pprimey+1./L.rymap.*ttprimey/mu0)*L.drx*L.dzx; % Recomputed plasma current

LY.Ip(tt) = sum(sum(Iy));
LY.FA(tt) = a.equilibrium.time_slice{tt}.global_quantities.psi_axis;
LY.FB(tt) = a.equilibrium.time_slice{tt}.global_quantities.psi_boundary;
LY.Wk(tt) = a.equilibrium.time_slice{tt}.global_quantities.w_mhd;
LY.PpQ(:,tt) = pprime;
LY.TTpQ(:,tt) = ttprime;
LY.psiN(:,tt) = psiN;
LY.PQ(:,tt) = a.equilibrium.time_slice{tt}.profiles_1d.pressure;
end
LY.t = LX.t;

% Put all vector in raw form
for field = fieldnames(LX).'
    LX.(field{1}) = v2rw(LX.(field{1}));
end

% Put all vector in raw form
for field = fieldnames(LY).'
    LY.(field{1}) = v2rw(LY.(field{1}));
end


end


% Put all vectors in raw form
function vout = v2rw(vin)
if numel(size(vin))==2 && sum(ismember(1,size(vin)))==1 && ~ischar(vin)
   vout = reshape(vin, numel(vin),1);
else
    vout = vin;
end
end