Sources Esope | Cast3M

neart3
C NEART3    SOURCE    CHAT      05/01/13    01:56:33     5004
      SUBROUTINE NEART3(XGAUSS,NBNN,XE1,XE2,QQ)
CCCCC
C     CALCUL DE LA MUTUELLE INDUCTANCE
C     ENTRE DEUX TRIANGLES PROCHES
C     COPIE DE RALI2 DE CORFOU
CCCCC
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 XGAUSS(3)
      REAL*8 XE1(3,NBNN),XE2(3,NBNN)
      REAL*8 U(5),V(5),W(5),UU(5),VV(5),WW(5),II(3)
*
      U0=XGAUSS(1)
      V0=XGAUSS(2)
      W0=XGAUSS(3)
      DO 10 I=1,NBNN
        U(I)=XE1(1,I)
        V(I)=XE1(2,I)
        W(I)=XE1(3,I)
        UU(I)=XE2(1,I)
        VV(I)=XE2(2,I)
        WW(I)=XE2(3,I)
   10 CONTINUE
      DO 11 I=1,2
        UU(NBNN+I)=UU(I)
        VV(NBNN+I)=VV(I)
        WW(NBNN+I)=WW(I)
   11 CONTINUE
*
      AN1=(UU(1)-U0)*(UU(2)-UU(1))+(VV(1)-V0)*(VV(2)-VV(1))
     .   +(WW(1)-W0)*(WW(2)-WW(1))
      AN2=(UU(1)-U0)*(UU(3)-UU(1))+(VV(1)-V0)*(VV(3)-VV(1))
     .   +(WW(1)-W0)*(WW(3)-WW(1))
      AA=(UU(2)-UU(1))*(UU(3)-UU(1))+(VV(2)-VV(1))*(VV(3)-VV(1))
     .  +(WW(2)-WW(1))*(WW(3)-WW(1))
      BB=(UU(2)-UU(1))**2+(VV(2)-VV(1))**2+(WW(2)-WW(1))**2
      CC=(UU(3)-UU(1))**2+(VV(3)-VV(1))**2+(WW(3)-WW(1))**2
      ALO=-CC*AN1+AA*AN2
      ALO=ALO/(BB*CC-AA*AA)
      AMO=-BB*AN2+AA*AN1
      AMO=AMO/(BB*CC-AA*AA)
      U1=UU(1)+ALO*(UU(2)-UU(1))+AMO*(UU(3)-UU(1))
      V1=VV(1)+ALO*(VV(2)-VV(1))+AMO*(VV(3)-VV(1))
      W1=WW(1)+ALO*(WW(2)-WW(1))+AMO*(WW(3)-WW(1))
      H0=SQRT((U0-U1)**2+(V0-V1)**2+(W0-W1)**2)
      DO 101 J=1,3
        K1=J
        K2=J+1
        K3=J+2
        XD=UU(K3)-UU(K2)
        YD=VV(K3)-VV(K2)
        ZD=WW(K3)-WW(K2)
        XV=(VV(K2)-VV(K1))*(WW(K3)-WW(K1))
     &    -(VV(K3)-VV(K1))*(WW(K2)-WW(K1))
        YV=(WW(K2)-WW(K1))*(UU(K3)-UU(K1))
     &    -(WW(K3)-WW(K1))*(UU(K2)-UU(K1))
        ZV=(UU(K2)-UU(K1))*(VV(K3)-VV(K1))
     &    -(UU(K3)-UU(K1))*(VV(K2)-VV(K1))
        R1=YD*ZV-YV*ZD
        R2=ZD*XV-ZV*XD
        R3=XD*YV-XV*YD
        P1=(UU(K1)-UU(K2))*R1+(VV(K1)-VV(K2))*R2+(WW(K1)-WW(K2))*R3
        P2=(U1-UU(K2))*R1+(V1-VV(K2))*R2+(W1-WW(K2))*R3
        II(J)=0
        IF(P1*P2.GT.0.) II(J)=1
        IF(P1*P2.LT.0.) II(J)=2
  101 CONTINUE
      QQ=0.
      DO 20 J=1,3
        K1=J
        K2=J+1
        K3=J+2
        IF(II(J).EQ.0) GOTO 20
        ALA=(UU(K2)-U1)*(UU(K3)-UU(K2))+(VV(K2)-V1)*(VV(K3)-VV(K2))
     .     +(WW(K2)-W1)*(WW(K3)-WW(K2))
        ALA=ALA/((UU(K3)-UU(K2))**2+(VV(K3)-VV(K2))**2
     &          +(WW(K3)-WW(K2))**2)
        XH=UU(K2)-ALA*(UU(K3)-UU(K2))
        YH=VV(K2)-ALA*(VV(K3)-VV(K2))
        ZH=WW(K2)-ALA*(WW(K3)-WW(K2))
        T1=(VV(K2)-V1)*(WW(K3)-W1)-(VV(K3)-V1)*(WW(K2)-W1)
        T2=(WW(K2)-W1)*(UU(K3)-U1)-(WW(K3)-W1)*(UU(K2)-U1)
        T3=(UU(K2)-U1)*(VV(K3)-V1)-(UU(K3)-U1)*(VV(K2)-V1)
        IF(II(J).EQ.1) GOTO 21
        T1=-T1
        T2=-T2
        T3=-T3
  21    CONTINUE
        R1=(YH-V1)*(WW(K2)-W1)-(ZH-W1)*(VV(K2)-V1)
        R2=(ZH-W1)*(UU(K2)-U1)-(XH-U1)*(WW(K2)-W1)
        R3=(XH-U1)*(VV(K2)-V1)-(YH-V1)*(UU(K2)-U1)
        PSCA=T1*R1+T2*R2+T3*R3
        THETA1=1.
        IF(PSCA.LT.0.) THETA1=-1.
        R1=(YH-V1)*(WW(K3)-W1)-(ZH-W1)*(VV(K3)-V1)
        R2=(ZH-W1)*(UU(K3)-U1)-(XH-U1)*(WW(K3)-W1)
        R3=(XH-U1)*(VV(K3)-V1)-(YH-V1)*(UU(K3)-U1)
        PSCA=T1*R1+T2*R2+T3*R3
        THETA2=1.
        IF(PSCA.LT.0.) THETA2=-1.
        PSCA=(XH-U1)*(UU(K2)-U1)+(YH-V1)*(VV(K2)-V1)
     &      +(ZH-W1)*(WW(K2)-W1)
        D0=SQRT((XH-U1)**2+(YH-V1)**2+(ZH-W1)**2)
        IF(D0.LT.1.E-5) GOTO 20
        D1=SQRT((UU(K2)-U1)**2+(VV(K2)-V1)**2+(WW(K2)-W1)**2)
        COSTH1=PSCA/(D0*D1)
        IF(ABS(COSTH1).GE.1.) THEN
          THETA1=0.
        ELSE
          THETA1=THETA1*ACOS(COSTH1)
        ENDIF
        PSCA=(XH-U1)*(UU(K3)-U1)+(YH-V1)*(VV(K3)-V1)
     &      +(ZH-W1)*(WW(K3)-W1)
        D2=SQRT((UU(K3)-U1)**2+(VV(K3)-V1)**2+(WW(K3)-W1)**2)
        COSTH2=PSCA/(D0*D2)
        IF(ABS(COSTH2).GE.1.) THEN
          THETA2=0.
        ELSE
          THETA2=THETA2*ACOS(COSTH2)
        ENDIF
        H1=D0
        IF(ABS(THETA1-THETA2).LE.1.E-4) GOTO 20
        FOR1=(SQRT((H0*COS(THETA2))**2+H1*H1)+H1*SIN(THETA2))
     .      *(SQRT((H0*COS(THETA1))**2+H1*H1)-H1*SIN(THETA1))
        FOR2=(SQRT((H0*COS(THETA2))**2+H1*H1)-H1*SIN(THETA2))
     .      *(SQRT((H0*COS(THETA1))**2+H1*H1)+H1*SIN(THETA1))
        FORA=H1*LOG(FOR1/FOR2)/2.
        FOR1=H0*SIN(THETA2)/SQRT(H0*H0+H1*H1)
        FOR2=H0*SIN(THETA1)/SQRT(H0*H0+H1*H1)
        IF(FOR1.LT.-1.)FOR1=-1.
        IF(FOR2.LT.-1.)FOR2=-1.
        IF(FOR1.GT.1.)FOR1=1.
        IF(FOR2.GT.1.)FOR2=1.
        FORB=H0*(ASIN(FOR1)-ASIN(FOR2))
        FORR=FORA+FORB-H0*(THETA2-THETA1)
        QQ=QQ+FORR
  20  CONTINUE
      RETURN
      END