Doppler models

From NebarnixWiki
Revision as of 08:51, 16 October 2015 by NebarnixWikiSysop (Talk | contribs)

Jump to: navigation, search

Understanding Doppler Shifts

page in work

Doppler shifting is caused by compression of the waves emitted from a source. For audio waves it is the physical objecting 'catching up' with waves emitted at a constant velocity (speed of sound) in front and 'leaving behind' the waves emitted backwards. Wikipedia has several articles to make this more understandable. For electromagnetic waves, it is actually time compression of the source, but it can be modeled as the same effect with the same equations.

We need to model the Doppler shifts of a radio signal in order to match real world data to an actual orbit (or position along a ground track of a known orbit).


Linear Case

Let's start simple. We have a road and there is an observer standing at some distance d from this road. A car is traveling at a constant velocity v. Let us assume that it is emitting engine noise at frequency f.


Circular Case

It is clear that this picture does not work for a satellite. A satellite is not traveling along a linear path, instead, a circle. Let's follow the same idea but close the track from a country road into a racetrack where the cars are traveling in a circle.

We can use the same variables, but now we can reference everyone against the center of the track instead of using the observer as the origin. The same car is traveling at the same constant velocity v on a track with a radius r while our observer is at distance rp from the center of the track, and an angle theta.

%% 2D circular Relative Velocity 
%using kilometers and seconds
clear all;

ObsAlta = 6378;       %sea level
SatAlt = 6378+850;   %850km orbit
OrbitalPeriod = 6120;       %102 minutes
SatAngularVel = 2*pi/OrbitalPeriod; 

t=0:OrbitalPeriod/100:2*OrbitalPeriod; %1000 seconds
InitialAngle = 0;

Vt = -(ObsAlta*SatAlt*SatAngularVel.*sin(InitialAngle+SatAngularVel.*t))./(sqrt(ObsAlta^2+SatAlt^2-2*ObsAlta*SatAlt.*cos(InitialAngle+SatAngularVel.*t)));
subplot(3,1,1);
plot(t,Vt_den);
title(['Distance to Observer: ' num2str(min(Vt_den)) 'km at ' num2str(TcrossingInterp2) 's']);
subplot(3,1,2);
plot(t,Vt);
title('Velocity to Observer');


Polar Spherical Case

The circular case models what happens when a satellite passes DIRECTLY over you, but cannot account for what happens when there is an offset to one side of the pass. We can only model this case by moving forward to spherical coordinates. Things get a little messy here, but we will use the exact same approaches we used before -- the equations are just a little longer because we have to account for the new unit vector Phi. We will use the standard where phi represents rotation along the equator and theta represents the angle between the north pole and the vector while r is again the radius.

This case works quite well, and in fact is almost enough to solve our problem. We can model a 90 degree inclination orbit which is representative of any other circular orbit. We can probably even use this model to solve for geo-locations within about 250km which isn't bad for a 'rough guess'.

What is that 250km all about? Well unfortunately this is the distance that a point will move along the equator during a 15 minute satellite pass (assuming a sat altitude of 850km). The effect will get better at the pole and worse at the equator, but we really can't ignore this.

%% 3D Relative Velocity
%using kilometers and seconds
%ObsPhi = 20*pi/180;
%ObsTheta = pi/2;

%ObsRad = 6378; %sea level
ObsRad = 6368;
ObsTheta =     1.74;
ObsPhi =     0.135;

SatRad = ObsRad + 854; %sea level + altitude
SatPhi = 0; %inclination of 90
SatTheta = pi/2;; %Initial position of 0 (equator)
OrbitalPeriod = 6127.2;       %102 minutes
SatAngularVel = 2*pi/OrbitalPeriod; 

Fscale = 4.9;
Foffset = 2048;

t=-.07*OrbitalPeriod:OrbitalPeriod/300:.07*OrbitalPeriod; %1000 seconds

Vt_num = (ObsRad*SatRad*SatAngularVel*(sin(ObsTheta)*cos(SatPhi-ObsPhi)*cos(SatAngularVel*t+SatTheta)-cos(ObsTheta)*sin(SatAngularVel*t+SatTheta)));
Vt_den = sqrt(ObsRad^2+SatRad^2-2*ObsRad*SatRad*(sin(SatTheta+SatAngularVel*t)*sin(ObsTheta)*cos(SatPhi-ObsPhi)+cos(SatTheta+SatAngularVel*t)*cos(ObsTheta)));
Vt = Vt_num ./ Vt_den;

subplot(3,1,1);
plot(t,Vt_den);
title(['Distance to Observer: ' num2str(min(Vt_den)) 'km at ' num2str(TcrossingInterp2) 's']);

subplot(3,1,2);
plot(t,Vt);
title('Velocity to Observer');

%%Doppler
c = 299792; %in km per second
Ft = 401.65e6; %401.65Mhz uplink freq
F =(2*Vt*Ft)./(c-Vt);
subplot(3,1,3);
plot(t,F,xdata,Fscale*ydata+Foffset,'-o');
title('Doppler Shift of Observer (at 401Mhz)');
max(F)

Polar Spherical Case with Rotating Earth

This is a simple modification. We can just make phi of the observer also time dependent then re-calculate the derivative. The curve changes shape noticeably, which is expected because 250 miles is no small distance.

%% 3D Relative Velocity with Earth's Rotation
%derivative of sqrt(A1^2+B1^2-2*A1*B1*(sin (v*t+A2)*sin (B2)*cos (A3-(v2*t+B3))+cos (v*t+A2)*cos (B2))) with respect to t
%(d)/(dt)(sqrt(A1^2+B1^2-2 A1 B1 (sin(v t+A2) sin(B2) cos(A3-(v2 t+B3))+cos(v t+A2) cos(B2)))) = -(A1 B1 (v2 sin(B2) sin(A2+t v) sin(A3-B3-t v2)+v sin(B2) cos(A2+t v) cos(A3-B3-t v2)-v cos(B2) sin(A2+t v)))/sqrt(A1^2-2 A1 B1 (sin(B2) sin(A2+t v) cos(A3-B3-t v2)+cos(B2) cos(A2+t v))+B1^2)
%A1 is sat radius => SatRad
%v is sat angular velocity = SatAngularVel
%A2 is sat theta start = SatTheta
%A3 is sat phi => SatPhi

%B1 is obs radius => ObsRad
%B2 is obs theta => ObsTheta
%B3 is obs phi start => ObsPhi
%v2 is obs angular velocity => ObsAngularVel

ydata = [1506, 1022, 493, -61, -618, -1154, -1648, -2089, -2472, -2798, -3070, -3296, -3482];
xdata = [42.311,72.211,102.4552,132.2255,162.215,192.2178,222.314,252.2178,282.2255,312.2255,342.2255,372.2255,402.2255];

ObsRad = 6368; %center of European pass (50 degrees)
ObsTheta = 1.74;
ObsPhi = 0.135;
SidrealPeriod = 86164; %23 hours 56 minutes and 4 seconds => seconds
ObsAngularVel = 2*pi/SidrealPeriod;  

SatRad = ObsRad + 854; %sea level + altitude
SatPhi = 0; %inclination of 90 (really its 98, but we don't realy have a term for this
SatTheta = pi/2; %Initial position of 0 (equator)
OrbitalPeriod = 6127.2; %102.12 minutes => seconds
SatAngularVel = 2*pi/OrbitalPeriod; 

Fscale = 4.9;
Foffset = 2048;

t = -0.07*OrbitalPeriod : OrbitalPeriod/1000 : 0.07*OrbitalPeriod; %1000 seconds
Vt_num = (SatRad * ObsRad * (ObsAngularVel*sin(ObsTheta).*sin(SatAngularVel*t+SatTheta).*sin(SatPhi-ObsPhi-ObsAngularVel*t)+SatAngularVel*sin(ObsTheta).*cos(SatTheta+SatAngularVel*t).*cos(SatPhi-ObsPhi-t*ObsAngularVel)-SatAngularVel*cos(ObsTheta).*sin(SatTheta+t*SatAngularVel)));
Vt_den = sqrt(ObsRad^2 + SatRad^2 - 2*SatRad*ObsRad*(sin(ObsTheta).*sin(SatAngularVel*t+SatTheta).*cos(SatPhi-ObsPhi-ObsAngularVel*t)+cos(ObsTheta).*cos(SatTheta+SatAngularVel*t)));
Vt = Vt_num ./ Vt_den;

subplot(3,1,1);
plot(t,Vt_den);

TcrossingInterp2= interp1(Vt,t,0);

title(['Distance to Observer: ' num2str(min(Vt_den)) 'km at ' num2str(TcrossingInterp2) 's']);

subplot(3,1,2);
plot(t,Vt);
title('Velocity to Observer');

%%Doppler
c = 299792; %in km per second
Ft = 401.65e6; %401.65Mhz uplink freq
F =(2*Ft*Vt)./(c-Vt);
subplot(3,1,3);
plot(t,F,xdata,Fscale*ydata+Foffset,'-o');
title('Doppler Shift of Observer (at 401Mhz)');
max(F);

Variable Inclination Spherical Case with Rotating Earth

Whittling away errors, we have now exposed the next problem. Orbits have inclinations, which means that the earth's angular velocity vector and the satellite's motion vector are no longer orthagonal. During ascending nodes and descending nodes the velocity of the earth slightly adds and subtracts from the relative velocity, and thus each node will have a slightly different shape based on this larger or smaller velocity.

Luckily this just means adding a time dependence to the satellite's phi variable, while also remembering that theta is 90 degrees at the equator but inclination is 0 degrees at the equator. We can model this by adding i as inclination and setting phi(t) = sin(i)(omega*t+phi0) and theta(t) = cos(i)(omega*t+theta0). Notice that theta0 and phi0 are now scaled by the inclination. This means that they are now relative to the position along the orbit and NOT absolute in terms of earth coordinates (but could be transformed to absolute by the inclination factor, limiting them to valid possible values)

I believe this THIS model finally has enough fidelity to use for geolocation based on orbital data and doppler shift!