      program string
      implicit real*8(a-h,o-z)
      parameter(np=161)
      dimension sklok(2,4),
     x     ska(np,np),ska0(np,np),skb(np,np),skb0(np,np),
     x     f(np),u(2,np),v(2,np),e(np),ipom(np)
      data alpha,gamma,  xl,  xf,  ro,  a,    q0,  xm0
     x     /1.d0, .01d0, 1d0, 1d0, 1d0, 1d0, -1d0, 1d0/


c-------------------------------------------------------------
c     string with the moving mass - space-time element method
c-------------------------------------------------------------
c	np - number of points
c	ne - number of elements (ne=np-1)
c	alpha, gamma - parameters of the space-time element method
c                    0.707 < alpha <= 1.0
c	               gamma=0 --- no damping
c	vm - velocity of the mass
c	h -time step
c	xl - length of the string
c	xf - tensile force
c	ro - mass density
c	a - cross section area
c	q0 - massless fpoint force
c	xm0 - moving mass
c	ia, ip - indices: ia=1 and ip=2 or ia=2 and ip=1
c

      vm=0.5
      h=xl/vm/(np-1)/40
      
      ip=1
      ne=np-1
      dx=xl/ne

      do 701 i3=1,np
         e(i3)=0d0
         f(i3)=0d0
         u(1,i3)=0d0
         v(1,i3)=0d0
         u(2,i3)=0d0
         v(2,i3)=0d0
         do 701 j3=1,np
            ska0(i3,j3)=0d0
            skb0(i3,j3)=0d0
 701  continue
 
c     .................................... element matrix
      call elem(sklok,xf,ro,a,dx,h,alpha,1)
      
c     .................................... global matrix assembly
      do 10 i=1,ne
         call dod(i,ska0,skb0,sklok,np)
 10   continue

c	tmax - final time of the experiment
c	nkrok - number of time steps

      tmax=xl/vm*1.0
      nkrok=int(tmax/h)
      
c.................................................... time step loop
      do 100 it=1,nkrok

      ip=3-ip
      ia=3-ip
      czas=it*h
      
c............................... nodal forces from the previous step
      do 11 i=1,np
 11      f(i)=-e(i)
 
c.................................. the load located in element 'iel'
      iel=int(czas*vm/dx-1d-10)+1
      if(iel.gt.ne)iel=0
      x0=(czas-h)*vm-iel*dx+dx
      xkappa=(x0+vm*alpha*h)/dx

      if(iel.gt.0)then
c........................................... moving load is added
         sk2=(czas*vm-iel*dx+dx)/dx
         sk1=1d0-sk2
         xm0s1=sk1*xm0
         xm0s2=sk2*xm0
         if(czas*vm.gt.0.)then
            f(iel  )=f(iel  )+sk1*q0
            f(iel+1)=f(iel+1)+sk2*q0
         endif
      endif
      
c................................ recuperation of matrix ska, skb
      do 72 ii=1,np
         do 72 jj=1,np
            ska(ii,jj)=ska0(ii,jj)
 72         skb(ii,jj)=skb0(ii,jj)
 
c.......................... mass matrix
      call maska(sklok,xm0,alpha,vm,h,x0,xf,a,dx,xkappa)      
      if(it.gt.0)call dod(iel,ska,skb,sklok,np)
      
c......................................... multiplication: ska*v1
c                            f:=F-A*v1  (f - right-hand side vector)
      do 31 i=1,np
        do 31 j=1,np
 31        f(i)=f(i)-ska(i,j)*v(ia,j)

      call warbrz(skb,f,np)

      call DGESV(np,1,skb,np,ipom,f,np,info)
      if(info.ne.0)stop

      do 66 ii=1,np
 66      v(ip,ii)=f(ii)

      beta=1.-alpha/(1.+gamma)

      do 51 ii=1,np
 51      u(ip,ii)=u(ia,ii)+h*((1-beta)*v(ia,ii)+beta*v(ip,ii))

      call elem(sklok,xf,ro,a,dx,h,alpha,2)
      do 83 ii=1,np
 83      e(ii)=0d0

      do 84 n=1,ne
         do 84 ii=1,2
            do 84 jj=1,2
 84            e(n-1+ii)=e(n-1+ii)+sklok(ii,jj)*u(ip,n-1+jj)

      xmnoz=xm0*vm**2/dx

      ul=u(ip,iel)
      up=u(ip,iel+1)
      vl=v(ip,iel)
      vp=v(ip,iel+1)
      
c...................... nodal forces
      e(iel  )=e(iel  )+xmnoz*( ul-up)/dx
      e(iel+1)=e(iel+1)+xmnoz*(-ul+up)/dx

      if(iel.eq.0)then
           ylok=0d0
      else
           ylok=sk1*u(ip,iel)+sk2*u(ip,iel+1)
      endif

      write(*,1001)czas*vm/xl,ylok,u(ip,(np-1)/2)
1001  format(3(1pe15.4))
100   continue
      end



c.................................................... mass marix
      subroutine maska(ste,xm0,alpha,v,h,x0,xf,a,b,xkappa)
      implicit real*8(a-h,o-z)
      dimension ste(2,4)

      do 9 i=1,2
         do 9 j=1,4
 9          ste(i,j)=0d0
 
c----------------stiffness
      xmnoz=xm0*v*v*h/b/b
      al=1-alpha
      ste(1,1)=ste(1,1)+xmnoz*al
      ste(1,2)=ste(1,2)-xmnoz*al
      ste(2,1)=ste(2,1)-xmnoz*al
      ste(2,2)=ste(2,2)+xmnoz*al

      ste(1,3)=ste(1,3)+xmnoz*(1.-al)
      ste(1,4)=ste(1,4)-xmnoz*(1.-al)
      ste(2,3)=ste(2,3)-xmnoz*(1.-al)
      ste(2,4)=ste(2,4)+xmnoz*(1.-al)
      
c-----------------damping
      xmnoz=xm0*v/b
      xkappa2=xkappa
      ste(1,1)=ste(1,1)-xmnoz*(1.-xkappa2)*(1-al)
      ste(1,2)=ste(1,2)+xmnoz*(1.-xkappa2)*(1-al)
      ste(2,1)=ste(2,1)-xmnoz*(xkappa2)*(1-al)
      ste(2,2)=ste(2,2)+xmnoz*(xkappa2)*(1-al)

      ste(1,3)=ste(1,3)-xmnoz*(1.-xkappa2)*al
      ste(1,4)=ste(1,4)+xmnoz*(1.-xkappa2)*al
      ste(2,3)=ste(2,3)-xmnoz*(xkappa2)*al
      ste(2,4)=ste(2,4)+xmnoz*(xkappa2)*al
      
c-----------------inertia
      xmnoz=xm0/h
      ste(1,1)=ste(1,1)-xmnoz*(1-xkappa)**2
      ste(1,2)=ste(1,2)-xmnoz*xkappa*(1-xkappa)
      ste(2,1)=ste(2,1)-xmnoz*xkappa*(1-xkappa)
      ste(2,2)=ste(2,2)-xmnoz*xkappa**2

      ste(1,3)=ste(1,3)+xmnoz*(1-xkappa)**2
      ste(1,4)=ste(1,4)+xmnoz*xkappa*(1-xkappa)
      ste(2,3)=ste(2,3)+xmnoz*xkappa*(1-xkappa)
      ste(2,4)=ste(2,4)+xmnoz*xkappa**2
      end
      
      
      
c....................................... string element
      subroutine elem(ste,f,ro,a,dx,h,alpha,itryb)
      implicit real*8(a-h,o-z)
      dimension ste(2,4)

      xr=f*a/dx

      ste(1,1)= xr
      ste(1,2)=-xr
      ste(2,1)=-xr
      ste(2,2)= xr

      if(itryb.eq.2)return

      xmnoz=h*alpha**2/2.d0

      do 16 i=1,2
         do 16 j=3,4
 16         ste(i,j)=ste(i,j-2)*xmnoz

      xmnoz=(1.d0-alpha/2.d0)*alpha*h

      do 17 i=1,2
         do 17 j=1,2
 17         ste(i,j)=ste(i,j)*xmnoz

c     inertia matrix...............................................
      xmas=ro*a*dx
      ste(1,1)=ste(1,1)-xmas/3d0/h
      ste(1,2)=ste(1,2)-xmas/6d0/h
      ste(2,1)=ste(2,1)-xmas/6d0/h
      ste(2,2)=ste(2,2)-xmas/3d0/h

      ste(1,3)=ste(1,3)+xmas/3d0/h
      ste(1,4)=ste(1,4)+xmas/6d0/h
      ste(2,3)=ste(2,3)+xmas/6d0/h
      ste(2,4)=ste(2,4)+xmas/3d0/h

      end
      
      
      
c........................................... boundary conditions
      subroutine warbrz(sk,r,np)
      implicit real*8(a-h,o-z)
      dimension sk(np,1),r(1),izam(2)
      data izam/1,0/
      izam(2)=np

      do 111 iii=1,2
      nrowb=izam(iii)

      do 430 j=1,np
      sk(nrowb,j)=0.d0
 430  sk(j,nrowb)=0.d0

      r(nrowb)=0.d0
 111  sk(nrowb,nrowb)=1.d0
      end




      subroutine dod(n,ska,skb,sklok,lukl)
      implicit real*8(a-h,o-z)
      dimension ska(lukl,1),skb(lukl,1),sklok(2,4)

      do 1 i=1,2
         do 1 j=1,2
            ska(n-1+i,n-1+j)=ska(n-1+i,n-1+j)+sklok(i,j)
 1          skb(n-1+i,n-1+j)=skb(n-1+i,n-1+j)+sklok(i,j+2)
      end
