1 : C CQ2KP SOURCE CHAT 05/01/12 22:26:39 5004 2 : SUBROUTINE CQ2KP(COO,PRESSI,AN,ALF11,ALF12, 3 : 1 ALF21,ALF22,RAYON,RE,CGAUS,TGAUS,NGAUS,A1,A4,S3,AM,A) 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 : CCCCCCCCCCCCCC CALCUL DE LA MATRICE KP POUR LES COQ2 TIREE INCA 10 : C 11 : DIMENSION A4(8,8),S3(8,8),AM(8,8),A1(8,8) 12 : DIMENSION COO(3,2),A(8,8) 13 : DIMENSION ALF11(*),ALF12(*),ALF21(*),ALF22(*),RAYON(*),RE(*) 14 : DIMENSION CGAUS(*),TGAUS(*) 14 : CALL ZDANUL(RE,64) 15 : RAYON(1)=COO(1,1) 16 : RAYON(3)=COO(1,2) 17 : RAYON(2)=(RAYON(1)+RAYON(3))*0.5D0 18 : DR=RAYON(3)-RAYON(1) 19 : DZ=COO(2,2)-COO(2,1) 20 : DS=SQRT(DR*DR+DZ*DZ) 21 : COSP=DR/DS 22 : SINP=DZ/DS 23 : XS2=0.5D0*DS 23 : CALL RESOLV(A4,DS,COSP,SINP,S3,A) 23 : CALL ZDANUL(ALF12,5) 23 : CALL ZDANUL(ALF22,7) 24 : DO 9 IN=1,NGAUS 25 : RAY=RAYON(2)+XS2*COSP*TGAUS(IN) 26 : X=TGAUS(IN)*XS2 27 : X2=X*X 28 : X3=X2*X 29 : X4=X3*X 30 : ANS=0.D0 31 : ANT=0.D0 32 : ALF11(1)=1.D0 33 : ALF11(2)=X 34 : ALF11(3)=X2 35 : ALF11(4)=X3 36 : ALF11(5)=X4 37 : DO 12 NC=1,5 38 : ALF21(NC)=ALF11(NC) 39 : 12 CONTINUE 40 : ALF21(6)=X4*X 41 : ALF21(7)=X4*X2 42 : ANS=ANS*RAY*XS2*CGAUS(IN) 42 : CALL AMULX(ALF11,5,ANS,ALF12) 43 : ANT=ANT*XS2*CGAUS(IN)/RAY 43 : CALL AMULX(ALF21,7,ANT,ALF22) 44 : 9 CONTINUE 44 : CALL SHIFTD(ALF12,ALF11,5) 44 : CALL SHIFTD(ALF22,ALF21,7) 45 : DO 15 NC=1,8 46 : DO 13 NCC=1,8 47 : A1(NCC,NC)=A4(NC,NCC) 48 : 13 CONTINUE 49 : 15 CONTINUE 49 : CALL ZDANUL(A,64) 50 : PROD1=COSP*COSP 51 : PROD2=SINP*SINP 52 : PROD3=SINP*COSP 53 : A(1,1)=PROD1*ALF21(1) 54 : A(1,2)=PROD1*ALF21(2) 55 : A(1,5)=PROD3*ALF21(1) 56 : A(1,6)=PROD3*ALF21(2) 57 : A(1,7)=PROD3*ALF21(3) 58 : A(1,8)=PROD3*ALF21(4) 59 : A(2,1)=A(1,2) 60 : A(2,2)=ALF11(1)+PROD1*ALF21(3) 61 : A(2,5)=PROD3*ALF21(2) 62 : A(2,6)=PROD3*ALF21(3) 63 : A(2,7)=PROD3*ALF21(4) 64 : A(2,8)=PROD3*ALF21(5) 65 : A(3,3)=ALF21(1) 66 : A(3,4)=ALF21(2) 67 : A(4,3)=A(3,4) 68 : A(4,4)=ALF11(1)+ALF21(3) 69 : A(5,1)=A(1,5) 70 : A(5,2)=A(2,5) 71 : A(5,5)=PROD2*ALF21(1) 72 : A(5,6)=PROD2*ALF21(2) 73 : A(5,7)=PROD2*ALF21(3) 74 : A(5,8)=PROD2*ALF21(4) 75 : A(6,1)=A(1,6) 76 : A(6,2)=A(2,6) 77 : A(6,5)=A(5,6) 78 : A(6,6)=ALF11(1)+A(5,7) 79 : A(6,7)=ALF11(2)*2.D0+A(5,8) 80 : A(6,8)=ALF11(3)*3.D0+PROD2*ALF21(5) 81 : A(7,1)=A(1,7) 82 : A(7,2)=A(2,7) 83 : A(7,5)=A(5,7) 84 : A(7,6)=A(6,7) 85 : A(7,7)=ALF11(3)*4.D0+PROD2*ALF21(5) 86 : A(7,8)=ALF11(4)*6.D0+PROD2*ALF21(6) 87 : A(8,1)=A(1,8) 88 : A(8,2)=A(2,8) 89 : A(8,5)=A(5,8) 90 : A(8,6)=A(6,8) 91 : A(8,7)=A(7,8) 92 : A(8,8)=ALF11(5)*9.D0+PROD2*ALF21(7) 93 : DSS=DS 94 : PDS=DSS*PRESSI*0.5D0 95 : DS2=DSS*DSS 96 : PS2DS=PDS*DS2 97 : PS4DS=PS2DS*DS2/80.D0 98 : PS2DS=PS2DS/12.D0 99 : PR=(RAYON(1)+4.D0*RAYON(2)+RAYON(3))*PDS/6.D0 100 : PSR=(-RAYON(1)+RAYON(3))*DSS*PDS/12.D0 101 : PS2R=(RAYON(1)+RAYON(3))*PDS*DS2/24.D0 102 : PS3R=PSR*DS2/4.D0 103 : A(1,5)=A(1,5)-PDS*0.5D0*COSP 104 : A(1,6)=A(1,6)+PR*0.5D0 105 : A(1,7)=A(1,7)-COSP*PS2DS*0.5D0+PSR 106 : A(1,8)=A(1,8)+1.5D0*PS2R 107 : A(5,1)=A(1,5) 108 : A(6,1)=A(1,6) 109 : A(7,1)=A(1,7) 110 : A(8,1)=A(1,8) 111 : A(2,5)=A(2,5)-PR*0.5D0 112 : A(2,6)=A(2,6)-PS2DS*0.5D0*COSP 113 : A(2,7)=A(2,7)+PS2R*0.5D0 114 : A(2,8)=A(2,8)-PS4DS*COSP*0.5D0+PS3R 115 : A(5,2)=A(2,5) 116 : A(6,2)=A(2,6) 117 : A(7,2)=A(2,7) 118 : A(8,2)=A(2,8) 119 : PROD1=PDS*SINP 120 : A(3,3)=A(3,3)-PROD1 121 : PROD2=PS2DS*SINP 122 : A(4,4)=A(4,4)-PROD2 123 : A(5,5)=A(5,5)-PROD1 124 : A(5,7)=A(5,7)-PROD2 125 : A(7,5)=A(5,7) 126 : A(6,6)=A(6,6)-PROD2 127 : PROD1=PS4DS*SINP 128 : A(6,8)=A(6,8)-PROD1 129 : A(8,6)=A(6,8) 130 : A(7,7)=A(7,7)-PROD1 131 : A(8,8)=A(8,8)-PROD1*DS2/5.6D0 131 : C CALL MULMAT (8,8,8,A,A4,S3) 131 : CALL MULMAT (S3,A,A4,8,8,8) 131 : C CALL MULMAT (8,8,8,A1,S3,A) 131 : CALL MULMAT (A,A1,S3,8,8,8) 131 : CALL ZDANUL(AM,64) 132 : CCCCCCCCCC CCC CCC 133 : AM(1,1)=ALF21(1) 134 : AM(3,3)=ALF21(1) 135 : AM(5,5)=ALF21(1) 136 : AM(1,2)=ALF21(2) 137 : AM(2,1)=ALF21(2) 138 : AM(3,4)=ALF21(2) 139 : AM(4,3)=ALF21(2) 140 : AM(5,6)=ALF21(2) 141 : AM(6,5)=ALF21(2) 142 : AM(2,2)=ALF21(3) 143 : AM(4,4)=ALF21(3) 144 : AM(6,6)=ALF21(3) 145 : AM(5,7)=ALF21(3) 146 : AM(7,5)=ALF21(3) 147 : AM(5,8)=ALF21(4) 148 : AM(8,5)=ALF21(4) 149 : AM(6,7)=ALF21(4) 150 : AM(7,6)=ALF21(4) 151 : AM(7,7)=ALF21(5) 152 : AM(6,8)=ALF21(5) 153 : AM(8,6)=ALF21(5) 154 : AM(7,8)=ALF21(6) 155 : AM(8,7)=ALF21(6) 156 : AM(8,8)=ALF21(7) 156 : C CALL MULMAT (8,8,8,AM,A4,S3) 156 : CALL MULMAT (S3,AM,A4,8,8,8) 156 : C CALL MULMAT (8,8,8,A1,S3,AM) 156 : CALL MULMAT (AM,A1,S3,8,8,8) 157 : AN2=AN*AN 158 : DO 30 NC=1,8 159 : DO 31 NCC=1,8 160 : A(NCC,NC)=A(NCC,NC)+AN2*AM(NCC,NC) 161 : 31 CONTINUE 162 : 30 CONTINUE 163 : PROD1=2.D0*COSP 164 : PROD2=2.D0*SINP 164 : CALL ZDANUL(AM,64) 165 : AM(1,3)=PROD1*ALF21(1) 166 : AM(1,4)=PROD1*ALF21(2) 167 : AM(2,3)=AM(1,4) 168 : AM(2,4)=PROD1*ALF21(3) 169 : AM(3,1)=AM(1,3) 170 : AM(3,2)=AM(2,3) 171 : AM(3,5)=PROD2*ALF21(1) 172 : AM(3,6)=PROD2*ALF21(2) 173 : AM(3,7)=PROD2*ALF21(3) 174 : AM(3,8)=PROD2*ALF21(4) 175 : AM(4,1)=AM(1,4) 176 : AM(4,2)=AM(2,4) 177 : AM(4,5)=AM(3,6) 178 : AM(4,6)=AM(3,7) 179 : AM(4,7)=AM(3,8) 180 : AM(4,8)=PROD2*ALF21(5) 181 : AM(5,3)=AM(3,5) 182 : AM(5,4)=AM(4,5) 183 : AM(6,3)=AM(3,6) 184 : AM(6,4)=AM(4,6) 185 : AM(7,3)=AM(3,7) 186 : AM(7,4)=AM(4,7) 187 : AM(8,3)=AM(3,8) 188 : AM(8,4)=AM(4,8) 189 : AM(3,5)=AM(3,5)-PDS 190 : AM(5,3)=AM(3,5) 191 : AM(3,7)=AM(3,7)-PS2DS 192 : AM(4,6)=AM(4,6)-PS2DS 193 : AM(7,3)=AM(3,7) 194 : AM(6,4)=AM(4,6) 195 : AM(4,8)=AM(4,8)-PS4DS 196 : AM(8,4)=AM(4,8) 196 : C CALL MULMAT (8,8,8,AM,A4,S3) 196 : CALL MULMAT (S3,AM,A4,8,8,8) 196 : C CALL MULMAT (8,8,8,A1,S3,AM) 196 : CALL MULMAT (AM,A1,S3,8,8,8) 197 : IF(AN.EQ.0) THEN 198 : COEF=2*XPI 199 : ELSE 200 : COEF=XPI 201 : ENDIF 202 : IROT=1 203 : DO 32 NC=1,8 204 : DO 33 NCC=1,8 205 : RE(IROT) =(A(NCC,NC)+AN*AM(NCC,NC))*2.D0+RE(IROT) 206 : RE(IROT) = RE(IROT)*COEF 207 : IROT=IROT+1 208 : 33 CONTINUE 209 : 32 CONTINUE 210 : RETURN 211 : END 212 : 213 : 214 : 215 : 216 : 217 :
© Cast3M 2003 - All rights reserved.
Disclaimer