C ANGSOL    SOURCE    CB215821  18/06/19    21:15:04     9862           
C   CALCUL DE L'ANGLE SOLIDE D'UN TRIEDRE
C
C    PIERRE VERPEAUX     MAI 1979
C
C
C     VERSION EN SIMPLE PRECISION   FORTRAN PUR
C
      SUBROUTINE ANGSOL(XO,XA,XB,XC,EPS,IFLAG,IFLIG)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)

      LOGICAL BNUM,BDEN
-INC CCREEL
      DIMENSION XO(3),XA(3),XB(3),XC(3),XU(3),XV(3),XW(3)
C  CALCUL DES TROIS VECTEURS COMPOSANT LE TRIEDRE
      IFLIG=0
      ZERO =REAL(XZERO)
      UNM  =REAL(-1.D0)
      UN   =REAL( 1.D0)
      DEUX =REAL( 2.D0)
      BNUM =.FALSE.
      BDEN =.FALSE.

      EPS  =ZERO

      DO 5 I=1,3
         A    =XO(I)
         XU(I)=XA(I)-A
         XV(I)=XB(I)-A
         XW(I)=XC(I)-A
 5    CONTINUE

C   CALCUL DES NORMES DES VECTEURS
      U=ZERO
      V=ZERO
      W=ZERO
      DO 6 I=1,3
         U=U+XU(I)**2
         V=V+XV(I)**2
         W=W+XW(I)**2
 6    CONTINUE

      IF(IFLAG.EQ.1) THEN
        IF (U.LT.XPETIT) GOTO 11
        IF (V.LT.XPETIT) GOTO 11
        IF (W.LT.XPETIT) GOTO 11
      ELSE
        IF (U.LT.XPETIT) RETURN
        IF (V.LT.XPETIT) RETURN
        IF (W.LT.XPETIT) RETURN
      ENDIF

      U=SQRT(U)
      V=SQRT(V)
      W=SQRT(W)

C   CALCUL DU PRODUIT MIXTE DES TROIS VECTEURS
      VOL=XU(1)*(XV(2)*XW(3)-XV(3)*XW(2))+
     #    XU(2)*(XV(3)*XW(1)-XV(1)*XW(3))+
     #    XU(3)*(XV(1)*XW(2)-XV(2)*XW(1))
      IF(IFLAG.EQ.1) THEN
        IF (ABS(VOL).LT.REAL(1.D-4)) CALL COPLAN(XO,XA,XB,XC,IFLIG)
        IF (IFLIG.EQ.1) GOTO 12
      ENDIF

C CALCUL DES ANGLES
      A1=(XV(1)*XW(1) + XV(2)*XW(2) + XV(3)*XW(3))/(V*W)
      B1=(XW(1)*XU(1) + XW(2)*XU(2) + XW(3)*XU(3))/(W*U)
      C1=(XU(1)*XV(1) + XU(2)*XV(2) + XU(3)*XV(3))/(U*V)

      A= ACOS(MAX(MIN(A1,UN),UNM))
      B= ACOS(MAX(MIN(B1,UN),UNM))
      C= ACOS(MAX(MIN(C1,UN),UNM))

C    CALCUL DES NUMERATEURS ET DENOMINATEURS DANS LE
C    SECOND MEMBRE DE LA FORMULE DE L'HUILIER
      P1=(A+B+C)/REAL(4.D0)
      P2=(B+C)  /DEUX
      P3=(A+C)  /DEUX
      P4=(A+B)  /DEUX

      DNUM= SIN(P1)* SIN(P2)* SIN(P3)* SIN(P4)
      IF (DNUM.LT.ZERO) THEN
        DNUM=ZERO
        BNUM=.TRUE.
      ELSE
        DNUM= SQRT(DNUM)
      ENDIF

      DDEN= COS(P1)* COS(P2)* COS(P3)* COS(P4)
      IF (DDEN.LT.ZERO) THEN
        DDEN=ZERO
        BDEN=.TRUE.
      ELSE
        DDEN= SQRT(DDEN)
      ENDIF

      IF (BNUM .AND. BDEN) THEN
        CALL ERREUR(21)
        RETURN
      ENDIF

C   CALCUL DE L'ANGLE SOLIDE
      EPS=SIGN(UN,VOL)*REAL(4.D0)*ATAN2(DNUM,DDEN)
      RETURN

C                        UN DES VECTEURS EST NUL
 11   IFLIG=1
 12   RETURN
      END