lenssurfaceupper[x_,a_,b_] = Piecewise[{{Sqrt[a^2-(x-b)^2],x>0 && x ≤ a+b}},0];
lenssurfacelower[x_,a_,b_] = Piecewise[{{-Sqrt[a^2-(x-b)^2],x>0 && x ≤ a+b}},0];

plotlens[a_,b_,nlens_] :=Plot[{lenssurfaceupper[x,a,b],lenssurfacelower[x,a,b]},{x,0,a+b},AspectRatio->1,PlotRange->{{-nlens*a,nlens*a},{-nlens*a,nlens*a}},ImageSize->600,PlotStyle->Directive[Thickness[0],Black],Filling->Axis,FillingStyle->Opacity[0.5]];

xint[θ_,a_,b_] = Piecewise[{{((b+(Sqrt[b^2+(a^2-b^2)*(1+Tan[θ]^2)]))/(1+Tan[θ]^2)),θ≠ Pi/2 && θ ≠ -Pi/2},{0,θ== Pi/2 || θ == -Pi/2}}]

deriv[θ_,a_,b_] = ((b-xint[θ,a,b])/(Sqrt[a^2-(xint[θ,a,b]-b)^2]));

β1[θ_,a_,b_] = ArcTan[-deriv[θ,a,b]];
β2[θ_,a_,b_] = ArcCos[b/a*Sin[θ]]-θ;

α[θ_,a_,b_] = Pi/2 - ArcCos[b/a*Sin[θ]];

δ[θ_,a_,b_,nlens_,nair_] = ArcSin[nlens/nair*Sin[α[θ,a,b]]];
δ2[θ_,a_,b_,nlens_,nair_] = ArcSin[nlens/nair*b/a*Sin[θ]];

γ[θ_,a_,b_,nlens_,nair_] = (Pi/2 - β2[θ,a,b]- δ[θ,a,b,nlens,nair])

dγ[θ_,a_,b_,nlens_,nair_] = D[(Pi/2 - β2[θ,a,b]- δ[θ,a,b,nlens,nair]),θ]

chyp[θ_,a_,b_] = Sqrt[a^2+b^2+2*a*b*Sin[β2[θ,a,b]]]

ray1[x_,θ_,a_,b_,nlens_,nair_] =  Piecewise[{
{x*Tan[θ],x<Re[xint[θ,a,b]]},
{(x-xint[θ,a,b])*Tan[γ[θ,a,b,nlens,nair]]+chyp[θ,a,b]*Sin[θ],( (x>xint[θ,a,b]) && (Pi/2≥Re[δ[θ,a,b,nlens,nair]]≥-Pi/2))}}];

ray2[x_,θ_,a_,b_,nlens_,nair_] = Piecewise[{{(x-xint[θ,a,b])*Tan[γ[θ,a,b,nlens,nair]]+chyp[θ,a,b]*Sin[θ],(Pi/2≥Re[δ[θ,a,b,nlens,nair]]≥-Pi/2)}}];

plotray[θ_,a_,b_,nlens_,nair_] := Piecewise[{
{Piecewise[{
{Plot[ray1[x,θ,a,b,nlens,nair],{x,0,nlens*a},PlotRange->{{-nlens*a,nlens*a},{-3*a,3*a}},PlotStyle->Directive[Thickness[0.003],Blue]],-Pi/2 ≤ Re[γ[θ,a,b,nlens,nair]] ≤ Pi/2},
{Show[{{Plot[ray1[x,θ,a,b,nlens,nair],{x,0,xint[θ,a,b]},PlotRange->{{-nlens*a,nlens*a},{-3*a,3*a}},PlotStyle->Directive[Thickness[0.003],Blue]],Plot[ray2[x,θ,a,b,nlens,nair],{x,-nlens*a,xint[θ,a,b]},PlotRange->{{-nlens*a,nlens*a},{-3*a,3*a}},PlotStyle->Directive[Thickness[0.003],Blue]]}}], Abs[Re[γ[θ,a,b,nlens,nair]]] > Pi/2}
}],Abs[θ]≠Pi/2},
{Piecewise[{
{Show[{{
Plot[ray1[x,θ,a,b,nlens,nair],{x,0,nlens*a},PlotRange->{{-nlens*a,nlens*a},{-3*a,3*a}},PlotStyle->Directive[Thickness[0.003],Blue]],
Graphics[{Thickness[0.003],Blue,Line[{{0,0},{0,Sign[θ]*Sqrt[a^2-b^2]}}]}]
}}],-Pi/2 ≤ Re[γ[θ,a,b,nlens,nair]] ≤ Pi/2},
{Show[{{Graphics[{Thickness[0.003],Blue,Line[{{0,0},{0,Sign[θ]*Sqrt[a^2-b^2]}}]}],Plot[ray2[x,θ,a,b,nlens,nair],{x,-nlens*a,xint[θ,a,b]},PlotRange->{{-nlens*a,nlens*a},{-3*a,3*a}},PlotStyle->Directive[Thickness[0.003],Blue]]}}], Abs[Re[γ[θ,a,b,nlens,nair]]] > Pi/2}
}],Abs[θ]==Pi/2}
}]

focusline[x_,θ_,a_,b_,nlens_,nair_]= Piecewise[{{(x-xint[θ,a,b])*Tan[γ[θ,a,b,nlens,nair]]+chyp[θ,a,b]*Sin[θ],θ≠0 && Re[δ[θ,a,b,nlens,nair]]≤ Pi/2 && Re[δ[θ,a,b,nlens,nair]] ≥ -Pi/2}}];

xfoc[θ_,a_,b_,nlens_,nair_] = Piecewise[{{chyp[θ,a,b]*(Cos[θ]-((Sin[θ])/(Tan[γ[θ,a,b,nlens,nair]]))),θ≠ 0 && θ ≥ -Pi/2 && θ ≤  Pi/2}}];

dxfoc[θ_,a_,b_,nlens_,nair_] = D[chyp[θ,a,b]*(Cos[θ]-((Sin[θ])/(Tan[γ[θ,a,b,nlens,nair]]))),θ];

plotfocus[θ_,a_,b_,nlens_,nair_] := Piecewise[{{Plot[focusline[x,θ,a,b,nlens,nair],{x,xfoc[θ,a,b,nlens,nair],xint[θ,a,b]},PlotStyle->Directive[Thickness[0.0025],Dashed,Red]], (Element[xfoc[θ,a,b,nlens,nair],Reals] && Re[xfoc[θ,a,b,nlens,nair]]<=0) || (b≤ 0 && Element[xfoc[θ,a,b,nlens,nair],Reals])}},Plot[focusline[x,θ,a,b,nlens,nair],{x,-nlens*a,xint[θ,a,b]},PlotStyle->Directive[Thickness[0.0025],Dashed,Red]]];

ani = Manipulate[Show[{
plotlens[1,b,3.38],
Table[plotray[θ,1,b,3.38,1],{θ,-0.99*Pi/2,0.99*Pi/2,0.99*Pi/20}],Table[plotfocus[θ,1,b,3.38,1],{θ,-0.99*Pi/2,0.99*Pi/2,0.99*Pi/20}]}],{b,-0.98,0.99,0.01}]

Export["C:\\Users\\Jake\\Documents\\GFET\\Modelling\\hyperhemisphere.gif",ani]

Show[{
plotlens[1,-1/3.38,3.38],
Table[plotray[θ,1,-1/3.38,3.38,1],{θ,-Pi/2,Pi/2,Pi/20}],
Table[plotfocus[θ,1,-1/3.38,3.38,1],{θ,-Pi/2,Pi/2,Pi/20}]}]

Show[{
plotlens[1,1/3.38,3.38],
Table[plotray[θ,1,1/3.38,3.38,1],{θ,-Pi/2,Pi/2,Pi/20}],
Table[plotfocus[θ,1,1/3.38,3.38,1],{θ,-Pi/2,Pi/2,Pi/20}]}]

Plot[{180/Pi*γ[Pi/180*θ,1,1/3.38,3.38,1],180/Pi*δ[Pi/180*θ,1,1/3.38,3.38,1],180/Pi*β2[Pi/180*θ,1,1/3.38]-90},{θ,-90,90},GridLines->Automatic,PlotRange->{-120,120}]

δ[Pi/4,1,1/3.38,3.38,1]

Plot[180/Pi*γ[Pi/2,1,1/nlens,nlens,1],{nlens,1,5},PlotRange->{0,90}]

180/Pi*γ[Pi/2,1,1/3.38,3.38,1]

Manipulate[Plot[xfoc[θ,1,b,3.42,1],{θ,-Pi/2,Pi/2},PlotRange->{-5,1}],{b,0,1/3.42}]

Manipulate[Plot[Abs[1/dγ[Pi/180*θ,1,b,3.38,1]],{θ,-90,90},PlotRange->{0,1}],{b,-1,1}]

Manipulate[Plot[180/Pi*γ[t*Pi/180,1,b,3.38,1],{t,-90,90},PlotRange->{{-90,90},{-180,180}}],{b,-1,1}]

1/dγ[0,1,1/5,5,1]

Tan[Pi/2-ArcCos[b/a*Sin[θ]]]

Sin[ArcCos[b/a*Sin[θ]]-θ]

TrigExpand[%]

xfocideal[d_,t_]=Sqrt[1+d^2-2*d^2*Sin[t]^2+2*d*Cos[t]*Sqrt[1-d^2*Sin[t]^2]]*(Cos[t]-1/d*Sqrt[1-d^2*Sin[t]^2])

xfocideal[d,t]^2

dxfocideal[d_,t_]=D[Sqrt[1+d^2-2*d^2*Sin[t]^2+2*d*Cos[t]*Sqrt[1-d^2*Sin[t]^2]]*(Cos[t]-1/d*Sqrt[1-d^2*Sin[t]^2]),t]

FullSimplify[%]

Plot[{xfocideal[d,0],-1/d+d},{d,0,1}]