function [dIdB,dIdx] = gradObjectiveFunc(xpl,xpt,x,e,ep,B,Breq)

% [dIdB,dIdx] = gradObjectiveFunc(xpl,xpt,x,e,ep,B,Breq)
%        calculates gradient of the following cost functional:
%        I(B,x) = (1/meas(Omega_p)) * \int_{Omega_p} |B(x)-Breq|^2 dx
%        w.r.t. magnetic field and grid displacements, respectively

    
dIdB = zeros(2*max(size(e)),1);
dIdx = zeros(size(x,1)*size(x,2),1);

Breq = [0.025;0];
area = abs(xpl)*abs(xpt);
np = max(size(ep));

for i=1:np
    ei = e(ep(i),:);
    
    x1 = x(:,ei(1));
    x2 = x(:,ei(2));
    x3 = x(:,ei(3));
    
    R = [x2-x1 x3-x1];
    detR = det(R);
    areai = abs(detR)/2;
    
    dIdB(2*ep(i)-1) = (2*(B(1,ep(i)) - Breq(1))*areai) / (area*power(norm(Breq),2));
    dIdB(2*ep(i)) = (2*(B(2,ep(i)) - Breq(2))*areai) / (area*power(norm(Breq),2));


    graddetR = [x2(2)-x3(2) x3(2)-x1(2) x1(2)-x2(2);x3(1)-x2(1) x1(1)-x3(1) x2(1)-x1(1)]';

    dIdx(2*ei(1)-1) = dIdx(2*ei(1)-1) + (power(norm(B(:,ep(i))-Breq),2)*sign(detR)*graddetR(1,1)) / (2*area*power(norm(Breq),2));
    dIdx(2*ei(1)) = dIdx(2*ei(1)) + (power(norm(B(:,ep(i))-Breq),2)*sign(detR)*graddetR(1,2)) / (2*area*power(norm(Breq),2));
    dIdx(2*ei(2)-1) = dIdx(2*ei(2)-1) + (power(norm(B(:,ep(i))-Breq),2)*sign(detR)*graddetR(2,1)) / (2*area*power(norm(Breq),2));
    dIdx(2*ei(2)) = dIdx(2*ei(2)) + (power(norm(B(:,ep(i))-Breq),2)*sign(detR)*graddetR(2,2)) / (2*area*power(norm(Breq),2));
    dIdx(2*ei(3)-1) = dIdx(2*ei(3)-1) + (power(norm(B(:,ep(i))-Breq),2)*sign(detR)*graddetR(3,1)) / (2*area*power(norm(Breq),2));
    dIdx(2*ei(3)) = dIdx(2*ei(3)) + (power(norm(B(:,ep(i))-Breq),2)*sign(detR)*graddetR(3,2)) / (2*area*power(norm(Breq),2));
end
