function err = calcContactHeelCOPError(var, icIdx, ic0, scaling, vToe, vParams, expX, expGRF, expCOP)

err = 10e8;

Fx = 0;
Fy = 0;
Fz = 0;
Mx = 0;
My = 0;
Mz = 0;
TK1cK2a = 0;

%%
%Populate the initial conditions
%%
ic = ic0;

for i=1:1:length(var)
   ic(icIdx(i)) = ic0(icIdx(i)) + var(i)/scaling; 
end


%%
% Compute passive toe torque
%%

TK1cK2a = calcToeTorque(ic,vToe);



%%
%Compute the GRF and COP error terms
%%
tmp = xDotMex(0, ic, vParams, [0 0 0 0 0 0 TK1cK2a]');
xdot = tmp(1:1:length(ic));
contactInfo = tmp(length(ic)+1:1:length(tmp));

%%
% Compute location and rotation of every relevant frame
%%

%Parameters
buildParamVariableList;

%State Information
x    = ic(8);
y    = ic(9);
z    = ic(10);
th   = ic(11);
zeta = ic(12);
eta  = ic(13);
xi   = ic(14);

%Contact Information
HFx = contactInfo(1);
HFy = contactInfo(2);
HFz = contactInfo(3);
Hcopx  = contactInfo(4);
Hcopy  = contactInfo(5);
Hcopz  = contactInfo(6);

FFx = contactInfo(7);
FFy = contactInfo(8);
FFz = contactInfo(9);
Fcopx  = contactInfo(10);
Fcopy  = contactInfo(11);
Fcopz  = contactInfo(12);
    
%Compute all kinematic positions and frame locations.
calcPosVecsRotMatrices;


%%
% Calculate the vector from the ankle to the heel contact point
% resolved in an ankle fixed frame.
%%

%EXP Data: Compute RExp01
xExp    = expX(8);
yExp    = expX(9);
zExp    = expX(10);
thExp   = expX(11);
zetaExp = expX(12);
etaExp  = expX(13);
xiExp   = expX(14);

t1 = cos(xiExp);
t2 = cos(etaExp);
t3 = sin(etaExp);
t4 = sin(zetaExp);
t5 = sin(xiExp);
t6 = cos(zetaExp);
t7 = t1 * t4;
t8 = t5 * t6;
t9 = t1 * t6;
t10 = t5 * t4;
RExp10 = [t1 * t2 t3 * t7 + t8 -t3 * t9 + t10; -t5 * t2 -t10 * t3 + t9 t3 * t8 + t7; t3 -t2 * t4 t2 * t6;];

%EXP Data: Vector from the ankle to the COP resolved in the inertial frame
rx1P0   = expCOP'-[xExp,yExp,zExp]';
rx1P0M  = sum(rx1P0.*rx1P0).^0.5;
ex1P0   = rx1P0./rx1P0M;
ex1P1   = (RExp10)*ex1P0;

%Sim Data
r1P0    = [Hcopx Hcopy Hcopz]'-[x,y,z]';
r1P0M   = sum(r1P0.*r1P0).^0.5;
e1P0    = r1P0./r1P0M;
e1P1    = (RMK1com)*e1P0;

%Error


err = 1000*sum((e1P1-ex1P1).^2);



