% An elastic torsion bar fixed at the left end and free at the right end 
% with initial deformation proportional to the distance from the fixed end
% and no initial velocity. Rotational displacement at the free end is pi/2
% radians. Solution via the Fourier series solution to the wave equation.
% 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
set(0,'DefaultAxesFontSize',12);
set(0,'DefaultTextFontSize',12);
xoverL=linspace(0,1,101);
omega1t=linspace(0,4*pi,401);
cs=[];
nterm=40;%40 modes
omegaoveromega1=[1];
for i=1:nterm
    fac=((2*i)-1)*pi/2;
    sine(i,:)=sin(fac*xoverL);
    B(1,i)=2*(((1/fac)^2)*sin(fac));
end
for i=2:40
    omegaoveromega1=[omegaoveromega1 2*i-1];
end
for i=1:nterm
  cs1=cos(omegaoveromega1(1,i)*omega1t); 
  cs=[cs;cs1];
end
framedata=[];

for k=1:401
    sum=zeros(1,101);
    for i=1:nterm
        sum=sum+B(1,i)*cs(i,k)*sine(i,:);
    end
    framedata=[framedata;(pi/2)*sum];
end

figure(1);clf;
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Torsional Displacement,  \theta(x,t)')
S=plot(xoverL,framedata(1,:),'k','EraseMode','xor');
text(-.1,1.1,'Press Enter to Continue')
axis([-.2 1.2 -2 2])
box on
pause
for k=11:10:101
    set(S,'Ydata',framedata(k,:));
    axis([-.2 1.2 -2 2])
    text(-.1,1.1,'Press Enter to Continue')
    pause
end

figure(2);clf;
axis([0 4*pi -2 2])
hold on
text(5.2 ,1.57,'x/L=')
box on
for j=1:10:101
    plot(omega1t,framedata(:,j))
    hold on
    text(6,framedata(1,j),num2str(xoverL(1,j)))
    hold on
end
xlabel('Dimensionless Time, \omega_1t')
ylabel('Torsional Displacement,  \theta(x,t)')
text(4,-1,'Press Enter to Continue')
pause

%now do animation
figure(3);clf;
axis([-.2 1.2 -2 2])
box on
xlabel('Dimensionless Distance, x/L')
hold on
ylabel('Torsional Displacement,  \theta(x,t)')
hold on
L=plot(xoverL, zeros(1,101),'k','EraseMode','xor');
hold off
texthandl1=text(-.1,1.1,'Press Enter to Set Initial Condition');
pause
set(L,'Ydata',framedata(1,:))
set(texthandl1,'String','Press Enter to Animate');
pause
set(texthandl1,'String','                       ');
for i=2:401
    set(L,'Ydata',framedata(i,:))
    pause(.04)
end
hold on
set(texthandl1,'String','Press Enter to Continue');
pause

figure(4);clf;
R=.2;R1=.16;
phi=(pi/15);
xcir=[];ycir=[];
for i=1:61
    gamma=2*pi*i/60;
    xcir=[xcir R*cos(gamma)];
    ycir=[ycir R*sin(gamma)];
end
xellipse=xcir*sin(phi);
yellipse=ycir;
xfront=xellipse(1,16:46);
xback=[xellipse(1,46:61) xellipse(1,2:16)];
yfront=yellipse(1,16:46);
yback =[yellipse(1,46:61) yellipse(1,2:16)];
axis([-.2 1.2 -1 1]);
ypatch=[-.5 [-.2 fliplr(yfront) .2] .5 .3 -.1 -.35 -.5];
xpatch=[0 [0 fliplr(xfront) 0] 0 -.15 -.08 -.17 0];
patch(xpatch, ypatch,'r');
hold on
x1=linspace(0.1, 1,10); 
for i=1:9
    plot(xfront+x1(1,i),yfront,'-');
    hold on
    plot(xback+x1(1,i),yback,'--')
    hold on
end
plot(xfront+1,yfront,'-');
hold on
plot(xback+1,yback,'-')
hold on
xbar1=[0 1];xbar2=[1 0];
ybar1=[.2 .2];ybar2=[ -.2 -.2];
hold on
box on
plot(xbar1,ybar1);plot(xbar2,ybar2)
R1sp=R1*sin(phi);R1cp=R1*cos(phi);
for i=1:10
    xtic=x1(1,i);ytic=R1;
    tic(1,i)=plot([xtic xtic],[ytic,-ytic],'EraseMode','xor','LineWidth',[2]);
end
%pause
xtic1=zeros(401,10);
ytic1=zeros(401,10);
for i=1:10;
    for j=1:401;
        xtic1(j,i)=R1sp*sin(framedata(j,10*i+1));
        ytic1(j,i)=R1*cos(framedata(j,10*i+1));
    end
end
t1=text(.3,.8,'Press Enter to Set Initial Condition');
xlabel('Distance Along the Bar, x/L')
pause
for i=1:10
    set(tic(1,i), 'Xdata', x1(1,i)*ones(1,2)+[xtic1(1,i) -xtic1(1,i)],'Ydata',[ytic1(1,i),-ytic1(1,i)]);
end
set(t1,'String','Press Enter to Animate               ');
pause
set(t1,'String','                                      ');

for k=2:401
   for i=1:10
        set(tic(1,i), 'Xdata', x1(1,i)*ones(1,2)+[xtic1(k,i) -xtic1(k,i)],'Ydata',[ytic1(k,i),-ytic1(k,i)]);
    end
    pause(.04)
end
set(t1,'String', 'Press Enter to Continue')
pause