cuidazhe提供的Response Spectrum程序
% *********************************************** % *** Response Spectrum By Duhamel Integral & *** % *** Newmark Linear Acceleration Method *** % *** Author:Yahong DENG (PHD) *** % *** Email:hoverdyh@zju.edu.cn *** % *** Chang'an University *** % *** 2008.12.11 Xi'an *** % *********************************************** clc; clear all; close all hidden; disp('***********************************************') disp('*** Response Spectrum By Duhamel Integral & ***') disp('*** Linear Acceleration Method ***') disp('*** Author:Yahong DENG (PHD) ***') disp('*** Email:hoverdyh@zju.edu.cn ***') disp('*** Chang An University ***') disp('*** 2008.12.11 XiAn ***') disp('***********************************************') disp(' ') disp('-------------------- ') disp('** 振动信号输入 ** ') disp('-------------------- ') disp(' ') filenin=input('请输入数据文件名(*.txt/*.dat/*.*)：','s'); fileid=fopen(filenin,'r'); sampfre=input('请输入数据采样频率：'); linsign=input('请输入数据起始行：'); for line=1:linsign-1 tline=fgetl(fileid); end filesty=input('请选择数据列类型（1：全部(单列/多列) 2：多列中一列）：'); if filesty==1 vibsign=fscanf(fileid,'%f',inf); vibsign=vibsign'; else colnumb=input('请输入文件总列数：'); colsign=input('请输入数据所在列数：'); inpsign=fscanf(fileid,'%f',[colnumb,inf]); vibsign=inpsign(colsign,:); end signumb=length(vibsign); time=(0:1/sampfre:(signumb-1)/sampfre); maxsign=0; for n=1:signumb if abs(vibsign(n))>maxsign maxsign=abs(vibsign(n)); maxtime=time(n); end end disp(' '); fprintf('---* 输入信号基本信息 *---\n'); fprintf('输入信号长度为：%d\n',signumb); fprintf('输入信号持时为：%.3f (s)\n',(signumb-1)/sampfre); fprintf('输入信号最大值（绝对值）：%f\n',maxsign); fprintf('输入信号最大值时间点：%.3f (s)\n',maxtime); fprintf('--- 恭喜，成功读入数据！！！ ---\n'); filesta=fclose(fileid); disp(' ') disp('---------------------- ') disp('** 反应谱参数设置 ** ') disp('---------------------- ') disp(' ') damnumb=input('请输入反应谱阻尼比个数（1-3）：'); switch damnumb case 1 damprat(1)=input('请输入阻尼比值（0-1）：'); case 2 damprat(1)=input('请输入第一个阻尼比值（0-1）：'); damprat(2)=input('请输入第二个阻尼比值（0-1）：'); case 3 damprat(1)=input('请输入第一个阻尼比值（0-1）：'); damprat(2)=input('请输入第二个阻尼比值（0-1）：'); damprat(3)=input('请输入第三个阻尼比值（0-1）：'); end haxisty=input('请选择反应谱计算点类型（1：等间距频率计算点 2：等间距周期计算点）：'); switch haxisty case 1 minfreq=input('请输入最小计算频率（一般可取0.1，且 > 0）：'); maxfreq=input('请输入最大计算频率（一般可取10-15）：'); delfreq=input('请输入计算频率间隔（一般可取0.01/0.02）：'); spenumb=round((maxfreq-minfreq)/delfreq)+1; freoper=(minfreq:delfreq:(spenumb-1)*delfreq+minfreq); case 2 minperi=input('请输入最小计算周期（一般可取0.1，且 > 0）：'); maxperi=input('请输入最大计算周期（一般可取2-6）：'); delperi=input('请输入计算周期间隔（一般可取0.01/0.02）：'); spenumb=round((maxperi-minperi)/delperi)+1; freoper=(minperi:delperi:(spenumb-1)*delperi+minperi); end disp(' ') disp('-------------------- ') disp('** 反应谱计算 ** ') disp('-------------------- ') disp(' ') compmet=input('请选择反应谱计算方法（1：杜哈姆积分卷积法(速度慢) 2：杜哈姆积分FFT法 3：纽马克线性加速度法(速度快,推荐)）：'); oritime=cputime; switch compmet case 1 if haxisty==1 for n=1:spenumb angfreq=2*pi*freoper(n); for m=1:damnumb angfred(m)=angfreq*sqrt(1-damprat(m)^2); termexp=exp(-damprat(m)*angfreq*time); termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2)); termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2))); uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre; uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre; uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre; dispspe(m,n)=max(abs(conv(vibsign,uniresd))); velospe(m,n)=max(abs(conv(vibsign,uniresv))); accespe(m,n)=max(abs(conv(vibsign,uniresa))); accspeb(m,n)=accespe(m,n)/maxsign; fftnumb=2^nextpow2(2*signumb); uniafft=fft(uniresa,fftnumb); signfft=fft(vibsign,fftnumb); accespe1=real(ifft(uniafft.*signfft,fftnumb)); accespe2=accespe1(1:signumb); accespe3=max(accespe2); end fprintf('共 %4d 个频率计算点，正在计算第 %4d 个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100); end disp(' ') elatime=cputime-oritime; fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime); else for n=1:spenumb angfreq=2*pi/freoper(n); for m=1:damnumb angfred(m)=angfreq*sqrt(1-damprat(m)^2); termexp=exp(-damprat(m)*angfreq*time); termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2)); termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2))); uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre; uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre; uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre; dispspe(m,n)=max(abs(conv(vibsign,uniresd))); velospe(m,n)=max(abs(conv(vibsign,uniresv))); accespe(m,n)=max(abs(conv(vibsign,uniresa))); accspeb(m,n)=accespe(m,n)/maxsign; end fprintf('共 %4d 个周期计算点，正在计算第 %4d 个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100); end disp(' ') elatime=cputime-oritime; fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime); end case 2 fftnumb=2^nextpow2(2*signumb-1); if haxisty==1 for n=1:spenumb angfreq=2*pi*freoper(n); for m=1:damnumb angfred(m)=angfreq*sqrt(1-damprat(m)^2); termexp=exp(-damprat(m)*angfreq*time); termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2)); termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2))); uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre; uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre; uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre; fftunid=fft(uniresd,fftnumb); fftuniv=fft(uniresv,fftnumb); fftunia=fft(uniresa,fftnumb); fftsign=fft(vibsign,fftnumb); dispspt=real(ifft(fftunid.*fftsign,fftnumb)); accespe(m,n)=max(dispspt(1:signumb)); velospt=real(ifft(fftuniv.*fftsign,fftnumb)); velospe(m,n)=max(velospt(1:signumb)); accespt=real(ifft(fftunia.*fftsign,fftnumb)); accespe(m,n)=max(accespt(1:signumb)); accspeb(m,n)=accespe(m,n)/maxsign; end fprintf('共 %4d 个频率计算点，正在计算第 %4d 个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100); end disp(' ') elatime=cputime-oritime; fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime); else for n=1:spenumb angfreq=2*pi/freoper(n); for m=1:damnumb angfred(m)=angfreq*sqrt(1-damprat(m)^2); termexp=exp(-damprat(m)*angfreq*time); termsin=sin(time.*angfred(m)).*((1-2*damprat(m)^2)/(1-damprat(m)^2)); termcos=cos(time.*angfred(m)).*((2*damprat(m))/(sqrt(1-damprat(m)^2))); uniresd=-1/angfred(m)*termexp.*sin(time.*angfred(m))/sampfre; uniresv=-termexp.*(cos(time.*angfred(m))-damprat(m)*sin(time.*angfred(m))/sqrt(1-damprat(m)^2))/sampfre; uniresa=angfred(m)*termexp.*(termsin+termcos)/sampfre; fftunid=fft(uniresd,fftnumb); fftuniv=fft(uniresv,fftnumb); fftunia=fft(uniresa,fftnumb); fftsign=fft(vibsign,fftnumb); dispspt=real(ifft(fftunid.*fftsign,fftnumb)); accespe(m,n)=max(dispspt(1:signumb)); velospt=real(ifft(fftuniv.*fftsign,fftnumb)); velospe(m,n)=max(velospt(1:signumb)); accespt=real(ifft(fftunia.*fftsign,fftnumb)); accespe(m,n)=max(accespt(1:signumb)); accspeb(m,n)=accespe(m,n)/maxsign; end fprintf('共 %4d 个周期计算点，正在计算第 %4d 个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100); end disp(' ') elatime=cputime-oritime; fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime); end case 3 if haxisty==1 for n=1:spenumb angfreq=2*pi*freoper(n); deltime=1/sampfre; for m=1:damnumb a(m)=angfreq*sqrt(1-damprat(m)^2)*deltime; b(m)=damprat(m)*angfreq*deltime; d(m)=sqrt(1-damprat(m)^2); e(m)=exp(-b(m)); a11(m)=e(m)*(damprat(m)*sin(a(m))/d(m)+cos(a(m))); a12(m)=e(m)*sin(a(m))/angfreq/d(m); a21(m)=e(m)*(-angfreq)*sin(a(m))/d(m); a22(m)=e(m)*(-damprat(m)*sin(a(m))/d(m)+cos(a(m))); b11(m)=e(m)*((2*damprat(m)^2-1+b(m))*sin(a(m))/angfreq^2/a(m)+(2*damprat(m)+angfreq*deltime)*cos(a(m))/angfreq^3/deltime)... -2*damprat(m)/angfreq^3/deltime; b12(m)=e(m)*((1-2*damprat(m)^2)*sin(a(m))/angfreq^2/a(m)-2*damprat(m)*cos(a(m))/angfreq^3/deltime)... -1/angfreq^2+2*damprat(m)/angfreq^3/deltime; b21(m)=e(m)*((-damprat(m)-angfreq*deltime)*sin(a(m))/angfreq/a(m)-cos(a(m))/angfreq^2/deltime)+1/angfreq^2/deltime; b22(m)=e(m)*(damprat(m)*sin(a(m))/angfreq/a(m)+cos(a(m))/angfreq^2/deltime)-1/angfreq^2/deltime; tempdis(m,1)=0; tempvel(m,1)=0; tempacc(m,1)=0; for k=2:signumb tempdis(m,k)=a11(m)*tempdis(m,k-1)+a12(m)*tempvel(m,k-1)+b11(m)*vibsign(k-1)+b12(m)*vibsign(k); tempvel(m,k)=a21(m)*tempdis(m,k-1)+a22(m)*tempvel(m,k-1)+b21(m)*vibsign(k-1)+b22(m)*vibsign(k); tempacc(m,k)=-(2*damprat(m)*angfreq.*tempvel(m,k)+angfreq^2.*tempdis(m,k)); end dispspe(m,n)=max(abs(tempdis(m,:))); velospe(m,n)=max(abs(tempvel(m,:))); accespe(m,n)=max(abs(tempacc(m,:))); accspeb(m,n)=accespe(m,n)/maxsign; end fprintf('共 %4d 个频率计算点，正在计算第 %4d 个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100); end disp(' ') elatime=cputime-oritime; fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime); else for n=1:spenumb angfreq=2*pi/freoper(n); deltime=1/sampfre; for m=1:damnumb a(m)=angfreq*sqrt(1-damprat(m)^2)*deltime; b(m)=damprat(m)*angfreq*deltime; d(m)=sqrt(1-damprat(m)^2); e(m)=exp(-b(m)); a11(m)=e(m)*(damprat(m)*sin(a(m))/d(m)+cos(a(m))); a12(m)=e(m)*sin(a(m))/angfreq/d(m); a21(m)=e(m)*(-angfreq)*sin(a(m))/d(m); a22(m)=e(m)*(-damprat(m)*sin(a(m))/d(m)+cos(a(m))); b11(m)=e(m)*((2*damprat(m)^2-1+b(m))*sin(a(m))/angfreq^2/a(m)+(2*damprat(m)+angfreq*deltime)*cos(a(m))/angfreq^3/deltime)... -2*damprat(m)/angfreq^3/deltime; b12(m)=e(m)*((1-2*damprat(m)^2)*sin(a(m))/angfreq^2/a(m)-2*damprat(m)*cos(a(m))/angfreq^3/deltime)... -1/angfreq^2+2*damprat(m)/angfreq^3/deltime; b21(m)=e(m)*((-damprat(m)-angfreq*deltime)*sin(a(m))/angfreq/a(m)-cos(a(m))/angfreq^2/deltime)+1/angfreq^2/deltime; b22(m)=e(m)*(damprat(m)*sin(a(m))/angfreq/a(m)+cos(a(m))/angfreq^2/deltime)-1/angfreq^2/deltime; tempdis(m,1)=0; tempvel(m,1)=0; tempacc(m,1)=0; for k=2:signumb tempdis(m,k)=a11(m)*tempdis(m,k-1)+a12(m)*tempvel(m,k-1)+b11(m)*vibsign(k-1)+b12(m)*vibsign(k); tempvel(m,k)=a21(m)*tempdis(m,k-1)+a22(m)*tempvel(m,k-1)+b21(m)*vibsign(k-1)+b22(m)*vibsign(k); tempacc(m,k)=-(2*damprat(m)*angfreq*tempvel(m,k)+angfreq^2*tempdis(m,k)); end dispspe(m,n)=max(abs(tempdis(m,:))); velospe(m,n)=max(abs(tempvel(m,:))); accespe(m,n)=max(abs(tempacc(m,:))); accspeb(m,n)=accespe(m,n)/maxsign; end fprintf('共 %4d 个周期计算点，正在计算第 %4d 个点，任务完成率： %6.2f %%，请稍等！！！\n',spenumb,n,n/spenumb*100); end disp(' ') elatime=cputime-oritime; fprintf('--- 恭喜，位移/速度/加速度反应谱计算完成，共耗时： %6.2f 秒 ！！！---\n',elatime); end end disp(' ') disp('------------------------ ') disp('** 数据 & 图形 输出 ** ') disp('------------------------ ') disp(' ') tempdat=date; datyorn=input('请选择数据输出选项（0：不输出 1：仅输出位移谱 2：仅输出速度谱 3：仅输出加速度谱(含β谱) 4：输出所有谱）：'); if datyorn==0 disp('--- 没有生成数据文件！！！---') elseif datyorn==1 fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s'); fileid=fopen(fnout,'wt'); switch damnumb case 1 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'%8s\t%15s\n','# F/T','# D_Spectrum'); fprintf(fileid,'* ------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\n',freoper(n),dispspe(1,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 2 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* --------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2'); fprintf(fileid,'* --------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 3 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* -------------------------------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2','# D_Spectrum_3'); fprintf(fileid,'* -------------------------------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n),dispspe(3,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); end elseif datyorn==2 fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s'); fileid=fopen(fnout,'wt'); switch damnumb case 1 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'%8s\t%15s\n','# F/T','# V_Spectrum'); fprintf(fileid,'* ------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\n',freoper(n),velospe(1,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 2 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* ---------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# V_Spectrum_1','# V_Spectrum_2'); fprintf(fileid,'* ---------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),velospe(1,n),velospe(2,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 3 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* --------------------------------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# V_Spectrum_1','# V_Spectrum_2','# V_Spectrum_3'); fprintf(fileid,'* --------------------------------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),velospe(1,n),velospe(2,n),velospe(3,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); end elseif datyorn==3 fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s'); fileid=fopen(fnout,'wt'); switch damnumb case 1 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* -----------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\n','# F/T','# A_Spectrum','# A_Spectrum_β'); fprintf(fileid,'* -----------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\n',freoper(n),accespe(1,n),accspeb(1,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 2 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* ------------------------------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# A_Spectrum_1','# A_Spectrum_1β','# A_Spectrum_2','# A_Spectrum_2β'); fprintf(fileid,'* ------------------------------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),accespe(1,n),accspeb(1,n),accespe(2,n),accspeb(2,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 3 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* ---------------------------------------------------------------------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# A_Spectrum_1','# A_Spectrum_1β','# A_Spectrum_2','# A_Spectrum_2β',... '# A_Spectrum_3','# A_Spectrum_3β'); fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),accespe(1,n),accspeb(1,n),... accespe(2,n),accspeb(2,n),accespe(3,n),accspeb(3,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); end else fnout=input('请输入输出数据文件名(*.txt/*.dat/*.*)：','s'); fileid=fopen(fnout,'wt'); switch damnumb case 1 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* --------------------------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum','# V_Spectrum','# A_Spectrum'); fprintf(fileid,'* --------------------------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),velospe(1,n),accespe(1,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 2 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2','# V_Spectrum_1','# V_Spectrum_2','# A_Spectrum_1','# A_Spectrum_2'); fprintf(fileid,'* ----------------------------------------------------------------------------------------------------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n),velospe(1,n),velospe(2,n),accespe(1,n),accespe(2,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); case 3 fprintf(fileid,'* ------------------------------------\n'); fprintf(fileid,'* OUTPUT FILE BY "RESPONSE_SPECTRUM.M"\n'); fprintf(fileid,'* YAHONG DENG (PHD)\n'); fprintf(fileid,'* CHANG_AN UNIVERSITY\n'); fprintf(fileid,'* %s\n',tempdat); fprintf(fileid,'* DATA\n'); fprintf(fileid,'* -------------------------------------------------------------------------------------------------------------------------------------------------------------\n'); fprintf(fileid,'%8s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\t%15s\n','# F/T','# D_Spectrum_1','# D_Spectrum_2','# D_Spectrum_3',... '# V_Spectrum_1','# V_Spectrum_2','# V_Spectrum_3','# A_Spectrum_1','# A_Spectrum_2','# A_Spectrum_3'); fprintf(fileid,'* -------------------------------------------------------------------------------------------------------------------------------------------------------------\n'); for n=1:spenumb fprintf(fileid,'%8.4f\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\t%15.6E\n',freoper(n),dispspe(1,n),dispspe(2,n),dispspe(3,n),... velospe(1,n),velospe(2,n),velospe(3,n),accespe(1,n),accespe(2,n),accespe(3,n)); end fprintf(fileid,'*DATAEND '); disp('--- 成功生成数据文件！！！---'); filesta=fclose(fileid); end end disp('') figyorn=input('请选择图形显示选项（0：不输出 1：仅输出位移谱 2：仅输出速度谱 3：仅输出加速度谱(含β谱) 4：输出所有谱）：'); if figyorn==0 disp(' ') disp('*******************************') disp('*** 程序结束，谢谢使用！ ***') disp('*******************************') disp(' ') elseif figyorn==1 for m=1:damnumb plot(freoper,dispspe(m,:)); grid on; hold on; end xlabel('频率(Hz)/周期(Sec)'); ylabel('幅值'); title('位移反应谱曲线'); disp('--- 成功生成曲线图形！！！---'); disp(' ') disp('*******************************') disp('*** 程序结束，谢谢使用！ ***') disp('*******************************') disp(' ') elseif figyorn==2 for m=1:damnumb plot(freoper,velospe(m,:)); grid on; hold on; end xlabel('频率(Hz)/周期(Sec)'); ylabel('幅值'); title('速度反应谱曲线'); disp('--- 成功生成曲线图形！！！---'); disp(' ') disp('*******************************') disp('*** 程序结束，谢谢使用！ ***') disp('*******************************') disp(' ') elseif figyorn==3 subplot(2,1,1) for m=1:damnumb plot(freoper,accespe(m,:),'r'); grid on; hold on; end xlabel('频率(Hz)/周期(Sec)'); ylabel('幅值'); title('加速度反应谱曲线'); legend('幅值谱 '); subplot(2,1,2) for m=1:damnumb plot(freoper,accspeb(m,:)); grid on; hold on; end xlabel('频率(Hz)/周期(Sec)'); ylabel('动力系数(β)'); legend('动力放大系数(β)谱 '); disp('--- 成功生成曲线图形！！！---'); disp(' ') disp('*******************************') disp('*** 程序结束，谢谢使用！ ***') disp('*******************************') disp(' ') else subplot(3,1,1); for m=1:damnumb plot(freoper,dispspe(m,:)); grid on; hold on; end xlabel('频率(Hz)/周期(Sec)'); ylabel('位移谱值'); title('位移/速度/加速度反应谱曲线'); subplot(3,1,2) for m=1:damnumb plot(freoper,velospe(m,:)); grid on; hold on; end xlabel('频率(Hz)/周期(Sec)'); ylabel('速度谱值'); subplot(3,1,3) for m=1:damnumb plot(freoper,accespe(m,:),'r'); grid on; hold on; end xlabel('频率(Hz)/周期(Sec)'); ylabel('加速度谱值'); disp('--- 成功生成曲线图形！！！---'); disp(' ') disp('*******************************') disp('*** 程序结束，谢谢使用！ ***') disp('*******************************') disp(' ') end |
GMT+8, 2025-3-13 09:27 , Processed in 0.037603 second(s), 17 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.
