th=[0:.01:pi];
f=th-atan((1-cos(th))./sin(th));

plot(th,f)
