      program beam
      implicit real*8(a-h,o-z)
      parameter(nbd=4,nbd2p1=2*nbd+1,ne=100,nszf=(ne+1)*2)
      dimension sklok(4,4),smlok(4,4),sk(nszf,nbd2p1),
     x     sm(nszf,nbd2p1),skkon(nszf,nbd2p1),t(nszf),f(nszf),
     x     u1(nszf),v1(nszf),p1(nszf)
c  input data [centymeter, gram microsekond]
c  Young modulus, cross sec. moment, Poissone coeff., cross sec. area
c  height of the cross section, mass density, length, travelling mass, 
c  gravity acceleration, shape coefficient
      data xxe, xxi, xni,  xxa,  xxro,  xxl0, xxm0,  graw,  xk
     x    /1d0,.2d0,.3d0, 2.5d0, 7.8d0, 200d0, 5d2, 10d-10, 1.2d0/

c   time step, velocity, Newmark beta, Newmark gamma
      data h, v, beta, gamma /50.d0, 1d-4, 0.25d0, 0.5d0/
c-----------------------------------------------------
c     Timoshenko with the moving mass - Newmark method
c-----------------------------------------------------
	data iel,ielmem/1,1/
      open(2,file='results.dat')

c     sklok - stiffness element matrix
c     smlok - mass element matrix
c     sk    - global stiffness matrix
c     sm    - global mass matrix

      dx=xxl0/ne

      do 701 i3=1,nszf
         f(i3)=0d0
         u1(i3)=0d0
         v1(i3)=0d0
         p1(i3)=0d0
         do 701 j3=1,2*nbd+1
            sk(i3,j3)=0d0
            sm(i3,j3)=0d0
 701     continue

c     elemental matrices: stiffness and inertia
      call elemk(sklok,xxe,xxa,xxi,dx,xni,xk)
      call elemm(smlok,xxro,xxa,xxi,dx)

c     global matrices assembly - loop over elements
      do 10 i=1,ne
         call dodkm(i,sk,sklok,sm,smlok,nszf,nbd)
 10   continue
 
      tmax=xxl0/v*1.2d0
c....................................... starting time loop
      do 100 it=1,int(tmax/h)
      time=it*h
      
      do 7 i=1,nszf
 7       f(i)=0.d0
c     'iel' -element that carry the mass
      iel=int(time*v/dx)+1
      if(iel.gt.ne)iel=0
	xkappa=time*v/dx-iel+1d0

      do 11 i=1,nszf
 11      f(i)=0.d0
      if(iel.gt.0)then
         sk2=(time*v-iel*dx+dx)/dx
         sk1=1d0-sk2
         xxm0s1=sk1*xxm0
         xxm0s2=sk2*xxm0
c........................................ external moving load
         f(2*iel-1)=f(2*iel-1)-xxm0s1*graw
         f(2*iel+1)=f(2*iel+1)-xxm0s2*graw
      endif
      
      call matk(beta,gamma,h,sk,sm,skkon,f,t,nbd,nszf,u1,v1,p1,
     x   xkappa,xxm0,dx,iel,ielmem,v)
	call warbrz(skkon,f,nszf,nbd)
      call solve(skkon,f,nszf,nbd)
      call weknas(beta,gamma,u1,v1,p1,h,f,nszf)
      if(iel.eq.0)then
           ylok=0d0
      else
           ylok=sk1*u1(2*iel-1)+sk2*u1(2*iel+1)
      endif

      write(*,7001)time*v/xxl0,ylok,u1(ne+1)
      write(2,7001)time*v/xxl0,ylok,u1(ne+1)
 7001 format(f10.7,1pe15.4,1pe15.4)
	ielmem=iel
100   continue
      end
      
c..................................................................
      subroutine elemk(sklok,xe,a,xi,xl,xni,xk)
      implicit real*8(a-h,o-z)
      dimension sklok(4,4)

      xg=xe/2.d0/(1d0+xni)
      s1=xe*xi/xl
      s2=xk*xg*a/xl
      sklok(1,1)=s2
      sklok(1,2)=s2*xl/2.d0
      sklok(1,3)=-s2
      sklok(1,4)=s2*xl/2.d0
      sklok(2,1)=sklok(1,2)
      sklok(2,2)=s2*xl**2/4.d0+s1
      sklok(2,3)=-s2*xl/2.d0
      sklok(2,4)=s2*xl**2/4.-s1
      sklok(3,1)=sklok(1,3)
      sklok(3,2)=sklok(2,3)
      sklok(3,3)=s2
      sklok(3,4)=-s2*xl/2.d0
      sklok(4,1)=sklok(1,4)
      sklok(4,2)=sklok(2,4)
      sklok(4,3)=sklok(3,4)
      sklok(4,4)=s2*xl**2/4.d0+s1
      end
      
c..................................................................
      subroutine elemm(smlok,ro,a,xi,xl)
      implicit real*8(a-h,o-z)
      dimension smlok(4,4)

      s=ro*a*xl
      ss=ro*xi*xl
      smlok(1,1)=s/3.d0
      smlok(1,2)=0.d0
      smlok(1,3)=s/6.d0
      smlok(1,4)=0.d0
      smlok(2,1)=0.d0
      smlok(2,2)=ss/3.d0
      smlok(2,3)=0.d0
      smlok(2,4)=ss/6.d0
      smlok(3,1)=s/6.d0
      smlok(3,2)=0.d0
      smlok(3,3)=s/3.d0
      smlok(3,4)=0.d0
      smlok(4,1)=0.d0
      smlok(4,2)=ss/6.d0
      smlok(4,3)=0.d0
      smlok(4,4)=ss/3.d0
 
      end

c....................................... global matrix assembly
      subroutine dodkm(n,sk,sklok,sm,smlok,nszf,nbd)
      implicit real*8(a-h,o-z)
      dimension sk(nszf,-nbd:nbd),sm(nszf,-nbd:nbd),
     x	sklok(4,4),smlok(4,4)

      do 1 i=1,4
         do 1 j=1,4
           sk(2*n-2+i,j-i)=sk(2*n-2+i,j-i)+sklok(i,j)
1          sm(2*n-2+i,j-i)=sm(2*n-2+i,j-i)+smlok(i,j)
      end
      
c............ solution of the system of equations (Gauss procedure)
      subroutine solve(s,f,n,nbd)
      implicit real*8(a-h,o-z)
      dimension s(n,-nbd:nbd),f(1),x(10,10)
 
	do 91 i=1,10
	do 91 j=1,10
	x(i,j)=0d0
91	if(j-i.ge.-nbd.and.j-i.le.nbd)x(i,j)=s(i,-i+j)
 
      do 3 i=1,n-1
        do 3 j=i+1,min0(i+nbd,n)
         do 4 k=1,nbd
           if(dabs(s(i,0)).lt.1e-20)then
             write(*,*)'solve - singular matrix=',i
             stop
           endif
4          s(j,k+i-j)=s(j,k+i-j)-s(j,i-j)*s(i,k)/s(i,0)
3        f(j)=f(j)-s(j,i-j)/s(i,0)*f(i)

      do 6 i=n,1,-1
        do 5 k=1,min0(nbd,n-i)
5         f(i)=f(i)-s(i,k)*f(i+k)
6       f(i)=f(i)/s(i,0)
      end
      
c........................................... boundary conditions
      subroutine warbrz(sk,r,nszf,nbd)
      implicit real*8(a-h,o-z)
      dimension sk(nszf,-nbd:nbd),r(1),izam(2)
      data izam/1,0/
   
      izam(2)=nszf-1

      do 111 iii=1,2
        nrowb=izam(iii)
        r(nrowb)=0.d0
        do 430 j=-nbd,nbd
          sk(nrowb,j)=0.d0
          if(nrowb-j.ge.1.and.nrowb-j.le.nszf)sk(nrowb-j,j)=0.d0
430     continue
        sk(nrowb,0)=1.d0
111   continue
      end
      
c................................. Newmark procedure - first stage
      subroutine matk(beta,gamma,h,sk,sm,skkon,f,t,nbd,nszf,
     x    u1,v1,p1,xkappa,xm0,b,iel,ielmem,v)
      implicit real*8(a-h,o-z)
      dimension sk(nszf,-nbd:nbd),sm(nszf,-nbd:nbd),
     x     skkon(nszf,-nbd:nbd),t(1),f(1),u1(1),v1(1),p1(1)

      a0=1.d0/beta/h/h
      a1=gamma/beta/h
      a2=1.d0/beta/h
      a3=1.d0/2.d0/beta-1.d0
      a4=gamma/beta-1.d0
      a5=h/2.d0*(gamma/beta-2.d0)
      a6=h*(1.d0-gamma)
      a7=gamma*h
      iel1=iel*2-1
      iel2=iel*2+1

      do 2 i=1,nszf
         do 2 j=-nbd,nbd
 2          skkon(i,j)=sk(i,j)+a0*sm(i,j)
c................................. point mass addition
      wsp11m= xm0*(1.0d0-xkappa)*(1.0d0-xkappa)
      wsp13m= xm0*xkappa*(1.0d0-xkappa)
      wsp31m= xm0*xkappa*(1.0d0-xkappa)
      wsp33m= xm0*xkappa*xkappa
      
      wsp11c=-xm0*v/b*(1.0d0-xkappa)
      wsp13c= xm0*v/b*(1.0d0-xkappa)
      wsp31c=-xm0*v/b*xkappa
      wsp33c= xm0*v/b*xkappa

      wsp11k=-xm0*v/b/h*(1.0d0-xkappa) 
      wsp13k= xm0*v/b/h*(1.0d0-xkappa)
      wsp31k=-xm0*v/b/h*xkappa
      wsp33k= xm0*v/b/h*xkappa

      wsp11 = wsp11m*a0 + wsp11c*a1 + wsp11k
      wsp13 = wsp13m*a0 + wsp13c*a1 + wsp13k
      wsp31 = wsp31m*a0 + wsp31c*a1 + wsp31k
      wsp33 = wsp33m*a0 + wsp33c*a1 + wsp33k

      skkon(iel1, 0)=skkon(iel1, 0) + wsp11
      skkon(iel1, 2)=skkon(iel1, 2) + wsp13
      skkon(iel2,-2)=skkon(iel2,-2) + wsp31
      skkon(iel2, 0)=skkon(iel2, 0) + wsp33

      do 3 i=1,nszf
 3      t(i)=a0*u1(i) + a2*v1(i) + a3*p1(i)

      do 31 i=1,nszf
        do 31 j=-nbd,nbd
        if(i+j.ge.1.and.i+j.le.nszf) f(i)=f(i) + sm(i,j)*t(i+j)
31      continue

      f(iel1)=f(iel1) + wsp11m*t(iel1) + wsp13m*t(iel2)
      f(iel2)=f(iel2) + wsp31m*t(iel1) + wsp33m*t(iel2)

      do 4 i=1,nszf
 4      t(i)=a1*u1(i) + a4*v1(i) + a5*p1(i)

      t_iel1 = a1*u1(iel1) + a4*v1(iel1) + a5*p1(iel1)
      t_iel2 = a1*u1(iel2) + a4*v1(iel2) + a5*p1(iel2)
c.................................. f:=f+C*t
      f(iel1)=f(iel1) + wsp11c*t_iel1 + wsp13c*t_iel2
      f(iel2)=f(iel2) + wsp31c*t_iel1 + wsp33c*t_iel2
c.................................. nodal forces of the mass
      eps=u1(2*ielmem+1)-u1(2*ielmem-1)
      f(iel1)=f(iel1)+(xm0*v/b/h)*(1.0d0-xkappa)*eps
      f(iel2)=f(iel2)+(xm0*v/b/h)*       xkappa *eps
      end
      
c.................................. Newmark procedure - second stage
      subroutine weknas(beta,gamma,u1,v1,p1,h,r,nszf)
      implicit real*8(a-h,o-z)
      dimension u1(1),v1(1),p1(1),r(1)

      a0=1.d0/beta/h/h
      a1=gamma/beta/h
      a2=1.d0/beta/h
      a3=1.d0/2.d0/beta-1.d0
      a4=gamma/beta-1.
      a5=h/2.d0*(gamma/beta-2.d0)
      a6=h*(1.d0-gamma)
      a7=gamma*h

      do 1 i=1,nszf
         u2=r(i)
         p2=a0*(u2-u1(i))-a2*v1(i)-a3*p1(i)
         v2=v1(i)+a6*p1(i)+a7*p2
         u1(i)=u2
         v1(i)=v2
         p1(i)=p2
 1    continue
      end

