clc
clear all
close all
dt = pi/1000;
theta=(0 : dt : 2*pi);
Nu=0.3;
f1=cos(theta/2).*(1-sin(theta/2).*sin((3*theta)/2));
f11=cos((theta+dt)/2).*(1-sin((theta+dt)/2).*sin((3*(theta+dt))/2));
f2=cos(theta/2).*(1+sin(theta/2).*sin((3*theta)/2));
f22=cos((theta+dt)/2).*(1+sin((theta+dt)/2).*sin((3*(theta+dt))/2));
f3=cos(theta/2).*sin(theta/2).*cos((3*theta)/2);
f33=cos((theta+dt)/2).*sin((theta+dt)/2).*cos((3*(theta+dt))/2);
r_pstress=(f1-f2).^2+f2.^2+f1.^2+6.*(f3.^2);
r_pstress1=(f11-f22).^2+f22.^2+f11.^2+6.*(f33.^2);
r_pstrain=(f1-f2).^2+(f2-Nu.*(f1+f2)).^2+(f1-Nu.*(f1+f2)).^2+6.*(f3.^2);
r_pstrain1=(f11-f22).^2+(f22-Nu.*(f11+f22)).^2+(f11-Nu.*(f11+f22)).^2+6.*(f33.^2);
polar(theta,r_pstress);
hold on
polar(theta,r_pstrain);
hold off

area_pstress = 0;
for theta = 0 : dt : (2*pi)-dt
    area_pstress = area_pstress + (0.5 * r_pstress[theta/dt+1] * r_pstress1[theta/dt+1] * sin(dt));
end

area_pstrain = 0;
for theta = 0 : dt : (2*pi)-dt
    area_pstrain = area_pstrain + (0.5 * r_pstrain[theta/dt+1] * r_pstrain1[theta/dt+1] * sin(dt));
end