A satellite navigation system consists of three main parts:
1. Space segment (a constellation of satellites orbiting Earth)
2. Control segment (a network of monitoring stations that track satellites and send commands to them)
3. User segment (radio equipment that recieves and processes signals from satellites)
The satellites broadcast messages that contain information on their coordinates. User GNSS receivers register time of signal arrival and measure time of signal transmission. This allows computing the so-called pseudo-ranges, i. e. distances from user antenna to visible satellites. Knowing distances to at least four of the satellites makes possible to determine position of the user. The fourth preudo-range is required because clocks of the receiver and the satellites are not synchronized and a shift between them needs to be estimated.

### GNSS observables

Typically, GNSS user equipment outputs three types of measurements: pseudo-ranges, pseudo-range rates, and carrier phases.
Pseudo-ranges are determined by receivers as
\begin{equation*} \tilde \rho_R = c \cdot (\tilde t_{sa} - \tilde t_{st}), \end{equation*}
where $$\tilde t_{sa}$$ is the time of signal arrival, measured by receiver clock, $$c$$ is the speed of light, and $$\tilde t_{st}$$ is the time of signal transmission. $$\tilde t_{st}$$ is computed using the code phase $$\tilde t^\prime_{st}$$ measured by the code tracking loop and the number of code repetitions determined from the navigation message.
Pseudo-range rates are computed as
\begin{equation*} \tilde{\dot{\rho}}_R = -\frac{c}{f_{ca}} \Delta \tilde{f}_{ca}. \end{equation*}
Herein $$f_{ca}$$ denotes the tramsmitted carrier frequency. The Doppler-shift $$\Delta \tilde{f}_{ca}$$ is obtained directly from the carrier tracking loop.
Carrier phase is obtained by sampling the phase of the reference oscillator of the carrier tracking loop. The measurement of the phase $$\phi(t_{sa})$$ is related to the pseudo-range as
\begin{equation*} \tilde{\rho}_{R} (t_{sa}) = \lambda \cdot \left( \frac{\tilde{\phi}(t_{sa})}{2\pi} + N \right), \end{equation*}
where $$N$$ is an unknown integer number of wavelengths contained in the pseudo-range. If the ambiguity is resolved, the carrier phase measurement allows to derive a pseudo-range measurement that is 10 to 100 times more accurate, than the one provided by the code tracking loop. Determination of the integer ambiguity $$N$$ requires additional measurements and can be achieved by applying LAMBDA method, for example.

### Calculation of position, velocity, and time

#### Measurement models

The range between the i-th satellite and the antenna measured by the navigation receiver is usually modeled as
\begin{eqnarray} \label{eq:prMeasEq} \tilde{\rho}^{i}_{R} &=& \left\Vert \mathbf{x}^{sat,i}_{i} - \mathbf{x}^{A}_{i} \right\Vert + c \cdot \left( \delta t_{r} - \delta t^{i} \right) + I^{i}_{R} + T^{i}_{R} + M^{i}_{R} + \nu^{i}_{R} = \\ \nonumber &=& \left\Vert \mathbf{x}^{sat,i}_{e} - \mathbf{x}^{A}_{e} \right\Vert + \delta \rho^{i}_{ie,R} + c \cdot \left( \delta t_{r} - \delta t^{i} \right) + I^{i}_{R} + T^{i}_{R} + M^{i}_{R} + \nu^{i}_{R} \end{eqnarray}
Herein
• $$\mathbf{x}^{sat,i}_{i}$$ is a vector of coordinates of the i-th satellite in inertial frame
• $$\mathbf{x}^{A}_{i}$$ is a vector of coordinates of the antenna in inertial frame
• $$\delta t_{r}$$ is the receiver clock error
• $$\delta t^{i}$$ is the clock error of the i-th satellite
• $$I^{i}_{R}$$ is the range error due to ionospheric delay
• $$T^{i}_{R}$$ is the range error due to tropospheric delay
• $$M^{i}_{R}$$ is the multipath error
• $$\nu^{i}_{R}$$ is the random measurement noise
• $$\mathbf{x}^{sat,i}_{e}$$ is a vector of coordinates of the i-th satellite in the ECEF coordinate frame
• $$\mathbf{x}^{A}_{e}$$ is a vector of coordinates of the antenna in the ECEF coordinate frame
• $$\delta \rho^{i}_{ie,R} \approx \frac{\omega_{ie}}{c} \left( x^{sat,i}_{e} y^{A}_{e} - y^{sat,i}_{e} x^{A}_{e} \right)$$ is the range Sagnac correction [1]
The measured pseudo-range is a function of the true satellite's coordinates, which are never known. Intoducing $$E^{i}_{R} = \left\Vert \mathbf{x}^{sat,i}_{e} - \mathbf{x}^{A}_{e} \right\Vert - \left\Vert \tilde{\mathbf{x}}^{sat,i}_{e} - \mathbf{x}^{A}_{e} \right\Vert$$  allows to rewrite the expression (\ref{eq:prMeasEq}) in terms of computed satellite coordinates as follows:
$$\label{eq:prMeasEq2} \tilde{\rho}^{i}_{R} = \left\Vert \tilde{\mathbf{x}}^{sat,i}_{e} - \mathbf{x}^{A}_{e} \right\Vert + \delta \rho^{i}_{ie,R} + E^{i}_{R} + c \cdot \left( \delta t_{r} - \delta t^{i} \right) + I^{i}_{R} + T^{i}_{R} + M^{i}_{R} + \nu^{i}_{R}$$
The term $$E^{i}_{R}$$ represents the error caused by inaccuracies in the broadcast ephemerides. The multipath error $$M^{i}_{R}$$ may be mitigated by antenna construction and installation. The ionospheric delay $$I^{i}_{R}$$ can be compenstated by dual-frequency GNSS receivers.
Measured pseudo-range rates can be modeled as
\begin{eqnarray} \label{eq:doMeasEq} \tilde{\dot{\rho}}_{R}^{i} &=&  ( \mathbf{e}^{i}_{i} )^{\mathrm{T}} \cdot \left( \mathbf{v}_{i}^{sat,i} - \mathbf{v}^{A}_{i} \right) + c \cdot \left( \delta \dot{t}_{r} - \delta \dot{t}^{i} \right) + \dot{I}^{i}_{r} + \dot{T}^{i}_{r} + \nu^{i}_{\dot{R}} = \\ \nonumber &=& ( \mathbf{e}^{i}_{e} )^{\mathrm{T}} \cdot \left( \mathbf{v}_{e}^{sat,i} - \mathbf{v}^{A}_{e} \right) + \delta \dot{\rho}^{i}_{ie,R} + c \cdot \left( \delta \dot{t}_{r} - \delta \dot{t}^{i} \right) + \dot{I}^{i}_{R} + \dot{T}^{i}_{R} + \nu^{i}_{\dot{R}} \end{eqnarray}
Herein
• $$\mathbf{v}^{sat,i}_{i}$$ is the velocity of the i-th satellite relative to inertial frame
• $$\mathbf{v}^{A}_{i}$$ is the velocity of the antenna relative to inertial frame
• $$\delta \dot{t}_{r}$$ is the error of the receiver clock drift
•  $$\delta \dot{t}^{i}$$ is the drift of the satellite clock
• $$\dot{I}^{i}_{R}$$ and $$\dot{T}^{i}_{R}$$ are the errors caused by the ionosphere and the troposhpere
• $$\nu^{i}_{\dot{R}}$$ is the random measurement error
• $$\mathbf{v}^{sat,i}_{e}$$ is the velocity of the i-th satellite relative to the ECEF coordinate frame
• $$\mathbf{v}^{A}_{e}$$ is the velocity of the antenna relative to the ECEF coordinate frame
• $$\delta \dot{\rho}^{i}_{ie,R} \approx \frac{\omega_{ie}}{c} \left( v^{sat,i}_{ey} x^{A}_{e} + y^{sat,i}_{e} v^{A}_{ex} - v^{sat,i}_{ex} y^{A}_{e} - x^{sat,i}_{e} v^{A}_{ey} \right)$$ is the range-rate Sagnac correction [1]
$$\mathbf{e}_{e}$$ is the line of sight vector that is determined by
$$\label{eq:lineOfSightDef} \mathbf{e}_{e} = \frac{\mathbf{x}^{sat,i}_{e} - \mathbf{x}^{A}_{e}}{\left\Vert \mathbf{x}^{sat,i}_{e} - \mathbf{x}^{A}_{e} \right\Vert}$$
Similar to pseudo-range, pseudo-range rate is a function of the true satellite's position and velocity, which are not known in practice. The equation (\ref{eq:doMeasEq}) can be rewritten in terms of the computed satellite's position and velocity as
$$\label{eq:doMeasEq2} \tilde{\dot{\rho}}_{R}^{i} = ( \tilde{\mathbf{e}}^{i}_{e} )^{\mathrm{T}} \cdot \left( \tilde{\mathbf{v}}_{e}^{sat,i} - \mathbf{v}^{A}_{e} \right) + c \cdot \left( \delta \dot{t}_{r} - \delta \dot{t}^{i} \right) + \delta \dot{\rho}^{i}_{ie,R} + E^{i}_{\dot{R}} + \dot{I}^{i}_{R} + \dot{T}^{i}_{R} + \nu^{i}_{\dot{R}}$$
where $$E^{i}_{\dot{R}} = ( \mathbf{e}^{i}_{e} )^{\mathrm{T}} \cdot \left( \mathbf{v}_{e}^{sat,i} - \mathbf{v}^{A}_{e} \right) - ( \tilde{\mathbf{e}}^{i}_{e} )^{\mathrm{T}} \cdot \left( \tilde{\mathbf{v}}_{e}^{sat,i} - \mathbf{v}^{A}_{e} \right)$$ is the error caused by inaccurate broadcast ephemerides.
The terms  $$E^{i}_{\dot{R}}$$, $$\dot{I}^{i}_{R}$$, and $$\dot{T}^{i}_{R}$$ are negligible, while the residual satellite clock drift $$\delta \dot{t}^{i}$$ is also small when a correction derived from the navigation message is applied. Thus the following simplified model is usually used in practice for velocity estimation
$$\label{eq:doMeasEq3} \tilde{\dot{\rho}}_{R}^{i} = ( \tilde{\mathbf{e}}^{i}_{e} )^{\mathrm{T}} \cdot \left( \tilde{\mathbf{v}}_{e}^{sat,i} - \mathbf{v}^{A}_{e} \right) + \delta \dot{\rho}^{i}_{ie,R} + c \cdot \delta \dot{t}_{r} + \nu^{i}_{\dot{R}}$$
Carrier phase measurement can be expressed similarly to pseudo-range as
$$\label{eq:crMeasEq} \tilde{\rho}_{\phi}^{i} = \lambda \cdot \left( \frac{\tilde{\phi}}{2\pi} + N \right) = \left\Vert \tilde{\mathbf{x}}^{sat,i}_{e} - \mathbf{x}^{A}_{e} \right\Vert + \delta \rho^{i}_{ie,R} + \lambda \cdot N + E^{i}_{R} + c \cdot \left( \delta t_{r} - \delta t^{i} \right) + I^{i}_{\phi} + T^{i}_{\phi} + M^{i}_{\phi} + \nu^{i}_{\phi}$$
Herein  $$I^{i}_{\phi}$$,  $$T^{i}_{\phi}$$, and $$M^{i}_{\phi}$$ are the errors caused by ionosphere, troposhere, and multipath, respectively. $$\nu^{i}_{\phi}$$ is the random measurement noise, $$\lambda$$ is the wavelength of the carrier, and $$N$$ is an unknown integer number, commonly called an ambiguity.

#### Single point solution

Assume that some initial approximation $$\tilde{\mathbf{x}}^{A,-}_{e}$$ of the antenna coordinates is provided. Then the following residual can be calculated for each measured pseudo-range
$$\label{eq:prResEq} \Delta \rho^{i}_{R} = \tilde{\rho}^{i}_{R} - \left\Vert \tilde{\mathbf{x}}^{sat,i}_{e} - \tilde{\mathbf{x}}^{A,-}_{e} \right\Vert - \delta \rho^{i}_{ie,R}$$
Introduce antenna position error as $$\delta \mathbf{x}^{A}_{e} = \mathbf{x}^{A}_{e} - \tilde{\mathbf{x}}^{A,-}_{e}$$, and let $$\delta \rho^{i}_{R}$$ denote the residual error of the pseudo-range measurement. Then in linear approximation
$$\Delta \rho^{i}_{R} = - ( \tilde{\mathbf{e}}^{i}_{e} )^{\mathrm{T}} \cdot \delta \mathbf{x}^{A}_{e} + c \cdot \delta t_{r} + \delta \rho^{i}_{R}$$
When $$N_{sat} \ge 4$$ satellites are visible, the linear system
\begin{eqnarray} \label{eq:prSys} \Delta \rho^{1}_{R} &=& - ( \tilde{\mathbf{e}}^{1}_{e} )^{\mathrm{T}} \cdot \delta \mathbf{x}^{A}_{e} + c \cdot \delta t_{r} + \delta \rho^{1}_{R} \\ \nonumber &\vdots& \\ \nonumber \Delta \rho^{N_{sat}}_{R} &=& - ( \tilde{\mathbf{e}}^{N_{sat}}_{e} )^{\mathrm{T}} \cdot \delta \mathbf{x}^{A}_{e} + c \cdot \delta t_{r} + \delta \rho^{N_{sat}}_{R} \end{eqnarray}
appears overdetermined and the quantities $$\delta \mathbf{x}^{A}_{e}$$ and $$\delta t_{r}$$ can be estimated  applying the least-squares method:
$$\label{eq:singlePosPointSol} \begin{pmatrix} \delta \hat{\mathbf{x}}^{A}_{e} \\ \delta \hat{t}_{r} \end{pmatrix} = \left( \mathbf{H}^{\mathrm{T}} \mathbf{W}^{-1}_{R} \mathbf{H} \right)^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{W}^{-1}_{R} \Delta \mathbf{z}$$
Herein $$\mathbf{H} = \begin{pmatrix} - ( \tilde{\mathbf{e}}^{1}_{e} )^{\mathrm{T}} & c \\ \vdots & \vdots \\ - ( \tilde{\mathbf{e}}^{N_{sat}}_{e} )^{\mathrm{T}} & c \end{pmatrix}, \; \Delta \mathbf{z} = \begin{pmatrix} \Delta \rho^{1}_{R} \\ \vdots \\ \Delta \rho^{N_{sat}}_{R} \end{pmatrix},$$ and $$\mathbf{W}_{R}$$ is a covariance matrix of the residual measurement error. $$\mathbf{W}_{R}$$ is often assumed to be a diagonal matrix, whose elements are inversely proportional to the elevation angle of the corresponding satellite.
When position error estimation is done, the initial approximation of the antenna coordinates could be corrected as follows
$$\tilde{\mathbf{x}}^{A,+}_{e} = \tilde{\mathbf{x}}^{A,-}_{e} + \delta \hat{\mathbf{x}}^{A}_{e}$$
Sometimes the initial approximation may be very rough, and error estimation needs to be repeated iteratively several times until reasonable accuracy is achieved.
Velocity can first be determined when position is known. Unlike in case of position, velocity is linearly related to pseudo-range rate measurements, and thus it can be estimated directly.
Consider the vectors
$$\label{eq:rrMeasSinglePointSol} \mathbf{z}_{\dot{R}} = \begin{pmatrix} \tilde{\dot{\rho}}_{R}^{1} - ( \tilde{\mathbf{e}}^{1}_{e} )^{\mathrm{T}} \cdot \tilde{\mathbf{v}}_{e}^{sat,1} - \delta \dot{\rho}^{1}_{ie,R} \\ \vdots \\ \tilde{\dot{\rho}}_{R}^{N_{sat}} - ( \tilde{\mathbf{e}}^{N_{sat}}_{e} )^{\mathrm{T}} \cdot \tilde{\mathbf{v}}_{e}^{sat,N_{sat}} - \delta \dot{\rho}^{N_{sat}}_{ie,R} \end{pmatrix}, \; \delta \dot{\pmb{\rho}}_{R} = \begin{pmatrix} \delta \dot{\rho}^{1}_{R} \\ \vdots \\ \delta \dot{\rho}^{N_{sat}}_{R} \end{pmatrix},$$
where $$\delta \dot{\rho}^{i}_{R}$$ denotes a residual error of the i-th pseudo-range rate measurement.
Then
$$\mathbf{z}_{\dot{R}} = \mathbf{H} \cdot \begin{pmatrix} \mathbf{v}^{A}_{e} \\ \delta \dot{t}_{r} \end{pmatrix} + \delta \dot{\pmb{\rho}}_{R},$$
and the estimates of the $$\mathbf{v}^{A}_{e}$$ and $$\delta \dot{t}_{r}$$ can be obtained using
$$\label{eq:singleVelPointSol} \begin{pmatrix} \hat{\mathbf{v}}^{A}_{e} \\ \delta \hat{\dot{t}}_{r} \end{pmatrix} = \left( \mathbf{H}^{\mathrm{T}} \mathbf{W}^{-1}_{\dot{R}} \mathbf{H} \right)^{-1} \mathbf{H}^{\mathrm{T}} \mathbf{W}^{-1}_{\dot{R}} \mathbf{z}_{\dot{R}}$$
Herein $$\mathbf{W}_{\dot{R}}$$ denotes a covariance matrix of the residual pseudo-range rate measurement error.
Note that the estimates of the receiver's clock bias and drift could be used to correct receiver's local time and provide the user with highly accurate global time reference.

#### Filtered solution

In practice GNSS receivers rarely exploit the described single point solution for position determination. Instead of it coordinates are estimated by means of Kalman filter. This approach has several advantages. First of all, the estimates of position and velocity become smoother. Secondly, exploiting Kalman filter allows the receiver to continue position calculation when less than four satellites are visible. Despite these advantages, the single point approach is still important, because it is used for getting the first position fix and initializing Kalman filter based algorithm. Moreover, single point solution is widely exploited in receiver autonomous integrity monitoring. A reader is encouraged to learn more about filtered solution from [1].