function [Fint,KT]=Beam2d2(Qa,Qb)
% function [Fint,KT]=Beam2d2(Qa,Qb)
% 
% 先行変位Qb,原稿変位Qaを与えてはりの要素力Fintと剛性マトリックスKTを計算する
% 
% Copyright (C) 2007 清水 信行, 曽我部 潔, 藤川 猛 All Rights Reserved.

% =================================================================================
% プログラムの使用条件
% 1. 本プログラムおよびデータを商用で，頒布してはならない．
% 2. 本プログラムおよびデータを著作者の許可なく営利目的で利用してはならない．
% 3. 本プログラムおよびデータの改変は，営利目的でない限り自由である．
% 4. 本プログラムおよびデータを使用することにより生じた損害などについては，著者，
%    プロクラム著作者，(社)日本機械学会，および(株)コロナ社は一切の責任を負わない．
% 5. 本プログラムの著作権は，本書編著者ならびにプログラム著作者に帰す．
% =================================================================================

global Ndof Numelm  Fint KT epsB LamdB L0 EA EI ;
    fi=zeros(3,1);fj=zeros(3,1); Fint=zeros(Ndof,1) ;KT=zeros(Ndof);
    Kg=zeros(6);Ko=zeros(6);KB=zeros(6);
    Lam6=zeros(6);Lamb=zeros(3);Lama=zeros(3);dLam=zeros(3,3);

dQ=Qa-Qb;                                                    %全体枠変位増分
for ke=1:Numelm                                              %要素数だけ反復
  Eps=sprintf(' beam No.= %d ',ke);
  is=(ke-1)*3;    js=ke*3;                                   %はり両端節点開始自由度
    Lamb=LamdB(:,:,ke);                                     %
  Qij=Qb(is+1:is+2)-Qb(js+1:js+2);  Lb=norm(Qij);            %先行長さ
  Qij=Qa(is+1:is+2)-Qa(js+1:js+2);  La=norm(Qij);            %現行長さ
  
    dQi=dQ(is+1:is+3);     dQj=dQ(js+1:js+3);
    dqi=Lamb*dQi ;         dqj=Lamb*dQj ;
    du =dqj(1)-dqi(1);     dw=dqj(2)-dqi(2);
%---------変形量--------        
        epsL=La-L0(ke)      ;                               
        epsi=dqi(3)-du/Lb +epsB(ke,1);
        epsj=dqj(3)-du/Lb +epsB(ke,2);
    
    Eps=sprintf(' %f  %f  %f  ',epsL,epsi,epsj);
    epsB(ke, 1)=epsi;     epsB(ke, 2)=epsj;     epsB(ke, 3)=epsL; %save
% -------座標変換マトリックス----------------------------
    [dLam]=dLmake(du,dw,Lb);
    Lama=dLam*Lamb;
    LamdB(:,:,ke)=Lama;   
    Lam6=[Lama zeros(3);zeros(3) Lama];
%----------要素力計算------------------------------------
  EAL=EA(ke)/L0(ke);
  EIL=EI(ke)/L0(ke);EIL2=EIL/Lb;EIL3=EIL2/La;
  
    fi(1)= (6*EIL2)*(epsi+epsj)  ;fj(1)=-fi(1);
    fi(2)=-(EAL)*epsL            ;fj(2)=-fi(2);
    fi(3)= (2*EIL)*(2*epsi+epsj) ;fj(3)=(2*EIL)*(epsi+2*epsj);
    Fi=Lama'*fi;    Fj=Lama'*fj;
% ---------要素力合成------------------------------------
    Fint(is+1:is+3)=Fint(is+1:is+3)+Fi;
    Fint(js+1:js+3)=Fint(js+1:js+3)+Fj;
    
% ----------剛性マトリックス計算--------------------------
        T0= EAL*epsL               ; TL=T0/La;  
        V0= 6*EIL2*(epsi+epsj)     ; VL=V0/La; 
    Kg(1,1)= TL; Kg(1,2)= VL; Kg(1,3)=0; Kg(1,4)=-TL; Kg(1,5)=-VL; Kg(1,6)=0;
    Kg(2,1)= VL; Kg(2,2)=  0; Kg(2,3)=0; Kg(2,4)=-VL; Kg(2,5)=  0; Kg(2,6)=0;
    Kg(3,1)=  0; Kg(3,2)=  0; Kg(3,3)=0; Kg(3,4)=  0; Kg(3,5)=  0; Kg(3,6)=0;
    Kg(4,1)=-TL; Kg(4,2)=-VL; Kg(4,3)=0; Kg(4,4)= TL; Kg(4,5)= VL; Kg(4,6)=0;
    Kg(5,1)=-VL; Kg(5,2)=  0; Kg(5,3)=0; Kg(5,4)= VL; Kg(5,5)=  0; Kg(5,6)=0;
    Kg(6,1)=  0; Kg(6,2)=  0; Kg(6,3)=0; Kg(6,4)=  0; Kg(6,5)=  0; Kg(6,6)=0;

Ko(1,1)= 12*EIL3;Ko(1,2)=   0; Ko(1,3)= 6*EIL2; Ko(1,4)=-12*EIL3; Ko(1,5)= 0;   Ko(1,6)=  6*EIL2;
Ko(2,1)=  0     ;Ko(2,2)= EAL; Ko(2,3)= 0     ; Ko(2,4)= 0      ; Ko(2,5)=-EAL; Ko(2,6)=  0     ;
Ko(3,1)=  6*EIL2;Ko(3,2)=   0; Ko(3,3)= 4*EIL ; Ko(3,4)=-6*EIL2 ; Ko(3,5)= 0;   Ko(3,6)=  2*EIL ;
Ko(4,1)=-12*EIL3;Ko(4,2)=   0; Ko(4,3)=-6*EIL2; Ko(4,4)= 12*EIL3; Ko(4,5)= 0;   Ko(4,6)= -6*EIL2;
Ko(5,1)=  0     ;Ko(5,2)=-EAL; Ko(5,3)= 0     ; Ko(5,4)= 0      ; Ko(5,5)= EAL; Ko(5,6)=  0     ;
Ko(6,1)=  6*EIL2;Ko(6,2)=   0; Ko(6,3)= 2*EIL ; Ko(6,4)=-6*EIL2 ; Ko(6,5)= 0;   Ko(6,6)=  4*EIL ;
    
% ---------剛性マトリックス合成---------------------------
    KB=Lam6'*(Ko+Kg)*Lam6;
    is1=is+1;is3=is+3;js1=js+1;js3=js+3;
    KT(is1:is3,is1:is3)=KT(is1:is3,is1:is3)+KB(1:3,1:3);
    KT(is1:is3,js1:js3)=KT(is1:is3,js1:js3)+KB(1:3,4:6);
    KT(js1:js3,is1:is3)=KT(js1:js3,is1:is3)+KB(4:6,1:3);
    KT(js1:js3,js1:js3)=KT(js1:js3,js1:js3)+KB(4:6,4:6);
end
%===================================================================================
