Télécharger pdirec.eso

Retour à la liste

Numérotation des lignes :

  1. C PDIREC SOURCE CHAT 05/01/13 02:11:29 5004
  2. c---------------------------------------------------------------------
  3. c
  4. SUBROUTINE PDIREC (NDIME, NSTR1, SIGMA, SGPRI, BETAM)
  5. c
  6. c=====================================================================
  7. c =
  8. c This routine computes the eigenvectors of a structural tensor. =
  9. c =
  10. c Note: sigma = (Sxx, Syy, Szz, Sxy, Sxz, Syz) =
  11. c sgpri = (S11, S22, S33) for NSTR1=6 =
  12. c sgpri = (S11, S22, Szz) for NSTR1=4 =
  13. c =
  14. c=====================================================================
  15. IMPLICIT INTEGER(I-N)
  16. integer*4 ndime, nstr1
  17. real*8 betam (ndime,*), sigma (*), sgpri (*)
  18. c
  19. integer*4 i
  20. real*8 twopi, error, str11, str22, str12, denta,
  21. . angle, ca , ca2 , sa , sa2 , value,
  22. . dummy (3,3) , absnor
  23. c
  24. data twopi, error / 6.283185307179586 d0, 1.0 d-12 /
  25. c
  26. if (ndime .eq. 2) then
  27. str11 = sigma (1)
  28. str22 = sigma (2)
  29. str12 = sigma (4)
  30. c
  31. absnor=sqrt(str11**2 + str22**2 + str12**2)
  32. c
  33. denta = str11 - str22
  34. c
  35. if (abs (str12)/absnor .le. error) then
  36. angle = 0.0 d0
  37. if (denta .lt. 0.0 d0) angle = twopi * 0.25 d0
  38. else
  39. if (abs (denta)/absnor .le. error) then
  40. angle = twopi * 0.125 d0
  41. if (str12 .lt. 0.0 d0) angle = -angle
  42. else
  43. angle = 0.5 d0 * atan (2.0d0*str12/abs(denta))
  44. end if
  45. end if
  46. c
  47. ca = cos (angle)
  48. sa = sin (angle)
  49. ca2 = ca **2
  50. sa2 = sa **2
  51. value = 2.0 d0 * sa * ca * str12
  52. c
  53. sgpri (1) = ca2 * str11 + sa2 * str22 + value
  54. sgpri (2) = sa2 * str11 + ca2 * str22 - value
  55. sgpri (3) = sigma (3)
  56. c
  57. if (sgpri (1) .lt. sgpri (2)) then
  58. c
  59. value = ca
  60. ca = -sa
  61. sa = value
  62. c
  63. value = sgpri (1)
  64. sgpri (1) = sgpri (2)
  65. sgpri (2) = value
  66. end if
  67. c=====----------------------------------------------------------------
  68. c build transformation matrix |
  69. c=====----------------------------------------------------------------
  70. betam (1,1) = ca
  71. betam (2,2) = ca
  72. betam (2,1) = -sa
  73. betam (1,2) = sa
  74. else
  75. if (sgpri (1) .eq. 0.0 d0 .and. sgpri (2) .eq. 0.0 d0
  76. . .and. sgpri (3) .eq. 0.0 d0) then
  77. call zero (betam, 3, 3)
  78. do 10 i = 1, 3
  79. betam (i,i) = 1.0 d0
  80. 10 continue
  81. return
  82. end if
  83. c
  84. call pdire3 (sigma, sgpri, dummy)
  85. c
  86. betam (1,1) = dummy (1,1)
  87. betam (1,2) = dummy (2,1)
  88. betam (1,3) = dummy (3,1)
  89. betam (2,1) = dummy (1,2)
  90. betam (2,2) = dummy (2,2)
  91. betam (2,3) = dummy (3,2)
  92. betam (3,1) = dummy (1,3)
  93. betam (3,2) = dummy (2,3)
  94. betam (3,3) = dummy (3,3)
  95. end if
  96. c
  97. return
  98. end
  99.  
  100.  
  101.  

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