Sources Esope | Cast3M

kateta
C KATETA    SOURCE    CB215821  17/11/30    21:16:38     9639           
      SUBROUTINE KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,GG,KIMP
     &     ,EXTINC,RAD)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
C*********************************************************************
C SP appele par KAXK
C    INTEGRATION ANGULAIRE POUR 2 POINTS PRIS SUR LES FACES 1 ET 2
C     entree:
C      RI,ZI   :COORDONNEES DU POINT I SUR LE SEGMENT 1
C      RJ,ZJ   :COORDONNEES DU POINT J SUR LE SEGMENT 2
C      DRI,DZI :COORDONNEES DU VECTEUR SEGMENT 1
C      DRI,DZI :COORDONNEES DU VECTEUR SEGMENT 2
C      AL      :INTERVALLES D INTEGRATION
C      NAL     :NOMBRE D INTERVALLES
C      EXTINC  :coefficient d'extinction de la cavite si absorbante
C               (multiplié par RAD, dimension du pb ds FACAXI)
C      RAD     : dimension du pb(le calcul st fait en coor. reduites)
C     resultat:
C      GG
C      si transparent: la primitive existe
C      si absorbant  : calcul de l'integrale par la methode des
C                      trapèzes
C
C*********************************************************************
      DIMENSION AL(2,NM)
 
      DMIN=1.D-5
      DMIN=DMIN*DMIN
 
      H=ZJ-ZI
 
      C1= RI*RJ*DZI*DZJ
      C1= RI*RJ*C1
      C4=  RI*RI + RJ*RJ + H*H
      C5= -2.D0*RI*RJ
 
*     milieu transparent
 
      IF(EXTINC.LT.1D-3) THEN
 
         C6=C4/C5
         C52=C5*C5
 
         ST=0.D0
         DO 100 I =1,NAL
            ST=ST+AL(2,I)-AL(1,I)
 100     CONTINUE
         G = (C1/C52)* ST
 
CC6>
         IF(ABS(C6+1.D0).GE.DMIN) THEN
 
C>>   LES POINTS SONT DIFFERENTS
 
            C2= - DZI*DZJ*(RI*RI+RJ*RJ) + H*( RJ*DRJ*DZI - RI*DRI*DZJ)
            C3= (RI*DZI + DRI*H)*(RJ*DZJ-DRJ*H)
            C2= C2*RI*RJ
            C3= C3*RI*RJ
            C7=C2-2*C1*C6
            C8=C3-C1*C6*C6
            C9=C6*C7-C8
            C10=(C7-C6*C8)*2.D0
            C11=SQRT(C6*C6-1.D0)
            C12=(C6-1.D0)/C11
 
            G2=0.D0
            G3=0.D0
 
            DO 110 IL=1,NAL
               T1=AL(1,IL)
               T2=AL(2,IL)
               GG2 =  (SIN(T2)/(C6+COS(T2)))-(SIN(T1)/(C6+COS(T1)))
               GG2=C9*GG2/C52/C11/C11
 
               IF (ABS(T2-XPI).LE.DMIN) THEN
                  Y2 = -XPI/2.D0
               ELSE
                  Y2 = ATAN(C12*TAN(T2/2.D0))
               ENDIF
               IF (ABS(T1-XPI).LE.DMIN) THEN
                  Y1 = -XPI/2.D0
               ELSE
                  Y1 = ATAN(C12*TAN(T1/2.D0))
               ENDIF
 
               GG3= (Y1-Y2)*C10 /C52/C11/C11/C11
 
               G2 = G2 + GG2
               G3 = G3 + GG3
 
 110        CONTINUE
            IF(KIMP.GE.4) WRITE(6,*) ' G1 G2 G3 '
            G  = G  + G2 + G3
CC6-
         ELSE
C       WRITE(6,*) ' POINTS PROCHES ',IA,IB
 
         ENDIF
CC6&lt;
         GG=G
 
*     milieu absorbant
 
      ELSE
 
         C2= - DZI*DZJ*(RI*RI+RJ*RJ) + H*( RJ*DRJ*DZI - RI*DRI*DZJ)
         C3= (RI*DZI + DRI*H)*(RJ*DZJ-DRJ*H)
         C2= C2*RI*RJ
         C3= C3*RI*RJ
 
*      arbitraire
         NT=50
 
         G  = 0.D0
         DO 200 IL=1,NAL
            G1=0.D0
            DO 201 K=1,NT
               DTETA = (AL(2,IL)-AL(1,IL))/NT
               TETA = DTETA/2.D0 + (K-1)*DTETA
               XNUM=C1*COS(TETA)*COS(TETA)+C2*COS(TETA)+C3
               XDEN=C4+C5*COS(TETA)
               DIST=SQRT(XDEN)
C         epaisseur optique
               EPAIS = EXTINC*RAD*DIST
               G1= G1+(XNUM/XDEN/XDEN)*(EXP(-EPAIS))*DTETA
 201        CONTINUE
            G = G + G1
 200     CONTINUE
 
         GG=G
 
      ENDIF
 
      IF(KIMP.GE.4)WRITE(6,*) ' KATETA  G ',GG
      RETURN
      END