Esope source files

   1 : C SELFT3    SOURCE    CHAT      05/01/13    03:13:30     5004           
   2 :       SUBROUTINE SELFT3(XGAUSS,NBNN,XE,QQ)
   3 : CCCCC
   4 : C     CALCUL DE LA SELF INDUCTANCE D'UN TRIANGLE
   5 : C     COPIE DE KALI1 DE CORFOU
   6 : CCCCC
   7 :       IMPLICIT INTEGER(I-N)
   8 :       IMPLICIT REAL*8(A-H,O-Z)
   9 : -INC CCOPTIO
  10 :       REAL*8 XGAUSS(3)
  11 :       REAL*8 XE(3,NBNN)
  12 :       REAL*8 U(5),V(5),W(5) 
  13 : *
  14 :       NDIM=IDIM
  15 :       NDIM2=IDIM-1
  16 : *
  17 :       U1=XGAUSS(1)
  18 :       V1=XGAUSS(2)
  19 :       W1=XGAUSS(3)
  20 :       DO 10 I=1,NBNN
  21 :         U(I)=XE(1,I)
  22 :         V(I)=XE(2,I)
  23 :         W(I)=XE(3,I)
  24 :    10 CONTINUE
  25 :       DO 11 I=1,2
  26 :         U(NBNN+I)=U(I)
  27 :         V(NBNN+I)=V(I)
  28 :         W(NBNN+I)=W(I)
  29 :    11 CONTINUE
  30 : *
  31 :       QQ=0.
  32 :       DO 20 J=1,3
  33 :         K1=J
  34 :         K2=J+1
  35 :         K3=J+2
  36 :         ALA=(U(K2)-U1)*(U(K3)-U(K2))+(V(K2)-V1)*(V(K3)-V(K2))
  37 :      &                              +(W(K2)-W1)*(W(K3)-W(K2))
  38 :         ALA=ALA/((U(K3)-U(K2))**2+(V(K3)-V(K2))**2+(W(K3)-W(K2))**2)
  39 :         XH=U(K2)-ALA*(U(K3)-U(K2))
  40 :         YH=V(K2)-ALA*(V(K3)-V(K2))
  41 :         ZH=W(K2)-ALA*(W(K3)-W(K2))
  42 :         T1=(V(K2)-V1)*(W(K3)-W1)-(V(K3)-V1)*(W(K2)-W1)
  43 :         T2=(W(K2)-W1)*(U(K3)-U1)-(W(K3)-W1)*(U(K2)-U1)
  44 :         T3=(U(K2)-U1)*(V(K3)-V1)-(U(K3)-U1)*(V(K2)-V1)
  45 :         R1=(YH-V1)*(W(K2)-W1)-(ZH-W1)*(V(K2)-V1)
  46 :         R2=(ZH-W1)*(U(K2)-U1)-(XH-U1)*(W(K2)-W1)
  47 :         R3=(XH-U1)*(V(K2)-V1)-(YH-V1)*(U(K2)-U1)
  48 :         PSCA=T1*R1+T2*R2+T3*R3
  49 :         THETA1=1.
  50 :         IF(PSCA.LT.0.) THETA1=-1.
  51 :         R1=(YH-V1)*(W(K3)-W1)-(ZH-W1)*(V(K3)-V1)
  52 :         R2=(ZH-W1)*(U(K3)-U1)-(XH-U1)*(W(K3)-W1)
  53 :         R3=(XH-U1)*(V(K3)-V1)-(YH-V1)*(U(K3)-U1)
  54 :         PSCA=T1*R1+T2*R2+T3*R3
  55 :         THETA2=1.
  56 :         IF(PSCA.LT.0.)THETA2=-1.
  57 :         PSCA=(XH-U1)*(U(K2)-U1)+(YH-V1)*(V(K2)-V1)+(ZH-W1)*(W(K2)-W1)
  58 :         D0=SQRT((XH-U1)**2+(YH-V1)**2+(ZH-W1)**2)
  59 :         D1=SQRT((U(K2)-U1)**2+(V(K2)-V1)**2+(W(K2)-W1)**2)
  60 :         COSTH1=PSCA/(D0*D1)
  61 :         THETA1=THETA1*ACOS(COSTH1)
  62 :         PSCA=(XH-U1)*(U(K3)-U1)+(YH-V1)*(V(K3)-V1)+(ZH-W1)*(W(K3)-W1)
  63 :         D2=SQRT((U(K3)-U1)**2+(V(K3)-V1)**2+(W(K3)-W1)**2)
  64 :         COSTH2=PSCA/(D0*D2)
  65 :         THETA2=THETA2*ACOS(COSTH2)
  66 :         ATA1=TAN(THETA1/2.)
  67 :         ATA2=TAN(THETA2/2.)
  68 :         ATA=((ATA2+1.)*(1.-ATA1))/((1.-ATA2)*(ATA1+1.))
  69 :         RATA=LOG(ATA)
  70 :         QQ=QQ+D0*RATA
  71 :   20  CONTINUE
  72 :       RETURN
  73 :       END
  74 :  
  75 :