% movingload2.m  Last modified 7/29/2008
% Animates the response of a simply-supported Bernoulli-Euler beam 
% with constant moving load (force only) Mg starting at the left end at t=0. 
% Input is the ratio of transit time to the first natural period of
% the beam. The load is a constant force only without mass.
% This m-file was written at the University of Wyoming in the Electrical
% and Computer Engineering Department and is to be distributed without
% cost.
clear all
string1='The ratio of transit time Tt to first natural period T1, Tt/T1 = ';
nmode=5;
Tratio=input(string1);
sine=[];
omegaoveromega1=0.5/(Tratio);
omegajoveromega1=zeros(1,10);
xoverL=linspace(0,1,101);
xoverL1=linspace(0,4,401);
factor=96/(pi^4);
for j=1:nmode
    sine(j,:)=sin(j*pi*xoverL);
    omegajoveromega1(1,j)=j^2;
end
xarrow=[0 0 .02 -.02 0];
yarrow=[-.5 -.01 -.1 -.1 -.01];

yoveryst=zeros(1,101);
omegat=linspace(0,4*pi,401);
for k=1:101%time loop
    sum=zeros(1,101);
    for j=1:nmode% series loop
        if j*omegaoveromega1==omegajoveromega1(1,j)
            term=0.5*(1/(j*omegaoveromega1))*((1/(j*j))*sin(j*omegat(1,k))-(1/omegaoveromega1)*omegat(1,k).*cos(j*omegat(1,k)))*sine(j,:);
        else
            num=(sin(j*omegat(1,k))-(j*omegaoveromega1/omegajoveromega1(1,j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*omegat(1,k)));
            term=(num/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2))*sine(j,:);
        end
        sum=sum+term;
    end
    yoveryst(k,:)=factor*sum;
end
y1=3*xoverL(1,1:51)-(4*(xoverL(1,1:51).^3));
yst=[y1 fliplr(y1(1,1:50))];
% evaluate initial conditions for unforced motion q and qprime are column
% vectors
omegat1=pi;
for j=1:nmode
    if j*omegaoveromega1==omegajoveromega1(1,j);
        q(j,1)=0.5*(1/(j*omegaoveromega1))*((1/(j*j))*sin(j*omegat1)-(1/omegaoveromega1)*omegat1*cos(j*omegat1));
        qprime(j,1)=0.5*(1/(j*omegaoveromega1))*(((1/j)-(1/omegaoveromega1))*cos(omegat1)+(j/(omegaoveromega1))*sin(omegat1));
    else
        num1=(sin(j*omegat1)-(j*omegaoveromega1/omegajoveromega1(1,j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*omegat1));
        q(j,1)=factor*(num1/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2));
        num2=j*(cos(j*omegat1)-cos((omegajoveromega1(1,j)/omegaoveromega1)*omegat1));
        qprime(j,1)=factor*num2/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2);
    end
    
end
for k=102:401
    sum=zeros(1,101);
    for j=1:nmode
        sum=sum+(q(j,1)*cos((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi))+qprime(j,1)*(omegaoveromega1/(j*j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi)))*sine(j,:);
    end
    yoveryst=[yoveryst;sum];
end
ytip=[];
for k=1:101
    ytip=[ytip yoveryst(k,k)];
end
ytip=[ytip zeros(1,300)];        
      
    
figure(1);clf;%Animate the beam deflection as load moves across the beam
h1=axes('Ydir','reverse');
axis([-.1 1.1 -2  2])
hold on
plot(xoverL, yst);
hold on
plot([0 1],[0 0],'o');
hold on
plot([-.1 0],[0 0]);
hold on;
plot([1 1.1],[0 0]);
hold on
text(-.05,1.7,['T_t/T_1 = ' num2str(Tratio)]);
hold on
A=plot(xarrow,yarrow,'LineWidth',[2],'EraseMode','xor');
Y=plot(xoverL,yoveryst(1,:),'Color','k','LineWidth',[3],'EraseMode','xor');
xlabel('Dimensionless Position, x/L')
hold on
ylabel('Deflection, 48EIy(x,t)/MgL^3')
hold on
box on
Maximum=max(max(yoveryst));
text(.5,1.7,['y_m_a_x/y_s_t_a_t_i_c_m_a_x = ' num2str(Maximum)]);
hold on
t1=text(-.05,-1.7,'Press Enter to Continue');
pause
set(t1,'String','                        ');
for k=1:401
    set(A,'Xdata',xarrow+xoverL1(1,k)*ones(1,5),'Ydata',yarrow+ytip(1,k)*ones(1,5));
    set(Y,'Ydata',yoveryst(k,:));
    pause(.06)
end
hold off
set(t1,'String','Press Enter to Continue');
pause

MoverMst=[];
for k=1:101%time loop
    sum1=zeros(1,101);
    for j=1:5% series loop
        if j*omegaoveromega1==omegajoveromega1(1,j)
            term1=0.5*j*j(1/(j*omegaoveromega1))*((1/(j*j))*sin(j*omegat(1,k))-(1/omegaoveromega1)*omegat(1,k).*cos(j*omegat(1,k)))*sine(j,:);
        else
            num=(sin(j*omegat(1,k))-(j*omegaoveromega1/omegajoveromega1(1,j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*omegat(1,k)));
            term1=(j*j*num/((omegajoveromega1(1,j)^2)-(j*omegaoveromega1)^2))*sine(j,:);
        end
        sum1=sum1+term1;
    end
    MoverMst=[MoverMst;(8/(pi*pi))*sum1];
end

for k=102:401
    sum1=zeros(1,101);
    for j=1:5
        sum1=sum1+(q(j,1)*cos((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi))+qprime(j,1)*(omegaoveromega1/(j*j))*sin((omegajoveromega1(1,j)/omegaoveromega1)*(omegat(1,k)-pi)))*sine(j,:);
    end
    MoverMst=[MoverMst;sum1*(8/(pi*pi))];
end
aoverL=xoverL;
Mstatic=[zeros(1,101)];
for k=2:100
     m1=4*(1-aoverL(1,k))*xoverL(1,1:k);
     m2=4*(aoverL(1,k)*ones(1,101-k)-aoverL(1,k)*xoverL(1,k+1:101));
     Mstatic=[Mstatic;m1 m2];
end
Mstatic=[Mstatic;zeros(1,101)];
MaxMoment=max(max(MoverMst));



figure(4);clf;
axis([-.1 1.1 -2  2])
hold on
plot([-.1 0],[0 0]);
hold on;
plot([1 1.1],[0 0]);
hold on
text(-.05,1.7,['T_t/T_1 = ' num2str(Tratio)]);
hold on
box on
text(.5,1.7,['M_m_a_x/M_s_t_a_t_i_c_m_a_x = ' num2str(MaxMoment)])
B=plot([0 0],[-.5 -.1],'LineWidth',[1],'EraseMode','xor');
Y=plot(xoverL,MoverMst(1,:),'Color','k','LineWidth',[1],'EraseMode','xor');

MM=plot(xoverL,Mstatic(1,:),'EraseMode','xor');
xlabel('Dimensionless Position, x/L')
hold on
ylabel('Bending Moment, 4M(x,t)/MgL')
hold on
t1=text(-.05,-1.7,'Press Enter to Continue');
pause
set(t1,'String','                        ');
for k=1:401
    set(B,'Xdata',xoverL1(1,k)*ones(1,2));
    set(Y,'Ydata',MoverMst(k,:));
    if k>101
        set(MM,'Ydata',zeros(1,101))
    else
        set(MM,'Ydata',Mstatic(k,:))
    end
    pause(.04)
end
hold off
set(t1,'String','Press Enter to Continue');
pause;