Esope source files

   1 : C LINPLA    SOURCE    CHAT      05/01/13    01:18:52     5004           
   2 :       SUBROUTINE LINPLA(P1,P2,P3,S1,S2,PP,IDE)
   3 : C=============================================================
   4 : C                                                            
   5 : C    INTERSECTION D UNN SEGMENT S1 S2 AVEC UN TRIANGLE P1P2 P3
   6 : C                    =
   7 : C     P1,P2 ET P3 SONT TROIS POINTS DANS UN PLAN             =
   8 : C     S1 ET S2 SONT DEUX POINTS DANS UNE LIGNE DROITE        =
   9 : C                                                            =
  10 : C     OUTPUT:                                                =
  11 : C      PP: POINT D'INTERSECTION ENTRE LA LIGNE ET LE PLAN 
  12 : C      ide =-2  pas de solution   
  13 : C      IDE =-1  intersection  en dehors de P1 P2 P3
  14 : C      ide = 0   intersection de la droite  dans p1 p2 p2   =
  15 : C      IDE = 1, intersection  du segment avec       P1 P2 P2 =
  16 : C                                                            =
  17 : C     G. L. JIANG  / MODIFIE BAZE                            =
  18 : C                                                            =
  19 : C=============================================================
  20 :       IMPLICIT INTEGER(I-N)
  21 :       IMPLICIT REAL*8(A-H,O-Z)
  22 :       DIMENSION P1(*),P2(*),P3(*),S1(*),S2(*),PP(*)
  23 :       DATA EPSI/0.0001/
  24 :       DO 1 I=1,3
  25 :  1    PP(I)=0.0
  26 :       IDE=-10
  27 :       A=(P2(2)-P1(2))*(P3(3)-P1(3))-(P3(2)-P1(2))*(P2(3)-P1(3))
  28 :       B=(P3(1)-P1(1))*(P2(3)-P1(3))-(P2(1)-P1(1))*(P3(3)-P1(3))
  29 :       C=(P2(1)-P1(1))*(P3(2)-P1(2))-(P3(1)-P1(1))*(P2(2)-P1(2))
  30 :       D=-A*P1(1)-B*P1(2)-C*P1(3)
  31 :       P=S2(1)-S1(1)
  32 :       Q=S2(2)-S1(2)
  33 :       R=S2(3)-S1(3)
  34 :       CR=P*A+Q*B+R*C
  35 :       P9=ABS(P)
  36 :       Q9=ABS(Q)
  37 :       R9=ABS(R)
  38 : c
  39 :       IF(ABS(CR).LE.0.00001) THEN
  40 :         IDE=-2
  41 :         GOTO 1000
  42 :       END IF
  43 :       DO 20 I=1,3
  44 :   20  PP(I)=S1(I)
  45 :       IF(P9.GT.EPSI) THEN
  46 :         C0=(B*Q+C*R)/P
  47 :         C1=A+C0
  48 :         C2=D-C0*S1(1)+B*S1(2)+C*S1(3)
  49 :         PP(1)=-C2/C1
  50 :         PM=(PP(1)-S1(1))/P
  51 :         PP(2)=PM*Q+S1(2)
  52 :         PP(3)=PM*R+S1(3)
  53 :         GOTO 999
  54 :       END IF
  55 :   11  CONTINUE
  56 :       IF(Q9.GT.EPSI) THEN
  57 :         C0=(A*P+C*R)/Q
  58 :         C1=B+C0
  59 :         C2=D-C0*S1(2)+A*S1(1)+C*S1(3)
  60 :         PP(2)=-C2/C1
  61 :         PM=(PP(2)-S1(2))/Q
  62 :         PP(1)=PM*P+S1(1)
  63 :         PP(3)=PM*R+S1(3)
  64 :         GOTO 999
  65 :       END IF
  66 :       IF(R9.GT.EPSI) THEN
  67 :         C0=(A*P+B*Q)/R
  68 :         C1=C+C0
  69 :         C2=D-C0*S1(3)+A*S1(1)+B*S1(2)
  70 :         PP(3)=-C2/C1
  71 :         PM=(PP(3)-S1(3))/R
  72 :         PP(1)=PM*P+S1(1)
  73 :         PP(2)=PM*Q+S1(2)
  74 :         GOTO 999
  75 :       END IF      
  76 : 999   continue       
  76 :       CALL SURF3P(P1,P2,P3,AT)
  76 :       CALL SURF3P(P1,P2,PP,A1)
  76 :       CALL SURF3P(P2,P3,pp,A2)      
  76 :       CALL SURF3P(P3,P1,pp,A3)
  77 :       PS=0.D0
  78 :       do 30 i=1,3
  79 :       ps = ps + (s1(i)-pp(i))*(s2(i)-pp(i))
  80 : 30    CONTINUE 
  81 : C   le segment traverse le triangle 
  82 : C   la direction traverse le triangle 
  83 : C   la intersection en dehors du triangle        
  84 :       
  85 :       adif = abs(A1+A2+A3-AT)
  86 : *      write(6,*) 'linpl2' , adif,ps            
  87 :       if(adif.LT.0.0001) then 
  88 :          if(ps.lt.0.D0.or.abs(ps).LT.1.D-3) then 
  89 :          ide = 1
  90 :          else
  91 :          ide = 0 
  92 :          endif 
  93 :       else
  94 :       ide=-1
  95 :       endif      
  96 :  1000 CONTINUE
  97 :       END
  98 :  
  99 :  
 100 :