c===============================================================================
c==========================================================Programa VigaContinua
c===============================================================================
      implicit real*8 (a-h,o-z)
      parameter(MXN=100,MXKP=4*MXN+2)
      integer*4 n
      dimension work(MXKP)

    1 write(6,'(a,$)') ' Number of spans: '
      read(5,*) n
      if (n.gt.mxn) call ErrorFatal('Excessive number of spans (n>MXN)')

      kplibre=1
      call DimensionamientoDinamico(MXKP,kplibre,n  ,kpk)
      call DimensionamientoDinamico(MXKP,kplibre,n  ,kpl)
      call DimensionamientoDinamico(MXKP,kplibre,n+1,kpd)
      call DimensionamientoDinamico(MXKP,kplibre,n+1,kpw)

      call LeerDatos(n,work(kpk),work(kpw))
      call ResolverSistema(n,work(kpk),work(kpl),work(kpd),work(kpw))
      call EscribirResultados(n,work(kpw))

      end

c______________________________________________________________________LeerDatos
      Subroutine LeerDatos(n,rk,w)
      implicit real*8 (a-h,o-z)
      dimension rk(n), w(0:n)

c.....Rigideces
      do i=1,n
        write(6,100) i
  100   format(' Span ',i5,' -----> Stiffness EI/L = ',$)
        read(5,*) rk(i)
      enddo

c.....Momentos
      do i=0,n
        write(6,110) i
  110   format(' Node ',i5,' -----> Moment      = ',$)
        read(5,*) w(i)
      enddo

      return
      end

c________________________________________________________________ResolverSistema
      Subroutine ResolverSistema(n,rk,rl,rd,w)
      implicit real*8 (a-h,o-z)
      dimension rk(n),rl(n),rd(0:n),w(0:n)

c.....Factorizacion
      i=0
        rd(i)=4.d+00*(      rk(i+1))
      do i=1,n-1
        rl(i)=2.d+00*rk(i)/rd(i-1)
        rd(i)=4.d+00*(rk(i)+rk(i+1))-rl(i)*rd(i-1)*rl(i)
      enddo
      i=n
        rl(i)=2.d+00*rk(i)/rd(i-1)
        rd(i)=4.d+00*(rk(i)        )-rl(i)*rd(i-1)*rl(i)

c.....Sistemas
      do i=1,n
        w(i)=w(i)-rl(i  )*w(i-1)
      enddo
      do i=0,n
        w(i)=w(i)/rd(i)
      enddo
      do i=n-1,0,-1
        w(i)=w(i)-rl(i+1)*w(i+1)
      enddo

      return
      end

c_____________________________________________________________EscribirResultados
      Subroutine EscribirResultados(n,w)
      implicit real*8 (a-h,o-z)
      dimension w(0:n)

      write(6,200)
 200  format( ' NODAL ROTATIONS',/
     .        ' ===== == === =====',/
     .        ' Node      Rotation (radians)',/
     .        ' -----     ---------------',/)

      do i=0,n
        write(6,210) i,w(i)
 210    format(1x,i10,d15.6)
      enddo

      return
      end

c_______________________________________________________DimensionamientoDinamico
      Subroutine DimensionamientoDinamico(MXKP,kplibre,ncomponentes,kp)

      kp=kplibre
      kplibre=kplibre+ncomponentes
      if (kplibre.gt.MXKP+1) call ErrorFatal('Insufficient memory')

      return
      end

c_____________________________________________________________________ErrorFatal
      Subroutine ErrorFatal(mensaje)
      character*(*) mensaje

      write(6,100) mensaje
  100 format(' Error: ',a)

      stop
      end
