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,2), TIME(2), U(6,2)
C               
      PARAMETER (NCRD = 3, PRECIS=1.0D-16, ASMALL = 1.0D-36)
      DIMENSION XINTRSCT(NCRD,2), DIST(2), CENTER(NCRD)
      DIMENSION CONESTARTCOORD(NCRD), UNORMAL(NCRD)
      DIMENSION ICROS(3,2)
C                    
      DATA ICROS / 2, 3, 1,
     $             3, 1, 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=5.0D0
      ZQ= Z0 + U(3,2)
      CENTER(1) = U(1,2)
      CENTER(2) = U(2,2)
      CENTER(3) = ZQ
      CONESTARTCOORD(1) = A*COSA
      CONESTARTCOORD(2) = 0.0D0
      CONESTARTCOORD(3) = 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
      DNDS(3,1)=0.0D0
      DNDS(1,2)=0.0D0
      DNDS(2,2)=0.0D0
      DNDS(3,2)=0.0D0
C
C     FIND INTERSECTION WITH SPHERE
C                    
      CALL SPHERELINEINTERSECTION(X,X(1,3),CENTER,A,NCRD,NINTRSCT,
     $     XINTRSCT,NCRD,DIST)
C
      IF (NINTRSCT .GT. 1) THEN
         DO KPT = 1, 2
            R = SQRT(XINTRSCT(1,KPT)**2+XINTRSCT(2,KPT)**2)
            IF (R .LE. CONESTARTCOORD(1)+PRECIS .AND.
     $           XINTRSCT(3,KPT) .LE. CONESTARTCOORD(3)+PRECIS) THEN
               H = -DIST(KPT)
               P(1) = XINTRSCT(1,KPT)
               P(2) = XINTRSCT(2,KPT)
               P(3) = XINTRSCT(3,KPT)
C
C              CALCULATE THE TANGENTS
C
               UNORMAL(1) = P(1) - CENTER(1)
               UNORMAL(2) = P(2) - CENTER(2)
               UNORMAL(3) = P(3) - CENTER(3)
               UNORMALMAG = SQRT(UNORMAL(1)**2 + UNORMAL(2)**2 +
     $              UNORMAL(3)**2)
               UNORMAL(1) = UNORMAL(1)/UNORMALMAG
               UNORMAL(2) = UNORMAL(2)/UNORMALMAG
               UNORMAL(3) = UNORMAL(3)/UNORMALMAG
               TGT(1,1) = -P(1)
               TGT(2,1) = -P(2)
               TGT(3,1) =  0.0D0
               IF (SQRT(P(1)**2+P(2)**2) .LT. ASMALL) THEN
                  TGT(1,1) = -1.0D0
                  TGT(2,1) =  0.0D0
               END IF
               TGTDOTUNORMAL = TGT(1,1)*UNORMAL(1) +
     $              TGT(2,1)*UNORMAL(2) + TGT(3,1)*UNORMAL(3)
               TGT(1,1) = TGT(1,1) - TGTDOTUNORMAL*UNORMAL(1)
               TGT(2,1) = TGT(2,1) - TGTDOTUNORMAL*UNORMAL(2)
               TGT(3,1) = TGT(3,1) - TGTDOTUNORMAL*UNORMAL(3)
               TGT1MAG = SQRT(TGT(1,1)**2 + TGT(2,1)**2 + TGT(3,1)**2)
               TGT(1,1) = TGT(1,1)/TGT1MAG
               TGT(2,1) = TGT(2,1)/TGT1MAG
               TGT(3,1) = TGT(3,1)/TGT1MAG
               TGT(1,2) = UNORMAL(2)*TGT(3,1) - UNORMAL(3)*TGT(2,1)
               TGT(2,2) = UNORMAL(3)*TGT(1,1) - UNORMAL(1)*TGT(3,1)
               TGT(3,2) = UNORMAL(1)*TGT(2,1) - UNORMAL(2)*TGT(1,1)
               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,0,COSA). THIS CASE IS NOT CONSIDERED 
C     BECAUSE THE SLAVE SURFACE WILL NOT COME INTO
C     CONTACT WITH THE CONE.
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