function [btotg, detJt, Tra,detra2,BK,R,R_]=Bmatrix(coordor,Nsm,Ntm,Num)

b1=zeros(3,30,8);
b2=zeros(3,24,8);
b2iso=zeros(3,24,8);
detJt=zeros(8,1);
BK=zeros(6,30);

% passaggio al sistema di riferimento locale
[coord,Tra,tradef,detra2,traspo,R,R_]=local_sys(coordor);

xn=coord(:,1);
yn=coord(:,2);
zn=coord(:,3);

pos=sqrt(3)/3;

GaussP=[-pos,-pos,-pos;
    pos,-pos,-pos;
    pos,pos,-pos;
    -pos,pos,-pos;
    -pos,-pos,pos;
    pos,-pos,pos;
    pos,pos,pos;
    -pos,pos,pos];

for ii=1:8
    
    s=GaussP(ii,1);
    t=GaussP(ii,2);
    u=GaussP(ii,3);
    
    % shape functions derivatives
    Ns=Nsm(ii,:);
    Nt=Ntm(ii,:);
    Nu=Num(ii,:);
    
    xs=Ns*xn;
    
    xu=Nu*xn;
    
    xt=Nt*xn;
    
    ys=Ns*yn;
    
    yu=Nu*yn;
    
    yt=Nt*yn;
    
    zs=Ns*zn;
    
    zu=Nu*zn;
    
    zt=Nt*zn;
    
    
    J=[xs, ys, zs
        xt, yt, zt
        xu, yu, zu];
    
    detJ=xs*yt*zu+ys*zt*xu+zs*xt*yu-zs*yt*xu-ys*xt*zu-xs*zt*yu;
    % derivate funzioni di forma rispetto a x,y,z
    Nstu=[Ns;Nt;Nu];
    Nxyz=J\Nstu;
    
    Nx=Nxyz(1,:);
    Ny=Nxyz(2,:);
    Nz=Nxyz(3,:);
    
    b1isot=[Nx,zeros(1,16)
        zeros(1,8),Ny,zeros(1,8)
        zeros(1,16),Nz];
    
    % Enhanced Assumed strain field
    b0MacNeal=1/detJ*[s 0 0 s*t 0 s*u;
        0 t 0 s*t t*u 0;
        0 0 u 0 t*u s*u];
    
    bn=[b1isot,b0MacNeal];
    
    bs=[Ny,Nx,zeros(1,8)
        zeros(1,8),Nz,Ny
        Nz,zeros(1,8),Nx];
    
    b1(:,:,ii)=bn;
    b2iso(:,:,ii)=bs;
    detJt(ii,1)=detJ;
    
end

%calcolo della b2 con procedura di MacNeal (pag. 294)

pint=[1,2,3,4; 5,6,7,8;
    2,3,7,6; 1,4,8,5;
    1,2,6,5; 4,3,7,8];

for jj=1:3
    for kk=1:2
        % calcolo della componente di b2 come media pesata con il detJ dei
        % 4 punti attorno al punto di integrazione 'virtuale' per il taglio
        totdetJ=0;
        bs=zeros(3,24);
        for ii=1:4
            detJ=detJt(pint(jj*2-1+kk-1,ii),1);
            bs(jj,:)=bs(jj,:)+b2iso(jj,:,pint(jj*2-1+kk-1,ii))*detJ;
            totdetJ=totdetJ+detJ;
        end
        for ii=1:4
            b2(jj,:,pint(jj*2-1+kk-1,ii))=bs(jj,:)/totdetJ;
        end
    end
end

indexdof=zeros(1,24);
for ii=0:7
    indexdof(1,ii*3+1:ii*3+3)=[1+ii,9+ii,17+ii];
end

b2(:,:,:)=b2(:,indexdof,:);
b1(:,1:24,:)=b1(:,indexdof,:);

Tr=eye(30);
for ii=1:8
    Tr(ii*3-2:ii*3,ii*3-2:ii*3)=traspo;
end


% applicazione delle trasformazioni di spostamento e coordinate nello spazio
% cartesiano distorto
btotg=zeros(6,30,8);
for ii=1:8
    btotg(:,:,ii)=tradef*[b1(:,:,ii);b2(:,:,ii),zeros(3,6)]*Tr;
    BK=BK+btotg(:,:,ii)*detJt(ii);
end

end % fucntion end