Esope source files

   1 : C KATETA    SOURCE    CHAT      05/01/13    00:51:42     5004           
   2 :       SUBROUTINE KATETA(RI,ZI,RJ,ZJ,DRI,DZI,DRJ,DZJ,NM,NAL,AL,GG,KIMP
   3 :      &     ,EXTINC,RAD)
   4 :       IMPLICIT INTEGER(I-N)
   5 :       IMPLICIT REAL*8 (A-H,O-Z)
   6 : C    Include contenant quelques constantes dont XPI :
   7 : -INC CCREEL
   8 : C*********************************************************************
   9 : C SP appele par KAXK
  10 : C    INTEGRATION ANGULAIRE POUR 2 POINTS PRIS SUR LES FACES 1 ET 2
  11 : C     entree:
  12 : C      RI,ZI   :COORDONNEES DU POINT I SUR LE SEGMENT 1
  13 : C      RJ,ZJ   :COORDONNEES DU POINT J SUR LE SEGMENT 2
  14 : C      DRI,DZI :COORDONNEES DU VECTEUR SEGMENT 1
  15 : C      DRI,DZI :COORDONNEES DU VECTEUR SEGMENT 2
  16 : C      AL      :INTERVALLES D INTEGRATION
  17 : C      NAL     :NOMBRE D INTERVALLES
  18 : C      EXTINC  :coefficient d'extinction de la cavite si absorbante
  19 : C               (multiplié par RAD, dimension du pb ds FACAXI)
  20 : C      RAD     : dimension du pb(le calcul st fait en coor. reduites)
  21 : C     resultat:
  22 : C      GG
  23 : C      si transparent: la primitive existe
  24 : C      si absorbant  : calcul de l'integrale par la methode des
  25 : C                      trapèzes
  26 : C     
  27 : C*********************************************************************
  28 :       DIMENSION AL(2,NM)
  29 : 
  30 :       DMIN=1.D-5
  31 :       DMIN=DMIN*DMIN
  32 : 
  33 :       H=ZJ-ZI
  34 : 
  35 :       C1= RI*RJ*DZI*DZJ
  36 :       C1= RI*RJ*C1
  37 :       C4=  RI*RI + RJ*RJ + H*H
  38 :       C5= -2.D0*RI*RJ
  39 : 
  40 : *     milieu transparent
  41 : 
  42 :       IF(EXTINC.LT.1D-3) THEN
  43 : 
  44 :          C6=C4/C5
  45 :          C52=C5*C5
  46 : 
  47 :          ST=0.D0
  48 :          DO 100 I =1,NAL
  49 :             ST=ST+AL(2,I)-AL(1,I)
  50 :  100     CONTINUE
  51 :          G = (C1/C52)* ST
  52 : 
  53 : CC6>  
  54 :          IF(DABS(C6+1.D0).GE.DMIN) THEN
  55 : 
  56 : C>>   LES POINTS SONT DIFFERENTS
  57 : 
  58 :             C2= - DZI*DZJ*(RI*RI+RJ*RJ) + H*( RJ*DRJ*DZI - RI*DRI*DZJ)
  59 :             C3= (RI*DZI + DRI*H)*(RJ*DZJ-DRJ*H)
  60 :             C2= C2*RI*RJ
  61 :             C3= C3*RI*RJ
  62 :             C7=C2-2*C1*C6
  63 :             C8=C3-C1*C6*C6
  64 :             C9=C6*C7-C8
  65 :             C10=(C7-C6*C8)*2.D0
  66 :             C11=DSQRT(C6*C6-1.D0)
  67 :             C12=(C6-1.D0)/C11
  68 : 
  69 :             G2=0.D0
  70 :             G3=0.D0
  71 : 
  72 :             DO 110 IL=1,NAL
  73 :                T1=AL(1,IL)
  74 :                T2=AL(2,IL)
  75 :                GG2 =  (DSIN(T2)/(C6+DCOS(T2)))-(DSIN(T1)/(C6+DCOS(T1)))
  76 :                GG2=C9*GG2/C52/C11/C11
  77 : 
  78 :                IF (DABS(T2-XPI).LE.DMIN) THEN
  79 :                   Y2 = -XPI/2.D0
  80 :                ELSE
  81 :                   Y2 = DATAN(C12*DTAN(T2/2.D0))
  82 :                ENDIF
  83 :                IF (DABS(T1-XPI).LE.DMIN) THEN
  84 :                   Y1 = -XPI/2.D0
  85 :                ELSE
  86 :                   Y1 = DATAN(C12*DTAN(T1/2.D0))
  87 :                ENDIF
  88 : 
  89 :                GG3= (Y1-Y2)*C10 /C52/C11/C11/C11
  90 : 
  91 :                G2 = G2 + GG2
  92 :                G3 = G3 + GG3
  93 : 
  94 :  110        CONTINUE
  95 :             IF(KIMP.GE.4) WRITE(6,*) ' G1 G2 G3 '
  96 :             G  = G  + G2 + G3
  97 : CC6-  
  98 :          ELSE
  99 : C       WRITE(6,*) ' POINTS PROCHES ',IA,IB
 100 : 
 101 :          ENDIF
 102 : CC6<  
 103 :          GG=G
 104 : 
 105 : *     milieu absorbant
 106 : 
 107 :       ELSE
 108 : 
 109 :          C2= - DZI*DZJ*(RI*RI+RJ*RJ) + H*( RJ*DRJ*DZI - RI*DRI*DZJ)
 110 :          C3= (RI*DZI + DRI*H)*(RJ*DZJ-DRJ*H)
 111 :          C2= C2*RI*RJ
 112 :          C3= C3*RI*RJ
 113 : 
 114 : *      arbitraire
 115 :          NT=50
 116 : 
 117 :          G  = 0.D0
 118 :          DO 200 IL=1,NAL
 119 :             G1=0.D0
 120 :             DO 201 K=1,NT
 121 :                DTETA = (AL(2,IL)-AL(1,IL))/NT
 122 :                TETA = DTETA/2.D0 + (K-1)*DTETA
 123 :                XNUM=C1*DCOS(TETA)*DCOS(TETA)+C2*DCOS(TETA)+C3
 124 :                XDEN=C4+C5*DCOS(TETA)
 125 :                DIST=DSQRT(XDEN)
 126 : C         epaisseur optique
 127 :                EPAIS = EXTINC*RAD*DIST
 128 :                G1= G1+(XNUM/XDEN/XDEN)*(DEXP(-EPAIS))*DTETA
 129 :  201        CONTINUE
 130 :             G = G + G1
 131 :  200     CONTINUE
 132 : 
 133 :          GG=G
 134 : 
 135 :       ENDIF
 136 : 
 137 :       IF(KIMP.GE.4)WRITE(6,*) ' KATETA  G ',GG
 138 :       RETURN
 139 :       END
 140 : 
 141 : 
 142 : 
 143 :  
 144 :  
 145 :