split get_geocentric_angle in two

[fgnievinski]

the calculation of the two geocentric angles is completely unrelated, they should not be part of the same function.

and we don't need to calculate both when the user only needs one.

https://github.com/vitorhjr/geo-alhazen/blob/main/geo-alhazen-aux/get_geocentric_angle.m

function [geo_ang,geo_ang_at] = get_geocentric_angle (Ha,Ht,e,g,Rs)

% Returns geocentric angles, where phi is the geocentric angle between 
% antenna and transmitter and phi1 is the geocentric angle between antenna 
% and reflection point
%
% INPUT:
% - Ha: Antenna height (in meters)
% - Ht: Transmitter/satelitte height (in meters)
% - e: Elevation angle of transmitter (in degree)
% - g: Grazing angle (in degree)
% - Rs: Earth radius (in meters)
%
% Obs: 
% - Rt and is only for phi_at, otherwise use h2=0 and e=0;
% - g is only for phi_as, otherwise use g=0.

geo_ang = get_geocentric_angle_sfc (g,Ha,Rs);

if (nargin < 2),  return;  end

geo_ang,geo_ang_at = get_geocentric_angle_trans (e,Ha,Ht,Rs)

end

function geo_ang_sfc = get_geocentric_angle_sfc (Ha,g,Rs)

% Returns geocentric angle between antenna and reflection point
%
% OUTPUT:
% - geo_ang_sfc: Geocentric angle antenna-reflection point
% 
% INPUT:
% - Ha: Antenna height (in meters)
% - g: Grazing angle (in degrees)
% - Rs: Earth surface radius (in meters)

%% Radius:
Ra = Rs+Ha;

%% Geocentric Angle:
geo_ang_sfc = acosd((Rs./Ra).*cosd(g))-g;

end

function [geo_ang_at, D] = get_geocentric_angle_trans (Ha,e,Ht,Rs)

% Returns geocentric angle between antenna and transmitter
%
% INPUT:
% - Ha: Antenna height (in meters)
% - Ht: Transmitter/satellite height (in meters)
% - e: Elevation angle of transmitter (in degrees)
% - Rs: Earth surface radius (in meters)

%% Radii:
Ra = Rs+Ha;
Rt = Rs+Ht;

%% Normalize radii to values closer to 1 
% (to avoid overflow when squaring, R^2):
R0 = Rs;
Ra = Ra./R0;
Rt = Rt./R0;
Rs = Rs./R0;

% Coefficients of quadratic polynomial
a=1;
b=-2*Ra.*cosd(90+e);
c=Ra.^2-Rt.^2;

% Delta of quadratic
delta=b.^2-4.*a.*c;

% Roots
r1=(-b+delta.^0.5)./(2*a);
r2=(-b-delta.^0.5)./(2*a);

% Choose the correct root
D=max(r1,r2);

% Geocentric angle
temp = (D.^2 -Ra.^2 -Rt.^2)./(-2*Ra.*Rt);
geo_ang_at=acos(temp);

% Denormalization
D = D*R0;

end

[fgnievinski]

I have also removed the definitions of default constants because it should be done in only one place to avoid bugs due to inconsistency -- indeed, here we had Ht = 20e6 while in get_satellite_height.m we had a different value: Ht = 20.2e6; https://github.com/vitorhjr/geo-alhazen/blob/main/geo-alhazen-aux/get_satellite_height.m

[fgnievinski]

I've tried to normalize radii to values closer to 1 to avoid overflow when squaring, R^2:

R0 = Rs;
Ra = Ra./R0;
Rt = Rt./R0;
Rs = Rs./R0;  %=1

[vitorhjr]

I split the function get_geocentric_angle into two new functions following the suggestions above:

• get_geocentric_angle_trans.m
• get_geocentric_angle_sfc.m

I updated the usages of the get_geocentric_angle inside the functions.

[vitorhjr]

The two new functions were tested and presented the expected results.

It was interesting that the normalized radii method generated better results than the older version when compared to the vectorial method.