program Maw
	
C
C This program calculates the collision cycle according to the
C theory of N.Maw, J.R.Barber & J.N.Fawcett: "The Oblique Impact
C of Elastic Spheres", Wear 38(1), 101-114 (1976).
C Other useful experimental data and information is also found in
C Maw, Barber & Fawcett: "The Role of Elastic Tangential Compliance
C in Oblique Impact", ASME J. of Lubrication Technology (103), 74-80 (1981).
C
C M.Y. Louge, May 3, 1993.
C
C Modified to permit impacts between inhomogeneous spheres. Version 1.
C M.Louge, March 19, 1997
C

	integer Ns,N

	Parameter (Ns=100)

	real nu1,nu2,E1,E2,D,k,G1,G2,Pi
C	real RR,r1,r2,RM,m1,m2,diam1,diam2,dens1,dens2
	real RR,r1,r2,RM,m1,m2,diam1,diam2
	real gt1,gn1,gt2,gn2
	real psi1,psi2,ee,psi1min,psi1max
	real duration,dt,time
	real rt,rn,gt,gn,bb
	real Tforce,Nforce,dgt,dgn
	real mu,zero
	real aa,h,Nforcemax
	real srel

	real eta
	
	integer jb,iter,inc,number,stick,mark,Npsi
	
	real m,a,b,u,f
	dimension M(Ns,Ns),A(Ns,Ns),B(Ns,Ns),U(Ns),F(Ns)
	real sca,scb,sum,test
	dimension rt(Ns),srel(Ns)
	
	integer i,j
	dimension INDX(Ns),mark(Ns)
	real ND
		
	logical SS(Ns),Flag,choice

C Input variables
*
*	N	=	number of strips
*
*	nu1,nu2	=	Poisson's ratio
*	E1,E2	=	Young's modulus (N/m2)
*	G1,G2	=	modulus of rigidity (N/m2)
*
*	m1,m2	=	particle mass (kg)
*	diam1,2	=	particle diameter
*
*	mu	=	coefficient of Coulomb friction
*
*	RR	=	reduced radius (m)
*	RM	=	reduced mass (kg)
*	D	=	constant on p.30 of Landau & Lifshitz (m sec2/kg)
*	k	= 	constant in Landau & Lifshitz p. 31 (kg/Ăm.sec2)
*
*	gt1	=	tangential relative velocity before impact (³0) (m/sec)
*	gn1	=	normal relative velocity before impact (>0)
*	gt2	=	tangential relative velocity after impact
*	gn2	=	normal relative velocity after impact (²0)
*
C	Note: the convention for gn & gt is that they are positive before impact.
*
*	dt	=	time increment (sec)
*	inc	=	duration/dt; number of time steps to resolve the impact
C
*	Note that psi1 and psi2 are not defined here as Maw et al:
*
*	psi1	=	gt1/gn1
*	psi2 	= 	gt2/gn1
*
*	choice	=	if true (T) then only one run
*		=	if false (F) then Npsi multiple runs
*							with minimum psi1 = psi1min
*							and maximum psi1= psi1max.
*
C
*	eta	=	1/K^2 in Maw et al; = 5/2 for homogeneous spheres.
C

C Input from data file "Maw.dat"

C	Open(1, File='Maw.dat', Form='Formatted', Status='Old')
C
C Note that in version 1 of the program, the mass of each sphere is input, rather
C than the density.  MYL 3/19/97.
C
	read (*,*) diam1
	read (*,*) m1
	read (*,*) E1
	read (*,*) nu1

	read (*,*) diam2
	read (*,*) m2
	read (*,*) E2
	read (*,*) nu2
	
	read (*,*) N
	read (*,*) mu
	read (*,*) inc
	read (*,*) choice
	  if (choice) then
	   read (*,*) psi1
	  else
	   read (*,*) psi1min,psi1max,Npsi
	   psi1 = psi1min
	   Open(4, File='Maw.psi', Form='Formatted', Status='New') 
	  end if
C
C New in version 1: enter eta from input file:
C
	read (*,*) eta
	
C	Close (1)
	
C By definition in this program, we make gn1 = 1.0 m/sec without
C loss of generality. Then the only independent input variable is psi1.
	
	gn1 = 1.0
	
	Pi = 3.14159265
	
C	m1 = dens1*Pi/6.*diam1**3
C	m2 = dens2*Pi/6.*diam2**3
	
	r1 = diam1/2.
	r2 = diam2/2.

	D = 3./4.*( (1-nu1*nu1)/E1 + (1-nu2*nu2)/E2 )
	
	RR = 1./( (1./r1) + (1./r2) )
	RM = 1./( (1./m1) + (1./m2) )
	
	k = 4./5./D * sqrt(RR)
	
	G1 = E1/2./(1.+nu1)
	G2 = E2/2./(1.+nu2)
	
	zero = 0.0

C Calculation of the Hertzian contact. These quantities will be
C used as reference for the output.
*
*	duration	= total duration of impact (sec)
*	aa		= maximum impact radius (m)
*	h		= maximum approach (m)
*	Nforcemax	= maximum normal force (N)
*	

	duration = 2.943275184*
     &  RM**(2./5.)/k**(2./5.)/ (abs(gn1)**(1./5.))
     
     	h = (RM/k)**(2./5.) * (abs(gn1))**(4./5.)
     	aa = sqrt(h*RR)
	
	Nforcemax = aa**3/D/RR
	
	dt = duration/float(inc)
	
C Initial displacements & velocities
*	
*	rt(j)	=	tangential displacement of strip j (m)
*	rn	=	normal approach (m)
*	
*	gt	=	tangential relative velocity (m/sec)
*	gn	=	normal relative velocity
*
*	time	=	current time/total duration of the impact
*
C Initial matrices
*
*	M(j,i)	=	stiffness matrix (dim-)
*	A(j,i)	=	stiffness matrix for stick (dim-)
*	B(j,i)	=	stiffness matrix for slip (dim-)
*
*	SS(j)	=	stick/slip state
*
*				SS = True	(stick)
*				SS = False	(slip)
*
C	Note: A(j,i) is written Aij in Maw et al.
*
*	sca	=	stiffness coefficient for stick (m2 sec2/kg)
*	scb	=	stiffness coefficient for slip (N/m2)
*

	sca = (1./2.)*( (2.-nu1)/G1 + (2.-nu2)/G2 ) * aa/4.
	
	do i=1,Ns
	
		do j=1,Ns
		
			if (j.le.i) then
			
	A(j,i) = Pi/2./float(i)/float(N)*
     &	(2.*float(i*i)-(float(j)-0.5)*(float(j)-0.5))
     
			else
			
	A(j,i) = ( 2.*float(i)*float(i) -
     &	(float(j)-0.5)*(float(j)-0.5)     )*
     &	( asin(float(i)/(float(j)-0.5)) ) / float(i)/float(N) +
     &	sqrt( (float(j)-0.5)*(float(j)-0.5)-float(i)*float(i) )/
     &	float(N)
	
			end if
			
		end do
		
	end do
	
	scb = 3./2./Pi/D/RR*aa
	
	do i=1,Ns
	
		do j=1,Ns
		
			if (j.le.i) then
			
	B(j,i) = sqrt(1.- ( (float(j)-0.5)/(float(i)) )**2)
     
			else
			
	B(j,i) = 0.0
	
			end if
			
		end do
		
	end do
	
C In case multiple values of psi1 are needed, then loop over all these:

5	continue

	gt1 = psi1*gn1

	do j=1,Ns
		rt(j) = 0.0
	end do
	
	rn = 0.0
	
	time = 0.0
	Tforce = 0.0
	Nforce = 0.0
	
	gt = gt1
	gn = gn1
	
C First guess of the sign of the relative motion:
	
	do j=1,N
		srel(j) = sign(mu,gt1)
	end do
	
	do i=1,Ns
		do j=1,Ns
			M(j,i) = 0.0
		end do
	end do
	
		if (choice) then
	Open(2, File='Maw.hist', Form='Formatted', Status='New')
	Open(3, File='Maw.slip', Form='Formatted', Status='New')
		end if
	
	if (Choice) write (2,'(f5.4,1HT,
     &	e12.6,1HT,f8.6,1HT,
     &	e12.6,1HT,f9.6,1HT,
     &	f8.6,1HT,
     &	f6.4,1HT,
     &	f6.4)')
     &	time,
     &	Tforce/Nforcemax,Nforce/Nforcemax,
     &	gt/gn1,gn/gn1,
     &	rn/h,
     &	zero,
     &	bb/aa

     
C Begin impact cycle.

C Error control*******
	number = 0
C Error control*******

C New time & displacements:

1	time = time + dt/duration

C The following is an estimate of the tangential displacement after the
C time interval dt.  This estimate is correct for all strips
C that stick.  However, the actual tangential displacement is generally
C different for the strips that slip.  After the time interval
C dt has converged, the actual tangential displacement will be computed
C for all strips, based on the matrix A(j,i) and the elementary
C tractions F(i).

	do j=1,N
		rt(j) = rt(j) + gt*dt
	end do
	
	rn = rn + gn*dt
	
C Stop when the penetration goes to zero:
	
	if (rn.lt.0.0) go to 4
	
C Error control: stop if too many time increments ************
*
	number = number + 1
*	
	if (choice) write (*,'(a,i3)') 'step number ',number

		if (number.gt.2*inc) then
			write (*,'(a)')
     &			'Exceeded number of cycle increments'
     			stop
		end if
*		
C Error control: **********************************************
	
	bb = sqrt(rn*RR)
	
C jb is the index of the first outermost strip to lose contact.
C A strip is deemed to lose contact if its mean radius a(j-0.5)/n
C is greater or equal to b.

	jb = int(bb/aa*float(N)+0.5) + 1
	
	iter = 0
	
C Consider the non-contact regions as slipping in any case:
C => Initial distribution of contacts before impact: All slip
C because jb=1.
	
2	do j=jb,N
	
		SS(j)=.False.
		
	end do
	
C Construct the stiffness matrix M(j,i)

	do j=1,n
	
			if (SS(j)) then
			
		do i=1,n
		M(j,i) = A(j,i)
		end do
		
			else
		
		do i=1,n
		M(j,i) = B(j,i)
		end do
		
			end if
			
	end do
	
C Calculate the vector of displacements/tractions for each strip:

	do j=1,n
	
			if (SS(j)) then
			
*
C 	Positive tangential displacement of sticking strips results in
C 	positive tractions exerted by the test particle onto the stand.
*
	    U(j) = + rt(j)/sca/scb
		
			else
*			
C	For slipping strips, a relative displacement in the positive
C	direction (srel>0) imparts a positive traction from that strip
C	onto the stand.  Because strips that do not
C	contact are treated as slipping, a statement ensures
C	that the traction that they experience is zero.
*
				if (j.lt.jb) then
				
	    U(j) = + srel(j)*sqrt( (bb/aa)**2 -
     &               ((float(j)-0.5)/float(n))**2 )
     
				else
	    U(j) = 0.0
				end if
			
			end if
			
	end do
	
C Calculate the tractions; invert the matrix M:

	Call LUDCMP(M,N,Ns,INDX,ND)	
	Call LUBKSB(M,N,Ns,INDX,U)

C Save the traction components; the vector F is in fact F(i)/scb

	do i=1,N
	F(i) = U(i)
	end do
	
C Test each contacting strip for possible changes of status stick/slip:

* Flag is normally "False".  It turns to "True" if there is
* one or more changes in the status of a strip.

	Flag = .False.
	
*****
* No test if there is not a single contact strip:
*	
	if (jb.eq.1) go to 3
*
*****
	do j=1,jb-1
	
		if (SS(j)) then
		
C Test of sticking regions.  The test makes sure that the
C absolute tangential traction on the strip does not exceed
C mu times the normal traction from the Hertzian solution:

		sum = 0.0
		do i=j,N
		sum = sum + B(j,i)*F(i)
		end do
		
		test = mu*sqrt( (bb/aa)**2 -
     &          ((float(j)-0.5)/float(n))**2 ) - abs(sum)
     
     			if (test.lt.0.0) then
		
		sum = 0.0
		do i=1,N
		sum = sum + A(j,i)*F(i)
		end do
		
		
			SS(j) = .False.
			Flag  = .True.
			
			end if
			
		else
		
C Test of slipping regions.  The test calculates the relative
C tangential displacement of the slipping strips.  The relative
C displacement is the true displacement of the slipping strip,
C evaluated using the matrix A(j,i) and the elementary tractions
C from the numerical solution F(i), minus the displacement of the
C particle at the point of contact rt(j).
C For consistency, the relative displacement thus calculated
C must be in a direction opposite to the direction of the
C local velocity, -srel(j)*[sca*scb*AjiFi - rt(j)] > 0.  This condition
C is exactly equivalent to the test of Eq. (18) in Maw et al (1976).
C There, Maw et al. write srel(j)*[rt(j)-sca*scb*AjiFi].

		sum = 0.0
		do i=1,N
		sum = sum + A(j,i)*F(i)
		end do
		
	test = srel(j)*(rt(j)/sca/scb - sum)
		
		srel(j) = sign(mu,rt(j)/sca/scb - sum)
		
			if (test.lt.0.0) then
			
			SS(j) = .True.
			Flag  = .True.
			
			end if
		
		end if

	end do
		
C Error control: stop if too many iterations are required
C to find the stick/slip state:

		if (Flag) then
		iter=iter+1
			if (iter.ge.100) then
				if (choice) then
			write (*,'(a)')
     &			'Maximum number of iterations'
     			stop
				else
			write (4,'(/)')
	psi1 = psi1 + (psi1max-psi1min)/float(Npsi-1)
			go to 5
				end if
			end if
			
	if (choice) write (*,*) (SS(i),i=1,N)

		go to 2
		end if

3	continue

C Calculate the tangential force:

	Tforce = 0.0
	do i=1,N
	Tforce = Tforce + float(i)*float(i)*F(i)
	end do
	
	Tforce = Tforce*scb/float(n)/float(n)*2.*Pi*aa*aa/3.
	
C Calculate the normal force:

	Nforce = bb*bb*bb/D/RR
	
C Calculate the new displacements of each strip from the
C matrix Aji and the solutions for Fi:

	do j=1,N
	
		sum = 0.0
		
			do i=1,N
			sum = sum + A(j,i)*F(i)
			end do
		
		rt(j) = sum*sca*scb
	
	end do
	
C Calculate the new velocities:

C	dgt = - Tforce*dt*(7./2.)/RM
	dgt = - Tforce*dt* (1. + eta ) / RM
	gt  =   gt + dgt
		
	dgn = - Nforce*dt/RM
	gn  =   gn + dgn
	
C Calculate the proportion of sticking strips in the contact area:

	stick = 0
	do i = 1,N
		if (SS(i)) stick = stick + 1
	end do

C Output the results at the end of the current time step:
	
		if (jb.ne.1) then
	
	if (choice) write (2,'(f5.4,1HT,
     &	e12.6,1HT,f8.6,1HT,
     &	e12.6,1HT,f9.6,1HT,
     &	f8.6,1HT,
     &	f6.4,1HT,
     &	f6.4)')
     &	time,
     &	Tforce/Nforcemax,Nforce/Nforcemax,
     &	gt/gn1,gn/gn1,
     &	rn/h,
     &	float(stick)/float(jb-1)*bb/aa,
     &	bb/aa
     
     		else

	if (choice) write (2,'(f5.4,1HT,
     &	e12.6,1HT,f8.6,1HT,
     &	e12.6,1HT,f9.6,1HT,
     &	f8.6,1HT,
     &	f6.4,1HT,
     &	f6.4)')
     &	time,
     &	Tforce/Nforcemax,Nforce/Nforcemax,
     &	gt/gn1,gn/gn1,
     &	rn/h,
     &	zero,
     &	bb/aa
     
		end if
		
C Write the distribution of stick/slip contacts in "Maw.slip":
	
			if (choice) then
		if (jb.ne.1) then
	do j=1,jb-1
		  if (SS(j)) then
			mark(j)=1
		  else
			mark(j)=0
		  end if	
	end do
	
	write (3,'(100(i1))') (mark(j),j=1,jb-1)
	
		end if
			end if

C Proceed to next time step:

	go to 1
	
C When the impact is over, calculate the overall results:
	
4	continue

	gt2 = gt
	gn2 = gn
	
	psi1 = gt1/gn1
	psi2 = gt2/gn1
	
	ee   = -gn2/gn1
	
		if (choice) then
	
	write (2,'(///)')
	
      write (2,'(1x,A,1x,e12.6)') 'time     ',time
      write (2,'(1x,A,1x,e12.6)') 'psi1     ',psi1
      write (2,'(1x,A,1x,e12.6)') 'psi2     ',psi2
      write (2,'(1x,A,1x,e12.6)') 'e        ',ee
      
      		else
	
	write (4,'(e12.6,1HT,e12.6)') psi1,psi2
	psi1 = psi1 + (psi1max-psi1min)/float(Npsi-1)
	
		if (psi1.le.psi1max) go to 5
	
		end if
		
	if (choice) then
				close (2)
				close (3)
	else
				close (4)
	end if
	
	stop
	end
	
C Subroutines from Numerical Recipes for the solution of
C a set of linear equations:
	
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      PARAMETER (NMAX=100,TINY=1.0E-20)
      DIMENSION A(NP,NP),INDX(N),VV(NMAX)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11      CONTINUE
        IF (AAMAX.EQ.0.) then
		write (*,'(a)') 'Singular matrix.'
		stop
	end if
        VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,N
        IF (J.GT.1) THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF (I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
              A(I,J)=SUM
            ENDIF
14        CONTINUE
        ENDIF
        AAMAX=0.
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
            A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
16      CONTINUE
        IF (J.NE.IMAX)THEN
          DO 17 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
17        CONTINUE
          D=-D
          VV(IMAX)=VV(J)
        ENDIF
        INDX(J)=IMAX
        IF(J.NE.N)THEN
          IF(A(J,J).EQ.0.)A(J,J)=TINY
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
18        CONTINUE
        ENDIF
19    CONTINUE
      IF(A(N,N).EQ.0.)A(N,N)=TINY
      RETURN
      END

      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
      DIMENSION A(NP,NP),INDX(N),B(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0)THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
11        CONTINUE
        ELSE IF (SUM.NE.0.) THEN
          II=I
        ENDIF
        B(I)=SUM
12    CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF(I.LT.N)THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
13        CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END