<jupyter_start><jupyter_text>DescriptionIn this notebook we explore the attitude update alone, assuming we know the the astrometric parameters of the sources.<jupyter_code># Module import
from IPython.display import Image
import sys
# append to path the folder that contains the analytic scanner
sys.path.append('../GaiaLab/scan/analytic_scanner')
import copy
from tqdm import trange
import quaternion 

# local imports
import frame_transformations as ft
from scanner import Scanner
from satellite import Satellite
from source import Source

import constants as const

from agis import Agis
from agis import Calc_source
from agis_functions import *
from analytic_plots import *

# Ipython magics
%load_ext autoreload
%autoreload 2
%reload_ext autoreload<jupyter_output><empty_output><jupyter_text>For scaling: Table resuming the scaling that has to be done (not yet implemented). However it can be kept in mind while playing with the parameters.<jupyter_code>Image('./figures/scaling_factors.png')<jupyter_output><empty_output><jupyter_text>1.- **Initializing objects:**  1.1 satellite<jupyter_code>t_init = 0 
t_end =  t_init + 1
my_dt = 1/24/4 # [days]

# objects
spline_degree = 3

gaia = Satellite(ti=t_init, tf=t_end, dt= my_dt, k=spline_degree)
print('#knots:', len(gaia.s_w.get_knots()))<jupyter_output>#knots: 95<jupyter_text>1.2 scanner<jupyter_code>double_telescope = True
my_scanner = Scanner(zeta_limit=np.radians(50) ,double_telescope=double_telescope)<jupyter_output><empty_output><jupyter_text>End scanner<jupyter_code>num_times_for_sources=100 # We will create 3 source per given time
times_for_source = np.linspace(t_init, t_end, num=num_times_for_sources, endpoint=True)

alphas, deltas = generate_angles_of_sources(times_for_source, sat=gaia, noise_factor=0)

sources = []
for n in range(len(alphas)):
    sources.append(Source(str(n), np.degrees(alphas[n]), np.degrees(deltas[n]), 0, 0, 0, 0))<jupyter_output><empty_output><jupyter_text>2- Scan and results  2.1 Scan<jupyter_code>obs_times = []
calc_sources = []
real_sources = []

# scan the sources and keep as real sources only the ones that are observed
for i in trange(len(sources)):
    s = sources[i]
    my_scanner.scan(gaia, s, ti=t_init, tf=t_end)
    my_scanner.compute_angles_eta_zeta(gaia, s)
    if len(my_scanner.obs_times)>0:
        obs_times += my_scanner.obs_times
        real_sources.append(s)
        calc_sources.append(Calc_source(obs_times=my_scanner.obs_times.copy(), source=s))
sources = real_sources
obs_times = list(np.sort(obs_times))
# check for consistency
if list(np.array(obs_times)[np.where(np.array(obs_times)>t_end)]):
    raise ValueError('observed time greater than end time')<jupyter_output><empty_output><jupyter_text>2.2 Visualise results of scan<jupyter_code># Plot the distribution of sources and observation along the time axis
plt.figure()
plt.title('Created_sources:'+str(len(alphas))+'  ||  scanner-found sources:'+str(len(calc_sources)))
plt.plot(obs_times,np.zeros(len(obs_times)), 'r+', label='obs_times')
plt.plot(times_for_source, np.zeros(len(times_for_source)), 'b.', label='times of sources', alpha=0.5)
sampled_times = np.linspace(t_init, t_end, num=500)
plt.plot(sampled_times, np.zeros(len(sampled_times)), 'k,', label='original times')
plt.xlabel('time [days]'), plt.ylabel('unit value') # , plt.ylim((-1e-6, 1e-6))

plt.grid(), plt.legend(), plt.show();
# plot the distibution of observations as a histogram with bin-size approximatively equal 
# to the one the attitude knot interval
plt.figure()
plt.title('Observation per (approximative) knot interval')
knots = gaia.s_x.get_knots()
bins = int(len(knots))
print('#bins:', bins)
plt.hist(obs_times, bins=bins, align='mid');
plt.xlabel('time [days]'), plt.ylabel('number of observations')
plt.grid(), plt.show();<jupyter_output>#bins: 95<jupyter_text>3.- *Update Attitude*  3.1 Create Solver<jupyter_code>spline_degree = 3
Solver = Agis(sat=gaia, calc_sources=calc_sources, real_sources=sources, updating='attitude',
              double_telescope=double_telescope,
              attitude_splines=[gaia.s_w, gaia.s_x, gaia.s_y, gaia.s_z], 
              spline_degree=spline_degree,
              attitude_regularisation_factor=5e-1)<jupyter_output><empty_output><jupyter_text>3.1.1 add noise to attitudeTo ignore if we want to start directly with the correct attitude<jupyter_code># Ignore this cell if you don't want to modify the initial attitude
# Can be used to check that when recreating the splines in the solver we (almost) do not create additional errors
my_times = np.array(obs_times)
print('Error before Noise: ', Solver.error_function())
c_noise = Solver.att_coeffs * np.random.rand(Solver.att_coeffs.shape[0], Solver.att_coeffs.shape[1]) * 1e-4
last_coef = 4
Solver.att_coeffs[:last_coef] = Solver.att_coeffs[:last_coef] + c_noise[:last_coef]
Solver.actualise_splines()
print('Error after Noise: ', Solver.error_function())<jupyter_output>Error before Noise:  0.0
Error after Noise:  6.291176341250841e-10<jupyter_text>3.2 Visualize Solver data  3.2.1 Visualize attitude<jupyter_code># compare initial difference between satellite and solver attitudes
fig = compare_attitudes(gaia, Solver, obs_times)<jupyter_output><empty_output><jupyter_text>3.2.2 Visualize Matrices<jupyter_code># Compute Normal matrix
N_aa = Solver.compute_attitude_LHS()
A=N_aa.copy()
plt.figure(figsize=(8,8))
threshold = 0
A[np.where(A==threshold)] = A.max()
plt.imshow(A, vmin=None, vmax=None)
plt.colorbar()
plt.show();<jupyter_output><empty_output><jupyter_text>3.2.3 Visualize eigen properties of Normal matrix<jupyter_code>eig_vals, eig_vecs = np.linalg.eigh(N_aa)
plt.figure()

plt.title('Ordered Eigenvalues.')
plt.semilogy(eig_vals, label='eigenvalues')
plt.legend()
plt.grid(), plt.show();<jupyter_output><empty_output><jupyter_text>3.3 **Iterate**<jupyter_code>errors = []  # list to store the errors
errors.append(Solver.error_function())
Solver.iterate(1, verbosity=1)
errors.append(Solver.error_function())<jupyter_output>***** Iteration: 1 *****
Error after iteration: 3.096401558604955e-09<jupyter_text>3.3.1 Visualized residuals<jupyter_code>fig3 = multi_compare_attitudes_errors(gaia, Solver, obs_times)<jupyter_output><empty_output><jupyter_text>3.3.2 Iterate more<jupyter_code>for i in range(5):
    Solver.iterate(1, verbosity=1)
    errors.append(Solver.error_function())<jupyter_output>***** Iteration: 2 *****
Error after iteration: 2.8037851933738734e-09
***** Iteration: 3 *****
Error after iteration: 2.828735925713584e-09
***** Iteration: 4 *****
Error after iteration: 2.8262694393481134e-09
***** Iteration: 5 *****
Error after iteration: 2.8265239336384923e-09
***** Iteration: 6 *****
Error after iteration: 2.826494620090904e-09<jupyter_text>3.3.3 Visualize convergence<jupyter_code>plt.figure()
plt.plot(errors[:], label='error')
plt.grid(), plt.legend(), plt.xlabel('iterations'), plt.ylabel('error')
plt.show()
fig3 = multi_compare_attitudes_errors(gaia, Solver, obs_times);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1284246575342466, 'fraction_numerical': 0.032248858447488586, 'mean_word_length': 6.753318584070796, 'pattern_counts': {'":': 0, '<': 60, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '498140', 'license': 'MIT', 'n_tokens_mistral': 2704, 'n_tokens_neox': 2526, 'n_words': 536, 'path': 'notebooks/03-Attitude_update.ipynb', 'repo_name': 'bombrun/GaiaLab', 'size': '7008'}