y,a=var('y,a')
phase(y,a) = (180/pi*(2*arccot((1 - 4*a^2)*csc((pi*y)/180)^2/a^3)/a + (1 - cos((pi*y)/180)))).full_simplify()
p=plot(phase(y,2/3),(y,0,90),thickness=3,legend_label='a=2/3',linestyle='--')+plot(phase(y,4/3),(y,0,90),thickness=3,legend_label='a=4/3',linestyle=':',color='black')+plot(phase(y,1),(y,0,90),thickness=3,legend_label='a=1',color='red')
p.set_legend_options(title='Phase Error (in degrees)')
p.save(filename='junk_sage.pdf',frame=True,gridlines=True,axes_labels=['dip angle',''])