function [res,SQ,dresdFQdot,dresdFtPQdot] = meqcde(L,LX,LY,F0,F1,Opy,Iy,Fx,Fy0,Ip,iD,IpQ,FtPQ,TQ)
% Current diffusion equation for one Domain.
% (0D, stationary-state, 1D in progress)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
assert(numel(iD)==1,'only one domain supported for CDE')

form = L.P.cde;
dt = (LX.t-LY.t); if dt==0, idt=0; else, idt=1/dt; end % avoid division by 0
i2pi = 1./(2*pi);
iTsQ = LX.signeo(:,iD)./TQ.^2; % sigma/F^2 profile appearing in various places
SQ =  (LX.rBt*i2pi) * cumtrapz(FtPQ, iTsQ);
switch form
  case {'cde0D','cde1D','cde0Dss'}
    
    FN = L.pQ.^2';
    [F0t,F1t] = meqpdomlim(LY.FA,LY.FB,L.nD); % flux limits for previous step
    FBdot = (F1(iD)-F1t(iD)) * idt;
    FAdot = (F0(iD)-F0t(iD)) * idt;
    FQ    = ( FN*(F1(iD)-F0(iD)) + F0(iD));
    FABdot = (FBdot + FAdot)/2;
    
    switch form
      case {'cde0D','cde1D'}
        error('%s does not converge and needs to be checked',form)
        FtPQdotFU = (FtPQ - LY.FtPQ(:,iD))*idt;  % dPhi/dt
        FQdotFU   = ( FN*(FBdot-FAdot) +  FAdot); % dPsi/dt
        % integral on pQ grid
        A1QFU =  (LX.rBt*i2pi) * cumtrapz(FtPQ, iTsQ.*FQdotFU  );
        A2QFU = -(LX.rBt*i2pi) * cumtrapz( FQ , iTsQ.*FtPQdotFU);
        if L.P.debug>3, fprintf('use FU\n\n'); end
        FtPQdot = FtPQdotFU;
        A1Q = A1QFU;
        A2Q = A2QFU;
      case 'cde0Dss'
        FtPQdotSS = zeros(size(FtPQ));
        A1QSS = SQ*FABdot;
        A2QSS = zeros(L.nQ,1); %Assume Phidot = 0 for stationary state;
        if L.P.debug>3, fprintf('use SS\n\n'); end
        FtPQdot = FtPQdotSS;
        A1Q = A1QSS;
        A2Q = A2QSS;
      otherwise
        error('invalid CDE form');
    end
    
    scal = 1e-6; % general scale of Ip, uesd for normalization
    % MEQ's COCOS 17 has sigmaBp=-1, hence ->  -1*(A1+A2) = (Ip-Ini) -> A1+A2 = Ini-Ip
    resIpQ = scal*(IpQ - LX.IniD(iD) + (A1Q + A2Q)); % residual of current diffusion equation on Q grid
    
    % select indices to keep
    switch form
      case '1D'
        % take a number of points equal to ng-1 (one should be left for pressure constraint?)
        iQ = ceil(linspace(0,1,L.ng)*(L.nQ)); iQ = iQ(2:end);
      otherwise
        iQ = L.nQ; % just take last one (integral 0 to 1)
    end
    res = resIpQ(iQ);
    
    
    if nargout>3
      dresdFQdot   = scal *  (LX.rBt*i2pi) * trapz(FtPQ, iTsQ.*eye(L.nQ));
      dresdFtPQdot = scal * -(LX.rBt*i2pi) * trapz( FQ , iTsQ.*eye(L.nQ));
    end
    
    if L.P.debug>3
      %  very detailed debug
      fprintf('Phibdot=%+4.2f, FAdot=%+4.2f, FBdot=%+4.2f, FA=%+4.2f, Phib=%+4.2f, FB=%+4.2f, Ip=%+4.2e, nO=%d, A1=%+4.2e A2=%+4.2e delta=%+4.2e\n',...
        FtPQdot(end),FAdot,FBdot,F0,F1,FtPQ(end),Ip,sum(Opy(:)),A1Q(end),A2Q(end),IpQ(end)+A1Q(end)+A2Q(end));
    end
  case 'OhmTor_rigid'
    res = meqcdeOhmTor_rigid(L,LX,Ip,Iy,Fy0,LY.Ia,LY.Iu,LY.Ip,idt);
  case 'OhmTor'
    res = meqcdeOhmTor(L,LX,Fx,Iy,Ip,LY.Fx,idt);
  otherwise
    error('invalid CDE case %s',form);
end

end

function [res] = meqcdeOhmTor(L,LX,Fx,Iy,Ip,zFx,idt)
% Toroidal projection Oms's law without rigid assuumption
% Iy*dFydt/Ip = - Rp (Ip - Ini)
dFx = (Fx - zFx);
res = (dFx(L.lxy)'*Iy(:)/Ip)*idt + LX.Rp*(Ip - LX.Ini);
end

function [res] = meqcdeOhmTor_rigid(L,LX,Ip,Iy,Fy0,zIa,zIu,zIp,idt)
% Toroidal projection of Ohm's law for rigid plasma
% Lp dIp/dt = -Rp Ip - Mpe dIe/dt
% Where by definition
% Lp = (Myy Iy(:))'Iy(:)/Ip^2 or provided externally
% Mpe = (Mye )'Iy(:)/Ip
% Rp = Iy (Ry Iy) /Ip^2 provided externally
%
% Myy Iy(:) is substituted with Fy - Mey *Ie for computational cost
% See F.Felici, Magnetic Modeling and Control of Tokamak, EPFL PHD course, 2020

zFy0 = [L.G.Mxa(L.lxy(:),:),L.G.Mxu(L.lxy(:),:)]*[zIa;zIu];  % Previous time step. Better to recompute instead of propagating this big matrix

dFy0dt = (Fy0-zFy0)*idt; % Derivative of flux from external current
dIpdt = (Ip - zIp)*idt; % Derivative of total plasma current

res = LX.Lp*dIpdt + LX.Rp*(Ip-LX.Ini) + dFy0dt'*Iy(:)/Ip;
end
