clearconsole()
using Plots
using LinearAlgebra
using DifferentialEquations
using Measures

#equation(d) = -7.65 * sind(d) + 9.87 * sind(2d + 206);
#plot(equation, 1:365)
Fact005(b)=1/sqrt((1-b^2)^2+(2*0.05*b)^2)
Fact010(b)=1/sqrt((1-b^2)^2+(2*0.10*b)^2)
Fact030(b)=1/sqrt((1-b^2)^2+(2*0.30*b)^2)
Fact0707(b)=1/sqrt((1-b^2)^2+(2*sqrt(2)/2*b)^2)

theta005(b)=atan(2*0.05*b,(1-b^2))
theta010(b)=atan(2*0.10*b,(1-b^2))
theta030(b)=atan(2*0.30*b,(1-b^2))
theta0707(b)=atan(2*sqrt(2)/2*b,(1-b^2))

SFact(b,x)=1/sqrt((1-b^2)^2+(2*x*b)^2)



plotFact=plot(Fact005,0:0.01:4,xtickfontsize=14,
ytickfontsize=14,linewidth=1.5,label="ζ=0.05",legendfontsize=12,
framestyle=:box,legend=:best)
plot!(plotFact,Fact010,0:0.01:4,linewidth=1.5,label="ζ=0.10")
plot!(plotFact,Fact030,0:0.01:4,linewidth=1.5,label="ζ=0.30")
plot!(plotFact,Fact0707,0:0.01:4,linewidth=1.5,label="ζ=0.707")
ylabel!("D",yguidefontsize=18)
xlabel!("β",xguidefontsize=18)
display(plotFact)
NameFig="DBeta.pdf"
savefig(NameFig)


plottheta=plot(theta005,0:0.01:4,xtickfontsize=14,
ytickfontsize=14,linewidth=1.5,label="ζ=0.05",legendfontsize=12,
framestyle=:box,legend=:bottomright)
plot!(plottheta,theta010,0:0.01:4,linewidth=1.5,label="ζ=0.10")
plot!(plottheta,theta030,0:0.01:4,linewidth=1.5,label="ζ=0.30")
plot!(plottheta,theta0707,0:0.01:4,linewidth=1.5,label="ζ=0.707")
ylabel!("̅θ",xguidefontsize=18)
xlabel!("β",xguidefontsize=18)
display(plottheta)
NameFig="thetaBeta.pdf"
savefig(NameFig)


plotWire=wireframe(0.01:0.05:4.0,0:0.05:1,SFact,zlims=(0,10),
clims=(0,10),cam=(45,45),zlabel=("D"),right_margin=3mm)
ylabel!("ζ",xguidefontsize=18)
xlabel!("β",xguidefontsize=18)
display(plotWire)
NameFig="WireFrame.pdf"
savefig(NameFig)