%
%  vibrationdata_local_stiffness_rotating.m  ver 1.1   by Tom Irvine
%
function[klocal]=...
              vibrationdata_local_stiffness_rotating(E,I,ne,dx,omega,L,rho)
%
jflag=0;
%
h=L/ne;   % assume constant length elements
%
klocal = zeros(ne, 4, 4);

for k=1:ne
%
    klocal(k,1,1)=12.;
    klocal(k,1,2)=6.;
    klocal(k,1,3)=-12.;
    klocal(k,1,4)=6.;
    klocal(k,2,2)=4.;
    klocal(k,2,3)=-6.;
    klocal(k,2,4)=2.;
    klocal(k,3,3)=12.;
    klocal(k,3,4)=-6.;
	klocal(k,4,4)=4.;
%
	klocal(k,1,1)=klocal(k,1,1)*1.;
	klocal(k,1,2)=klocal(k,1,2)*dx(k);
	klocal(k,1,3)=klocal(k,1,3)*1.;
	klocal(k,1,4)=klocal(k,1,4)*dx(k);
	klocal(k,2,2)=klocal(k,2,2)*dx(k)^2;
	klocal(k,2,3)=klocal(k,2,3)*dx(k);
	klocal(k,2,4)=klocal(k,2,4)*dx(k)^2;
	klocal(k,3,3)=klocal(k,3,3)*1.;
	klocal(k,3,4)=klocal(k,3,4)*dx(k);
	klocal(k,4,4)=klocal(k,4,4)*dx(k)^2;
%
% symmetry
%
	for i=1:4
		for j=i:4
            if( dx(k) < 1.0e-10 )
                disp(' dx error ');
                jflag=1;
                break;
            end
			klocal(k,i,j)=klocal(k,i,j)*E*I/dx(k)^3;            
        end
        if(jflag==1)
            return;
        end
    end
        
%
    j=k;
%

    r0=(j-1)*h/L;
    
    rbar  = r0 + h/2;  % mid-radius of the element
    
    Ne = 0.5*rho*omega^2*(L^2 - rbar^2);   % mid-span axial tension
    
    kr = (Ne/(30*h)) * [ 36,  3*h, -36,  3*h;
                          3*h, 4*h^2, -3*h, -h^2;
                         -36, -3*h,  36, -3*h;
                          3*h, -h^2, -3*h, 4*h^2 ];

%%%
%
	for i=1:4
		for j=1:4
            klocal(k,i,j)=klocal(k,i,j)+kr(i,j);
        end
    end
%
	for i=1:4
		for j=i:4
			klocal(k,j,i)=klocal(k,i,j);
        end
    end    
 %
end