diff --git a/src/neptune.f90 b/src/neptune.f90
index e05cbc8..3953ffe 100644
--- a/src/neptune.f90
+++ b/src/neptune.f90
@@ -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
@@ -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
@@ -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
@@ -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
diff --git a/src/neptuneClass.f90 b/src/neptuneClass.f90
index bfecd0e..d7f6070 100644
--- a/src/neptuneClass.f90
+++ b/src/neptuneClass.f90
@@ -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, &
@@ -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
@@ -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()
@@ -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.)
@@ -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, &
@@ -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.)
@@ -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
@@ -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
@@ -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()
@@ -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.)
diff --git a/src/neptuneParameters.f90 b/src/neptuneParameters.f90
index 41dd234..3724b27 100644
--- a/src/neptuneParameters.f90
+++ b/src/neptuneParameters.f90
@@ -89,6 +89,7 @@ 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"
@@ -96,6 +97,7 @@ module neptuneParameters
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"
diff --git a/src/numint.f90 b/src/numint.f90
index 15ee7ce..b89b5b8 100644
--- a/src/numint.f90
+++ b/src/numint.f90
@@ -126,6 +126,8 @@ module numint
integer, private :: lk ! k in last integration step for interpolator
integer, private :: fail ! number of consecutive failures
logical, private :: flag_srp ! becomes true as soon as SRP perturbation has been activated and varstormcow is initialized
+ real(dp), private :: cd_cov ! number of consecutive failures
+ logical, private :: cd_cov_flag ! becomes true as soon as SRP perturbation has been activated and varstormcow is initialized
@@ -141,6 +143,8 @@ module numint
procedure,public :: get_covariance_integration_method
procedure,public :: getStateTransitionMatrix
procedure,public :: get_current_step_size
+ procedure,public :: getCdCovFlag
+ procedure,public :: getCdCov
!** set
procedure,public :: resetCountIntegrator
@@ -154,6 +158,8 @@ module numint
procedure,public :: setCovariancePropagationFlag
procedure,public :: setSrpCorrect
procedure,public :: set_start_epoch
+ procedure,public :: setCdCovFlag
+ procedure,public :: setCdCov
!** others
procedure,public :: integrateStep
procedure,private :: varstormcow
@@ -291,6 +297,44 @@ subroutine setSrpCorrect(this,val)
return
end subroutine
+ !========================================================================
+ !
+ !> @anchor setCdCovFlag
+ !!
+ !> @brief Set the flag for the C_d consider covariance
+ !> @author Daniel Lück
+ !!
+ !> @date
+ !! - 15.11.2024 (initial design)
+ !!
+ !!
+ !------------------------------------------------------------------------
+ subroutine setCdCovFlag(this,val)
+ class(Numint_class),intent(inout) :: this
+ logical, intent(in) :: val
+ this%cd_cov_flag = val
+ return
+ end subroutine
+
+ !========================================================================
+ !
+ !> @anchor setCdCovFlag
+ !!
+ !> @brief Set the flag for the C_d consider covariance
+ !> @author Daniel Lück
+ !!
+ !> @date
+ !! - 15.11.2024 (initial design)
+ !!
+ !!
+ !------------------------------------------------------------------------
+ subroutine setCdCov(this,val)
+ class(Numint_class),intent(inout) :: this
+ real(dp), intent(in) :: val
+ this%cd_cov = val
+ return
+ end subroutine
+
!========================================================================
!
!> @anchor getIntegrationLogfileChannel
@@ -823,6 +867,7 @@ subroutine resetCountSetMatrix(this)
!! VB 24.01.2016 (fixed issue with call counter: now via modul variable and reset function, which is called by NEPTUNE)
!! VB 14.07.2017 (changed to time offset in parameter list and adding to start epoch)
!! DLU 02.03.2021 (added RK8 Integration)
+!! DLU 08.04.2025 (added consider cd approach)
!!
!!
!------------------------------------------------------------------------
@@ -840,7 +885,8 @@ subroutine getStateTransitionMatrix( this, &
r, &
v, &
reqt, &
- set &
+ set, &
+ sensitivity_matrix &
)
class(Numint_class),intent(inout) :: this
@@ -853,13 +899,19 @@ subroutine getStateTransitionMatrix( this, &
type(Thirdbody_class),intent(inout) :: thirdbody_model
type(Tides_class),intent(inout) :: tides_model
type(Derivatives_class),intent(inout) :: derivatives_model
- type(Reduction_type),intent(inout) :: reduction
+ type(Reduction_type),intent(inout) :: reduction
real(dp), dimension(3), intent(in) :: r
real(dp), dimension(3), intent(in) :: v
real(dp), intent(in) :: reqt
real(dp), dimension(setdim,setdim), intent(out) :: set
+ real(dp), dimension(setdim+1), intent(out) :: sensitivity_matrix
+ real(dp), dimension(3) :: r_itrf
+ real(dp), dimension(3) :: v_itrf
+ real(dp), dimension(3) :: acc_atmosphere
+ real(dp), dimension(6) :: pnm
real(dp), dimension(setdim,setdim) :: k1,k2,k3,k4,k5,k6,k7,k8,k9,k10
+ real(dp), dimension(3) :: k1_atm, k2_atm, k3_atm, k4_atm, k5_atm, k6_atm, k7_atm, k9_atm, k10_atm
real(dp), dimension(setdim,setdim) :: pdm
integer :: reset
@@ -887,6 +939,23 @@ subroutine getStateTransitionMatrix( this, &
!** reset state error transition matrix
call identity_matrix(set)
+ call reduction%inertial2earthFixed(r, v, time_mjd, r_itrf, v_itrf)
+
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ r, & ! <-- DBL(3) radius vector in GCRF
+ v, & ! <-- DBL(3) velocity vector in GCRF
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ time_mjd, & ! <-- DBL current time (MJD)
+ acc_atmosphere & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ pnm = 0.d0
+ pnm(4:6) = acc_atmosphere
+
!** get partial derivatives matrix
call derivatives_model%deriv_cov( &
gravity_model, &
@@ -908,6 +977,13 @@ subroutine getStateTransitionMatrix( this, &
+ 0.5d0*this%covIntegrationStep**2.d0*matmul(pdm,pdm) &
+ 1.d0/6.d0*this%covIntegrationStep**3.d0*matmul(matmul(pdm,pdm),pdm)
+ sensitivity_matrix(1:6) = this%covIntegrationStep*pnm &
+ + 0.5d0*this%covIntegrationStep**2.d0*matmul(pdm,pnm) &
+ + 1.d0/6.d0*this%covIntegrationStep**3.d0*matmul(matmul(pdm,pdm),pnm) &
+ + 1.d0/24.d0*this%covIntegrationStep**4.d0*matmul(matmul(matmul(pdm,pdm),pdm),pnm)
+
+ sensitivity_matrix(7) = 1.d0
+
!** update covariance matrix
!cov = matmul(matmul(set,cov),transpose(set))
@@ -1026,7 +1102,75 @@ subroutine getStateTransitionMatrix( this, &
if(hasFailed()) return
! RK4:
- set = set + this%covIntegrationStep/6.d0*(k1 + 2.0*k2 + 2.d0*k3 + k4)
+ set = set + this%covIntegrationStep/6.d0*(k1 + 2.d0*k2 + 2.d0*k3 + k4)
+ sensitivity_matrix = 0.
+
+ if(this%cd_cov_flag) then
+
+ call reduction%inertial2earthFixed(state(1)%r, state(1)%v, t1d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state(1)%r, &
+ state(1)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t1d, & ! <-- DBL current time (MJD)
+ k1_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call reduction%inertial2earthFixed(state(2)%r, state(2)%v, t2d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state(2)%r, &
+ state(2)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t2d, & ! <-- DBL current time (MJD)
+ k2_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state(2)%r, &
+ state(2)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t2d, & ! <-- DBL current time (MJD)
+ k3_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call reduction%inertial2earthFixed(state(3)%r, state(3)%v, t3d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state(3)%r, &
+ state(3)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t3d, & ! <-- DBL current time (MJD)
+ k4_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+
+ sensitivity_matrix(4:6) = this%covIntegrationStep/6.d0*(k1_atm + 2.0*k2_atm + 2.d0*k3_atm + k4_atm)
+ sensitivity_matrix(7) = 1.d0
+
+ end if
!write(52,'(36(e14.7e2,x))') set
!write(*,*)
@@ -1097,8 +1241,8 @@ subroutine getStateTransitionMatrix( this, &
t2, & ! <-- DBL requested time (s)
.false., &
reqt_loc, & ! <-- DBL current time (s)
- state_rk8(2)%r, & ! <--> DBL() radius vector (km)
- state_rk8(2)%v, & ! <--> DBL() velocity vector (km/s)
+ state_rk8(2)%r, & ! <--> DBL() radius vector (km)
+ state_rk8(2)%v, & ! <--> DBL() velocity vector (km/s)
reset, & ! <-- INT reset flag
dt & ! --> DBL propagated time (s)
)
@@ -1126,8 +1270,8 @@ subroutine getStateTransitionMatrix( this, &
t3, & ! <-- DBL requested time (s)
.false., &
reqt_loc, & ! <-- DBL current time (s)
- state_rk8(3)%r, & ! <--> DBL() radius vector (km)
- state_rk8(3)%v, & ! <--> DBL() velocity vector (km/s)
+ state_rk8(3)%r, & ! <--> DBL() radius vector (km)
+ state_rk8(3)%v, & ! <--> DBL() velocity vector (km/s)
reset, & ! <-- INT reset flag
dt & ! --> DBL propagated time (s)
)
@@ -1155,8 +1299,8 @@ subroutine getStateTransitionMatrix( this, &
t4, & ! <-- DBL requested time (s)
.false., &
reqt_loc, & ! <-- DBL current time (s)
- state_rk8(4)%r, & ! <--> DBL() radius vector (km)
- state_rk8(4)%v, & ! <--> DBL() velocity vector (km/s)
+ state_rk8(4)%r, & ! <--> DBL() radius vector (km)
+ state_rk8(4)%v, & ! <--> DBL() velocity vector (km/s)
reset, & ! <-- INT reset flag
dt & ! --> DBL propagated time (s)
)
@@ -1184,8 +1328,8 @@ subroutine getStateTransitionMatrix( this, &
t5, & ! <-- DBL requested time (s)
.false., &
reqt_loc, & ! <-- DBL current time (s)
- state_rk8(5)%r, & ! <--> DBL() radius vector (km)
- state_rk8(5)%v, & ! <--> DBL() velocity vector (km/s)
+ state_rk8(5)%r, & ! <--> DBL() radius vector (km)
+ state_rk8(5)%v, & ! <--> DBL() velocity vector (km/s)
reset, & ! <-- INT reset flag
dt & ! --> DBL propagated time (s)
)
@@ -1213,8 +1357,8 @@ subroutine getStateTransitionMatrix( this, &
t6, & ! <-- DBL requested time (s)
.false., &
reqt_loc, & ! <-- DBL current time (s)
- state_rk8(6)%r, & ! <--> DBL() radius vector (km)
- state_rk8(6)%v, & ! <--> DBL() velocity vector (km/s)
+ state_rk8(6)%r, & ! <--> DBL() radius vector (km)
+ state_rk8(6)%v, & ! <--> DBL() velocity vector (km/s)
reset, & ! <-- INT reset flag
dt & ! --> DBL propagated time (s)
)
@@ -1242,8 +1386,8 @@ subroutine getStateTransitionMatrix( this, &
t7, & ! <-- DBL requested time (s)
.false., &
reqt_loc, & ! <-- DBL current time (s)
- state_rk8(7)%r, & ! <--> DBL() radius vector (km)
- state_rk8(7)%v, & ! <--> DBL() velocity vector (km/s)
+ state_rk8(7)%r, & ! <--> DBL() radius vector (km)
+ state_rk8(7)%v, & ! <--> DBL() velocity vector (km/s)
reset, & ! <-- INT reset flag
dt & ! --> DBL propagated time (s)
)
@@ -1271,8 +1415,8 @@ subroutine getStateTransitionMatrix( this, &
t9, & ! <-- DBL requested time (s)
.false., &
reqt_loc, & ! <-- DBL current time (s)
- state_rk8(9)%r, & ! <--> DBL() radius vector (km)
- state_rk8(9)%v, & ! <--> DBL() velocity vector (km/s)
+ state_rk8(9)%r, & ! <--> DBL() radius vector (km)
+ state_rk8(9)%v, & ! <--> DBL() velocity vector (km/s)
reset, & ! <-- INT reset flag
dt & ! --> DBL propagated time (s)
)
@@ -1291,8 +1435,8 @@ subroutine getStateTransitionMatrix( this, &
solarsystem_model, &
thirdbody_model, &
reduction, &
- state_rk8(1)%r, &
- state_rk8(1)%v, &
+ state_rk8(1)%r, &
+ state_rk8(1)%v, &
set, &
t1d, &
k1)
@@ -1421,7 +1565,120 @@ subroutine getStateTransitionMatrix( this, &
! RK8:
set = set + (this%covIntegrationStep/840.d0)*(41.d0*k1 + 27.d0*k4 + 272.d0*k5 + 27.d0*k6 + 216.d0*k7 + 216.d0*k9 + 41.d0*k10)
+ sensitivity_matrix = 0.
+
+ if(this%cd_cov_flag) then
+
+ call reduction%inertial2earthFixed(state_rk8(1)%r, state_rk8(1)%v, t1d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state_rk8(1)%r, &
+ state_rk8(1)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t1d, & ! <-- DBL current time (MJD)
+ k1_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call reduction%inertial2earthFixed(state_rk8(4)%r, state_rk8(4)%v, t4d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state_rk8(4)%r, &
+ state_rk8(4)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t4d, & ! <-- DBL current time (MJD)
+ k4_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call reduction%inertial2earthFixed(state_rk8(5)%r, state_rk8(5)%v, t5d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state_rk8(5)%r, &
+ state_rk8(5)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t5d, & ! <-- DBL current time (MJD)
+ k5_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call reduction%inertial2earthFixed(state_rk8(6)%r, state_rk8(6)%v, t6d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state_rk8(6)%r, &
+ state_rk8(6)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t6d, & ! <-- DBL current time (MJD)
+ k6_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+ call reduction%inertial2earthFixed(state_rk8(7)%r, state_rk8(7)%v, t7d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state_rk8(7)%r, &
+ state_rk8(7)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t7d, & ! <-- DBL current time (MJD)
+ k7_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call reduction%inertial2earthFixed(state_rk8(9)%r, state_rk8(9)%v, t9d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state_rk8(9)%r, &
+ state_rk8(9)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t9d, & ! <-- DBL current time (MJD)
+ k9_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+ call reduction%inertial2earthFixed(state_rk8(10)%r, state_rk8(10)%v, t10d, r_itrf, v_itrf)
+ call atmosphere_model%getAtmosphereAcceleration( &
+ gravity_model, & ! <-> TYPE gravity model
+ satellite_model, & ! <-> TYPE satellite model
+ solarsystem_model, & ! <-> TYPE solarsystem model
+ reduction, &
+ state_rk8(10)%r, &
+ state_rk8(10)%v, &
+ r_itrf, & ! <-- DBL(3) radius vector in ITRF
+ v_itrf, & ! <-- DBL(3) velocity vector in ITRF
+ t10d, & ! <-- DBL current time (MJD)
+ k10_atm & ! --> DBL(3) acceleration vector in inertial frame
+ )
+ if(hasFailed()) return
+
+
+ sensitivity_matrix(4:6) = (this%covIntegrationStep/840.d0)*(41.d0*k1_atm + 27.d0*k4_atm + 272.d0*k5_atm + 27.d0*k6_atm + 216.d0*k7_atm + 216.d0*k9_atm + 41.d0*k10_atm)
+ sensitivity_matrix(7) = 1.d0
+
+ end if
!** save last state
this%lastState = state_rk8(8)
@@ -1504,6 +1761,46 @@ real(dp) function get_current_step_size(this)
return
end function
+!========================================================================
+!
+!> @anchor getCdCov
+!!
+!> @brief Returns the step size of the integrator
+!> @author Christopher Kebschull
+!!
+!> @date
+!! - 29.10.2020 (initial design)
+!!
+!!
+!------------------------------------------------------------------------
+real(dp) function getCdCov(this)
+ class(Numint_class),intent(inout) :: this
+
+ getCdCov = this%cd_cov
+
+ return
+end function
+
+!========================================================================
+!
+!> @anchor getCdCovFlag
+!!
+!> @brief Returns the step size of the integrator
+!> @author Christopher Kebschull
+!!
+!> @date
+!! - 29.10.2020 (initial design)
+!!
+!!
+!------------------------------------------------------------------------
+logical function getCdCovFlag(this)
+ class(Numint_class),intent(inout) :: this
+
+ getCdCovFlag = this%cd_cov_flag
+
+ return
+end function
+
!========================================================================
!
diff --git a/src/rdinp.f90 b/src/rdinp.f90
index 554e0db..ba9f199 100644
--- a/src/rdinp.f90
+++ b/src/rdinp.f90
@@ -352,6 +352,10 @@ subroutine rdinp( &
covar%frame = REF_FRAME_GCRF
end if
+ !** Drag sigma
+ call nxtbuf('#', 0, ich_inp, cbuf)
+ read(cbuf,*) ctemp
+ ierr = neptune%setNeptuneVar("PAR_CDRAG_COV", ctemp)
!========================================================
!
! perturbation switches
@@ -599,6 +603,16 @@ subroutine rdinp( &
end if
ierr = neptune%setNeptuneVar("CORRELATION_MATRIX", ctemp)
+ !** correlation matrix computation switch
+ call nxtbuf('#', 0, ich_inp, cbuf)
+ read(cbuf,*) itemp
+ if(itemp == 0) then
+ ctemp = "OFF"
+ else
+ ctemp = "ON"
+ end if
+ ierr = neptune%setNeptuneVar("CONSIDER_CD_COV", ctemp)
+
!** further options
!----------------------------
diff --git a/work/input/neptune.inp b/work/input/neptune.inp
index b509f96..60bcf57 100644
--- a/work/input/neptune.inp
+++ b/work/input/neptune.inp
@@ -94,6 +94,8 @@
0.d0 0.d0 0.d0 covariance matrix 4th row
0.d0 0.d0 0.d0 0.d0 covariance matrix 5th row
0.d0 0.d0 0.d0 0.d0 0.d0 covariance matrix 6th row
+#
+ 0.1d0 c_d covariance
#
#------------------------------------------------------------------------------
#
@@ -158,6 +160,13 @@
#
0
#
+# Consider C_D Covariance
+#
+# 1 = on (only considered if covariance matrix propagation is on)
+# 0 = off
+#
+ 0
+#
#------------------------------------------------------------------------------
#
# Further options