% Default
% wireRad = 1.0e-3
% springRad = 5.0e-3
% pitch = 3.0e-3
% noOfLoops = 10
% Ex = 200e9
% nuxy = 0.3
% load = -100

data = [[1.0e-3 5.0e-3 3.0e-3 10 200e9 0.3 -100 -0.00680499];...
        [1.0e-3 5.0e-3 3.0e-3 10 200e9 0.3 100 0.00644066];...
        [1.0e-3 5.0e-3 3.0e-3 10 200e9 0.15 -100 -0.00608688];...
        [1.0e-3 5.0e-3 3.0e-3 10 200e9 0.45 -100 -0.00770514];...
        [1.0e-3 5.0e-3 3.0e-3 10 10e9 0.3 -1 -1.37747e-3];...
        [1.0e-3 5.0e-3 3.0e-3 10 100e9 0.3 -10 -1.37747e-3];...
        [1.0e-3 5.0e-3 3.0e-3  5 200e9 0.3 -10 -0.351917e-3];...
        [1.0e-3 5.0e-3 2.1e-3 10 200e9 0.3 10 0.683918e-3];...
        [1.0e-3 5.0e-3 4.0e-3 10 200e9 0.3 -10 -0.694533e-3];...
        [1.0e-3 3.0e-3 3.0e-3 10 200e9 0.3 -10 -0.149202e-3];...
        [1.0e-3 1.0e-2 3.0e-3 10 200e9 0.3 -10 -0.0055028];...
        [3.0e-4 5.0e-3 3.0e-3 10 200e9 0.3 -10 -7.71733e-3];...
        [1.4e-3 5.0e-3 3.0e-3 10 200e9 0.3 10 0.34779e-3]];

% Calculate spring stiffness (force /displacement), shear modulus,
% torsional area moment of inertia
[rows, cols] = size(data);
for ii=1:rows
  k(ii) = data(ii, 7)/data(ii, 8);
  nu = data(ii,6);
  G = data(ii, 5)/(2*(1+nu));
  d = 2*data(ii,1);
  D = 2*data(ii,2);
  n = data(ii,4);
  J = pi*d^4/32;
  p = data(ii,3);
  alpha = atan2(p, pi*D);
  c = D/d;
  f11 = (n*pi*D^3)/(4*G*J*cos(alpha))*((1+0.5/c^2)*cos(alpha)^2 + (1+0.25/c^2)/(1+nu)*sin(alpha)^2);
  kcalc(ii) = 1/f11;
  GG(ii) = G;
  JJ(ii) = J;
  Alpha(ii) = alpha*180/pi;
end
[data(:,1:7) GG' JJ' Alpha' data(:,8) k' kcalc']
loglog(k, kcalc, 'o', 'LineWidth', 3, 'MarkerSize', 10); hold on;
loglog([100 7.0e4],[100, 7.0e4], 'r-', 'LineWidth', 3);
xlabel('Stiffness (closed form) (N/m)', 'FontSize', 18);
ylabel('Stiffness (Abaqus) (N/m)', 'FontSize', 18);