---
title: "HW4"
subtitle: "Software Defined Radio"
output: pdf_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)
```

```{r select-date}
# TODO
selected_date <- as.Date("2000-1-1")
```


\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}
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_1 = (5e-9 + A_I * cos(X_I(tGPS))) * F
  }else{
    I_1 = 5e-9 * F
  }
I_1
}
```

```{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,ys,xlab="Time of Day [Minutes]",ylab="Delay [s]") 
  
```