Sources Esope | Cast3M

cq2kp
C CQ2KP     SOURCE    CHAT      05/01/12    22:26:39     5004
      SUBROUTINE CQ2KP(COO,PRESSI,AN,ALF11,ALF12,
     1     ALF21,ALF22,RAYON,RE,CGAUS,TGAUS,NGAUS,A1,A4,S3,AM,A)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
C
CCCCCCCCCCCCCC  CALCUL DE LA MATRICE KP POUR LES COQ2 TIREE INCA
C
      DIMENSION A4(8,8),S3(8,8),AM(8,8),A1(8,8)
      DIMENSION COO(3,2),A(8,8)
      DIMENSION ALF11(*),ALF12(*),ALF21(*),ALF22(*),RAYON(*),RE(*)
      DIMENSION CGAUS(*),TGAUS(*)
      CALL ZDANUL(RE,64)
      RAYON(1)=COO(1,1)
      RAYON(3)=COO(1,2)
      RAYON(2)=(RAYON(1)+RAYON(3))*0.5D0
      DR=RAYON(3)-RAYON(1)
      DZ=COO(2,2)-COO(2,1)
      DS=SQRT(DR*DR+DZ*DZ)
      COSP=DR/DS
      SINP=DZ/DS
      XS2=0.5D0*DS
      CALL RESOLV(A4,DS,COSP,SINP,S3,A)
      CALL ZDANUL(ALF12,5)
      CALL ZDANUL(ALF22,7)
      DO 9 IN=1,NGAUS
         RAY=RAYON(2)+XS2*COSP*TGAUS(IN)
         X=TGAUS(IN)*XS2
         X2=X*X
         X3=X2*X
         X4=X3*X
         ANS=0.D0
         ANT=0.D0
         ALF11(1)=1.D0
         ALF11(2)=X
         ALF11(3)=X2
         ALF11(4)=X3
         ALF11(5)=X4
         DO 12 NC=1,5
            ALF21(NC)=ALF11(NC)
 12      CONTINUE
         ALF21(6)=X4*X
         ALF21(7)=X4*X2
         ANS=ANS*RAY*XS2*CGAUS(IN)
         CALL AMULX(ALF11,5,ANS,ALF12)
         ANT=ANT*XS2*CGAUS(IN)/RAY
         CALL AMULX(ALF21,7,ANT,ALF22)
 9    CONTINUE
      CALL SHIFTD(ALF12,ALF11,5)
      CALL SHIFTD(ALF22,ALF21,7)
      DO 15 NC=1,8
         DO 13 NCC=1,8
            A1(NCC,NC)=A4(NC,NCC)
 13      CONTINUE
 15   CONTINUE
      CALL ZDANUL(A,64)
      PROD1=COSP*COSP
      PROD2=SINP*SINP
      PROD3=SINP*COSP
      A(1,1)=PROD1*ALF21(1)
      A(1,2)=PROD1*ALF21(2)
      A(1,5)=PROD3*ALF21(1)
      A(1,6)=PROD3*ALF21(2)
      A(1,7)=PROD3*ALF21(3)
      A(1,8)=PROD3*ALF21(4)
      A(2,1)=A(1,2)
      A(2,2)=ALF11(1)+PROD1*ALF21(3)
      A(2,5)=PROD3*ALF21(2)
      A(2,6)=PROD3*ALF21(3)
      A(2,7)=PROD3*ALF21(4)
      A(2,8)=PROD3*ALF21(5)
      A(3,3)=ALF21(1)
      A(3,4)=ALF21(2)
      A(4,3)=A(3,4)
      A(4,4)=ALF11(1)+ALF21(3)
      A(5,1)=A(1,5)
      A(5,2)=A(2,5)
      A(5,5)=PROD2*ALF21(1)
      A(5,6)=PROD2*ALF21(2)
      A(5,7)=PROD2*ALF21(3)
      A(5,8)=PROD2*ALF21(4)
      A(6,1)=A(1,6)
      A(6,2)=A(2,6)
      A(6,5)=A(5,6)
      A(6,6)=ALF11(1)+A(5,7)
      A(6,7)=ALF11(2)*2.D0+A(5,8)
      A(6,8)=ALF11(3)*3.D0+PROD2*ALF21(5)
      A(7,1)=A(1,7)
      A(7,2)=A(2,7)
      A(7,5)=A(5,7)
      A(7,6)=A(6,7)
      A(7,7)=ALF11(3)*4.D0+PROD2*ALF21(5)
      A(7,8)=ALF11(4)*6.D0+PROD2*ALF21(6)
      A(8,1)=A(1,8)
      A(8,2)=A(2,8)
      A(8,5)=A(5,8)
      A(8,6)=A(6,8)
      A(8,7)=A(7,8)
      A(8,8)=ALF11(5)*9.D0+PROD2*ALF21(7)
      DSS=DS
      PDS=DSS*PRESSI*0.5D0
      DS2=DSS*DSS
      PS2DS=PDS*DS2
      PS4DS=PS2DS*DS2/80.D0
      PS2DS=PS2DS/12.D0
      PR=(RAYON(1)+4.D0*RAYON(2)+RAYON(3))*PDS/6.D0
      PSR=(-RAYON(1)+RAYON(3))*DSS*PDS/12.D0
      PS2R=(RAYON(1)+RAYON(3))*PDS*DS2/24.D0
      PS3R=PSR*DS2/4.D0
      A(1,5)=A(1,5)-PDS*0.5D0*COSP
      A(1,6)=A(1,6)+PR*0.5D0
      A(1,7)=A(1,7)-COSP*PS2DS*0.5D0+PSR
      A(1,8)=A(1,8)+1.5D0*PS2R
      A(5,1)=A(1,5)
      A(6,1)=A(1,6)
      A(7,1)=A(1,7)
      A(8,1)=A(1,8)
      A(2,5)=A(2,5)-PR*0.5D0
      A(2,6)=A(2,6)-PS2DS*0.5D0*COSP
      A(2,7)=A(2,7)+PS2R*0.5D0
      A(2,8)=A(2,8)-PS4DS*COSP*0.5D0+PS3R
      A(5,2)=A(2,5)
      A(6,2)=A(2,6)
      A(7,2)=A(2,7)
      A(8,2)=A(2,8)
      PROD1=PDS*SINP
      A(3,3)=A(3,3)-PROD1
      PROD2=PS2DS*SINP
      A(4,4)=A(4,4)-PROD2
      A(5,5)=A(5,5)-PROD1
      A(5,7)=A(5,7)-PROD2
      A(7,5)=A(5,7)
      A(6,6)=A(6,6)-PROD2
      PROD1=PS4DS*SINP
      A(6,8)=A(6,8)-PROD1
      A(8,6)=A(6,8)
      A(7,7)=A(7,7)-PROD1
      A(8,8)=A(8,8)-PROD1*DS2/5.6D0
C     CALL MULMAT (8,8,8,A,A4,S3)
      CALL MULMAT (S3,A,A4,8,8,8)
C     CALL MULMAT (8,8,8,A1,S3,A)
      CALL MULMAT (A,A1,S3,8,8,8)
      CALL ZDANUL(AM,64)
CCCCCCCCCC  CCC CCC
      AM(1,1)=ALF21(1)
      AM(3,3)=ALF21(1)
      AM(5,5)=ALF21(1)
      AM(1,2)=ALF21(2)
      AM(2,1)=ALF21(2)
      AM(3,4)=ALF21(2)
      AM(4,3)=ALF21(2)
      AM(5,6)=ALF21(2)
      AM(6,5)=ALF21(2)
      AM(2,2)=ALF21(3)
      AM(4,4)=ALF21(3)
      AM(6,6)=ALF21(3)
      AM(5,7)=ALF21(3)
      AM(7,5)=ALF21(3)
      AM(5,8)=ALF21(4)
      AM(8,5)=ALF21(4)
      AM(6,7)=ALF21(4)
      AM(7,6)=ALF21(4)
      AM(7,7)=ALF21(5)
      AM(6,8)=ALF21(5)
      AM(8,6)=ALF21(5)
      AM(7,8)=ALF21(6)
      AM(8,7)=ALF21(6)
      AM(8,8)=ALF21(7)
C     CALL MULMAT (8,8,8,AM,A4,S3)
      CALL MULMAT (S3,AM,A4,8,8,8)
C     CALL MULMAT (8,8,8,A1,S3,AM)
      CALL MULMAT (AM,A1,S3,8,8,8)
      AN2=AN*AN
      DO 30 NC=1,8
         DO 31 NCC=1,8
            A(NCC,NC)=A(NCC,NC)+AN2*AM(NCC,NC)
 31      CONTINUE
 30   CONTINUE
      PROD1=2.D0*COSP
      PROD2=2.D0*SINP
      CALL ZDANUL(AM,64)
      AM(1,3)=PROD1*ALF21(1)
      AM(1,4)=PROD1*ALF21(2)
      AM(2,3)=AM(1,4)
      AM(2,4)=PROD1*ALF21(3)
      AM(3,1)=AM(1,3)
      AM(3,2)=AM(2,3)
      AM(3,5)=PROD2*ALF21(1)
      AM(3,6)=PROD2*ALF21(2)
      AM(3,7)=PROD2*ALF21(3)
      AM(3,8)=PROD2*ALF21(4)
      AM(4,1)=AM(1,4)
      AM(4,2)=AM(2,4)
      AM(4,5)=AM(3,6)
      AM(4,6)=AM(3,7)
      AM(4,7)=AM(3,8)
      AM(4,8)=PROD2*ALF21(5)
      AM(5,3)=AM(3,5)
      AM(5,4)=AM(4,5)
      AM(6,3)=AM(3,6)
      AM(6,4)=AM(4,6)
      AM(7,3)=AM(3,7)
      AM(7,4)=AM(4,7)
      AM(8,3)=AM(3,8)
      AM(8,4)=AM(4,8)
      AM(3,5)=AM(3,5)-PDS
      AM(5,3)=AM(3,5)
      AM(3,7)=AM(3,7)-PS2DS
      AM(4,6)=AM(4,6)-PS2DS
      AM(7,3)=AM(3,7)
      AM(6,4)=AM(4,6)
      AM(4,8)=AM(4,8)-PS4DS
      AM(8,4)=AM(4,8)
C      CALL MULMAT (8,8,8,AM,A4,S3)
      CALL MULMAT (S3,AM,A4,8,8,8)
C     CALL MULMAT (8,8,8,A1,S3,AM)
      CALL MULMAT (AM,A1,S3,8,8,8)
      IF(AN.EQ.0) THEN
         COEF=2*XPI
      ELSE
         COEF=XPI
      ENDIF
      IROT=1
      DO 32 NC=1,8
         DO 33 NCC=1,8
            RE(IROT) =(A(NCC,NC)+AN*AM(NCC,NC))*2.D0+RE(IROT)
            RE(IROT) = RE(IROT)*COEF
            IROT=IROT+1
 33      CONTINUE
 32   CONTINUE
      RETURN
      END