冲击响应谱计算的matlab程序

disp(' ') disp(' srs.m ver 2.0 July 3, 2006') disp(' by Tom Irvine Email: tomirvine@aol.com') 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 columns: time(sec) & acceleration') disp(' ') % clear t; clear y; clear yy; clear n; clear fn; clear a1; clear a2 clear b1; clear b2; clear jnum; clear THM; clear resp; clear x_pos; clear x_neg; % iunit=input(' Enter acceleration unit: 1= G 2= m/sec^2 '); % disp(' ') disp(' Select file 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,'%g %g',[2 inf]); THM=THM'; else THM = 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',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 0.05) '); else Q=input(' Enter the amplification 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 % disp(' ') disp(' Calculating response..... ') % % SRS engine % for j=1:1000 % omega=2.*pi*fn(j); omegad=omega*sqrt(1.-(damp^2)); 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(-domegadt); b2(j)=b2(j)*( (omega/omegad)*(1.-2.*(damp^2))*sind -2.*damp*cosd ); b3(j)=0; % else E=exp(-damp*omega*dt); K=omegad*dt; C=E*cos(K); S=E*sin(K); Sp=S/K; % a1(j)=2*C; a2(j)=-E^2; b1(j)=1.-Sp; b2(j)=2.*(Sp-C); b3(j)=E^2-Sp; end forward=[ b1(j), b2(j), b3(j) ]; back =[ 1, -a1(j), -a2(j) ]; % resp=filter(forward,back,yy); % x_pos(j)= max(resp); x_neg(j)= min(resp); % jnum=j; if fn(j) > sr/8. break end fn(j+1)=fn(1)*(2. ^ (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] = uiputfile('*','Save SRS data as'); writepfname = fullfile(writepname, writefname); writedata = [fn' x_pos' (abs(x_neg))' ]; fid = fopen(writepfname,'w'); fprintf(fid,' %g %g %g\n',writedata'); fclose(fid); %% % disp(' Enter output filename '); % SRS_filename = input(' ','s'); % % fid = fopen(SRS_filename,'w'); % for j=1:jnum % fprintf(fid,'%7.2f %10.3f %10.3f \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) ) > srs_max srs_max = max( abs(x_neg )); end srs_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==1 ylabel('Peak Accel (G)'); else ylabel('Peak Accel (m/sec^2)'); end xlabel('Natural Frequency (Hz)'); 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','negative',2); % ymax= 10^(round(log10(srs_max)+0.8)); ymin= 10^(round(log10(srs_min)-0.6)); % fmax=max(fn); fmin=fmax/10.; % fmax= 10^(round(log10(fmax)+0.5)); % if fn(1) >= 0.1 fmin=0.1; end if fn(1) >= 1 fmin=1; end if fn(1) >= 10 fmin=10; end if fn(1) >= 100 fmin=100; end axis([fmin,fmax,ymin,ymax]); % disp(' ') disp(' Plot pseudo velocity? '); vchoice=input(' 1=yes 2=no ' ); if(vchoice==1) figure(2); % % Convert to pseudo velocity % for j=1:jnum if iunit==1 x_pos(j)=386.*x_pos(j)/(2.*pi*fn(j)); x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j)); else x_pos(j)=x_pos(j)/(2.*pi*fn(j)); x_neg(j)=x_neg(j)/(2.*pi*fn(j)); end end % srs_max = max(x_pos); if max( abs(x_neg) ) > srs_max srs_max = max( abs(x_neg )); end srs_min = min(x_pos); if min( abs(x_neg) ) < srs_min srs_min = min( abs(x_neg )); end % plot(fn,x_pos,fn,abs(x_neg),'-.'); % if iunit==1 ylabel('Velocity (in/sec)'); else ylabel('Velocity (m/sec)'); end xlabel('Natural Frequency (Hz)'); Q=1./(2.*damp); out5 = sprintf(' Pseudo Velocity Shock Response Spectrum Q=%g ',Q); title(out5); grid; set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log'); legend ('positive','negative',2); % ymax= 10^(round(log10(srs_max)+0.8)); ymin= 10^(round(log10(srs_min)-0.6)); % axis([fmin,fmax,ymin,ymax]); end