Esope source files

   1 : C NEART3    SOURCE    CHAT      05/01/13    01:56:33     5004           
   2 :       SUBROUTINE NEART3(XGAUSS,NBNN,XE1,XE2,QQ)
   3 : CCCCC
   4 : C     CALCUL DE LA MUTUELLE INDUCTANCE
   5 : C     ENTRE DEUX TRIANGLES PROCHES
   6 : C     COPIE DE RALI2 DE CORFOU
   7 : CCCCC
   8 :       IMPLICIT INTEGER(I-N)
   9 :       IMPLICIT REAL*8(A-H,O-Z)
  10 :       REAL*8 XGAUSS(3)
  11 :       REAL*8 XE1(3,NBNN),XE2(3,NBNN)
  12 :       REAL*8 U(5),V(5),W(5),UU(5),VV(5),WW(5),II(3) 
  13 : *
  14 :       U0=XGAUSS(1)
  15 :       V0=XGAUSS(2)
  16 :       W0=XGAUSS(3)
  17 :       DO 10 I=1,NBNN
  18 :         U(I)=XE1(1,I)
  19 :         V(I)=XE1(2,I)
  20 :         W(I)=XE1(3,I)
  21 :         UU(I)=XE2(1,I)
  22 :         VV(I)=XE2(2,I)
  23 :         WW(I)=XE2(3,I)
  24 :    10 CONTINUE
  25 :       DO 11 I=1,2
  26 :         UU(NBNN+I)=UU(I)
  27 :         VV(NBNN+I)=VV(I)
  28 :         WW(NBNN+I)=WW(I)
  29 :    11 CONTINUE
  30 : *
  31 :       AN1=(UU(1)-U0)*(UU(2)-UU(1))+(VV(1)-V0)*(VV(2)-VV(1))
  32 :      .   +(WW(1)-W0)*(WW(2)-WW(1))
  33 :       AN2=(UU(1)-U0)*(UU(3)-UU(1))+(VV(1)-V0)*(VV(3)-VV(1))
  34 :      .   +(WW(1)-W0)*(WW(3)-WW(1))
  35 :       AA=(UU(2)-UU(1))*(UU(3)-UU(1))+(VV(2)-VV(1))*(VV(3)-VV(1))
  36 :      .  +(WW(2)-WW(1))*(WW(3)-WW(1))
  37 :       BB=(UU(2)-UU(1))**2+(VV(2)-VV(1))**2+(WW(2)-WW(1))**2
  38 :       CC=(UU(3)-UU(1))**2+(VV(3)-VV(1))**2+(WW(3)-WW(1))**2
  39 :       ALO=-CC*AN1+AA*AN2
  40 :       ALO=ALO/(BB*CC-AA*AA)
  41 :       AMO=-BB*AN2+AA*AN1
  42 :       AMO=AMO/(BB*CC-AA*AA)
  43 :       U1=UU(1)+ALO*(UU(2)-UU(1))+AMO*(UU(3)-UU(1))
  44 :       V1=VV(1)+ALO*(VV(2)-VV(1))+AMO*(VV(3)-VV(1))
  45 :       W1=WW(1)+ALO*(WW(2)-WW(1))+AMO*(WW(3)-WW(1))
  46 :       H0=SQRT((U0-U1)**2+(V0-V1)**2+(W0-W1)**2)
  47 :       DO 101 J=1,3
  48 :         K1=J
  49 :         K2=J+1
  50 :         K3=J+2
  51 :         XD=UU(K3)-UU(K2)
  52 :         YD=VV(K3)-VV(K2)
  53 :         ZD=WW(K3)-WW(K2)
  54 :         XV=(VV(K2)-VV(K1))*(WW(K3)-WW(K1))
  55 :      &    -(VV(K3)-VV(K1))*(WW(K2)-WW(K1))
  56 :         YV=(WW(K2)-WW(K1))*(UU(K3)-UU(K1))
  57 :      &    -(WW(K3)-WW(K1))*(UU(K2)-UU(K1))
  58 :         ZV=(UU(K2)-UU(K1))*(VV(K3)-VV(K1))
  59 :      &    -(UU(K3)-UU(K1))*(VV(K2)-VV(K1))
  60 :         R1=YD*ZV-YV*ZD
  61 :         R2=ZD*XV-ZV*XD
  62 :         R3=XD*YV-XV*YD
  63 :         P1=(UU(K1)-UU(K2))*R1+(VV(K1)-VV(K2))*R2+(WW(K1)-WW(K2))*R3
  64 :         P2=(U1-UU(K2))*R1+(V1-VV(K2))*R2+(W1-WW(K2))*R3
  65 :         II(J)=0
  66 :         IF(P1*P2.GT.0.) II(J)=1
  67 :         IF(P1*P2.LT.0.) II(J)=2
  68 :   101 CONTINUE
  69 :       QQ=0.
  70 :       DO 20 J=1,3
  71 :         K1=J
  72 :         K2=J+1
  73 :         K3=J+2
  74 :         IF(II(J).EQ.0) GOTO 20
  75 :         ALA=(UU(K2)-U1)*(UU(K3)-UU(K2))+(VV(K2)-V1)*(VV(K3)-VV(K2))
  76 :      .     +(WW(K2)-W1)*(WW(K3)-WW(K2))
  77 :         ALA=ALA/((UU(K3)-UU(K2))**2+(VV(K3)-VV(K2))**2
  78 :      &          +(WW(K3)-WW(K2))**2)
  79 :         XH=UU(K2)-ALA*(UU(K3)-UU(K2))
  80 :         YH=VV(K2)-ALA*(VV(K3)-VV(K2))
  81 :         ZH=WW(K2)-ALA*(WW(K3)-WW(K2))
  82 :         T1=(VV(K2)-V1)*(WW(K3)-W1)-(VV(K3)-V1)*(WW(K2)-W1)
  83 :         T2=(WW(K2)-W1)*(UU(K3)-U1)-(WW(K3)-W1)*(UU(K2)-U1)
  84 :         T3=(UU(K2)-U1)*(VV(K3)-V1)-(UU(K3)-U1)*(VV(K2)-V1)
  85 :         IF(II(J).EQ.1) GOTO 21
  86 :         T1=-T1
  87 :         T2=-T2
  88 :         T3=-T3
  89 :   21    CONTINUE
  90 :         R1=(YH-V1)*(WW(K2)-W1)-(ZH-W1)*(VV(K2)-V1)
  91 :         R2=(ZH-W1)*(UU(K2)-U1)-(XH-U1)*(WW(K2)-W1)
  92 :         R3=(XH-U1)*(VV(K2)-V1)-(YH-V1)*(UU(K2)-U1)
  93 :         PSCA=T1*R1+T2*R2+T3*R3
  94 :         THETA1=1.
  95 :         IF(PSCA.LT.0.) THETA1=-1.
  96 :         R1=(YH-V1)*(WW(K3)-W1)-(ZH-W1)*(VV(K3)-V1)
  97 :         R2=(ZH-W1)*(UU(K3)-U1)-(XH-U1)*(WW(K3)-W1)
  98 :         R3=(XH-U1)*(VV(K3)-V1)-(YH-V1)*(UU(K3)-U1)
  99 :         PSCA=T1*R1+T2*R2+T3*R3
 100 :         THETA2=1.
 101 :         IF(PSCA.LT.0.) THETA2=-1.
 102 :         PSCA=(XH-U1)*(UU(K2)-U1)+(YH-V1)*(VV(K2)-V1)
 103 :      &      +(ZH-W1)*(WW(K2)-W1)
 104 :         D0=SQRT((XH-U1)**2+(YH-V1)**2+(ZH-W1)**2)
 105 :         IF(D0.LT.1.E-5) GOTO 20
 106 :         D1=SQRT((UU(K2)-U1)**2+(VV(K2)-V1)**2+(WW(K2)-W1)**2)
 107 :         COSTH1=PSCA/(D0*D1)
 108 :         IF(ABS(COSTH1).GE.1.) THEN
 109 :           THETA1=0.
 110 :         ELSE
 111 :           THETA1=THETA1*ACOS(COSTH1)
 112 :         ENDIF
 113 :         PSCA=(XH-U1)*(UU(K3)-U1)+(YH-V1)*(VV(K3)-V1)
 114 :      &      +(ZH-W1)*(WW(K3)-W1)
 115 :         D2=SQRT((UU(K3)-U1)**2+(VV(K3)-V1)**2+(WW(K3)-W1)**2)
 116 :         COSTH2=PSCA/(D0*D2)
 117 :         IF(ABS(COSTH2).GE.1.) THEN
 118 :           THETA2=0.
 119 :         ELSE
 120 :           THETA2=THETA2*ACOS(COSTH2)
 121 :         ENDIF
 122 :         H1=D0
 123 :         IF(ABS(THETA1-THETA2).LE.1.E-4) GOTO 20
 124 :         FOR1=(SQRT((H0*COS(THETA2))**2+H1*H1)+H1*SIN(THETA2))
 125 :      .      *(SQRT((H0*COS(THETA1))**2+H1*H1)-H1*SIN(THETA1))
 126 :         FOR2=(SQRT((H0*COS(THETA2))**2+H1*H1)-H1*SIN(THETA2))
 127 :      .      *(SQRT((H0*COS(THETA1))**2+H1*H1)+H1*SIN(THETA1))
 128 :         FORA=H1*LOG(FOR1/FOR2)/2.
 129 :         FOR1=H0*SIN(THETA2)/SQRT(H0*H0+H1*H1)
 130 :         FOR2=H0*SIN(THETA1)/SQRT(H0*H0+H1*H1)
 131 :         IF(FOR1.LT.-1.)FOR1=-1.
 132 :         IF(FOR2.LT.-1.)FOR2=-1.
 133 :         IF(FOR1.GT.1.)FOR1=1.
 134 :         IF(FOR2.GT.1.)FOR2=1.
 135 :         FORB=H0*(ASIN(FOR1)-ASIN(FOR2))
 136 :         FORR=FORA+FORB-H0*(THETA2-THETA1)
 137 :         QQ=QQ+FORR
 138 :   20  CONTINUE
 139 :       RETURN
 140 :       END
 141 :  
 142 :