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

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

global Ndof Numelm Fint KT L0 EA0 EI0 alpha0;
Fint=zeros(Ndof,1) ;KT=zeros(Ndof,Ndof);

for ke=1:Numelm                              %要素数だけ反復
  Eps=sprintf(' %d  new itereation  ',ke);
  is=(ke-1)*3;  js=ke*3;                     %はり両端の開始自由度
     [alpha,co,si,La]=Lmake(Qa(is+1),Qa(is+2),Qa(js+1),Qa(js+2));
%---------変形量--------        
  epsL=La-L0(ke);
  epsi=alpha0(ke,1)+Qa(is+3)-alpha;
  epsj=alpha0(ke,2)+Qa(js+3)-alpha;     
%----------要素力計算------------------------------------
  EAL=EA0(ke)/L0(ke);
  EIL=EI0(ke)/L0(ke);EIL2=EIL/La;EIL3=EIL2/La;
  
    Fi(1)=-(EAL)*epsL*co -(6*EIL2)*(epsi+epsj)*si;
    Fi(2)=-(EAL)*epsL*si +(6*EIL2)*(epsi+epsj)*co;
    Fi(3)= 2*EIL*(2*epsi+epsj);
    Fj(1)= (EAL)*epsL*co +(6*EIL2)*(epsi+epsj)*si;
    Fj(2)= (EAL)*epsL*si -(6*EIL2)*(epsi+epsj)*co;
    Fj(3)= 2*EIL*(epsi+2*epsj);
% ---------要素力合成------------------------------------
for I=1:3    
    Fint(is+I)=Fint(is+I)+Fi(I);
    Fint(js+I)=Fint(js+I)+Fj(I);
end
% ----------剛性マトリックス計算--------------------------
cc=co*co;ss=si*si;cs=co*si;sc=si*co;L=La;
        T0= EAL*epsL               ; TL=T0/La;  
        V0=-6*EIL2*(epsi+epsj)     ; VL=V0/La; 
K1=[cc sc 0;cs ss 0;0 0 0];K2=[-cc -cs 0;-sc -ss 0;0 0 0];
    K01=EAL*[K1 K2;K2' K1];
K1=[12*ss -12*sc -6*La*si;-12*cs  12*cc  6*La*co;-6*La*si  6*La*co 4*La*La];
K2=[-12*ss 12*sc -6*La*si; 12*cs -12*cc  6*La*co; 6*La*si -6*La*co 2*La*La];
K3=[12*ss -12*sc  6*La*si;-12*sc  12*cc -6*La*co; 6*La*si -6*La*co 4*La*La];
    K02=EIL3*[K1 K2;K2' K3];
K1=[ss -sc 0; -sc cc 0;0 0 0];K2=[-ss sc 0; sc -cc 0;0 0 0]; 
    Kg1=TL*[K1 K2;K2' K1];
K1=[2*sc ss-cc 0; ss-cc -2*sc 0 ;0 0 0];
    Kg2=-VL*[K1 -K1;-K1' K1];
KB=K01+K02 +Kg1+Kg2;

% ---------剛性マトリックス合成---------------------------
for I=1:3;for J=1:3
    KT(is+I,is+J)=KT(is+I,is+J)+KB(I,J);
    KT(is+I,js+J)=KT(is+I,js+J)+KB(I,3+J);
    KT(js+I,is+J)=KT(js+I,is+J)+KB(3+I,J);
    KT(js+I,js+J)=KT(js+I,js+J)+KB(3+I,3+J);
end;end
end

