sx1 = 1
sy1 = 2
sz1 = 3
txy1 = 4
tyz1 = 5
tzx1 = 6

sx2 = 7
sy2 = 8
sz2 = 9
txy2 = 10
tyz2 = 11
tzx2 = 12

IF $1=1
sigmax = sx1
sigmay = sy1
sigmaz = sz1
tauxy = txy1
tauyz = tyz1
tauzx = tzx1
ENDIF
IF $1=2
sigmax = sx2
sigmay = sy2
sigmaz = sz2
tauxy = txy2
tauyz = tyz2
tauzx = tzx2
ENDIF
IF $1=3
sigmax = sx1 +  sx2
sigmay = sy1 +  sy2
sigmaz = sz1 +  sz2
tauxy = txy1 + txy2
tauyz = tyz1 + tyz2
tauzx = tzx1 + tzx2
ENDIF

I1 = sigmax + sigmay + sigmaz
I2 = sigmax*sigmay + sigmay*sigmaz + sigmaz*sigmax - tauxy^2 - tauyz^2 - tauzx^2
I3 = sigmax*sigmay*sigmaz - sigmax*tauyz^2 - sigmay*tauzx^2 - sigmaz*tauxy^2 + 2*tauxy*tauyz*tauzx

c1 = sqrt(abs(I1^2 - 3*I2))
phi = 1.0/3.0 * acos((2.0*I1^3 - 9.0*I1*I2 + 27.0*I3)/(2.0*c1^3))
c2 = I1/3.0
c3 = 2.0/3.0 * c1

sigma1 = c2 + c3 * cos(phi)
sigma2 = c2 + c3 * cos(phi - 2.0*pi/3.0)
sigma3 = c2 + c3 * cos(phi - 4.0*pi/3.0)

PRINT I1 I2 I3 phi sigma1 sigma2 sigma3