clear
syms wb x real
lb = 5.5e-2;
lf = 70e-3;
wf = 1e-3;
theta = deg2rad(60);
r_egg = 21.5e-3;
dx_beam = 1e-3;
dx_beam_eff = dx_beam * (1 + cos(theta / 2));
xrange = 7.088e-2;
egg_contact = r_egg * (1 + 1 / sin(theta / 2));
x_contact = xrange - egg_contact;
x_max = x_contact + dx_beam_eff;
t = 0.00635;
E = 1.9e9;
gamma = 0.85147;
ktheta = 2.65;
Ib = t * wb^3 / 12;
If = t * wf^3 / 12;
Lcp = (1 - gamma) * lf / 2;
B = asin(x / lf);
kbeam = 192 * E * Ib / lb^3;
kbeam2 = kbeam * sin(theta / 2) / 2;
kgripper = 1 / (1 / kbeam + 1 / (2 * kbeam2));
kt = gamma * ktheta * E * If / lf;
Ustored = 2 * kt * B^2;
Fstored = diff(Ustored, x);
Kstored = diff(Fstored, x);
x_eff = @(x_in) max(0, x_in - x_contact);
K_exert = @(wb_in) double(subs(kgripper, wb, wb_in));
F_exert = @(x_in, wb_in) K_exert(wb_in) * x_eff(x_in);
n = 50;
x_vals = linspace(x_contact, x_max, n);
wb_vals = linspace(0, 5e-3, n);
[WB, X] = meshgrid( wb_vals,x_vals);
Z = arrayfun(@(x, wb) F_exert(x, wb), X, WB);
figure;
surf(X,WB, Z);
xlabel('x (m)');
ylabel('wb (m)');
zlabel('F_{exert}(x, wb) (N)');
title('3D Plot of F_{exert}(x, wb)');
shading interp;
colormap winter;
colorbar;
view(3);