SUBROUTINE RSURFU(H,P,TGT,DNDS,X,TIME,U,CINAME,SLNAME,
     1   MSNAME,NOEL,NODE,LCLOSE)
C     
      INCLUDE 'ABA_PARAM.INC'
C     
      CHARACTER*80 CINAME,SLNAME,MSNAME
      DIMENSION P(3),TGT(3,2),DNDS(3,2),X(3,3),TIME(2),U(6,2)
C               
      PARAMETER (NCRD = 2, PRECIS=1.0D-16, ASMALL = 1.0D-36)
      DIMENSION XINTRSCT(NCRD,2), DIST(2), CENTER(NCRD)
      DIMENSION CONESTARTCOORD(NCRD), D(2,2), F(2)
C     
C     DEFINE THE FOLLOWING QUANTITIES:
C     A = RADIUS 'A' OF THE SPHERICAL HEAD
C     SINA = SINE (CONE ANGLE ALPHA)
C     COSA = COSINE (CONE ANGLE ALPHA)
C     Z0 = ORIGINAL 'Z' COORDINATE OF POINT 'Q'
C     
      A=5.0D0
      SINA=0.5D0
      COSA=0.866025403784438647D0
      Z0=6.0D0
      ZQ=Z0 + U(2,2)
      CENTER(1) = U(1,2)
      CENTER(2) = ZQ
      CONESTARTCOORD(1) = A*COSA
      CONESTARTCOORD(2) = ZQ-A*SINA
C
C     THE FOLLOWING QUANTITIES ARE NOT NEEDED
C     FOR SURFACE-TO-SURFACE CONTACT APPROACH
C
      DNDS(1,1)=0.0D0
      DNDS(2,1)=0.0D0
C
C     FIND INTERSECTION WITH CIRCLE
C                    
      CALL SPHERELINEINTERSECTION(X,X(1,3),CENTER,A,NCRD,NINTRSCT,
     $     XINTRSCT,NCRD,DIST)
C                    
      IF (NINTRSCT .GT. 1) THEN
         DO KPT = 1, 2
            IF (XINTRSCT(1,KPT) .LE. CONESTARTCOORD(1)+PRECIS .AND.
     $           XINTRSCT(1,KPT) .GE. -PRECIS .AND.
     $           XINTRSCT(2,KPT) .LE. CONESTARTCOORD(2)+PRECIS .AND.
     $           XINTRSCT(2,KPT) .GE. ZQ-A-PRECIS) THEN
               H = -DIST(KPT)
               P(1) = ABS(XINTRSCT(1,KPT))
               P(2) = XINTRSCT(2,KPT)
               TGT(1,1)= -(CENTER(2)-P(2))/A
               TGT(2,1)= -P(1)/A
               RETURN
            END IF
         END DO
      END IF
C     
C     FIND INTERSECTION WITH CONE DEFINED WITH
C     A STARTING POINT CONESTARTCOORD AND A VECTOR
C     (SINA,COSA)
C     
      F(1) = CONESTARTCOORD(1) - X(1,1)
      F(2) = CONESTARTCOORD(2) - X(2,1)
      D(1,1) = X(1,3)
      D(1,2) = -SINA
      D(2,1) = X(2,3)
      D(2,2) = -COSA
      DETERMINANT = D(1,1)*D(2,2) - D(1,2)*D(2,1)
      IF (DETERMINANT .GT. ASMALL) THEN
         GAMMA = -D(2,1)*F(1) + D(1,1)*F(2)
         IF (GAMMA .GE. -PRECIS) THEN
            H = -(D(2,2)*F(1) - D(1,2)*F(2))
            P(1) = CONESTARTCOORD(1) + GAMMA*SINA
            P(2) = CONESTARTCOORD(2) + GAMMA*COSA
            TGT(1,1)= -SINA
            TGT(2,1)= -COSA
         END IF
      END IF
C
      RETURN
      END
      
      SUBROUTINE SPHERELINEINTERSECTION(Q,V,CENTER,R,NCRD,NINTRSCT,
     $     XINTRSCT,NROWDIMOFXINTRSCT,DIST)
C     
C     INTERSECTION OF A SPHERE DEFINED BY CENTER AND RADIUS R AND 
C     A LINE DEFINED BY POINT Q AND A UNIT VECTOR V.
C     RETURNS THE FOLLOWING:
C        NINTRSCT = NUMBER OF POINTS INTERSECTED
C        XINTRSCT = COORDINATES OF THE POINTS ON THE SPHERE
C                   INTERSECTED BY THE LINE
C
      INCLUDE 'ABA_PARAM.INC'
C                    
      PARAMETER (TWO = 2.0D0, ASMALL = 1.0D-36)
      DIMENSION Q(NCRD), V(NCRD), CENTER(NCRD)
      DIMENSION XINTRSCT(NROWDIMOFXINTRSCT,*), DIST(2)
C
      A = V(1)**2 + V(2)**2
      B = TWO*(V(1)*(Q(1)-CENTER(1)) + V(2)*(Q(2)-CENTER(2)))
      C = CENTER(1)**2 + CENTER(2)**2 + Q(1)**2 + Q(2)**2 -
     $     TWO*(CENTER(1)*Q(1) + CENTER(2)*Q(2)) - R**2
      IF (NCRD .EQ. 3) THEN
         A = A + V(3)**2
         B = B + TWO*V(3)*(Q(3)-CENTER(3))
         C = C + CENTER(3)**2 + Q(3)**2 - TWO*CENTER(3)*Q(3)
      END IF
C                    
      DISCRIMINANT = B**2 - 4.0D0*A*C
C                    
      NINTRSCT = 0
      IF (DISCRIMINANT .LT. 0.0D0) THEN
      ELSE
         QTEMP = -0.5D0*(B+SIGN(SQRT(DISCRIMINANT),B))
         IF (ABS(QTEMP) .GT. ASMALL) THEN
            NINTRSCT = NINTRSCT + 1
            DIST(NINTRSCT) = C/QTEMP
            DO KCRD = 1, NCRD
               XINTRSCT(KCRD,NINTRSCT) = Q(KCRD) +
     $              DIST(NINTRSCT)*V(KCRD)
            END DO
         END IF
         IF (ABS(A) .GT. ASMALL) THEN
            NINTRSCT = NINTRSCT + 1
            DIST(NINTRSCT) = QTEMP/A
            DO KCRD = 1, NCRD
               XINTRSCT(KCRD,NINTRSCT) = Q(KCRD) +
     $              DIST(NINTRSCT)*V(KCRD)
            END DO
         END IF            
      END IF
C
      RETURN 
      END