# Estimate GNSS Receiver Position with Simulated Satellite Constellations

Track the position of a ground vehicle using a simulated Global Navigation Satellite System (GNSS) receiver. The satellites are simulated using the `satelliteScenario` object, the satellite signal processing of the receiver are simulated using the `lookangles` (Navigation Toolbox) and `pseudoranges` (Navigation Toolbox) functions, and the receiver position is estimated with the `receiverposition` (Navigation Toolbox) function.

### Overview

This example focuses on the space segment, or satellite constellations, and the GNSS sensor equipment for a statellite system. To obtain an estimate of the GNSS receiver position, the navigation processor requires the satellite positions from the space segment and the pseudoranges from the ranging processor in the receiver.

### Specify Simulation Parameters

Load the MAT-file that contains the ground-truth position and velocity of a ground vehicle travelling toward the Natick, MA campus of Mathworks, Inc.

Specify the start, stop, and sample time of the simulation. Also, specify the mask angle, or minimum elevation angle, of the GNSS receiver.

```% Load ground truth trajectory. load("routeNatickMA.mat","lat","lon","pos","vel","lla0"); recPos = pos; recVel = vel; % Specify simulation times. startTime = datetime(2020,10,28,8,0,0,"TimeZone","America/New_York"); simulationSteps = size(pos,1); dt = 1; stopTime = startTime + seconds((simulationSteps-1)*dt); % Specify mask angle. maskAngle = 5; % degrees % Convert receiver position from North-East-Down (NED) to geodetic % coordinates. Requires Navitation Toolbox(TM). receiverLLA = ned2lla(recPos,lla0,"ellipsoid"); % Set RNG seed to allow for repeatable results. rng("default");```

Visualize the `geoplot` for the ground truth trajectory.

```figure geoplot(lat,lon) geobasemap("topographic") title("Ground Truth Trajectory")``` ### Simulate Satellite Positions Over Time

The `satelliteScenario` object enables you to specify initial orbital parameters and visualize them using the `satelliteScenarioViewer` object. This example uses the `satelliteScenario` and the a MAT-file with initial orbital parameters to simulate the GPS constellations over time. Alternatively, you could use the `gnssconstellation` (Navigation Toolbox) object which simulates satellite positions using nominal orbital parameters, and only the current simulation time is needed to calculate the satellite positions.

```% Create scenario. sc = satelliteScenario(startTime, stopTime, dt); load("initialOrbitalParameters.mat","semiMajorAxis","eccentricity",... "inclination","rightAscensionOfAscendingNode",... "argumentOfPeriapsis","trueAnomaly","prnStr"); % Initialize satellites. satellite(sc,semiMajorAxis,eccentricity,inclination,... rightAscensionOfAscendingNode,argumentOfPeriapsis,trueAnomaly,... "Name",prnStr); % Preallocate results. numSats = numel(sc.Satellites); allSatPos = zeros(numSats,3,simulationSteps); allSatVel = zeros(numSats,3,simulationSteps); % Save satellite states over entire simulation. for i = 1:numel(sc.Satellites) [oneSatPos, oneSatVel] = states(sc.Satellites(i),"CoordinateFrame","ecef"); allSatPos(i,:,:) = permute(oneSatPos,[3 1 2]); allSatVel(i,:,:) = permute(oneSatVel,[3 1 2]); end```

### Calculate Pseudoranges

Use the satellite positions to calculate the pseudoranges and satellite visibilities throughout the simulation. The mask angle is used to determine the satellites that are visibile to the receiver. The pseudoranges are the distances between the satellites and the GNSS receiver. The term pseudorange is used because this distance value is calculated by multiplying the time difference between the current receiver clock time and the timestamped satellite signal by the speed of light.

```% Preallocate results. allP = zeros(numSats,simulationSteps); allPDot = zeros(numSats,simulationSteps); allIsSatVisible = false(numSats,simulationSteps); % Use the skyplot to visualize satellites in view. sp = skyplot([], []); for idx = 1:simulationSteps satPos = allSatPos(:,:,idx); satVel = allSatVel(:,:,idx); % Calculate satellite visibilities from receiver position. [satAz,satEl,allIsSatVisible(:,idx)] = lookangles(receiverLLA(idx,:),satPos,maskAngle); % Calculate pseudoranges and pseudorange rates using satellite and % receiver positions and velocities. [allP(:,idx),allPDot(:,idx)] = pseudoranges(receiverLLA(idx,:),satPos,recVel(idx,:),satVel); set(sp,"AzimuthData",satAz(allIsSatVisible(:,idx)), ... "ElevationData",satEl(allIsSatVisible(:,idx)), ... "LabelData",prnStr(allIsSatVisible(:,idx))) drawnow limitrate end``` ### Estimate Receiver Position from Pseudoranges and Satellite Positions

Finally, use the satellite positions and pseudoranges to estimate the receiver position with the `receiverposition` function.

```% Preallocate results. lla = zeros(simulationSteps,3); gnssVel = zeros(simulationSteps,3); hdop = zeros(simulationSteps,1); vdop = zeros(simulationSteps,1); for idx = 1:simulationSteps p = allP(:,idx); pdot = allPDot(:,idx); isSatVisible = allIsSatVisible(:,idx); satPos = allSatPos(:,:,idx); satVel = allSatVel(:,:,idx); % Estimate receiver position and velocity using pseudoranges, % pseudorange rates, and satellite positions and velocities. [lla(idx,:),gnssVel(idx,:),hdop(idx,:),vdop(idx,:)] = receiverposition(p(isSatVisible),... satPos(isSatVisible,:),pdot(isSatVisible),satVel(isSatVisible,:)); end```

### Visualize Results

Plot the ground truth position and the estimated receiver position on a `geoplot`.

```figure geoplot(lat,lon,lla(:,1),lla(:,2)) geobasemap("topographic") legend("Ground Truth","Estimate")``` Plot the absolute error in the position estimate. The error is smoothed by a moving median to make the plot more readable. The error in the x- and y-axis is smaller because there are satellites on either side of the receiver. The error is the z-axis is larger because there are only satellites above the receiver, not below it. The error changes over time as the receiver moves and some satellites come in and out of view.

```estPos = lla2ned(lla,lla0,"ellipsoid"); winSize = floor(size(estPos,1)/10); figure plot(smoothdata(abs(estPos-pos),"movmedian",winSize)) legend("x","y","z") xlabel("Time (s)") ylabel("Error (m)") title("Position (NED) Error")``` ## Related Examples 