%fixedguidedbeambending.m - This script finds the x- and y- forces for
%motion of the end of a guided beam along a line with angle gamma from the
%vertical. It models behavior of a bistable fixed-guided segment.

t = 1.0e-5;%m - in-plane thickness
l = 1.1e-3;%m - beam length
gamma = -3;%degrees - angle of the line of action
w = 2e-4;%m - out-of-plane width
E = 8e9;% Pa - Young's modulus

bmax = -.2*l*sind(gamma);%the largest value of vertical deflection
numbs = 131;%the number of steps required to get to the maximum deflection
% If numbs is large, resolution will be better but solve time may be long.

b = [bmax/numbs:bmax/numbs:bmax];%create a vector of desired vertical deflections

a = -b*tand(gamma);%create a vector of desired horizontal deflections

%Call the function to find the solution.
[Fx,Fy,errx,erry,s,M,sax,sbend,x,y,mode,k] = beambending(a,b,l,t,w,E,100);

figure(1)%This is a plot of desired and actual points, to verify that the solution converged.
clf
plot(errx+a,erry+b,'o',a,b)

xlabel('Horizontal Motion')
ylabel('Vertical Motion')

ft = -Fx*sind(gamma)+Fy*cosd(gamma);%Calculate the force along the line of action


del = sqrt((a).^2+b.^2);%Calculate the deflection along the line of action
figure(2)%Plot the force vs. deflection
clf
plot(del,ft)

xlabel('Displacement')
ylabel('Force')

equilib = find((ft(1:end-1).*ft(2:end))<0);
U = cumtrapz([0 del],[0 ft]);
if isempty(equilib)
    display('Not Bistable')
elseif numel(equilib)==2
    stable = equilib(2);
    unstable = equilib(1);
    dels = del(stable) - ft(stable)*(del(stable+1) - del(stable))/(ft(stable+1) - ft(stable));
    Us = U(stable+1) + (U(stable+2)-U(stable+1))*(dels - del(stable))/(del(stable+1) - del(stable));
    stablepos = dels
    delu = del(unstable) - ft(unstable)*(del(unstable+1) - del(unstable))/(ft(unstable+1) - ft(unstable));
    unstablepos = delu
end
    
smin = min(sax - abs(sbend));
smax = max(sax + abs(sbend));
figure(3)%Plot the maximum and minimum stresses
clf
plot(del,smax,del,-smin)
xlabel('Displacement')
ylabel('Maximum Tensile and Compressive Stress')
legend('Tensile','Compressive','Location','SouthEast')

figure(4)%Plot the energy stored vs. displacement
clf
plot([0 del],U)
xlabel('Displacement')
ylabel('Energy')
figure(5)%Plot the deflection beam shapes
clf
plot(x,y)
title('Deflected Beam Shapes')