沖擊響應(yīng)譜計(jì)算的matlab程序

1、disp('')disp(' ver July 3, 2006')disp(' by Tom Irvine Email')disp(' ')disp(' This program calculates the shock response spectrum')disp(' of an acceleration time history, which is pre-loaded into Matlab.') disp(' The time history must have two colum

2、ns: time(sec) & acceleration') disp(' ') %clear t;clear y;clear yy;clear n;clear fn;clear a1;clear a2clear b1;clear b2;clear jnum;clear THM;clear resp;clear x_pos;clear x_neg; %iunit=input(' Enter acceleration unit: 1= G 2= m/secA2 ');%disp(' ')disp(' Select file

3、input method ');disp(' 1=external ASCII file ');disp(' 2=file preloaded into Matlab ');file_choice = input('');%if(file_choice=1)filename, pathname = uigetfile('*.*');filename = fullfile(pathname, filename);%fid = fopen(filename,'r');THM = fscanf(fid,'

4、%g %g',2 inf);THM=THM'elseTHM = input(' Enter the matrix name: ');end %t=double(THM(:,1);y=double(THM(:,2);%tmx=max(t);tmi=min(t);n = length(y);%out1 = sprintf('n %d samples n',n);disp(out1)%dt=(tmx-tmi)/(n-1);sr=1./dt;%out1 = sprintf(' SR = %g samples/sec dt = %g sec n&#

5、39;,sr,dt);disp(out1)%fn(1)=input(' Enter the starting frequency (Hz)');if fn(1)>sr/30.fn(1)=sr/30.;end%idamp=input(' Enter damping format: 1= damping ratio 2= Q '); %disp(' ')if(idamp=1)damp=input(' Enter damping ratio (typically ');elseQ=input(' Enter the amp

6、lification factor (typically Q=10) ');damp=1./(2.*Q);end%disp(' ')disp(' Select algorithm: ')disp(' 1=Kelly-Richman 2=Smallwood ');ialgorithm=input(' ');%tmax=(tmx-tmi) + 1./fn(1);limit = round( tmax/dt );n=limit;yy=zeros(1,limit);for i=1:length(y)yy(i)=y(i);end%d

7、isp(' ')disp(' Calculating response.')% SRS engine%for j=1:1000 %omega=2.*pi*fn(j);omegad=omega*sqrt(dampA2);cosd=cos(omegad*dt);sind=sin(omegad*dt);domegadt=damp*omega*dt;%if(ialgorithm=1)a1(j)=2.*exp(-domegadt)*cosd;a2(j)=-exp(-2.*domegadt);b1(j)=2.*domegadt;b2(j)=omega*dt*exp(-dom

8、egadt);b2(j)=b2(j)*( (omega/omegad)*.*(dampA2)*sind -2.*damp*cosd ); b3(j)=0; % elseE=exp(-damp*omega*dt);K=omegad*dt;C=E*cos(K);S=E*sin(K);Sp=S/K; % a1(j)=2*C;a2(j)=-EA2;b1(j)=;b2(j)=2.*(Sp-C);b3(j)=EA2-Sp;endforward= b1(j), b2(j), b3(j);back =1, -a1(j), -a2(j);%resp=filter(forward,back,yy);%x_pos(

9、j)= max(resp);x_neg(j)= min(resp);%jnum=j;if fn(j) > sr/8.breakendfn(j+1)=fn(1)*(2. A (j*(1./12.);end% Output options%disp('')disp(' Select output option ');choice=input(' 1=plot only 2=plot & output text file ');disp('')%if choice = 2%writefname, writepname =

10、uipufile('*','Save SRS data as'); writepfname = fullfile(writepname, writefname);writedata = fn' x_pos' (abs(x_neg)'fid = fopen(writepfname,'w');fprintf(fid,' %g %g %gn',writedata'); fclose(fid);% disp(' Enter output filename ');% SRS_filename

11、= input(' ','s');% fid = fopen(SRS_filename,'w');% for j=1:jnum% fprintf(fid,'% % % n',fn(j),x_pos(j),abs(x_neg(j);% end% fclose(fid);end% Plot SRS%disp(' ')disp(' Plotting output')% Find limits for plot%srs_max = max(x_pos);if max( abs(x_neg) ) > s

12、rs_max srs_max = max( abs(x_neg );endsrs_min = min(x_pos);if min( abs(x_neg) ) < srs_min srs_min = min( abs(x_neg );end%figure(1);plot(fn,x_pos,fn,abs(x_neg),'-.');%if iunit=1ylabel('Peak Accel (G)');elseylabel('Peak Accel (m/secA2)');endxlabel('Natural Frequency (Hz)&

13、#39;);Q=1./(2.*damp);out5 = sprintf(' Acceleration Shock Response Spectrum Q=%g ',Q);title(out5);grid;set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');legend ('positive','negativ

14、e',2);%ymax= 10A(round(log10(srs_max)+);ymin= 10A(round(log10(srs_min);%fmax=max(fn);fmin=fmax/10.;%fmax= 10A(round(log10(fmax)+);%if fn(1) >=fmin=;endif fn(1) >= 1fmin=1;endif fn(1) >= 10fmin=10;endif fn(1) >= 100fmin=100;endaxis(fmin,fmax,ymin,ymax);%disp(' ')disp(' Plo

15、t pseudo velocity?');vchoice=input(' 1=yes 2=no ' );if(vchoice=1)figure(2);% Convert to pseudo velocity%for j=1:jnumif iunit=1x_pos(j)=386*x_pos(j)/(2.*pi*fn(j);x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j);elsex_pos(j)=x_pos(j)/(2.*pi*fn(j);x_neg(j)=x_neg(j)/(2.*pi*fn(j);endend%srs_max = max(x_po

16、s);if max( abs(x_neg) ) > srs_maxsrs_max = max( abs(x_neg );endsrs_min = min(x_pos);if min( abs(x_neg) ) < srs_minsrs_min = min( abs(x_neg );end%plot(fn,x_pos,fn,abs(x_neg),'-.');%if iunit=1ylabel('Velocity (in/sec)');elseylabel('Velocity (m/sec)');endxlabel('Natural Frequency (Hz)');Q=1./(2.*damp);out5 = sprintf(' Pseudo Velocity Shock Response Spectrum Q=%g ',Q);title(out5);grid;set(gca,'MinorGridLineStyle