Infinite distance (trigonometric): non-iterative, too?

[fgnievinski]

it seems the solution should not be restricted to the iterative or numerical case -- it'd be possible to obtain a trigonometric solution for infinite distance based on a preliminary solution for finite distance using analytical algorithms, no?

another question: does the infinite distance assumption keep the finite-distance specular point unchanged (position coordinates and geocentric angle and arc_len)?

[fgnievinski]

pedir ajuda para Fujimura, Vuorinen et al.

[fgnievinski]

rascunho:

Dear Dr. Vuorinen,

we'd like to respectfully ask if your research group could help us with the following mathematical problem.

• Would it be possible to simplify the quartic involved in the Ptolemy-Alhazen problem by assuming the source is at infinite distance, i.e., defined only by the direction angles?

In our application, the transmitting source is a GPS satellite at 20,000 km altitude and the receiving antenna is only a few meters above the Earth's surface, assumed spherical with 6,000 km radius.

[vitorhjr]

I would only change the information about Earth radius, that we use the average radius of approximately 6370 km. Then, the email draft is:

Dear Dr. Vuorinen,

We'd like to respectfully ask if your research group could help us with the following mathematical problem.

• Would it be possible to simplify the quartic involved in the Ptolemy-Alhazen problem by assuming the source is at infinite distance, i.e., defined only by the direction angles?

In our application, the transmitting source is a GPS satellite at 20,000 km altitude and the receiving antenna is only a few meters above the Earth's surface, assuming a sphere with a 6,370 km radius.

Regards,

Vitor Hugo de Almeida Junior
PhD student at Federal University of Rio Grande do Sul

Do you agree with this draft @fgnievinski? I will send it scheduled for Monday.

[fgnievinski]

parece bom

…

On Fri., May 5, 2023, 21:13 Vitor Hugo de Almeida Junior, < ***@***.***> wrote: I would only change the information about Earth radius, that we use the average radius of approximately 6370 km. Then, the email draft is: Dear Dr. Vuorinen, We'd like to respectfully ask if your research group could help us with the following mathematical problem. - Would it be possible to simplify the quartic involved in the Ptolemy-Alhazen problem by assuming the source is at infinite distance, i.e., defined only by the direction angles? In our application, the transmitting source is a GPS satellite at 20,000 km altitude and the receiving antenna is only a few meters above the Earth's surface, assuming a sphere with a 6,370 km radius. Regards, Vitor Hugo de Almeida Junior PhD student at Federal University of Rio Grande do Sul Do you agree with this draft @fgnievinski <https://github.com/fgnievinski>? I will send it scheduled for Monday. — Reply to this email directly, view it on GitHub <#31 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABGJHBAAIMGCOMYVYEO6P2DXEWJR7ANCNFSM6AAAAAATZ46NGY> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[vitorhjr]

Response from Dr. Vuorinen:

Dear Mr. Vitor Hugo De Almeida Junior,

Thank you for your email. We will think about this and inform you about our ideas, if any.

Best regards, M Vuorinen

[vitorhjr]

Response from Dr. Vuorinen

Dear Mr. Vitor Hugo De Almeida Junior,

The solution to the "source at infinite-distance" Ptolemy-Alhazen problem that you posed is attached. We are currently working on a manuscript on related questions will perhaps include
this result there.

Please keep me informed about your work on this topic and whether you will use
this result in your work.

Sincerely, Matti Vuorinen

cal20230511.pdf

[vitorhjr]

Dr. Vuorinen also send the following:

Dear Mr. Vitor Hugo De Almeida Junior,

Good to hear. I have a Mathematica notebook on this and I can send it to you
if you are interested.

Good luck for your research!

M Vuorinen

[vitorhjr]

Implementation of Dr. Vuorinen's mathematical model

The mathematical model is based on a complex unit circle. The incident ray reach the sphere from the right side parallel to the real axis.

• $f=r e^{i\theta}$ is the receiver, where r is the normalized radius of the receiver. In our case, we can assume r as $r=R_a/R_0$, which Ra is the radius of the antenna and R0 is the radius of the sphere.
• w is the reflection point. It's solved through the quartic polinomial $r e^{-i\theta}w^4-w^3+w-r e^{i\theta}$.
• $\tilde{w} = w+1$. It lies in the same line of w and is parallel to the real axis
• The polynomial has four roots. We can choose the one that satisfies $e^{i\phi}$
• $\phi$ is the phase angle of the reflection point and $0&lt; \phi &lt;\pi/2$
• As $\phi$ is very close $\theta$, we can define $\phi=\theta$. Thus, we define this as a condition to get the correct root of the polynomial
• The specular point is obtained through the equation: $pos_spec = (w/(e^{-i \pi/2-\theta}))*R_0-[0 R_0i]$

Matlab algorithm

R0 = get_earth_radius (); % sphere radius
Ra = R0 + Ha; % antenna radius
theta = deg2rad(90-e); % phase angle of the antenna
ra = Ra./R0; % normalized antenna radius

% Polynomial coefficients
c4 = ra.*exp(-1i.*theta);
c3 = -1;
c2 = 0;
c1 = 1;
c0 = -ra.*exp(1i.*theta);

ws = roots ([c4 c3 c2 c1 c0]); % roots

phis = angle(ws); % phase angle of the specular point
ind = argmin(abs(phis-theta)); % index of the correct root
w = ws (ind); % correct root
pos_spec_complex_local = w./exp(-1i*(pi/2-theta)).*R0-complex(0,R0); % specular position with local complex coordinates
x_complex_local = real(pos_spec_complex_local);
y_complex_local = imag(pos_spec_complex_local);

[vitorhjr]

I implement Vuorinen's model in the geo-alhazen repository in the following commit (ba83c99).

• get_reflection_spherical_vuorinen: function with the mathematical model
• get_spherical_reflection: implemented Vuorinen model in the global function to retrieve the specular reflection on a sphere
• demo_vuorinen: demonstration of the get_reflection_spherical_vuorinen and a comparison with Fujimura model.

[vitorhjr]

Questions about this model:

• Is the elevation zero the limit for this model?
• What do we have to do in case we have two roots with the same condition? Could it occur?

[vitorhjr]

I started discussions about Vuorinen's algorithm and results at #38.

[fgnievinski]

here's a revised implementation, with grazing angle and interferometric delay more consistent with the assumption of infinite distance; the calculation of the specular position was unaffected:

function [graz_ang, geo_ang_as, x_spec, y_spec, dx_trans, dy_trans] = get_reflection_spherical_vuorinen (e, Ha, Ht, Rs)
% GET_REFLECTION_SPHERICAL_VUORINEN Calculates reflection on spherical 
% surface based on Vuorinen (2023) equations (internal document - filename cal20230511.tex)
%
% On the Ptolemy-Alhazen problem: source at infinite distance. 2023, May
% 11. Internal document. 
% 
% INPUT:
% - Ha: antenna/receiver height (in meters)
% - e: elevation angle (in radians)
% - Rs: Earth radius (in meters)
% 
% NOTE: transmitter/satellite height is assumed infinite.
% 
% OUTPUT:
% - x_spec, y_spec: reflection point in local frame (vector, in meters)
% - dx_trans, dy_trans: transmitter direction in local frame (unit vector, in meters)
% - graz_ang: grazing angle of spherical reflection that satisfies Snell's Law (in degrees)
% - geo_ang_as: geocentric angle between antenna and reflection point (in degrees) 

% Antenna radius
Ra = Rs+Ha;

% Normalized antenna radius
ra = Ra./Rs;

% Phase angle
theta = deg2rad(90-e);

% Quartic polynomial coefficients
c4 = ra.*exp(-1i.*theta);
c3 = -1;
c2 = 0;
c1 = 1;
c0 = -ra.*exp(1i.*theta);

% Polynomial roots
ws = roots ([c4 c3 c2 c1 c0]);

% Candidate phase angles
phis = angle(ws);

% Phase angle at reflection point:
% (TODO: test interf. delay, not phase angle.)
ind = argmin(abs(phis-theta));
phi = phis(ind);
w = ws (ind);

% Geocentric angle at reflection point:
geo_ang_as = deg2rad (theta-phi);

% Reflection point in a quasigeocentric frame
pos_spec_complex = w./exp(-1i*(pi./2-theta)).*Rs-complex(0,Rs);
X_spec = real(pos_spec_complex);
Y_spec = imag(pos_spec_complex);
pos_spec_geo = [Y_spec Y_spec];

% Satellite direction (unit vector)
% in either local or quasigeocentric frames:
dx_trans = cosd(e);
dy_trans = sind(e);
dir_trans = [dx_trans, dy_trans];

% Grazing angle
graz_ang = get_grazing_angle_infinite (pos_ant_geo, pos_spec_geo, dir_trans);

function [g,ga,gt] = get_grazing_angle_infinite (pos_ant, pos_spec, dir_sat, pos_geocenter)

% Return grazing angle calculated by vector relation
% 
% INPUT:
% - pos_ant: position vector of antenna (x,y)
% - pos_spec: position vector of reflection point (x,y)
% - pos_sat: direction vector of satellite (dx,dy)
% - pos_geocenter: position vector of geocenter (x,y), optional
% 
% OUTPUT:
% - g: grazing angle (in degrees)
% - ga: left grazing angle (in degrees)
% - gt: right grazing angle (in degrees)

% relative position vector, from antenna to surface:
pos_sfc_ant = pos_spec - pos_ant;

% distance from antenna to surface:
dist_ant_sfc = norm_all(pos_ant_sfc);

% unit direction vectors, from surface to antenna or satellite:
dir_ant_sfc = pos_ant_sfc ./ norm_all(pos_ant_sfc);

% bistatic scattering angle at the surface, between antenna and satellite:
rst = acosd(dot_all(dir_ant_sfc, dir_sat));

% grazing angle:
g = 90 - 0.5*rst;

if (nargout < 2),  return;  end

% vertical direction at the surface reflection point:
pos_sfc_geocenter = pos_spec - pos_geocenter;
dir_vert = pos_sfc_geocenter ./ norm_all(pos_sfc_geocenter);

% left and right grazing angles
ga = 90 - acosd(dot_all(dir_ant_sfc, dir_vert));
gt = 90 - acosd(dot_all(dir_vert, dir_sat_sfc));

end

function [delay, graz_ang, arc_len, slant_dist, x_spec, y_spec, x_trans, y_trans, elev_spec, delay_direct] = ...
            get_spherical_reflection (e, Ha, Ht, Rs, algorithm, trajectory, frame)
% GET_SPHERICAL_REFLECTION  Calculates specular reflection on spherical surface.
%
% INPUT:
% - e: transmitter elevation angles (matrix; in degrees) 
% - Ha: receiver antenna height (matrix; in meters)
% Notes: 
% - transmitter elevation angle is defined 90� at zenith or zero at the tangent plane horizon, as seen from the receiver antenna; it may be negative.
% - transmitter elevation angle is defined 90° at zenith or zero at the tangent plane horizon, as seen from the receiver antenna; it may be negative.
% - matrix input may also be a vector or a scalar
% - non-scalar input must have the same size
% 
% OUTPUT:
% - delay: interferometric propagation delay on spherical surface (matrix, in meters)
% - graz_ang: grazing angle at reflection point to transmitter or to receiver (matrix, in degrees)
% - arc_len: arc lenght between subreceiver point and reflection point (matrix, in meters)
% - slant_dist: slant distance between receiver and reflection point (matrix, in meters)
% - x_spec, y_spec: reflection point coordinates (matrices, in meters)
% - x_trans, y_trans: transmitter point coordinates (matrices, in meters)
% - elev_spec: elevation angle from antenna to reflection point (matrix, in degrees)
% - delay_direct: direct distance from antenna to satellite
% 
% Notes: 
% - grazing angle is defined 90° at zenith or zero at the tangent plane horizon, as seen from the reflection point; 
%   it is the same to the transmitter or to the receiver, as it satisfies Snell's law; it is never negative.
% - output will be a matrix of the same size as input.
% 
% OPTIONAL INPUT
% - Ht: Transmitter/satelitte height (scalar, in meters)
% - Rs: Earth surface radius (scalar, in meters)
% - algorithm: (char), see source code for details
%       'fujimura'
%       'martin-neira'
%       'miller&vegh'
%       'helm'
%       'fermat'
%       'itu' (under research)
%       'line-sphere fermat' (under research)
%       'circle inversion' (under research)
%       'vuorinen' (finite distance)
% - trajectory: (char), see source code for details
%       'orbital'  % orbital trajectory (constant geocentric distance)
%       'horizontal'  % horizontal trajectory (constant y-axis coordinate)
%       'circular'  % circular trajectory around the receiving antenna (constant direct distance)
% - frame: (char) coordinate reference frame ('local' - default - or 'quasigeo')

    if (nargin < 1),  e = [];  end
    if (nargin < 2),  Ha = [];  end
    if (nargin < 3),  Ht = [];  end
    if (nargin < 4),  Rs = [];  end
    if (nargin < 5),  algorithm = [];  end
    if (nargin < 6),  trajectory = [];  end
    if (nargin < 7),  frame = [];  end

    if isempty(e),   e  = 45;  end
    if isempty(Ha),  Ha = 10;  end
    if isempty(Rs),  Rs = get_earth_radius();  end
    if isempty(Ht),  Ht = get_satellite_height();  end
    if isempty(algorithm),  algorithm = 'fujimura';  end
    if isempty(trajectory),  trajectory = 'orbital';  end
    if isempty(frame),  frame = 'local';  end

    %% Check input size compatibility:
    assert(isscalar(Ht))
    assert(isscalar(Rs))
    siz = size(e);  e = e(:);  n = prod(siz);
    siz2 = size(Ha);  Ha = Ha(:);  n2 = prod(siz2);
    assert(n2==1 || (n2==n && isequal(siz2, siz)))
    %assert(isscalar(Ha) || isequal(size(Ha), size(e)))

    %% Select algorithm
    is_finite = true;
    switch lower(algorithm)
        case {'fujimura'}
            f = @get_reflection_spherical_fujimura;
        case {'miller','miller&vegh','millerandvegh'}
            f = @get_reflection_spherical_miller;
        case {'martinneira','martin-neira'}
            f = @get_reflection_spherical_martinneira;
        case {'helm'}
            f = @get_reflection_spherical_helm;
         case {'fermat','numerical'}
            f = @get_reflection_spherical_fermat;
        case {'vuorinen'}
            %f = @get_reflection_spherical_vuorinen;
            f = @(ei, Hai, Htsi, Rs) get_reflection_spherical_vuorinen(ei, Hai, Rs);
            is_finite = false;
        otherwise
            error('Unknown algorithm "%s"', char(algorithm));
    end

    %% Get ancillary results:
    Hts = get_satellite_trajectory (e, Ha, Ht, Rs, trajectory);
    ehor = get_horizon_elevation_angle (Ha, Rs);
    
    %% Pre-allocate memory:
    temp = NaN(n,1);
    graz_ang  = temp;
    geo_ang   = temp;
    x_spec    = temp;
    y_spec    = temp;
    x_trans   = temp;
    y_trans   = temp;

    %% Run selected algorithm
    i2 = 1;
    for i=1:n
        if (n2>1),  i2 = i;  end
        if (e(i) < ehor(i2)),  continue;  end
        [graz_ang(i), geo_ang_as(i), x_spec(i), y_spec(i), x_trans(i), y_trans(i)] = f(...
            e(i), Ha(i2), Hts(i), Rs);
    end

    %% Additional parameters
    [delay, arc_len, slant_dist, elev_spec, delay_direct] = get_spherical_reflection_extra (...
        n2, Ha, Rs, geo_ang_as, x_spec, y_spec, x_trans, y_trans, is_finite);

    %% Reshape output matrices as in input matrices:
    delay = reshape(delay, siz);
    graz_ang = reshape(graz_ang, siz);
    arc_len = reshape(arc_len, siz);
    slant_dist = reshape(slant_dist, siz);
    x_spec = reshape(x_spec, siz);
    y_spec = reshape(y_spec, siz);
    x_trans = reshape(x_trans, siz);
    y_trans = reshape(y_trans, siz);
    elev_spec = reshape(elev_spec, siz);
    delay_direct = reshape(delay_direct, siz);

    %% Optionally, convert from reference frame:
    if strcmpi(frame, 'local'),  return;  end
    [x_spec, y_spec] = get_quasigeo_coord (x_spec, y_spec, Rs);
    [x_trans, y_trans] = get_quasigeo_coord (x_trans, y_trans, Rs);    
end

%%
function [delay, arc_len, slant_dist, elev_spec, delay_direct] = get_spherical_reflection_extra (n2, Ha, Rs, geo_ang_as, x_spec, y_spec, x_trans, y_trans, is_finite)

    % Arc Length from subreceiver point to reflection point:
    arc_len = deg2rad(geo_ang_as)*Rs;

    % Receiving antenna position vector in local frame:
    pos_ant = [zeros(n2,1) Ha]; 

    % Relative position vectors:
    pos_spec = [x_spec(:) y_spec(:)];
    pos_trans = [x_trans(:) y_trans(:)];
    pos_trans_spec = pos_trans - pos_spec;
    pos_ant_spec   = pos_ant   - pos_spec;
    pos_ant_trans  = pos_ant   - pos_trans;

    % Slant distance from reflection point to receiver:
    slant_dist = norm_all(pos_ant_spec);

    % Interferometric propagation delay:
    % (TODO: split into a separate function)
    if is_finite
        delay_reflect_out = slant_dist;
        delay_reflect_inc = norm_all(pos_trans_spec);
        delay_direct = norm_all(pos_ant_trans);
        delay_reflect = delay_reflect_inc + delay_reflect_out;
        delay = delay_reflect - delay_direct;
    else
        delay_reflect_out = slant_dist;
        [~, delay_direct_aux] = proj_pt_line (pos_spec, pos_ant, dir_sat)
        delay = pos_ant_spec - delay_direct_aux;
        delay_direct = inf(size(delay));
    end

    % Reflection elevation angle:
    dir_rec_spec = pos_ant_spec./delay_reflect_out;
    dir_vert = [0 1];
    elev_spec = -90 + acosd(dot_all(-dir_vert, -dir_rec_spec));
    debugit = 0;
    if debugit
      elev_specb = e - 2*graz_ang;
      elev_specb(e==90) = -90;
      fpk e erb-elev_spec  % DEBUG
    end

end

%%
function Hts = get_satellite_trajectory (e, Ha, Ht, R0, trajectory)
    assert(isvector(Ha))
    assert(isscalar(Ht))
    %n = numel(Ha);  % WRONG!
    n = numel(e);

    %% initial positions
    if isscalar(Ha)
        pos_ant  = [0 Ha]; %Antenna position
    else
        siz = size(Ha);
        Ha = Ha(:);
        zero = zeros(siz);
        pos_ant  = [zero Ha]; %Antenna position
    end
    pos_foot = [0 0];  %Antenna's foot in local frame (origin of local frame)
    pos_geo  = [0 -R0];  % Quasigeocentric origin
    pos_foot_geo = pos_foot - pos_geo; % Antenna's foot in quasigeocentric frame

    %% transmitting satellite's trajectory:
    switch lower(trajectory)
        case 'orbital'  
            % orbital trajectory (constant geocentric radius, variable direct distance)
            Hts = repmat(Ht, [n 1]);
        case 'circular'  
            % circular trajectory around the receiving antenna (constant direct distance, variable geocentric radius)
            %pos_trans = Ht.*[cosd(e), sind(e)];  % WRONG!
            direct_dist = Ht - Ha;
            pos_trans_ant = direct_dist.*[cosd(e(:)), sind(e(:))];
            pos_trans = pos_ant + pos_trans_ant;
            pos_trans_geo = pos_trans + pos_foot_geo;
            Hts = norm_all(pos_trans_geo);
        case 'horizontal'  
            % horizontal trajectory (constant y-axis coordinate, variable geocentric radius, variable direct distance)
            vert_dist = (Ht-Ha).*ones(n,1);
            horiz_dist = vert_dist.*cotd(e(:));
            pos_trans_ant = [horiz_dist vert_dist];
            pos_trans = pos_ant + pos_trans_ant;
            pos_trans_geo = pos_trans + pos_foot_geo;
            Hts = norm_all(pos_trans_geo);
        otherwise
            error('Unknown transmitting satellite trajectory "%s".', char(trajectory))
    end

end

[vitorhjr]

Updates based on your previous suggestions:

New m file:

• get_grazing_angle_infinite.m

Updates:

• get_reflection_spherical_vuorinen.m
• get_spherical_reflection.m

[vitorhjr]

I've made changes in your code to run it properly. However, it should be revised to confirm it is correct.

get_grazing_angle_infinite

Your suggestion
g = 90 - 0.5*rst;

My change:
g = 0.5*rst;

get_spherical_reflection

Your suggestions

delay_reflect_out = slant_dist;
[~, delay_direct_aux] = proj_pt_line (pos_spec, pos_ant, dir_sat)
delay = pos_ant_spec - delay_direct_aux;
delay_direct = inf(size(delay));

My change:

dir_trans = [x_trans, y_trans]; % In finite case, x_trans and y_trans are unit directions
delay_reflect_out = slant_dist;
[~, delay_direct_aux] = proj_pt_line (pos_spec, pos_ant, dir_trans);
% delay = pos_ant_spec - delay_direct_aux;
delay = delay_reflect_out - delay_direct_aux;
delay_direct = inf(size(delay));

[vitorhjr]

For e>0, these suggested changes are working well but for e<=0 they aren't. I continue the discussion about this issue in https://github.com/ufrgs-gnss-lab/atm-interf-sph-dev/issues/61

[vitorhjr]

I updated the function get_grazing_angle_infinite where I corrected the grazing angle computation

Also, I updated get_spherical_reflection_vuorinen where I corrected the geocentric angle between antenna and specular point.