% tDL2 is for L2 dissipation

function [tD,dD_dCP] = Incremental_drag_diss(Aux,C,C_1)
global Invariant ini_Val Material phi dt V CC

if isfield(Material,'eta_s')
    eta_s = Material.eta_s;
else
    eta_s = 0;
end

p = ini_Val.p;
gaussno = ini_Val.gaussno;

dr = C(:,1)-C_1(:,1);
dz = C(:,2)-C_1(:,2);


%[C_v,C1_v]=interpolate_func_2(Invariant,V);

%Cortex velocity on membrane

CC= C;

Invariant_mem_to_cortex=precal_generic(abs(C(:,1)/Material.R0));
[v_c_m]=interpolate_func_2(Invariant_mem_to_cortex,V);
tmp_s = Material.tmp_s;
h_cut_off = Material.h_cut_off;
tmp_z  = 0.5*(1-tanh(tmp_s*(C(:,2)-h_cut_off)));
%v_c_m = v_c_m .* tmp_z;
%v_c_m = zeros(size(dr));


tD = sum((tmp_z.*(dr - v_c_m.*dt).^2.*Aux.au.*C(:,1)).*Invariant.gweight);


%tD = sum(dr.^2.*Aux.au.*C(:,1).*Invariant.gweight);
tD = tD.*eta_s.*pi/dt;
tD = tD +  Material.eta_END * dr(end)^2/dt;

if nargout == 2
    %Derivatives of E (local)
    temp1 = (tmp_z.*(dr - v_c_m.*dt).^2.*Aux.au ...
         + 2*tmp_z.*(dr - v_c_m.*dt).*Aux.au.*C(:,1))*ones(1,p+1);
    temp2 = (tmp_z.*(dr - v_c_m.*dt).^2.*C(:,1))*ones(1,p+1);
    temp2b = - (0.5*tmp_s*sech(tmp_s*(C(:,2)-h_cut_off)).^2.*(dr - v_c_m.*dt).^2.*C(:,1).*Aux.au)*ones(1,p+1);
    
    dD_dr_l = (temp1.*Invariant.Nshape + temp2.*Aux.da_dr_l).*Aux.Mgweight;
    dD_dz_l = (temp2.*Aux.da_dz_l + temp2b.*Invariant.Nshape).*Aux.Mgweight;  

    
    ns = length(dD_dr_l)/gaussno;
    vs = [1:1:ns];
    dD_dCPr = reshape(sum(reshape(dD_dr_l,gaussno,ns,p+1),1),ns,p+1);
    dD_dCPz = reshape(sum(reshape(dD_dz_l,gaussno,ns,p+1),1),ns,p+1);
    dD_dCP = zeros(2,ini_Val.ncontrol);

    for i = 0:p
        dD_dCP(1,vs(:)+i) = dD_dCP(1,vs(:)+i) + dD_dCPr(:,i+1)';
        dD_dCP(2,vs(:)+i) = dD_dCP(2,vs(:)+i) + dD_dCPz(:,i+1)';
    end
            
    dD_dCP= dD_dCP * pi/dt*eta_s;
end
dD_dCP(1,end) = dD_dCP(1,end) + 2 * Material.eta_END* dr(end)/dt;