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 76 : call zero (betam, 3, 3) 77 : do 10 i = 1, 3 78 : betam (i,i) = 1.0 d0 79 : 10 continue 80 : return 81 : end if 82 : c 82 : call pdire3 (sigma, sgpri, dummy) 83 : c 84 : betam (1,1) = dummy (1,1) 85 : betam (1,2) = dummy (2,1) 86 : betam (1,3) = dummy (3,1) 87 : betam (2,1) = dummy (1,2) 88 : betam (2,2) = dummy (2,2) 89 : betam (2,3) = dummy (3,2) 90 : betam (3,1) = dummy (1,3) 91 : betam (3,2) = dummy (2,3) 92 : betam (3,3) = dummy (3,3) 93 : end if 94 : c 95 : return 96 : end 97 : 98 :
© Cast3M 2003 - All rights reserved.
Disclaimer