Skip to content
GitLab
Commit 3a990d70 authored by Ananya's avatar Ananya
Browse files

Corrected peak maximum calculation

parent 93a52dbd
Pipeline #1615 failed with stage
in 2 seconds
Hide whitespace changes
Inline Side-by-side
CHANGELOG.md
View file @ 3a990d70
......@@ -19,6 +19,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
## [Unreleased]
## [0.6.25] - 2022-09-14
### Added
+ Peak maximum in R19_Occ_Info_PeakMax
### Changed
+ Corrected the calculation of elctron density maximum
## [0.6.24] - 2022-09-10
### Changed
......
pyproject.toml
View file @ 3a990d70
[tool.poetry]
name = "radiocc"
version = "0.6.24"
version = "0.6.25"
description = "Radio occultations"
license = "Apache-2.0"
authors = ["Ananya Krishnan <ananyakrishnaniiserk@gmail.com>", "Greg Henry <gregoire.henry@oma.be>"]
......
radiocc/constants.py
View file @ 3a990d70
......@@ -15,7 +15,7 @@ CC = C1 * kB # m^3 "K. L. Cahoy(2005)" Refractive Volume
m_bar = 7.221e-26 # kg "S.Tellmann (2013)"
V = 3.4e3 # m.s^-1
# radius thresholds
Threshold_Cor = 3.75e6 # 4.0e6 # 4.2E6 #4128432 # Wang et al. (2018) ??
Threshold_Cor = 3.6e6 # 4.0e6 # 4.2E6 #4128432 # Wang et al. (2018) ??
B_Cor_Type = "Linear" # 'Linear' #'Quadratic' # 'Linear' or 'Quadratic'
Threshold_neut_upp = 3466000 # M.Patzold et al.(2014) + Peter et al.(2014) ==> 80 km
Threshold_Surface = 3386000 # K. L. Cahoy (2005)
......
radiocc/old/R12_Electron_Density.py
View file @ 3a990d70
......@@ -19,8 +19,8 @@ def Elec(N_data, refractivity, e, eps0, me,Bande,fsup,c):
Ne = np.full(N_data,np.nan)#np.array([np.nan for i in range( N_data )], float)
index_peakmax = []
Ne_peakmax = []
#index_peakmax = []
#Ne_peakmax = []
#ET_peakmax = []
for i in range(N_data):
......@@ -30,12 +30,13 @@ def Elec(N_data, refractivity, e, eps0, me,Bande,fsup,c):
#Ne_dn[i] = - refractivity[i]*(8*np.pi**2*me*eps0*f_TX**2)/(e**2)*1E-6
#print(f_TX)
"""
for i in range(N_data):
if (Ne[i]==np.nanmax(np.array(Ne))):
index_peakmax.append(i) #index at maximum electron density
Ne_peakmax.append(Ne[i]) # maximum electron density
return Ne, index_peakmax, Ne_peakmax
"""
return Ne #, index_peakmax, Ne_peakmax
def TEC_Calc(N_data, Ne, Corr_ind, i_Surface, bend_radius, R_Mars):
......
radiocc/old/R17_Run_Routines1.py
View file @ 3a990d70
......@@ -373,7 +373,7 @@ def run(
if item == "IONO":
Ne, index_peakmax, Ne_peakmax = R12_Electron_Density.Elec(
Ne = R12_Electron_Density.Elec(
N_data, refractivity, e, eps0, me, Bande, fsup, c
)
TEC = R12_Electron_Density.TEC_Calc(
......@@ -435,6 +435,6 @@ def run(
)
if item == "IONO":
return refractivity, Ne, index_peakmax, Ne_peakmax, TEC, bend_radius, imp_param
return refractivity, Ne, TEC, bend_radius, imp_param
if item == "ATMO":
return refractivity, Nn, i_neut_upp, P_low, P_med, P_upp, T_low, T_med, T_upp, bend_radius
radiocc/old/R19_Occ_Info_PeakMax.py
View file @ 3a990d70
......@@ -7,8 +7,28 @@ import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pudb import set_trace as bp # noqa:F401
def occultation_info(planet,ground,index_peakmax,epoch,frame,sc,DATA_ID,DATA_DIR,Threshold_Surface,EPHE_DIR,distance):
def occultation_info(planet,ground,N_data,epoch,frame,sc,DATA_ID,DATA_DIR,Threshold_Surface,EPHE_DIR,distance, ALTITUDE, Ne):
index_peakmax = []
Ne_peakmax = []
alt1 = 200000
alt2 = 100000
idx1 = []
for i in range(N_data):
if (ALTITUDE[i]<=alt1)&(ALTITUDE[i]>=alt2):
idx1.append(i)
a = idx1[0]
b = idx1[-1]
for i in range(a,b):
if Ne[i] > 0:
if (Ne[i]==np.nanmax(np.array(Ne[a:b]))):
index_peakmax.append(i)
Ne_peakmax.append(Ne[i])
Ne_max = Ne_peakmax[0]
#index at the Maximum electron density, sometimes there maybe two or more indices with same electron density values,
# then we always choose the first index
index = index_peakmax[0]
......@@ -291,7 +311,7 @@ def occultation_info(planet,ground,index_peakmax,epoch,frame,sc,DATA_ID,DATA_DIR
################### FICHIER INFO_MEX......
FILE_OUT = 'INFO_AfterProcess.txt'
FILE_OUT = 'INFO_At_PeakMax.txt'
output = []
......@@ -300,6 +320,8 @@ def occultation_info(planet,ground,index_peakmax,epoch,frame,sc,DATA_ID,DATA_DIR
output.append('UTC : %s' % utc)
#output.append('Epoch : %s' %epoch)
output.append('Electron density Maximum [10^6 cm-3]: %f' %Ne_max)
output.append('MEX Sun Distance [km]: %f' %Dist_sun_sc)
......
radiocc/old/R4_State_Vectors_twoway_v3.py
View file @ 3a990d70
......@@ -66,15 +66,15 @@ def ephemerides(
vel_GS_dn = np.full((N_data, 3), np.nan)
#for i in range(N_data):
for i in range(N_data):
#time at which SC send signal to mars via downlink =
#time at which SC received signal from Mars via uplink
[states1, lt1] = spy.spkezr(SC, ET, Ref_frame, 'CN+S', GS)
tMEX1 = ET - lt1
[states2, lt2] = spy.spkezr(GS, tMEX1, Ref_frame, 'CN+S', SC)
tA1 = tMEX1 - lt2
for i in range(N_data):
[states1, lt1] = spy.spkezr(SC, ET[i], Ref_frame, 'CN+S', GS)
tMEX1 = ET - lt1
[states2, lt2] = spy.spkezr(GS, tMEX1, Ref_frame, 'CN+S', SC)
tA1 = tMEX1 - lt2
#for i in range(N_data):
#Calculation of uplink occultation time
tA = tA1[i]
tMEX = tMEX1[i]
......@@ -128,7 +128,7 @@ def ephemerides(
dkoab = koab - oldkoab
itoab += 1
#Calculation of downlink occultation time
......@@ -156,12 +156,12 @@ def ephemerides(
#print(tMEX, pos_SC_ssb)
#position of GS at the time of occulation during uplink in the solar system barycenter frame
[states8, lt8] = spy.spkpos('SSB', ET, Ref_frame, 'LT', GS)
pos_GS_ssb = states8
[states8, lt8] = spy.spkpos('SSB', ET[i], Ref_frame, 'LT', GS)
pos_GS_ssb[i] = states8[:3]
#position of Mars relative to spacecraft in working frame
pos_Mars_SCC = pos_Mars_ssb[i] - pos_SC_ssb[i]
pos_Mars_SC[i] = pos_Mars_ssb[i] - pos_SC_ssb[i]
#position of GS relative to spacecraft in working frame
pos_GS_SC[i] = pos_GS_ssb[i] - pos_SC_ssb[i]
......@@ -192,25 +192,25 @@ def ephemerides(
[states9,lt9] = spy.spkezr(GS, tA, Ref_frame, 'NONE', 'SSB')
pos_GS_ssb_tA[i] = states9[:3] * 1e3
vel_GS_ssb_tA[i] = states9[3:] * 1e3
#the position and velocity of mars wrt SSB at occultation time toab
[states10,lt10] = spy.spkezr("mars", toab, Ref_frame, 'NONE', 'SSB')
pos_Mars_ssb_toab[i] = states10[:3] * 1e3
vel_Mars_ssb_toab[i] = states10[3:] * 1e3
#The position and velocity of SC wrt SSB (m) at reception time tMEX
[states11,lt11] = spy.spkezr(SC, tMEX, Ref_frame, 'NONE', 'SSB')
pos_SC_ssb_tMEX[i] = states11[:3] * 1e3
vel_SC_ssb_tMEX[i] = states11[3:] * 1e3
#The position and velocity of GS at transmission time A wrt the position of Mars at occultation time toab
pos_GS_Mars_toab = pos_GS_ssb_tA[i] - pos_Mars_ssb_toab
vel_GS_Mars_toab = vel_GS_ssb_tA[i] - vel_Mars_ssb_toab
#The position and velocity of SC at tMEX wrt the position of Mars at occultation time toab
pos_SC_Mars_toab = pos_SC_ssb_tMEX -pos_Mars_ssb_toab
vel_SC_Mars_toab = vel_SC_ssb_tMEX -vel_Mars_ssb_toab
#The position and velocity of Mars wrt SSB (m) at occultation time tocd
#print(tocd)
[states12, lt12] = spy.spkezr( "mars", tocd, Ref_frame, 'NONE', 'SSB')
......@@ -227,7 +227,7 @@ def ephemerides(
pos_GS_Mars_tocd = pos_GS_ssb_ET - pos_Mars_ssb_tocd
vel_GS_Mars_tocd = vel_GS_ssb_ET - vel_Mars_ssb_tocd
#The position of SC at transmission time tMEX wrt the position of Mars at occultation time tocd
pos_SC_Mars_tocd = pos_SC_ssb_tMEX - pos_Mars_ssb_tocd
vel_SC_Mars_tocd = vel_SC_ssb_tMEX - vel_Mars_ssb_tocd
......
radiocc/reconstruct.py
View file @ 3a990d70
......@@ -187,7 +187,7 @@ def run( # noqa: C901
re,
f,
) = R5_Foot_Print.long_lat(
N_data, pos_SC_dn, vel_SC_dn, ET, tMEX, Ref_frame, l2_data.FOLDER_TYPE
N_data, pos_SC_dn, vel_SC_dn, ET, tocd, Ref_frame, l2_data.FOLDER_TYPE
)
# R7 coordinate frames
......@@ -418,7 +418,7 @@ def run( # noqa: C901
radiocc.gui.button_next.set_sensitive(False)
else:
if item == "IONO":
refractivity, Ne, index_peakmax, Ne_peakmax, TEC, bend_radius, imp_param = R17_Run_Routines1.run(
refractivity, Ne, TEC, bend_radius, imp_param = R17_Run_Routines1.run(
Doppler,
Doppler_debias,
distance,
......@@ -594,7 +594,7 @@ def run( # noqa: C901
if radiocc.cfg.skip_correction_guess and radiocc.gui is not None:
radiocc.gui.button_next.set_sensitive(True)
if item == "IONO":
refractivity, Ne, index_peakmax, Ne_peakmax, TEC, bend_radius, imp_param = R17_Run_Routines1.run(
refractivity, Ne, TEC, bend_radius, imp_param = R17_Run_Routines1.run(
Doppler,
Doppler_debias2,
distance,
......@@ -723,23 +723,27 @@ def run( # noqa: C901
N_data, fsup, error, c, distance, Threshold_neut_upp, V, e, eps0, me, Bande, CC
)
R19_Occ_Info_PeakMax.occultation_info(
NAIF_P,
NAIF_GS,
index_peakmax,
ET,
Ref_frame,
NAIF_SC,
DATA_PRO,
DATA_ID_i_profile,
Threshold_Surface,
EPHE_DIR,
distance,
)
ALTITUDE = bend_radius - d_radius * 1e3
if item == "IONO":
R19_Occ_Info_PeakMax.occultation_info(
NAIF_P,
NAIF_GS,
N_data,
ET,
Ref_frame,
NAIF_SC,
DATA_PRO,
DATA_ID_i_profile,
Threshold_Surface,
EPHE_DIR,
distance,
ALTITUDE,
Ne,
)
ALTITUDE = bend_radius - d_radius * 1e3
#ALTITUDE = bend_radius - d_radius * 1e3
# Return variables to be exported.
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment