x_pla = Ha./tand(e); % Reflection point x axis on plane
if e~=0
[x_sph] = fminsearch(f, x_pla,opt); % Reflection point x axis on sphere by Fermat's Principle
elseif e==0
[x_sph] = fminsearch(f, 0, opt);
end
with:
x_pla = Ha./tand(e); % Reflection point x axis on plane
if (e==0), x_pla = 0; end
[x_sph] = fminsearch(f, x_pla,opt); % Reflection point x axis on sphere by Fermat's Principle
replace:
with: