Skip to content

Add covariance consider cd #57

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
danlueck wants to merge 9 commits into
base: master
Choose a base branch
from
Open

Add covariance consider cd #57

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
aa5e751
added code and inputs for considering cd uncertainty for covariance p…
danlueck Nov 15, 2024
52265ee
Merge branch 'develop' into add-covariance-consider-cd
danlueck Nov 15, 2024
6d8a64a
Merge branch 'master' into add-covariance-consider-cd
danlueck Feb 27, 2025
a4b0d47
Merge branch 'master' into add-covariance-consider-cd
danlueck Apr 8, 2025
6a5c35e
Merge branch 'master' into add-covariance-consider-cd
danlueck Apr 8, 2025
d382642
Add RK8 integration for consider c_d propagation
danlueck Apr 8, 2025
8b07d68
Merge branch 'master' into add-covariance-consider-cd
Feb 24, 2026
32e6d7b
Remove duplicated variable
Feb 25, 2026
8c7861d
Update function headers
Feb 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 24 additions & 2 deletions src/neptune.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -574,6 +574,11 @@ subroutine propagate_set( &
real(dp) :: request_time ! requested time in numerical integration loop
real(dp) :: propCounterAtReset ! propCounter at last reset - required to prevent infinite loops
real(dp),dimension(6,6) :: set ! state error transition matrix
real(dp),dimension(7) :: sensitivity_matrix ! sensitivity matrix
real(dp),dimension(7,7) :: set_ext ! state error transition matrix
real(dp),dimension(7,7) :: cumSet_ext ! state error transition matrix
real(dp),dimension(7,7) :: covar_in_ext ! state error transition matrix
real(dp),dimension(7,7) :: covar_out_ext ! state error transition matrix
real(dp) :: start_epoch_sec ! start epoch in seconds (MJD)
type(kepler_t) :: kep ! mean kepler elements for correlation matrix computation
type(state_t) :: last_state_out ! saving the last state vector which has been written to output
Expand Down Expand Up @@ -728,6 +733,9 @@ subroutine propagate_set( &
if(neptune%numerical_integrator%getCovariancePropagationFlag()) then

cumSet = set_in%elem
cumSet_ext = 0.d0
cumSet_ext(1:6,1:6) = cumSet
cumSet_ext(7,7) = 1.d0
!call identity_matrix(cumSet) ! initial state error transition matrix is the unity matrix
call neptune%numerical_integrator%resetCountSetMatrix() ! the counter for the number of calls to the getStateTransitionMatrix routine is being reset

Expand Down Expand Up @@ -1034,7 +1042,8 @@ subroutine propagate_set( &
state_out%r, & ! <-- DBL() radius vector (km)
state_out%v, & ! <-- DBL() velocity vector (km/s)
request_time, & ! <-- DBL requested time
set & ! <--> DBL() state error transition matrix
set, & ! <--> DBL() state error transition matrix
sensitivity_matrix &
)
if(hasFailed()) return

Expand All @@ -1043,7 +1052,20 @@ subroutine propagate_set( &
set_out%elem = cumSet

!** compute new covariance matrix for given time
covar_out%elem = matmul(matmul(cumSet,covar_in%elem),transpose(cumSet))
if(.not. neptune%numerical_integrator%getCdCovFlag()) then
covar_out%elem = matmul(matmul(cumSet,covar_in%elem),transpose(cumSet))
else
set_ext = 0.d0
set_ext(1:6,1:6) = set
set_ext(1:7,7) = sensitivity_matrix
cumSet_ext = matmul(set_ext,cumSet_ext)
covar_in_ext = 0.d0
covar_in_ext(1:6,1:6) = covar_in%elem
covar_in_ext(7,7) = neptune%numerical_integrator%getCdCov()
covar_out_ext = matmul(matmul(cumSet_ext,covar_in_ext),transpose(cumSet_ext))
covar_out%elem(1:6,1:6) = covar_out_ext(1:6,1:6)
end if

if(neptune%correlation_model%getNoisePropagationFlag()) then
corrMat = neptune%correlation_model%getCorrelationMatrix(request_time)
covar_out%elem(1:6,1:6) = covar_out%elem(1:6,1:6) + corrMat
Expand Down
30 changes: 26 additions & 4 deletions src/neptuneClass.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ module neptuneClass
C_JUPITER, C_SATURN, C_NEPTUNE, C_URANUS, C_SRP, PAR_INT_RELEPS, PAR_INT_ABSEPS, PAR_INT_COV_STEP, OPTION_OUTPUT, &
OPTION_PN_LOOKUP, OPTION_SRP_CORRECT, OPTION_INT_LOGFILE, OPTION_HARMONICS, OPTION_EOP, OPTION_CORRELATION, &
C_EPOCH_END_GD, C_EPOCH_START_GD, C_OPT_SOL_FORECAST, C_PAR_INT_COV_STEP, C_PAR_CDRAG, &
PAR_CDRAG, C_PAR_CROSS_SECTION, PAR_CROSS_SECTION, C_PAR_MASS, PAR_MASS, C_COV_GEOPOTENTIAL, &
PAR_CDRAG, C_PAR_CROSS_SECTION, PAR_CROSS_SECTION, PAR_CDRAG_COV, C_PAR_MASS, PAR_MASS, C_COV_GEOPOTENTIAL, &
C_PAR_INT_COV_METHOD, C_OPT_STORE_DATA, C_OPT_ATMOSPHERE_MODEL, C_PAR_INT_METHOD, C_OUTPUT_STEP, &
C_GEOPOTENTIAL, C_OUTPUT_COV_UVW, C_OUTPUT_COV_ECI, C_OUTPUT_VAR_ECI, C_OUTPUT_VAR_UVW, &
C_OUTPUT_AMA, C_WIND, C_ATMOSPHERE, &
Expand All @@ -51,7 +51,7 @@ module neptuneClass
C_OUTPUT_ASA, C_OUTPUT_ACD, C_OUTPUT_ACG, C_OUTPUT_ACN, &
C_OUTPUT_ACM, C_OUTPUT_ACS, C_OUTPUT_ACJ, C_OUTPUT_ACV, &
C_OUTPUT_AME, C_OUTPUT_ACC, C_OUTPUT_FILES, C_OPT_HARMONICS, C_OPT_SRP_CORRECT, C_OPT_INT_LOG, &
C_OPT_PN_LOOKUP, C_OPT_EOP, C_CORRELATION, C_COV_MOON, C_COV_SUN, C_COV_SRP, C_COV_DRAG, C_COV_PROP, &
C_OPT_PN_LOOKUP, C_OPT_EOP, C_CORRELATION, C_CONSIDER_CD_COV, C_COV_MOON, C_COV_SUN, C_COV_SRP, C_COV_DRAG, C_COV_PROP, &
C_MANEUVERS, C_OCEAN_TIDES, C_ALBEDO, C_RUN_ID, INPUT_UNDEFINED, &
C_FILE_DE_EPHEM, C_FILE_LEAP_SPICE, C_FILE_TXYS, C_FILE_PROGRESS, C_OPT_PROGRESS, C_BOUNDARY_CHECK
use numint, only: Numint_class
Expand Down Expand Up @@ -485,6 +485,7 @@ subroutine initialize_input_array(this)
call this%set_input(parName=C_OPT_SOL_FORECAST, valType='double', initFlag=.true.)
call this%set_input(parName=C_OUTPUT_STEP, valType='double', initFlag=.true.)
call this%set_input(parName=C_OPT_STORE_DATA, valType='double', initFlag=.true.)
call this%set_input(parName=PAR_CDRAG_COV, valType='double', initFlag=.true.)

! ON/OFF parameters (also set default values, if available)
do i = 1, this%derivatives_model%get_neptune_perturbation_number()
Expand All @@ -501,6 +502,7 @@ subroutine initialize_input_array(this)
end do

call this%set_input(parName=C_CORRELATION, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_CONSIDER_CD_COV, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_OPT_EOP, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_OPT_PN_LOOKUP, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_OPT_INT_LOG, valType='boolean', initFlag=.true.)
Expand Down Expand Up @@ -1415,7 +1417,8 @@ integer function setNeptuneVar_char( &
!---------------------------------
case(C_ATMOSPHERE, C_SUN, C_MOON, C_SRP, C_SOLID_TIDES, &
C_MERCURY, C_VENUS, C_JUPITER, C_MARS, C_SATURN, C_URANUS, C_NEPTUNE, &
C_OCEAN_TIDES, C_MANEUVERS, C_CORRELATION, C_WIND, &
C_OCEAN_TIDES, C_MANEUVERS, C_CORRELATION, C_CONSIDER_CD_COV, &
C_WIND, &
C_ALBEDO, C_OPT_EOP, C_OPT_PN_LOOKUP, C_OPT_INT_LOG, &
C_OPT_SRP_CORRECT, &
C_OUTPUT_FILES, C_OUTPUT_ACC, C_OUTPUT_ACG, C_OUTPUT_ACD, &
Expand Down Expand Up @@ -1578,6 +1581,19 @@ integer function setNeptuneVar_char( &
call this%correlation_model%setNoisePropagationFlag(.false.)
end if

!** CD Covariance propagation
case(C_CONSIDER_CD_COV)
call this%set_input(parName=C_CONSIDER_CD_COV, val=val, set=.true.)


if(itemp == SWITCHED_ON) then
write(*,*) "consider cd turned on"
call this%numerical_integrator%setCdCovFlag(.true.)
else
write(*,*) "consider cd turned off"
call this%numerical_integrator%setCdCovFlag(.false.)
end if

case(C_OPT_EOP)
call this%set_input(parName=C_OPT_EOP, val=val, set=.true.)

Expand Down Expand Up @@ -1791,7 +1807,7 @@ integer function setNeptuneVar_char( &
case(C_PAR_MASS, C_PAR_CROSS_SECTION, C_PAR_CDRAG, C_OUTPUT_STEP, &
C_PAR_CREFL, C_PAR_REENTRY, C_PAR_INT_RELEPS, C_PAR_INT_ABSEPS, &
C_PAR_INT_COV_STEP, C_PAR_EARTH_RADIUS, C_OPT_STORE_DATA, &
C_OPT_SOL_FORECAST)
C_OPT_SOL_FORECAST, PAR_CDRAG_COV)

read(val,*,iostat=ios) dtemp

Expand Down Expand Up @@ -1867,6 +1883,10 @@ integer function setNeptuneVar_char( &
call this%output%write_to_output(C_OUTPUT_STEP, dtemp)
this%output_step = dtemp

case(PAR_CDRAG_COV)
call this%set_input(parName=PAR_CDRAG_COV, val=val, set=.true.)
call this%numerical_integrator%setCdCov(dtemp)

end select

end if
Expand Down Expand Up @@ -2830,6 +2850,7 @@ subroutine initializeInputArray(this)
call this%set_input(parName=C_OPT_SOL_FORECAST, valType='double', initFlag=.true.)
call this%set_input(parName=C_OUTPUT_STEP, valType='double', initFlag=.true.)
call this%set_input(parName=C_OPT_STORE_DATA, valType='double', initFlag=.true.)
call this%set_input(parName=PAR_CDRAG_COV, valType='double', initFlag=.true.)

! ON/OFF parameters (also set default values, if available)
do i = 1, this%derivatives_model%get_neptune_perturbation_number()
Expand All @@ -2846,6 +2867,7 @@ subroutine initializeInputArray(this)
end do

call this%set_input(parName=C_CORRELATION, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_CONSIDER_CD_COV, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_OPT_EOP, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_OPT_PN_LOOKUP, valType='boolean', initFlag=.true.)
call this%set_input(parName=C_OPT_INT_LOG, valType='boolean', initFlag=.true.)
Expand Down
2 changes: 2 additions & 0 deletions src/neptuneParameters.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,13 +89,15 @@ module neptuneParameters
character(len=*), parameter :: C_COV_SUN = "COVARIANCE_SUN"
character(len=*), parameter :: C_COV_SRP = "COVARIANCE_SRP"
character(len=*), parameter :: C_CORRELATION = "CORRELATION_MATRIX"
character(len=*), parameter :: C_CONSIDER_CD_COV = "CONSIDER_CD_COV"

character(len=*), parameter :: C_HARMONIC_C = "HARMONIC_C"
character(len=*), parameter :: C_HARMONIC_SD_C = "HARMONIC_SD_C"
character(len=*), parameter :: C_HARMONIC_SD_S = "HARMONIC_SD_S"
character(len=*), parameter :: C_HARMONIC_S = "HARMONIC_S"
character(len=*), parameter :: C_HARMONICS = "HARMONICS"
character(len=*), parameter :: C_INITIAL_COVARIANCE = "INITIAL_COVARIANCE"
character(len=*), parameter :: PAR_CDRAG_COV = "PAR_CDRAG_COV"
character(len=*), parameter :: C_INITIAL_STATE = "INITIAL_STATE"
character(len=*), parameter :: C_OPT_CANONICAL = "OPT_CANONICAL"
character(len=*), parameter :: C_OPT_AP_FORECAST = "OPT_AP_FORECAST"
Expand Down
Loading