diff --git a/CertificateOfParticipation Poster.pdf b/CertificateOfParticipation Poster.pdf
new file mode 100644
index 0000000..d16c054
Binary files /dev/null and b/CertificateOfParticipation Poster.pdf differ
diff --git a/Flying_on_the_sphere0.pdf b/Flying_on_the_sphere0.pdf
new file mode 100644
index 0000000..6906a27
Binary files /dev/null and b/Flying_on_the_sphere0.pdf differ
diff --git a/GatherProcessData/Historical/GatherOpenSky.py b/GatherProcessData/Historical/GatherOpenSky.py
index 810096b..b050c3d 100644
--- a/GatherProcessData/Historical/GatherOpenSky.py
+++ b/GatherProcessData/Historical/GatherOpenSky.py
@@ -1,4 +1,4 @@
 import OpenSkyManager as mng
 
 
-osky = mng.OpenSkyManager('Lars_Johansson57', '0pen_Sky_N@|se' )
\ No newline at end of file
+osky = mng.OpenSkyManager('Lars_Johansson57', 'se' )
\ No newline at end of file
diff --git a/GatherProcessData/Historical/Open_Sky_Testing.py b/GatherProcessData/Historical/Open_Sky_Testing.py
index 08b932c..b7d0eae 100644
--- a/GatherProcessData/Historical/Open_Sky_Testing.py
+++ b/GatherProcessData/Historical/Open_Sky_Testing.py
@@ -30,7 +30,7 @@ def main():
     while continues == False:
         try:
             opensky.username = "Lars_Johansson57"
-            opensky.password = "0pen_Sky_N@|se"
+            opensky.password = "0se"
 
 
             #If the following works, well be able to get flight data for every 0.5 seconds
diff --git a/Planes/Planes1.R b/Planes/Planes1.R
new file mode 100644
index 0000000..ab95103
--- /dev/null
+++ b/Planes/Planes1.R
@@ -0,0 +1,369 @@
+#Writen by: Lars Daniel Johansson Niño
+#Created date: 08/08/2024
+#Purpose: Perform FPCA with flight data. 
+
+
+library(refund)
+library(tidyverse)
+library(RFPCA)
+library(abind)
+library(pracma)
+library(Rmpfr)
+library(ggplot2)
+
+options(digits = 22)  # or any number up to 22
+
+
+fd <- as.data.frame(read.csv('smoothed_less_p_data_flights_egll_esgg1d2.csv')  ) #Reads data on processed curves, latitude. 
+a <- pi/180
+fd <- fd%>%mutate(x_s0 = sin(a*x_org)*cos(a*y_org), y_s0 = sin(a*x_org)*sin(a*y_org), z_s0 = cos(a*x_org) )
+obs_xyz <- fd%>%select(n, x_s0, y_s0, z_s0)
+
+obs_t <- fd%>%select(n, t)
+
+no_t <- length(unique(obs_t$t))
+n_un <- max(unique(obs_xyz$n))
+obs_xyz_list <- list()
+obs_t_list <- list()
+
+for(i in 1:n_un){
+  curr_n <- t(obs_xyz%>%filter(n == i)%>%select(x_s0, y_s0, z_s0))
+  curr_t <- t(obs_t%>%filter(n == i)%>%select(t))
+  rownames(curr_n) <- NULL
+  rownames(curr_t) <- NULL
+  
+  obs_xyz_list[[i]] <- curr_n 
+  obs_t_list[[i]] <- c(curr_t)
+  
+}
+
+bw <- 00.1
+kern <- 'gauss'
+K <- 12
+nRegG = 100
+mfdSp <- structure(1, class='Sphere')
+resSp <- RFPCA(obs_xyz_list, obs_t_list, list(userBwMu=bw, userBwCov=bw * 2, npoly = 3,  nRegGrid = nRegG, kernel=kern, error = FALSE, maxK=K, mfd=mfdSp)) #Performs RFPCA on sparsified observations. 
+
+tol <-  0.0001
+
+mu_est <- t(resSp$muWork)
+mu_org_grid <- t(resSp$muObs) #Gets estimated mean function on matrix of dimension (# points in grid) x 3
+
+norm_mu_sphr <- sqrt(mu_org_grid[,1]^2 + mu_org_grid[,2]^2 + mu_org_grid[,3]^2)
+is_on_sphere_mu_sphr <- abs(norm_mu_sphr - 1) < tol
+prop_mu_on_sphere_sphr <- sum(is_on_sphere_mu_sphr)/length(is_on_sphere_mu_sphr)
+
+phi <- resSp$phi # nWorkGrid x D x K
+
+to_sphere <- function(mu, phi, l){
+  phi_fin <- phi
+  
+  for(i in 1:l){
+    v <- phi_fin[i, ]
+    curr_mu <- mu[i,]
+    v_n <- sqrt(sum(v^2) )
+    phi_fin[i,] <- cos(v_n)*curr_mu + sin(v_n)*v*(1/v_n)
+  }
+  phi_fin
+}
+
+
+
+sph_obs_df <- data.frame(x = mu_est[,1], y = mu_est[,2], z = mu_est[,3], lbl = rep('mu_est', length(mu_est[,1])    )    ) #Stores mean function in appropiate dataframe. 
+
+
+for (i in 1:resSp$K){
+  print("------------------------------------------------------------------------------------------------------")
+  c_phi <- phi[,,i]
+  s_phi <- to_sphere(mu_est, (1 )*c_phi, nRegG) #3*sqrt(resSp$lam[i])  
+  
+  aux <- s_phi[,1]^2 + s_phi[,2]^2 + s_phi[,3]^2
+  print(mean(aux))
+  print(var(aux))
+  
+  curr_df <- data.frame(x = s_phi[,1], y = s_phi[,2], z = s_phi[,3],lbl = rep(paste('phi',as.character(i)), length(s_phi[,1])))
+  curr_df_e <-data.frame(x = c_phi[,1], y = c_phi[,2], z = c_phi[,3],lbl = rep(paste('phi_eu',as.character(i)), length(s_phi[,1]))      )
+  sph_obs_df <- rbind(sph_obs_df, curr_df)
+  sph_obs_df <- rbind(sph_obs_df, curr_df_e)
+}
+
+
+write.csv(sph_obs_df, 'rfpca_spherical_flights_egll_esgg1d2.csv')
+
+
+xi <- resSp$xi #Matrix of dimensions n_un x k, i.e. its element [i, ] corresponds to the estimated scores for the i-th observed function. 
+phi_org_grid <- resSp$phiObsTrunc #Gets tensor of dimensions (# points in grid) x 3 x k, i.e. its element [,,i] corresponds to the i-th estimated eigenfunction. 
+mu_org_grid <- t(resSp$muObs) #Gets estimated mean function on matrix of dimension (# points in grid) x 3
+
+
+#-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+
+phi_acc <- lapply(rep(1, n_un), function(x){mu_org_grid}) #Creates a vector with the mean function repeated n_un times i.e. phi_acc_0[[i]] = mu_org_grid for i = 1,2,...., n_un
+phi_acc <- abind(phi_acc, along = 3) #Converts phi_acc to a tensor of dimensions (# points in grid) x 3 x n_un 
+phi_acc <- array(phi_acc, dim = c(as.vector(dim(phi_acc)),1)) #Converts phi_acc to a tensor of dimensions (# points in grid) x 3 x n_un x 1. i.e. to later store everything in the same tensor. 
+
+
+
+for(i in 1:K){
+  phi_acc_curr <- lapply(xi[,i], function(x){x*phi_org_grid[,,i]}) #Multiples individual observation scores by i-th estimated eigenfunction. 
+  phi_acc_curr <- abind(phi_acc_curr, along = 3) #Converts phi_acc_curr to a tensor of dimensions (# points in grid) x 3 x n_un
+  phi_acc_curr <- array(phi_acc_curr, dim = c(as.vector(dim(phi_acc_curr)),1))#Converts phi_acc_curr to a tensor of dimensions (# points in grid) x 3 x n_un x 1. i.e. to later store everything in the same tensor. 
+  phi_acc <- abind(phi_acc, phi_acc_curr, along = 4) #Adds phi_acc_curr to phi_acc as a new slice. 
+} #Obs: in the end, phi_acc will be a tensor of dimension (#points in grid) x 3 x n_un x K. 
+
+
+geodesic_dist_sphere <- function(p, q, l){ #p, and q are two matrices of the same dimensions
+
+  warn <- 0
+  prob_x <- 0
+  inn_prods <- rep(0,l)
+  toRet <- rep(0, l) #Creates vector of zeros to store output. 
+  for( i in 1:l){
+    inn_prod <- sum(p[i,]*q[i,]) #Takes inner product between p[i,] and q[i,]
+    inn_prods[i] <- inn_prod
+    toRet[i] <- acos(inn_prod) #Takes spherical distance between ith rows p[i,] and q[i,]
+    if(is.nan(toRet[i]) ){
+      
+      if(sqrt(sum(p[i,]^2))>1  | sqrt(sum(p[i,]^2))<1 ){
+        prob_x <- prob_x + 1
+      }
+      warn <- warn + 1
+      if(inn_prod > 1)
+        toRet[i] <- acos(1)
+      else 
+        toRet[i] <- acos(-1)
+    }
+      
+  }
+  
+  #May I be forgiven for this way of programming. Warning, warn_track and n_curr are modified within the geodesic_dist_sphere function. 
+  warn_track_curr <- data.frame(n = c(n_curr), warn = c(warn), prop = c(warn/l), k_n = c(k_curr), prob_x = c(prob_x), diff_warn_x = c(abs(warn-prob_x)) , max_inn_prod = c(max(inn_prods)), min_inn_prod = c(min(inn_prods)))
+  warn_track <<- rbind(warn_track, warn_track_curr)
+  n_curr <<- n_curr + 1
+  
+  return (toRet)
+
+ } #Asumes that both 
+
+
+warn_track <- data.frame(n = c(0), warn = c(0), prop = c(0), k_n = c(0), prob_x = c(0), diff_warn_x = c(0), max_inn_prod = c(NaN), min_inn_prod = c(NaN))
+k_curr <- 0
+n_curr <- 1
+
+obs_X <- abind(obs_xyz_list, along = 3)
+norm_X <- matrix(0, nrow = n_un, ncol = no_t)
+
+for(i in 1:n_un){
+  curr_X <- obs_X[,,i]
+  norm_X[i, ] <- sqrt( curr_X[1,]^2 + curr_X[2,]^2  +curr_X[3,]^2 )
+}
+
+
+hist(norm_X)
+max(norm_X)
+min(norm_X)
+
+
+geod_dist <- lapply(1:n_un, function(x){  geodesic_dist_sphere(t(obs_X[,,x]), mu_org_grid, no_t )         })
+geod_dist <- abind(geod_dist, along = 2)
+geod_dist <- array(geod_dist, dim = c(as.vector(dim(geod_dist)), 1))
+  
+aux_ten <-array(0, c(no_t, 3, n_un))#Creates auxiliary tensor to accumulate xi*phi(t) values up to the k-th eigenfunction. 
+
+for(k in 1:K){
+  
+  k_curr <- k
+  n_curr <- 1
+  
+  curr_acc_phi <- phi_acc[,,,k+1] #Obtains values xi_k phi_k(t) for the original grid. 
+  aux_ten <- aux_ten + curr_acc_phi #Obtains the previous + xi_k phi_k(t)
+  curr_acc_k_sphr <- lapply(1:n_un, function(x){  to_sphere(mu_org_grid, aux_ten[,,x], no_t)})#Passes sum(xi_k phi_k(t)) to the sphere. 
+  curr_acc_k_sphr <- abind(curr_acc_k_sphr, along = 3) #Accomodates the former as a tensor of dimension no_t x 3 x n_un for better management.
+  l <- sqrt(curr_acc_k_sphr[,1,]^2 + curr_acc_k_sphr[,2,]^2 + curr_acc_k_sphr[,3,]^2) #Computes pointwise norm to see whether points are on the unit sphere. 
+  lmean <-apply(l, 2, mean) #Calculates, by observation, the mean of the norm values.
+  lvar <- apply(l, 2, var) #Calculates, by observation, the variance of the norm values.
+  
+  curr_geod_dist <- lapply(1:n_un, function(x){geodesic_dist_sphere( t(obs_X[,,x])  , curr_acc_k_sphr[,,x] ,no_t)}) #Calculates, pointwise, the geodesic distance for exp(sum xi_k phi_k(t)) to the actual observation. 
+  curr_geod_dist <- abind(curr_geod_dist, along = 2) #Accomodates the formeer as a tensor of dimension no_t x n_un for better management. 
+  #phi_acc <- abind(phi_acc, phi_acc_curr, along = 4) #Adds phi_acc_curr to phi_acc as a new slice. 
+  geod_dist <- abind(geod_dist, curr_geod_dist, along = 3)#Adds curr_dist_sqrd to dist_sqrd as a new slice.
+  print("..........................................................................................................")
+  print("mean of means")
+  print(mean(lmean))
+  print("Mean of variances")
+  print(mean(lvar))
+  print("..........................................................................................................")
+}
+
+
+ggplot(warn_track, aes(x =warn)) +geom_histogram(bins = 7, fill = "steelblue") +facet_wrap(~ k_n) +theme_minimal()
+ggplot(warn_track, aes(x =prob_x)) +geom_histogram(bins = 7, fill = "blue") +facet_wrap(~ k_n) +theme_minimal()
+ggplot(warn_track, aes(x =diff_warn_x)) +geom_histogram(bins = 7, fill = "cyan") +facet_wrap(~ k_n) +theme_minimal()
+hist(warn_track$max_inn_prod)
+hist(warn_track$min_inn_prod)
+
+#--------------------------------------------------------------------------------------------------------------------------------------------------
+
+dist_sqrd <- matrix(0, nrow = n_un, ncol = K + 1)
+
+for(t in 1:no_t){
+  dist_sqrd <- dist_sqrd + geod_dist[t, ,]^2
+}
+U0 <- mean(dist_sqrd[,1])
+U_m0 <- colMeans(dist_sqrd[, 1:K+1])
+FVE <- ((U0 - U_m0)/U0)
+
+FVE_k_d <- formatC(FVE, digits = 4, format = "fg", flag = "#")
+FVE_to_95 <- FVE >= 0.95
+
+#-----------------------------------------------------------------------------------------------------------------------
+#A little sanity check, see if final calculations coincide with what they are supposed to be. 
+
+xi <- resSp$xi #Matrix of dimensions n_un x k, i.e. its element [i, ] corresponds to the estimated scores for the i-th observed function. 
+phi_org_grid <- resSp$phiObsTrunc #Gets tensor of dimensions (# points in grid) x 3 x k, i.e. its element [,,i] corresponds to the i-th estimated eigenfunction. 
+mu_org_grid <- t(resSp$muObs) #Gets estimated mean function on matrix of dimension (# points in grid) x 3
+
+
+u0g <- lapply(1:n_un, function(x){  geodesic_dist_sphere(t(obs_X[,,x]), mu_org_grid, no_t )         })
+
+u0 <- 0
+for( i in 1:n_un){
+  u0 <- u0 + sum( u0g[[i]]^2)
+}
+u0 <- u0/n_un
+
+
+phi1 <- phi_org_grid[,,1]
+phi2 <- phi_org_grid[,,2]
+
+
+xi1 <- xi[,1]
+xi2 <- xi[,2]
+u1 <- c()
+
+for( i in 1:n_un){
+  
+  curr_sphr <- to_sphere(mu_org_grid, xi1[i]*phi1 + xi2[i]*phi2, no_t)
+  curr <- geodesic_dist_sphere(t(obs_X[,,i]),  curr_sphr, no_t  )
+  u1 <- c(u1, sum(curr^2))
+}
+
+u1 <- mean(u1)
+
+curr_sphr <- to_sphere(mu_org_grid, xi1[1]*phi1, no_t)
+curr <- geodesic_dist_sphere(t(obs_X[,,1]),  curr_sphr, no_t  )
+l <- geod_dist[,,2]
+l <- l[,2]
+
+
+ll <- dist_sqrd[1,2]
+aa <- sum(curr^2)
+
+
+
+
+#----------------------------------------------------------------------------------------------------------------------------
+
+
+mu_est_work<- t(resSp$muWork)
+phi_work <- resSp$phiObsTrunc
+
+work_obs_df <- data.frame(x = mu_est_work[,1], y = mu_est_work[,2], z = mu_est_work[,3], lbl = rep('mu_est', length(mu_est_work[,1]))    )
+aux <- mu_est_work[,1]^2 + mu_est_work[,2]^2 + mu_est_work[,3]^2 
+
+
+print("Principal components")
+
+for(i in 1:resSp$K){
+  print("------------------------------------------------------------------------")
+  c_phi <- phi_work[,,i]
+  curr_df <- data.frame(x = c_phi[,1], y = c_phi[,2], z = c_phi[,3],    lbl =  rep(paste('phi',as.character(i)), length(c_phi[,1]))        )
+  
+  aux <- c_phi[,1]^2 + c_phi[,2]^2 + c_phi[,3]^2
+  print(mean(aux))
+  print(var(aux))
+  
+  work_obs_df <- rbind(work_obs_df, curr_df)
+}
+
+write.csv(work_obs_df, 'rfpca_work_flights_egll_esgg1d2.csv')
+
+
+
+#----------------------------------------------------------------------------------------------------------------------------------
+
+
+var_lam <- resSp$lam
+pve_ish <- var_lam/sum(var_lam)
+
+pve <- cumsum(pve_ish)
+pve_k_d <- formatC(pve, digits = 4, format = "fg", flag = "#")
+
+pve_to_95 <- pve >= 0.95
+
+mfd <- structure(1, class='Euclidean')
+resEu <- RFPCA(obs_xyz_list, obs_t_list, list(userBwMu=bw, userBwCov=bw * 2, npoly = 3,  nRegGrid = nRegG, kernel=kern, error = FALSE,  maxK=K, mfd=mfd)) #Performs RFPCA on sparsified observations. 
+
+mu_est_eu <- t(resEu$muWork)
+mu_est_eu_obs <- t(resEu$muObs)
+
+
+
+norm_eu_mu <- sqrt(  mu_est_eu_obs[,1]^2 + mu_est_eu_obs[,2]^2 +mu_est_eu_obs[,3]^2)
+is_on_sphere_mu_eu <- abs(norm_eu_mu - 1) < tol 
+
+
+prop_mu_on_sphere_eu <- 100*(sum(is_on_sphere_mu_eu)/length(is_on_sphere_mu_eu) )
+
+phi_eu <- resEu$phi
+
+phi_eu_obsTrunc <- resEu$phiObsTrunc
+aux <- phi_eu_obsTrunc[, 1, ]^2 + phi_eu_obsTrunc[,2,]^2 + phi_eu_obsTrunc[,3,]^2
+norm_phi_eu <- (sqrt(aux) - ones(no_t, K) ) < tol
+prop_norm_phi_eu <- 100*(colSums(norm_phi_eu)/no_t)
+
+
+##----------------------------------------------------------------------------------------------------
+##A little sanity check, 
+
+phi1_eu <- phi_eu_obsTrunc[,,1]
+aux2 <- phi1_eu[,1]^2 + phi1_eu[,2]^2  + phi1_eu[,3]^2 
+l <- aux[,1]
+l <- l - aux2
+
+##----------------------------------------------------------------------------------------------------
+
+
+
+
+
+eu_obs_df <- data.frame(x = mu_est_eu[,1], y = mu_est_eu[,2], z = mu_est_eu[,3], lbl = rep('mu_est', length(mu_est[,1]))    )
+aux <- mu_est_eu[,1]^2 + mu_est_eu[,2]^2 + mu_est_eu[,3]^2 
+
+print("Eucledian mean")
+print(mean(aux))
+print(var(aux))
+
+print("Principal components")
+
+for(i in 1:resEu$K){
+  print("------------------------------------------------------------------------")
+  c_phi <- phi_eu[,,i]
+  curr_df <- data.frame(x = c_phi[,1], y = c_phi[,2], z = c_phi[,3],    lbl =  rep(paste('phi',as.character(i)), length(c_phi[,1]))        )
+  
+  aux <- c_phi[,1]^2 + c_phi[,2]^2 + c_phi[,3]^2
+  print(mean(aux))
+  print(var(aux))
+  
+  eu_obs_df <- rbind(eu_obs_df, curr_df)
+}
+
+
+var_lam_eu <- resEu$lam
+pve_ish_eu <- var_lam_eu/sum(var_lam_eu)
+write.csv(eu_obs_df, 'rfpca_euclidean_flights_egll_esgg1d2.csv')
+
+
+matplot( t(resEu$muObs), type='l', lty=2)
+matplot(t(resSp$muObs), type='l', lty=1, add=TRUE)
diff --git a/Planes/Planes2.R b/Planes/Planes2.R
new file mode 100644
index 0000000..67bb78b
--- /dev/null
+++ b/Planes/Planes2.R
@@ -0,0 +1,75 @@
+#Writen by: Lars Daniel Johansson Niño
+#Created date: 08/08/2024
+#Purpose: Perform FPCA with flight data. 
+
+
+library(refund)
+library(dplyr)
+library(tidyverse)
+library(RFPCA)
+
+fd <- as.data.frame(read.csv('smoothed_less_p_data_flights_egll_esgg1d2.csv')  ) #Reads data on processed curves, latitude. 
+
+obs_xy <- fd %>% select(n, x_org, y_org)
+obs_t <- fd%>%select(n, t)
+
+n_un <- max(unique(obs_xy$n))
+obs_xy_list <- list()
+obs_t_list <- list()
+
+for(i in 1:n_un){
+  curr_n <- t(obs_xy%>%filter(n == i)%>%select(x_org, y_org))
+  curr_t <- t(obs_t%>%filter(n == i)%>%select(t))
+  rownames(curr_n) <- NULL
+  rownames(curr_t) <- NULL
+  
+  obs_xy_list[[i]] <- curr_n 
+  obs_t_list[[i]] <- c(curr_t)
+  
+}
+
+
+
+l <- dim(curr_t)[2]
+
+bw <- 00.1
+kern <- 'gauss'
+K <- 12
+nRegG = 80
+
+mfd <- structure(1, class='Euclidean')
+resEu <- RFPCA(obs_xy_list, obs_t_list, list(userBwMu=bw, userBwCov=bw * 2, npoly = 3,  nRegGrid = nRegG, kernel=kern, maxK=K, mfd=mfd)) #Performs RFPCA on sparsified observations. 
+
+
+
+mu_est_eu <- t(resEu$muWork)
+phi_eu <- resEu$phi
+
+eu_obs_df <- data.frame(x = mu_est_eu[,1], y = mu_est_eu[,2], lbl = rep('mu_est', length(mu_est_eu[,1]))    )
+
+
+print("Principal components")
+
+for(i in 1:resEu$K){
+  print("------------------------------------------------------------------------")
+  
+  c_phi <- phi_eu[,,i]
+  curr_df <- data.frame(x = c_phi[,1], y = c_phi[,2],    lbl =  rep(paste('phi',as.character(i)), length(c_phi[,1]))        )
+
+  
+  eu_obs_df <- rbind(eu_obs_df, curr_df)
+}
+
+
+var_lam_eu <- resEu$lam
+pve_ish_eu <- var_lam_eu/sum(var_lam_eu)
+pve2 <- cumsum(pve_ish_eu)
+
+pve2_k_d <- formatC(pve2, digits = 4, format = "fg", flag = "#")
+
+write.csv(eu_obs_df, 'fpca_org_flights_egll_esgg1d2.csv')
+
+
+matplot( resEu$muObs[1,], resEu$muObs[2,], type='l', lty=2)
+
+resEu$muReg
\ No newline at end of file
diff --git a/Planes/ves_res_rfpca.py b/Planes/ves_res_rfpca.py
new file mode 100644
index 0000000..ac7ea4c
--- /dev/null
+++ b/Planes/ves_res_rfpca.py
@@ -0,0 +1,149 @@
+import os
+scrpt_path = os.path.abspath(__file__) #Gets absolute path of the current script. 
+scrpt_dir = os.path.dirname(scrpt_path) # Get the directory name of the current script
+os.chdir(scrpt_dir) #Chenges working direct  ory to the python file´s path. 
+
+import matplotlib.pyplot as plt
+import numpy as np
+import polars as pl
+import pandas as pd
+from matplotlib import cm
+from mpl_toolkits.mplot3d import Axes3D
+
+
+open_file_name =  'rfpca_euclidean_flights_egll_esgg1d2.csv' 
+#open_file_name =  'rfpca_spherical_flights_egll_esgg1d2.csv'
+#open_file_name = 'rfpca_work_flights_egll_esgg1d2.csv' 
+#open_file_name =  'rfpca_flights_egll_esgg1d2.csv'  
+#open_file_name = 'fpca_org_flights_egll_esgg1d2.csv'
+
+subs = 's0'
+color_col = 'cs'
+
+K =12
+
+sphere_data = pd.read_csv(open_file_name)
+
+fig, ax = plt.subplots(4, 3, subplot_kw={'projection': '3d'}, figsize=(15, 9))
+ax = ax.flatten()
+
+
+cmap = cm.get_cmap("tab10", K) 
+
+#fig = plt.figure()
+#ax = fig.add_subplot(projection='3d')
+
+#P = np.zeros((450,3, K))
+
+for i in range(K):
+    phi = sphere_data[sphere_data['lbl'] == ('phi ' + str(i+1))]
+    color = cmap(i)  # Assign a unique color
+
+    
+    x = phi['x'].to_numpy()
+    y = phi['y'].to_numpy()
+    z = phi['z'].to_numpy()
+
+    norm = np.sqrt(x**2 + y**2 + z**2)
+    print(f"phi {i+1}  mean of norms: {np.mean(norm)} \nvariance of norm {np.var(norm)}\n")
+
+    #P[:, 0, i] = x
+    #P[:,1,i] = y 
+    #P[:, 2, i] = z
+    ax[i].plot(x,y,z,  linewidth = 2.2, alpha = 0.9, color = color,  label=f"EF {i+1}")
+    ax[i].legend()
+
+r = 2
+print(f"With respect to {r}")
+#for i in range(0,K):
+    #L = P[:,:, i] - P[:, :, r]
+    #l = np.sum(L)/L.shape[0]*L.shape[1]
+    #print(f"For {i} mean entry difference {l}")
+
+
+
+
+
+
+mu = sphere_data[sphere_data['lbl'] == 'mu_est']
+mu_x = mu['x'].to_numpy()
+mu_y = mu['y'].to_numpy()
+mu_z = mu['z'].to_numpy()
+
+
+phi_mu = np.arccos(mu_z)
+theta_mu = np.arcsin(mu_y/np.sin(phi_mu) )
+
+phi_min = np.min(phi_mu)
+phi_max = np.max(phi_mu)
+
+theta_min = np.min(theta_mu)
+theta_max = np.max(theta_mu)
+l = 0.01
+
+
+figS = plt.figure(figsize=(8, 6))
+axS = figS.add_subplot(111, projection='3d')
+
+
+# Plot the unit sphere
+#uS = np.linspace(0, 2 * np.pi, 100)
+#vS = np.linspace(0, np.pi, 100)
+
+uS = np.linspace(theta_min-l, theta_max+l, 100)
+vS = np.linspace(phi_min-l, phi_max+l, 100)
+
+xS = np.outer(np.cos(uS), np.sin(vS))
+yS = np.outer(np.sin(uS), np.sin(vS))
+zS = np.outer(np.ones_like(uS), np.cos(vS))
+
+axS.plot_surface(xS, yS, zS, color='lightgray', alpha=0.3, linewidth=0)
+
+axS.plot(mu_x, mu_y, mu_z, color = 'red')
+
+#ax.plot(mu_x ,mu_y, mu_z, color = "red", linewidth = 2)
+
+#--Create figure to plot data--------------------------------------------------------------------------------------------------------------
+
+flight_data_file =  'smoothed_less_p_data_flights_egll_esgg1d2.csv'# 'egll_esgg_clustered_0vals.csv'  #'less_p_data_flights_egll_esgg1d2.csv'
+
+data = pl.read_csv(  flight_data_file, has_header= True)
+flight_n= data['n'].unique().to_list()
+
+fig1, ax1 = plt.subplots(1, 2, subplot_kw={'projection': '3d'}, figsize=(15, 9))
+ax1 = ax1.flatten()
+
+ax1[1].plot(mu_x ,mu_y, mu_z, color = "red", linewidth = 2.5, label = 'Mu')
+ax1[1].legend()
+
+
+unique_cs = data[color_col].unique().to_list() #Creates list of unique callsigns for flights. 
+cmap = plt.cm.viridis  # You can choose a different colormap, e.g., 'plasma', 'cividis', etc.
+colors = cmap(np.linspace(0, 1, len(unique_cs))) #Creates a unique color for each callsign. 
+
+for f in flight_n: #[0:1]:
+    print(f"In flight {f}")
+    cf = data.filter(pl.col('n') == f) #Selects those values with number n. 
+    cf_call_s = cf.get_column(color_col)[0]
+    cf_color = colors[unique_cs.index(cf_call_s)]
+    x_cf = cf[f'x_{subs}'].to_list()
+    y_cf = cf[f'y_{subs}'].to_list()
+    z_cf = cf[f'z_{subs}'].to_list()
+    axS.plot(x_cf, y_cf, z_cf, color = cf_color, alpha = 0.4)
+
+    t = cf['t'].to_list()
+    #print(f"Number of points {len(x_cf)}")
+    ax1[0].plot(x_cf, y_cf, z_cf, color = cf_color)
+ax1[0].legend()
+
+
+plt.tight_layout()
+plt.show()
+
+
+
+
+
+plt.show()
+
+
diff --git a/Volando_en_la_esfera0.pdf b/Volando_en_la_esfera0.pdf
new file mode 100644
index 0000000..cf1f709
Binary files /dev/null and b/Volando_en_la_esfera0.pdf differ