hw4 calculcations

hw4.Rmd

@@ -27,14 +27,17 @@ 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)
@@ -66,6 +69,12 @@ extract_klobuchar_parameters <- function(rinex_file) {
 }
 
 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)
 ```
@@ -75,54 +84,145 @@ print(klobuchar_params)
 selected_date <- as.Date("2000-1-1")
 ```
 
+
+\newpage
+
 # Klobuchar calculation
-1. Calculate earth-centered angle
-What are R_E and h?
+- $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)
+### 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 )$$
-3. Compute the longitude of the IPP
+```{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}$$
-4. Find the geomagnetic latitude of the IPP
-What is phi_P?
+```{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)) $$
-5. Find the local time at the IPP
+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}
+tGPS = 1 # TODO
+t = 43200 * lambda_I / pi + tGPS 
+```
+```{r klob-5-demo}
+t
+```
 
-6. Compute the amplitude of ionospheric delay
+### 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
+### 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$.
 
-8. Compute the phase of ionospheric delay
+```{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}$$
 
-9. Compute the slant factor (ionospheric mapping function)
+```{r klob-8}
+X_I = (2 * pi * (t-50400))/P_I
+```
+```{r klob-8-demo, echo=FALSE}
+X_I
+```
+
+### 9. Compute the slant factor (ionospheric mapping function)
 $$ F=\left[1 - (\frac{R_E}{R_E + h} \cos E)^2\right]^{-1/2} $$
-10. Compute the ionospheric time delay
+
+```{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}
 $$
-The delay I1 is given in seconds and is referred to the GPS L1 or
+
+```{r klob-10}
+if (abs(X_I)< pi/2) {
+  I_1 = (5e-9 + A_I * cos(X_I)) * F
+}else{
+  I_1 = 5e-9 * F
+}
+```
+
+```{r klob-10-demo, echo=FALSE}
+I_1
+```
+The delay $I_1$ is given in seconds and is referred to the GPS L1 or
 Beidou B1 frequencies.
-
-
-```{r simulate}
-# input constants
-lat = 45.33709 # phi
-lon = 14.42496 # lambda
-E = 90
-A = 180
-
-```