##calculate the orientation of the principal stress and deterimine stress components and orientation of the principal axes of strain
## initialization of variables
import math
import numpy
## Material constants
Ex=14700. ## M Pa
Ey=1000. ## M Pa
Ez=735. ## M Pa
Gxy=941. ## M Pa
Gxz=1147. ## M Pa
Gyz=103. ## M pa
Vxy=0.292
Vxz=0.449
Vyz=0.39
## Stresses at a point
Sxx=7. ## M pa
Syy=2.1 ## M Pa
Szz=-2.8 ##M Pa
Sxy=1.4 ## M Pa
Sxz=0. ##M Pa
Syz=0. ## M Pa
## part (a)
th=1/2.*math.atan(2.*Sxy/(Sxx-Syy))*180./math.pi
I1=Sxx+Syy+Szz
I2=Sxx*Syy-Sxy**2+Szz*Sxx-Sxz**2+Szz*Syy-Syz**2
M=numpy.matrix([[Sxx, Sxy, Sxz],
[Sxy, Syy, Syz],
[Sxz, Syz, Szz]])
I3=numpy.linalg.det(M)
p=([1, -I1, I2 ,-I3])
S=numpy.roots(p)
## results
print('Part (a) \n')
print"%s %.2f %s"%('The maximum principal stress is S1 = ',S[0],' M Pa')
print"%s %.2f %s"%('\n and occurs in direction th = ',th,' degrees')
print"%s %.2f %s %.2f %s "%('\n and the intermediate principal stress S2 = ',S[2],' M Pa'and 'occurs in the direction th = ',th+90,' degrees \n')
print"%s %.2f %s"%(' The minimum principal stress is S3 = Szz = ',S[1],' M Pa')
Ex=Ex*10**6
Ey=Ey*10**6
Ez=Ez*10**6
Gxy=Gxy*10**6
Gxz=Gxz*10**6
Gyz=Gyz*10**6
## part (b) is to find strains
Exx=Sxx/Ex-Vxy*Syy/Ey-Vxz*Szz/Ez
Eyy=-Vxy*Sxx/Ex+Syy/Ey-Vyz*Szz/Ez
Ezz=-Vxz*Sxx/Ex-Vyz*Syy/Ey+Szz/Ez
Exy=Sxy/Gxy
Exz=Sxz/Gxz
Eyz=Syz/Gyz
print('\n Part (b)')
print('\n The srains are')
print"%s %.2e %s %.2e %s %.2e %s"%('\n Exx = ',Exx,'' and'Eyy =',Eyy,''and 'Ezz = ',Ezz,'')
print"%s %.2e %s %.2d %s %.2d %s"%('\n Exy = ',Exy,''and ' Exz = ',Exz,''and 'Eyz = ',Eyz,'')
## Wrong Exx value in the textbook
th=1/2.*math.atan(Exy/(Exx-Eyy))
th=th*180./math.pi
th2=th+90.
print('\n part (c)')
print"%s %.2f %s %.2f %s "%('\n theta = ',th,'' and 'or theta=',th2,'degrees')
## Wrong theta too since Example given in textbook is wrong