Tg=0.75; end if fz==3 Tg=0.90; end end %%%%%%%%%

ceita=zn; %%%%%阻尼比 lmt1=0.02+(0.05-ceita)/8; if lmt1<0 lmt1=0; end

lmt2=1+(0.05-ceita)/(0.06+1.7*ceita); if lmt2<0.55 lmt2=0.55; end

sjzs=0.9+(0.05-ceita)/(0.5+5*ceita); %%%%%分段位置 T1 T2 T3 T1=0.1; T2=Tg; T3=5*Tg; T_jg=2*pi./pl; %%%% 第一段 0～T1 if T_jg<=T1

arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg; end

%%%% 第二段 T1～T2 if T1

%%%% 第三段 T2~T3 if T2

arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax; end

%%%% 第四段 T3～6.0 if T3

arf_jg=(lmt2*0.2^sjzs-lmt1*(T_jg-5*Tg))*arfmax; end

%%%% 第五段 6.0～ if 6.0

arf_jg=(lmt2*0.2^sjzs-lmt1*(6.0-5*Tg))*arfmax; end

%%%%%%反应谱值 拟加速度值 rs_z=arf_jg*9.8; end

3 生成人造地震波主程序：

%%%主程序%%%%

%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M， 反应谱控制容差rc Tyz=[0.04:0.016:0.1,0.15:0.05:3.0,3.2:0.05:5.0]; rc=0.06; nTyz=length(Tyz);

ceita=0.035;%%%阻尼比：0.035 for i=1:nTyz

Syz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1); %%%%8度，2类场地，第1地震分组 end

%%%%%% 变换的频率差：2*pi*0.005(可以保证长周期项5s附近有5项三角级数)； %%%%频率变化范围 N1=30, 30*0.005*2*pi ;N2=3000, 5000*0.005*2*pi plc=2*pi*0.005;

pl=30*0.005*2*pi:0.005*2*pi:10000*0.005*2*pi; npl=length(pl); P=0.9; %%%保证率

%%%%%%人造地震动持续时间40s，时间间隔：0.02s Td=40; dt=0.02; t=0:0.02:40; nt=length(t); %%%%%%% 衰减包络函数 t1=8; %%%%上升段

t2=8+24; %%%%%平稳段; 下降段则为40－32＝8s c=0.6; %%%%衰减段参数 for i=1:nt if t(i)<=t1

f(i)=(t(i)/t1)^2; end

if t(i)>t1 & t(i)=t2

f(i)=exp(-c*(t(i)-t2)); end end

%%%%%%% 反应谱转换功率谱 for i=1:npl

Sw(i)=(2*ceita/(pi*pl(i)))*r_s_1(pl(i),ceita,8,2,1)^2/(-2*log(-1*pi*log(P)/(pl(i)*Td)));

Aw(i)=sqrt(4*Sw(i)*plc); end

%%%%%%%%%%%%%% 合成地震动 at=zeros(nt,1);atj=zeros(nt,1); for i=1:npl

fai(i)=rand(1)*2*pi; for j=1:nt

atj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i)))); end at=at+atj; end

%%%%%%% 计算反应谱验证是否满足rc在5%的要求,需要时程动力分析 %%%%%%%%%%%% response spectra of callidar %%%%%%% parameter g=9.8; m=1; x0=0; v0=0;

ww=2*pi./Tyz; %%%%%%%% load ag=at; %%%%%%%修改

%%%%%%% solution for y=1:nTyz z=0.037; w=ww(y); c=2*z*w; k=w^2; for i=1:nt-1

p(i)=-ag(i+1)+ag(i); a0=m\\(-ag(i)-c*v0-k*x0); kk=k+(dt^2)\\(6*m)+dt\\(3*c);

pp=p(i)+m*(dt\\(6*v0)+3*a0)+c*(3*v0+2\\(dt*a0)); dx=kk\\pp;

dv=dt\\(3*dx)-3*v0-2\\(dt*a0); x1=x0+dx; x0=x1; v1=v0+dv; v0=v1;

as(i)=a0;

as(i)=as(i)+ag(i); vs(i)=v0; xs(i)=x0; end

maxas(y)=max(as); maxvs(y)=max(vs); maxxs(y)=max(xs); end

for i=1:nTyz

rspa(i)=maxas(i); end

%%%%%%% 比较容差 for i=1:nTyz

rcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:)); end jsnum=1;

while max(rcrsp(:))>rc %%%%%循环体函数 blxs=Syz./rspa; for xsxs=1:npl

if 2*pi/pl(xsxs)

for sxsx=1:nTyz-1

if (2*pi/pl(xsxs)>=Tyz(sxsx)) & (2*pi/pl(xsxs)<=Tyz(sxsx+1))

blxs1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx))*(2*pi/pl(xsxs)-Tyz(sxsx))/(Tyz(sxsx+1)-Tyz(sxsx)); end end

if 2*pi/pl(xsxs)>Tyz(nTyz) blxs1(xsxs)=blxs(nTyz); end end

Aw=Aw.*blxs1;

%%%%%%%%%%%%%% 合成地震动 at=zeros(nt,1); atj=zeros(nt,1); for i=1:npl

for j=1:nt

atj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));