C LINPLA    SOURCE    CHAT      05/01/13    01:18:52     5004           
      SUBROUTINE LINPLA(P1,P2,P3,S1,S2,PP,IDE)
C=============================================================
C                                                            
C    INTERSECTION D UNN SEGMENT S1 S2 AVEC UN TRIANGLE P1P2 P3
C                    =
C     P1,P2 ET P3 SONT TROIS POINTS DANS UN PLAN             =
C     S1 ET S2 SONT DEUX POINTS DANS UNE LIGNE DROITE        =
C                                                            =
C     OUTPUT:                                                =
C      PP: POINT D'INTERSECTION ENTRE LA LIGNE ET LE PLAN 
C      ide =-2  pas de solution   
C      IDE =-1  intersection  en dehors de P1 P2 P3
C      ide = 0   intersection de la droite  dans p1 p2 p2   =
C      IDE = 1, intersection  du segment avec       P1 P2 P2 =
C                                                            =
C     G. L. JIANG  / MODIFIE BAZE                            =
C                                                            =
C=============================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION P1(*),P2(*),P3(*),S1(*),S2(*),PP(*)
      DATA EPSI/0.0001/
      DO 1 I=1,3
 1    PP(I)=0.0
      IDE=-10
      A=(P2(2)-P1(2))*(P3(3)-P1(3))-(P3(2)-P1(2))*(P2(3)-P1(3))
      B=(P3(1)-P1(1))*(P2(3)-P1(3))-(P2(1)-P1(1))*(P3(3)-P1(3))
      C=(P2(1)-P1(1))*(P3(2)-P1(2))-(P3(1)-P1(1))*(P2(2)-P1(2))
      D=-A*P1(1)-B*P1(2)-C*P1(3)
      P=S2(1)-S1(1)
      Q=S2(2)-S1(2)
      R=S2(3)-S1(3)
      CR=P*A+Q*B+R*C
      P9=ABS(P)
      Q9=ABS(Q)
      R9=ABS(R)
c
      IF(ABS(CR).LE.0.00001) THEN
        IDE=-2
        GOTO 1000
      END IF
      DO 20 I=1,3
  20  PP(I)=S1(I)
      IF(P9.GT.EPSI) THEN
        C0=(B*Q+C*R)/P
        C1=A+C0
        C2=D-C0*S1(1)+B*S1(2)+C*S1(3)
        PP(1)=-C2/C1
        PM=(PP(1)-S1(1))/P
        PP(2)=PM*Q+S1(2)
        PP(3)=PM*R+S1(3)
        GOTO 999
      END IF
  11  CONTINUE
      IF(Q9.GT.EPSI) THEN
        C0=(A*P+C*R)/Q
        C1=B+C0
        C2=D-C0*S1(2)+A*S1(1)+C*S1(3)
        PP(2)=-C2/C1
        PM=(PP(2)-S1(2))/Q
        PP(1)=PM*P+S1(1)
        PP(3)=PM*R+S1(3)
        GOTO 999
      END IF
      IF(R9.GT.EPSI) THEN
        C0=(A*P+B*Q)/R
        C1=C+C0
        C2=D-C0*S1(3)+A*S1(1)+B*S1(2)
        PP(3)=-C2/C1
        PM=(PP(3)-S1(3))/R
        PP(1)=PM*P+S1(1)
        PP(2)=PM*Q+S1(2)
        GOTO 999
      END IF      
999   continue       
      CALL SURF3P(P1,P2,P3,AT)
      CALL SURF3P(P1,P2,PP,A1)
      CALL SURF3P(P2,P3,pp,A2)      
      CALL SURF3P(P3,P1,pp,A3)
      PS=0.D0
      do 30 i=1,3
      ps = ps + (s1(i)-pp(i))*(s2(i)-pp(i))
30    CONTINUE 
C   le segment traverse le triangle 
C   la direction traverse le triangle 
C   la intersection en dehors du triangle        
      
      adif = abs(A1+A2+A3-AT)
*      write(6,*) 'linpl2' , adif,ps            
      if(adif.LT.0.0001) then 
         if(ps.lt.0.D0.or.abs(ps).LT.1.D-3) then 
         ide = 1
         else
         ide = 0 
         endif 
      else
      ide=-1
      endif      
 1000 CONTINUE
      END