clc disp(' This script calculates the undamped anddamped shock spectrum '); disp(' for the six most common pulse waveforms, namely a rectangular pulse,'); disp(' a half-sine pulse, a full cycle sine pulse, a symmetric triangular'); disp(' pulse and two nonsymmetric triangular pulses.'); disp(' '); disp(' The following force pulse waveforms can be treated:'); disp(' 1--Rectangular Pulse'); disp(' 2--Half Cycle Sine Pulse'); disp(' 3--Full Cycle Sine Pulse'); disp(' 4--Symmetric Triangular Pulse'); disp(' 5--Triangular Pulse, Ramp Rise, Step Fall'); disp(' 6--Triangular Pulse, Step Rise, Ramp Fall'); disp(' '); disp(' Input the number of the pulse you desire (1, 2, 3, 4, 5 or 6):'); disp(' '); ptype=[]; while isempty(ptype)|(ptype~=1 & ptype~=2 & ptype~=3 & ptype~=4 & ptype~=5 & ptype~=6) ptype=input(' Pulse Number = '); end dt=0.01;% scaled time increment wn=2*pi; p=zeros(1,500); tau=[0:0.01:5];%t/Tn taud=[0.1:0.1:5];%td/Tn zeta1=[0 .02 .1 .2 .4]; Rd=zeros(5,50); for k=1:5 z=zeta1(1,k); wD=wn*sqrt(1-z*z); Co=cos(wD*dt);Si=sin(wD*dt); E=exp(-z*wn*dt);H=z/sqrt(1-z*z);J=sqrt(1-z*z); PHI=zeros(2,2); % reserve space GAMMA=zeros(2,2);%reserve space A=E*(Co+H*Si);B=E*Si/wD;Ap=-E*((wn/J)*Si);Bp=E*(Co-H*Si); PHI=[A B;Ap Bp];%fill the transition matrix C=(((2*z)/(wn*dt))+E*((((1-2*z*z)/(wD*dt))-H)*Si-(1+((2*z)/(wn*dt)))*Co)); D=(1-((2*z)/(wn*dt))+E*(((2*z*z-1)/(wD*dt))*Si+((2*z/(wn*dt)))*Co)); Cp=(-(1/dt)+E*(((wn/J)+(H/dt))*Si+(1/dt)*Co)); Dp=(1/(dt))*(1-E*(H*Si+Co)); GAMMA=[C D;Cp Dp];%fill the input distribution matrix for i=1:50 id=fix((taud(1,i)+.001)/.01); tau1=[0:.01:id*.01]; if ptype==1;%rectangular pulse p=[ones(1,id),zeros(1,501-id)]; text1='Rectangular Pulse Shock Spectrum'; elseif ptype==2;%half cycle sine pulse p=[sin(pi*tau1/taud(1,i)) zeros(1,501-id-1)]; text1='Half Cycle Sine Pulse Shock Spectrum'; elseif ptype==3;%full cycle sine pulse p=[sin(2*pi*tau1/taud(1,i)) zeros(1,501-id-1)]; text1='One Cycle Sine Pulse Shock Spectrum'; elseif ptype==4;% symmetric triangular pulse ramp=[0:1:(id/2)]/(id/2); p=[ramp fliplr(ramp(1,1:((id/2)))) zeros(1,500-id)]; text1='Symmetric Triangular Pulse Shock Spectrum'; elseif ptype==5;%triangular pulse, ramp up, step down ramp=[0:1:id]/id; p=[ramp zeros(1,500-id)]; text1='Triangular Pulse,Ramp Rise, Step Fall Shock Spectrum'; elseif ptype==6;%triangular pulse, step up, ramp down ramp=ones(1,id+1)-([0:1:id]/id); p=[ramp zeros(1,500-id)]; text1='Triangular Pulse, Step Rise, Ramp Fall Shock Spectrum'; end %Solution generated by using Chopra's Interpolation of Excitation x=[0;0];X=[x]; for j=1:500 x=PHI*x+GAMMA*[p(1,j);p(1,j+1)]; X=[X x]; end Umax=max(max(abs(X(1,:)))); Rd(k,i+1)=Umax; if i==20&&k==1 figure(1);clf Pmax=max(p);Pmin=min(p); plot([-.2 0 tau],[0 0 p]) axis([-.2 5 Pmin-.1 Pmax+.1]) text(3,.8,'t_d/T_n = 2') xlabel('Time t/T_n') ylabel('Forcing Function p(t)/p_0') text(3,.6,'Press Enter to Continue') grid on pause end end end figure(2);clf plot([0 taud],Rd(1,:));hold on plot([0 taud],Rd(2,:),'-.');hold on plot([0 taud],Rd(3,:),'-+');hold on plot([0 taud],Rd(4,:),'-x');hold on plot([0 taud],Rd(5,:),'-o');hold on legend('\zeta = 0','\zeta = 0.02','\zeta = 0.1','\zeta = 0.2','\zeta = 0.4') xlabel('Dimensionless Pulse,Duration, t_d/T_n'); ylabel('R_d = u(t)_m_a_x/(u_s_t)_0') grid on text(1.2,.3,text1)