Esope source files

   1 : C ANGSOL    SOURCE    CHAT      05/01/12    21:22:26     5004           
   2 : C   CALCUL DE L'ANGLE SOLIDE D'UN TRIEDRE
   3 : C
   4 : C    PIERRE VERPEAUX     MAI 1979
   5 : C
   6 : C
   7 : C     VERSION EN SIMPLE PRECISION   FORTRAN PUR
   8 : C
   9 :       SUBROUTINE ANGSOL(XO,XA,XB,XC,EPS,IFLAG,IFLIG)
  10 :       IMPLICIT INTEGER(I-N)
  11 :       IMPLICIT REAL*8 (A-H,O-Z)
  12 : -INC CCREEL
  13 :       DIMENSION XO(3),XA(3),XB(3),XC(3),XU(3),XV(3),XW(3)
  14 : C  CALCUL DES TROIS VECTEURS COMPOSANT LE TRIEDRE
  15 :       IFLIG=0
  16 :       DO 5 I=1,3
  17 :          A=XO(I)
  18 :          XU(I)=XA(I)-A
  19 :          XV(I)=XB(I)-A
  20 :          XW(I)=XC(I)-A
  21 :  5    CONTINUE
  22 : C   CALCUL DES NORMES DES VECTEURS
  23 :       U=0.D0
  24 :       V=0.D0
  25 :       W=0.D0
  26 :       DO 6 I=1,3
  27 :          U=U+XU(I)**2
  28 :          V=V+XV(I)**2
  29 :          W=W+XW(I)**2
  30 :  6    CONTINUE
  31 :       U= SQRT(U)
  32 :       V= SQRT(V)
  33 :       W= SQRT(W)
  34 :       IF(IFLAG.EQ.1) THEN
  35 :          IF (U.LT.XPETIT) GOTO 11
  36 :          IF (V.LT.XPETIT) GOTO 11
  37 :          IF (W.LT.XPETIT) GOTO 11
  38 :       ENDIF
  39 :       IF (U.LT.XPETIT) GOTO 10
  40 :       IF (V.LT.XPETIT) GOTO 10
  41 :       IF (W.LT.XPETIT) GOTO 10
  42 : C   CALCUL DU PRODUIT MIXTE DES TROIS VECTEURS
  43 :       VOL=XU(1)*(XV(2)*XW(3)-XV(3)*XW(2))-XU(2)*(XV(1)*XW(3)-XV(3)*XW(1)
  44 :      #     )+XU(3)*(XV(1)*XW(2)-XV(2)*XW(1))
  45 :       IF(IFLAG.EQ.1) THEN
  45 :          IF (ABS(VOL).LT.0.0001D0) CALL COPLAN(XO,XA,XB,XC,IFLIG)
  46 :          IF (IFLIG.EQ.1) GOTO 12
  47 :       ENDIF
  48 : C CALCUL DES ANGLES
  49 :       A= ACOS((XV(1)*XW(1)+XV(2)*XW(2)+XV(3)*XW(3))/(V*W))
  50 :       B= ACOS((XW(1)*XU(1)+XW(2)*XU(2)+XW(3)*XU(3))/(W*U))
  51 :       C= ACOS((XU(1)*XV(1)+XU(2)*XV(2)+XU(3)*XV(3))/(U*V))
  52 :       P=(A+B+C)/2.D0
  53 : C    CALCUL DES NUMERATEURS ET DENOMINATEURS DANS LE
  54 : C    SECOND MEMBRE DE LA FORMULE DE L'HUILIER
  55 :       P1=P/2.D0
  56 :       P2=(P-A)/2.D0
  57 :       P3=(P-B)/2.D0
  58 :       P4=(P-C)/2.D0
  59 :       DNUM= SIN(P1)* SIN(P2)* SIN(P3)* SIN(P4)
  60 :       IF (DNUM.LT.0.D0) DNUM=0.D0
  61 :       DNUM= SQRT(DNUM)
  62 :       DDEN= COS(P1)* COS(P2)* COS(P3)* COS(P4)
  63 :       IF (DDEN.LT.0.D0) DDEN=0.D0
  64 :       DDEN= SQRT(DDEN)
  65 : C   CALCUL DE L'ANGLE SOLIDE
  66 :       EPS= ATAN2(DNUM,DDEN)
  67 :       EPS=4.D0*EPS
  68 :       IF (VOL.LE.0.D0) EPS=-EPS
  69 :       RETURN
  70 :  10   EPS=0.D0
  71 :       RETURN
  72 : C                        UN DES VECTEURS EST NUL
  73 :  11   IFLIG=1
  74 :  12   RETURN
  75 :       END
  76 : 
  77 : 
  78 : 
  79 : 
  80 : 
  81 :  
  82 :  
  83 :