254 lines
6.1 KiB
Plaintext
254 lines
6.1 KiB
Plaintext
---
|
|
title: "HW4"
|
|
subtitle: "Software Defined Radio"
|
|
output: html_document
|
|
date: "`r Sys.Date()`"
|
|
---
|
|
|
|
```{r setup, include=FALSE}
|
|
knitr::opts_chunk$set(echo = TRUE)
|
|
```
|
|
|
|
# Assignment
|
|
Write an R script to calculate the output of the Klobuchar model.
|
|
Provide a GUI means for entering input data, and the name and the
|
|
path to navigation message file with coefficients of the Klobuchar
|
|
model. Load the coefficients of the Klobuchar model from the given
|
|
RINEX n file using a dedicated function in the R script. Provide a
|
|
feature to calculate the value of the ionospheric delay of the GPS
|
|
signal during 24 hours of a given day, in 1-minute increments, with
|
|
the given coefficients of the ionospheric model. Verify the correctness
|
|
of the algorithm implementation with a given numerical example. After
|
|
that, conduct a 24-hour simulation of the GPS ionospheric delay in 1-
|
|
minute increments for the position $\varphi_u$ = 45.33709°, $\lambda_u$ = 14.42496°, $E$
|
|
= 90°, $A$ = 180°. Store the estimated 24-hour simulation data set in a
|
|
dedicated file on the local disc.
|
|
|
|
The problem description, method, R script and implementation results
|
|
should be presented in a single HTML document, which should be
|
|
sent to the subject teacher's e-mail by May 21, 2026 at 23:59 CEST.
|
|
\newpage
|
|
|
|
# Input
|
|
|
|
### Select the file
|
|
```{r select-file}
|
|
library(rstudioapi)
|
|
rinex_file <- rstudioapi::selectFile()
|
|
```
|
|
|
|
### Read the file
|
|
```{r read-file}
|
|
extract_klobuchar_parameters <- function(rinex_file) {
|
|
ion_alpha <- numeric(4)
|
|
ion_beta <- numeric(4)
|
|
|
|
# Read the file line by line
|
|
con <- file(rinex_file, "r")
|
|
while (length(line <- readLines(con, n = 1)) > 0) {
|
|
if (grepl("ION ALPHA", line, fixed = TRUE)) {
|
|
# Extract the 4 numeric values from the line
|
|
ion_alpha <- as.numeric(gsub("D", "e", strsplit(line, "\\s+")[[1]][2:5]))
|
|
}
|
|
if (grepl("ION BETA", line, fixed = TRUE)) {
|
|
# Extract the 4 numeric values from the line
|
|
ion_beta <- as.numeric(gsub("D", "e", strsplit(line, "\\s+")[[1]][2:5]))
|
|
}
|
|
# Stop reading if we've reached the end of the header (END OF HEADER)
|
|
if (grepl("END OF HEADER", line, fixed = TRUE)) {
|
|
break
|
|
}
|
|
}
|
|
close(con)
|
|
|
|
# Return the parameters as a named list
|
|
list(
|
|
ION_ALPHA = ion_alpha,
|
|
ION_BETA = ion_beta
|
|
)
|
|
}
|
|
|
|
klobuchar_params <- extract_klobuchar_parameters(rinex_file)
|
|
alpha <- function(index){
|
|
klobuchar_params$ION_ALPHA[index+1]
|
|
}
|
|
beta <- function(index){
|
|
klobuchar_params$ION_BETA[index+1]
|
|
}
|
|
|
|
print(klobuchar_params)
|
|
```
|
|
|
|
|
|
\newpage
|
|
|
|
# Klobuchar calculation
|
|
- $R_E$ is the earth's mean radius, $6371$ km.
|
|
- $h$ is the height of the ionosphere, which is $350$ km for GPS ($375$ for Beidou).
|
|
```{r klob-inputs}
|
|
# constants
|
|
R_E = 6371e3
|
|
h = 350e3
|
|
# inputs
|
|
lat = 45.33709 # phi
|
|
lon = 14.42496 # lambda
|
|
E = 90 # elevation
|
|
A = 180 # azimuth
|
|
```
|
|
|
|
### 1. Calculate earth-centered angle
|
|
$$ \psi = \pi/2 - E - \arcsin(\frac{R_E}{R_E + h} \cos(E)) $$
|
|
```{r klob-1}
|
|
psi = pi/2 - E - asin(R_E / (R_E+h) * cos(E))
|
|
```
|
|
|
|
```{r klob-1-demo, echo=FALSE}
|
|
psi
|
|
```
|
|
|
|
### 2. Compute the latitude of the IPP (ionospheric pierce point)
|
|
$$ \phi_I = \arcsin( sin \varphi_u \cos \psi + \cos\varphi_u \sin\psi \cos A )$$
|
|
```{r klob-2}
|
|
phi_I = asin(sin(lat) * cos(psi) + cos(lat) * sin(psi)* cos(A))
|
|
```
|
|
|
|
```{r klob-2-demo, echo=FALSE}
|
|
phi_I
|
|
```
|
|
|
|
### 3. Compute the longitude of the IPP
|
|
$$ \lambda_I= \lambda_u + \frac{\psi \sin A}{\cos \phi_I}$$
|
|
```{r klob-3}
|
|
lambda_I = lat + (psi * sin(A))/cos(phi_I)
|
|
```
|
|
```{r klob-3-demo, echo=FALSE}
|
|
lambda_I
|
|
```
|
|
|
|
### 4. Find the geomagnetic latitude of the IPP
|
|
$$ \phi_m = \arcsin(\sin\phi_I \sin\phi_P + \cos\phi_I \cos\phi_P \cos(\lambda_I - \lambda_P)) $$
|
|
with $\phi_P = 78.3$°, $\lambda_P = 291.0$° the coordinates of the geomagnetic pole.
|
|
```{r klob-4}
|
|
phi_P = 78.3
|
|
lambda_P = 291
|
|
phi_m = asin(sin(phi_I) * sin(phi_P) + cos(phi_I) * cos(phi_P) * cos(lambda_I - lambda_P))
|
|
```
|
|
|
|
```{r klob-4-demo, echo=FALSE}
|
|
phi_m
|
|
```
|
|
### 5. Find the local time at the IPP
|
|
$$ t = 43200 \lambda_I / \pi + t_{GPS} $$
|
|
With:
|
|
$\lambda_I$ in radians, $t$ in seconds.
|
|
where $0 \leq t \leq 86 400$. Therefore:
|
|
If $t \geq 86 400$, subtract $86 400$. If $t < 0$, add $86 400$.
|
|
```{r klob-5}
|
|
t <- function(tGPS){
|
|
43200 * lambda_I / pi + tGPS
|
|
}
|
|
```
|
|
```{r klob-5-demo}
|
|
t(1)
|
|
```
|
|
|
|
### 6. Compute the amplitude of ionospheric delay
|
|
$$ A_I = \sum_{n=0}^{3} \alpha_n(\phi_m/\pi)^n $$
|
|
If $A_I < 0$, then $A_I =0$.
|
|
```{r klob-6}
|
|
n = c(0,1,2,3)
|
|
A_I = sum(alpha(n) * (phi_m/pi) ^ n)
|
|
if (A_I < 0) {
|
|
A_I = 0
|
|
}
|
|
```
|
|
```{r klob-6-demo, echo=FALSE}
|
|
A_I
|
|
```
|
|
|
|
### 7. Compute the period of ionospheric delay
|
|
$$ P_I = \sum_{n=0}^{3} \beta_n(\phi_m/\pi)^n $$
|
|
If $P_I < 72 000$ then $P_I = 72000$.
|
|
|
|
```{r klob-7}
|
|
P_I = sum(beta(n) * (phi_m/pi) ^ n)
|
|
if (P_I < 72000) {
|
|
P_I = 72000
|
|
}
|
|
```
|
|
```{r klob-7-demo, echo=FALSE}
|
|
P_I
|
|
```
|
|
|
|
### 8. Compute the phase of ionospheric delay
|
|
$$ X_I = \frac{2\pi ( t-50400)}{P_I}$$
|
|
|
|
```{r klob-8}
|
|
X_I <- function(tGPS){
|
|
(2 * pi * (t(tGPS)-50400))/P_I
|
|
}
|
|
```
|
|
```{r klob-8-demo, echo=FALSE}
|
|
X_I(1)
|
|
```
|
|
|
|
### 9. Compute the slant factor (ionospheric mapping function)
|
|
$$ F=\left[1 - (\frac{R_E}{R_E + h} \cos E)^2\right]^{-1/2} $$
|
|
|
|
```{r klob-9}
|
|
# with the given input this becomes constant 1 due to FP inaccuracies
|
|
F = (1- (R_E/(R_E+h) * cos(E) )^2)^(-1/2)
|
|
```
|
|
|
|
```{r klob-9-demo, echo=FALSE}
|
|
F
|
|
```
|
|
|
|
### 10. Compute the ionospheric time delay
|
|
$$
|
|
I_1 = \begin{cases}
|
|
[5 \cdot 10^{-9} + A_I \cos X_I] \times F, & \quad |X_I| < \pi/2 \\
|
|
5 \cdot 10^{-9} \times F, & \quad |X_I| \geq \pi/2
|
|
\end{cases}
|
|
$$
|
|
|
|
```{r klob-10}
|
|
I_1 <- function(tGPS){
|
|
if (abs(X_I(tGPS))< pi/2) {
|
|
I = (5e-9 + A_I * cos(X_I(tGPS))) * F
|
|
}else{
|
|
I = 5e-9 * F
|
|
}
|
|
I
|
|
}
|
|
```
|
|
|
|
```{r klob-10-demo, echo=FALSE}
|
|
I_1(1)
|
|
```
|
|
The delay $I_1$ is given in seconds and is referred to the GPS L1 or
|
|
Beidou B1 frequencies.
|
|
|
|
# 24-hour simulation
|
|
```{r simulate}
|
|
xs <- seq(from = 0, to = 24*60, by = 1)
|
|
ys <- numeric(length(xs))
|
|
for (x in xs){
|
|
y <- I_1(x*60)
|
|
ys[x] <- y
|
|
}
|
|
```
|
|
|
|
```{r simulate-demo, echo=FALSE}
|
|
head(ys)
|
|
plot(xs[-length(xs)],ys[-length(xs)],xlab="Time of Day [Minutes]",ylab="Delay [s]")
|
|
|
|
```
|
|
|
|
# Output
|
|
|
|
```{r output}
|
|
write.csv(ys, './klobuchar.csv')
|
|
```
|
|
# Results
|
|
Due to the calculation in Section xxx the output is a constant $5 \cdot 10^{-9} \cdot 1 $. |