subroutine genmat(K,E,I,L,A,phi)

  ! Esta rutina genera la matriz de rigidez para una sección de barra
  ! reticulada en ejes globales.  Sus argumentos son:
  
  ! E: Módulo Elástico
  ! I: Momento de Inercia de la sección de la barra
  ! L: Longitud de la barra
  ! A: Área de la barra
  ! phi: Ángulo de la barra respecto a la horizontal medido contrahorariamente
  
  real*8, dimension(6,6), intent(out) :: K
  real*8, intent(in) :: E,I,L,A,phi
  
  K(1,1) = (12d0*E*I/L**3)*sin(phi)**2+(E*A/L)*cos(phi)**2
  K(1,2) = (E*A/L)*cos(phi)*sin(phi)-(12d0*E*I/L**3)*cos(phi)*sin(phi)
  K(1,3) = -(6d0*E*I/L**2)*sin(phi)
  K(1,4) = -(12d0*E*I/L**3)*sin(phi)**2-(E*A/L)*cos(phi)**2
  K(1,5) = (12d0*E*I/L**3)*cos(phi)*sin(phi)-(E*A/L)*cos(phi)*sin(phi)
  K(1,6) = -(6d0*E*I/L**2)*sin(phi)
  K(2,1) = (E*A/L)*cos(phi)*sin(phi)-(12d0*E*I/L**3)*cos(phi)*sin(phi)
  K(2,2) = (E*A/L)*sin(phi)**2+(12d0*E*I/L**3)*cos(phi)**2
  K(2,3) = (6d0*E*I/L**2)*cos(phi)
  K(2,4) = (12d0*E*I/L**3)*cos(phi)*sin(phi)-(E*A/L)*cos(phi)*sin(phi)
  K(2,5) = -(E*A/L)*sin(phi)**2-(12d0*E*I/L**3)*cos(phi)**2
  K(2,6) = (6d0*E*I/L**2)*cos(phi)
  K(3,1) = -(6d0*E*I/L**2)*sin(phi)
  K(3,2) = (6d0*E*I/L**2)*cos(phi)
  K(3,3) = (4d0*E*I/L)
  K(3,4) = (6d0*E*I/L**2)*sin(phi)
  K(3,5) = -(6d0*E*I/L**2)*cos(phi)
  K(3,6) = 2d0*E*I/L
  K(4,1) = -(12d0*E*I/L**3)*sin(phi)**2-(E*A/L)*cos(phi)**2
  K(4,2) = (12d0*E*I/L**3)*cos(phi)*sin(phi)-(E*A/L)*cos(phi)*sin(phi)
  K(4,3) = (6d0*E*I/L**2)*sin(phi)
  K(4,4) = (12d0*E*I/L**3)*sin(phi)**2+(E*A/L)*cos(phi)**2
  K(4,5) = (E*A/L)*cos(phi)*sin(phi)-(12d0*E*I/L**3)*cos(phi)*sin(phi)
  K(4,6) = (6d0*E*I/L**2)*sin(phi)
  K(5,1) = (12d0*E*I/L**3)*cos(phi)*sin(phi)-(E*A/L)*cos(phi)*sin(phi)
  K(5,2) = -(E*A/L)*sin(phi)**2-(12d0*E*I/L**3)*cos(phi)**2
  K(5,3) = -(6d0*E*I/L**2)*cos(phi)
  K(5,4) = (E*A/L)*cos(phi)*sin(phi)-(12d0*E*I/L**3)*cos(phi)*sin(phi)
  K(5,5) = (E*A/L)*sin(phi)**2+(12d0*E*I/L**3)*cos(phi)**2
  K(5,6) = -(6d0*E*I/L**2)*cos(phi)
  K(6,1) = -(6d0*E*I/L**2)*sin(phi)
  K(6,2) = (6d0*E*I/L**2)*cos(phi)
  K(6,3) = 2d0*E*I/L
  K(6,4) = (6d0*E*I/L**2)*sin(phi)
  K(6,5) = -(6d0*E*I/L**2)*cos(phi)
  K(6,6) = (4d0*E*I/L)

  ! Y llegados aquí todos se preguntan... Guillem. ¿Cómo haces para escribir
  ! estos cipotes y no equivocarte? Pues no escribiendolos

end subroutine genmat

subroutine kglobal(K,n,Klocal,nnode1,nnode2)

  ! Esta subrutina toma la matriz de rigidez en ejes globales calculada
  ! con la rutina genmat y la añade a la matriz de rigidez de la
  ! estructura completa.
  
  ! Los argumentos son:
  
  ! K: Matriz de rigidez de la estructura completa
  ! n: Número de nodos de la estructura
  ! Klocal: Matriz de rigidez de la barra (el sufijo local no significa 
  ! que esté en ejes locales
  ! nnode1: número del nodo de un extremo de la barra
  ! nnode2: número del nodo del otro extremo de la barra
  
  integer, intent(in) :: n
  integer, intent(in) :: nnode1,nnode2
  real*8, dimension(6,6), intent(in) :: Klocal
  real*8, dimension(3*n,3*n), intent(inout):: K
  
  K(3*(nnode1-1)+1:3*(nnode1-1)+3,3*(nnode1-1)+1:3*(nnode1-1)+3)=&
       K(3*(nnode1-1)+1:3*(nnode1-1)+3,3*(nnode1-1)+1:3*(nnode1-1)+3)&
       +Klocal(1:3,1:3)
  K(3*(nnode2-1)+1:3*(nnode2-1)+3,3*(nnode2-1)+1:3*(nnode2-1)+3)=&
       K(3*(nnode2-1)+1:3*(nnode2-1)+3,3*(nnode2-1)+1:3*(nnode2-1)+3)&
       +Klocal(4:6,4:6)
  K(3*(nnode2-1)+1:3*(nnode2-1)+3,3*(nnode1-1)+1:3*(nnode1-1)+3)=&
       K(3*(nnode2-1)+1:3*(nnode2-1)+3,3*(nnode1-1)+1:3*(nnode1-1)+3)&
       +Klocal(4:6,1:3)
  K(3*(nnode1-1)+1:3*(nnode1-1)+3,3*(nnode2-1)+1:3*(nnode2-1)+3)=&
       K(3*(nnode1-1)+1:3*(nnode1-1)+3,3*(nnode2-1)+1:3*(nnode2-1)+3)&
       +Klocal(1:3,4:6)

end subroutine kglobal


!*****************

subroutine lubksb(a,indx,b,n)
  implicit none
  integer, intent(in) :: n
  real*8, dimension(n,n), intent(in) ::a
  integer, dimension(n), intent(in) :: indx
  real*8, dimension(n), intent(inout) :: b
  
  integer :: i,ii,ll
  real*8:: summ
  
  ii=0
  do i=1,n
     ll=indx(i)
     summ=b(ll)
     b(ll)=b(i)
     if (ii /= 0) then
        summ = summ-dot_product(a(i,ii:i-1),b(ii:i-1))
     else if (summ /= 0d0) then
        ii=i
     end if
     b(i)=summ
  end do
  do i=n,1,-1
     b(i)=(b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
  end do
  
end subroutine lubksb

!***************

subroutine ludcmp(a,indx,d,n)
  
  implicit none
  
  integer, intent(in) :: n
  real*8, dimension(n,n), intent(inout) :: a
  integer, dimension(n), intent(out) :: indx
  real*8, intent(out) :: d
  
  real*8, dimension(n) :: vv
  integer :: j,imax
  
  d=1d0
  vv=maxval(abs(a),dim=2)
  do j=1,n
     imax=(j-1)+maxloc(vv(j:n)*abs(a(j:n,j)),1)
     
     if (j /= imax) then
        call swap(a(imax,:),a(j,:),n)
        d=-d
        vv(imax)=vv(j)
     end if
     indx(j)=imax
     a(j+1:n,j)=a(j+1:n,j)/a(j,j)
     a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-spread(a(j+1:n,j),dim=2,ncopies=n-j)*&
          spread(a(j,j+1:n),dim=1,ncopies=n-j)
  end do
  
end subroutine ludcmp

!**********************

subroutine swap(a,b,n)
  implicit none
  integer, intent(in) :: n
  real*8, dimension(n), intent(inout) :: a,b
  
  real*8 :: sav
  integer :: i
  
  do i=1,n
     sav=a(i)
     a(i)=b(i)
     b(i)=sav
  end do
  
end subroutine swap

subroutine dlusolve(a,b,n)
  integer, intent(in) :: n
  real*8, dimension(n,n), intent(inout) :: a
  real*8, dimension(n), intent(inout) :: b
  
  integer, dimension(n) :: indx
  real*8 :: d
  
  call ludcmp(a,indx,d,n)
  call lubksb(a,indx,b,n)
  
end subroutine dlusolve

subroutine carga(p,b,n,nnodo,componente)

  ! Esta subrutina modifica el vector de cargas dadas los siguientes
  ! argumentos:
  
  ! p: valor de la carga
  ! b: vector de cargas
  ! n: número de nodos de la estructura
  ! nnodo: número del nodo al que se aplica la carga
  ! componente: componente a la que se aplica la carga
  !   1: componente x
  !   2: componente y
  !   3: Momento en z
  
  integer, intent(in) :: n,nnodo,componente
  real*8, intent(in) :: p
  real*8, dimension(3*n), intent(inout) :: b
  
  b(3*(nnodo-1)+componente)=p

end subroutine carga

subroutine ccapoyo(K,b,n,nnodo,componente)

  ! Esta subrutina modifica la matriz de rigidez del sistema y el
  ! vector de cargas para introducir condiciones de contorno de
  ! apoyo en un nodo según una dirección.  Los argumentos son:
  
  ! K: Matriz de rigidez
  ! b: vector de cargas
  ! n: número de nodos de la estructura
  ! componente: componente sobre la que se ejecuta la carga
  !   1: dirección según x
  !   2: dirección según y
  !   3: momento según z
  
  integer, intent(in) :: n,nnodo,componente
  real*8, dimension(3*n,3*n), intent(inout) :: K
  real*8, dimension(3*n), intent(inout) :: b
  
  K(3*(nnodo-1)+componente,:)=0
  K(:,3*(nnodo-1)+componente)=0
  K(3*(nnodo-1)+componente,3*(nnodo-1)+componente)=1
  b(3*(nnodo-1)+componente)=0
  
end subroutine ccapoyo  
  


