function coupling
format short e
%
%-------------------------选择时间步长Δt、参数β和γ，计算积分常数
beta=0.25;                    %β=0.25
gamma=0.5;                    %γ=0.5
dt=0.2;                       %间步长Δt=0.2（秒）
a0=1/(beta*dt^2);
a1=gamma/(beta*dt);
a2=1/(beta*dt);
a3=1/(2*beta)-1;
a4=gamma/beta-1;
a5=0.5*dt*(gamma/beta-2);
a6=(1-gamma)*dt;
a7=gamma*dt;
%-------------------------计算车辆和桥梁的刚度矩阵、质量矩阵、阻尼矩阵
[Kv,Mv,Cv,Kt,Ct,W]=vehicle;
[Kb,Mb,Cb,LOC]=bridge;
%[Kb,Mb,Cb,LOC]=bridgeplane;
%-------------------------计算车辆和桥梁的等效刚度矩阵Kb和Kv
Kb_=Kb+a0*Mb+a1*Cb;
Kv_=Kv+a0*Mv+a1*Cv;
%-------------------------取初始值
Ub=ones(size(Kb,1),1);               %桥梁初始位移
Ub1=ones(size(Kb,1),1);              %桥梁初始速度
Ub2=ones(size(Kb,1),1);              %桥梁初始加速度
Uv=ones(size(Kv,1),1);               %车辆初始位移
Uv1=ones(size(Kv,1),1);              %车辆初始速度
Uv2=ones(size(Kv,1),1);              %车辆初始加速度
ub(:,1)=Ub;
ub1(:,1)=Ub1;
ub2(:,1)=Ub2;
uv(:,1)=Uv;
uv1(:,1)=Uv1;
uv2(:,1)=Uv2;
L=4;%桥的长度
v=20*1000/3600;        %车速km/h
x=2;

i=1;
t(1)=0;
I=1;
%t(i)*v<=L
%--------------------开始一个newmark时间步的迭代----------------
while I<=2;
    k=0;%迭代次数
    delta=0.03;
    %--------------计算中间变量
    Sb=Mb*(a0*ub(:,i)+a2*ub1(:,i)+a3*ub2(:,i))+Cb*(a1*ub(:,i)+a4*ub1(:,i)+a5*ub2(:,i));
    Sv=Mv*(a0*uv(:,i)+a2*uv1(:,i)+a3*uv2(:,i))+Cv*(a1*uv(:,i)+a4*uv1(:,i)+a5*uv2(:,i));
    %--------------计算路面不平度
    for j=1:2
        [r,r1]=roadsurface(x,16);
        rb(j)=r;
        rb1(j)=r1;
        x=x+v*dt;        
    end
   
    %-----------------------开始迭代
    while delta>0.02             %判断是否满足精度要求
          %-------------桥梁在t+Δt时刻的加速度、速度
          if(k==0)
              ub2(:,i+1)=a0*(ub(:,i)-ub(:,i))-a2*ub1(:,i)-a3*ub2(:,i)
              ub1(:,i+1)=ub1(:,i)+a6*ub2(:,i)+a7*ub2(:,i+1)
             
          else
              ub2(:,i+1)=a0*(ub(:,i+1)-ub(:,i))-a2*ub1(:,i)-a3*ub2(:,i);
              ub1(:,i+1)=ub1(:,i)+a6*ub2(:,i)+a7*ub2(:,i+1);
             
          end
  
          %--------------形成车辆方程右端荷载项
          %za=zeros(4,1);
          za(1)=ub(3*LOC(6,I),i)+rb(2);
          za(2)=ub(3*LOC(7,I),i)+rb(2);
          za(3)=ub(3*LOC(5,I),i)+rb(1);
          za(4)=ub(3*LOC(8,I),i)+rb(1);
          za=reshape(za,4,1);
          %za1=zeros(4,1);
          za1(1)=ub1(3*LOC(6,I),i)+rb1(2);
          za1(2)=ub1(3*LOC(7,I),i)+rb1(2);
          za1(3)=ub1(3*LOC(5,I),i)+rb1(1);
          za1(4)=ub1(3*LOC(8,I),i)+rb1(1);
          za1=reshape(za1,4,1);
         
          %----------------车辆等效荷载列阵
          Pv=Kt*za+Ct*za1;
          Pv_=Pv+Sv;               
         
          %-------------求解t+Δt时刻的车辆位移列阵
         
          uv(:,i+1)=inv(Kv_)*Pv_;
         
          %-------------车辆在t+Δt时刻的加速度、速度
          uv2(:,i+1)=a0*(uv(:,i+1)-uv(:,i))-a2*uv1(:,i)-a3*uv2(:,i);
          uv1(:,i+1)=uv1(:,i)+a6*uv2(:,i)+a7*uv2(:,i+1);
          %-------------计算桥梁方程右端荷载项
          Pb=zeros(size(Kb,1),1);
          for j=i:4
              Pb0(j)=W(j)+Kt(3+j,j)*(uv(j+3,i+1)-za(j,1))+Ct(3+j,j)*(uv1(j+3,i+1)-za1(j,1));
          end
          Pb(3*LOC(6,I),1)=Pb0(1);
          Pb(3*LOC(7,I),1)=Pb0(2);
          Pb(3*LOC(5,I),1)=Pb0(3);
          Pb(3*LOC(8,I),1)=Pb0(4);
          Pb_=Pb+Sb;               %桥梁等效荷载列阵
         
          %--------------求解t+Δt时刻的桥梁位移列阵
          ub(:,i+1)=inv(Kb_)*Pb_;
         
         
          %--------------桥梁在t+Δt时刻的加速度、速度
          ub2(:,i+1)=a0*(ub(:,i+1)-ub(:,i))-a2*ub1(:,i)-a3*ub2(:,i);
          ub1(:,i+1)=ub1(:,i)+a6*ub2(:,i)+a7*ub2(:,i+1);
         
          %---------------------判断桥梁位移是否收敛-----------------
          u1=ub(:,i);
          u2=ub(:,i+1);
         
          delta=norm((u2-u1)/(u2+eps))
          k=k+1
    end
   
    %--------------------进入下一时刻

    i=i+1;
    t(i)=t(i-1)+dt;
    x=v*dt;
    I=I+1;
end