Télécharger cq2kp.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales