Skip to content

Commit

Permalink
Merge pull request #37 from Space-Systems/chk/gravity_test_issue
Browse files Browse the repository at this point in the history 
Improved maneuver handling and fixes
  • Loading branch information
@ckebschull
ckebschull authored Feb 13, 2024
2 parents 8104164 + 31963da commit 8dbe5f5
Show file tree
Hide file tree
Showing 21 changed files with 6,191 additions and 5,777 deletions.
18 changes: 9 additions & 9 deletions build.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,14 +20,14 @@ git submodule update --init --recursive
# #
################################################################################
# Build OPI
cd OPI || exit
cd OPI || exit || exit
# Create the build directory if it does not exist
if [[ ! -d "build" ]]; then
mkdir build
else
rm -rf build/*
fi
cd build || exit
cd build || exit || exit
echo "Updating cmake"
cmake -DCMAKE_BUILD_TYPE=$BUILD_TYPE -DCMAKE_INSTALL_PREFIX=../ -DCMAKE_Fortran_COMPILER=$Fortran_COMPILER -DENABLE_CXX11=ON -DENABLE_CL_SUPPORT=OFF -DENABLE_CUDA_SUPPORT=OFF -DENABLE_PYTHON=OFF ../
echo "Building OPI"
Expand All @@ -44,14 +44,14 @@ cd ../../
# #
################################################################################
# Build libslam
cd libslam || exit
cd libslam || exit || exit
# Create the build directory if it does not exist
if [[ ! -d "build" ]]; then
mkdir build
else
rm -rf build/*
fi
cd build || exit
cd build || exit || exit
echo "Updating cmake"
cmake -DCMAKE_BUILD_TYPE=$BUILD_TYPE -DCMAKE_Fortran_COMPILER=$Fortran_COMPILER -DENABLE_OpenMP_SUPPORT=OFF -DENABLE_POSTGRESQL_SUPPORT=OFF ../
echo "Building libslam"
Expand Down Expand Up @@ -98,7 +98,7 @@ if [[ ! -d "lib" ]]; then
mkdir lib
fi
# Create the links to libraries needed by CAMP
cd lib || exit
cd lib || exit || exit
ln -sf ../libslam/lib/libslam-Fortran.$LIBSUFFIX .
ln -sf ../OPI/lib/libOPI-Fortran.$LIBSUFFIX .
ln -sf ../OPI/lib/libOPI.$LIBSUFFIX .
Expand All @@ -108,7 +108,7 @@ if [[ ! -d "include" ]]; then
mkdir include
fi
# Create the links to includes needed by NEPTUNE
cd include || exit
cd include || exit || exit
ln -sf ../libslam/include/SLAM .
ln -sf ../OPI/include/OPI .
cd ..
Expand All @@ -118,7 +118,7 @@ if [[ ! -d "build" ]]; then
else
rm -rf build/*
fi
cd build || exit
cd build || exit || exit
echo "Updating cmake"
export PFUNIT_DIR=..//pFUnit/build/installed
cmake -DCMAKE_BUILD_TYPE=$BUILD_TYPE -DCMAKE_Fortran_COMPILER=$Fortran_COMPILER -DENABLE_OPI_SUPPORT=ON -DSKIP_MSIS_2=ON ../
Expand All @@ -129,7 +129,7 @@ if [[ $? -ne 0 ]]; then
exit $?
fi
echo "Leaving NEPTUNE"
cd ../work || exit
ln -sf ../bin/neptune-sa .
cd ../work || exit || exit
ln -sf ../bin/neptune-sa . .
cd ..
echo "Done"
3 changes: 2 additions & 1 deletion src/gravity.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ module gravity
use slam_error_handling, only: isControlled, hasToReturn, hasFailed, FATAL, WARNING, &
E_SPECIAL, E_UNKNOWN_PARAMETER, E_OUT_OF_BOUNDS, checkIn, checkOut
use slam_io, only: openFile, closeFile, SEQUENTIAL, IN_FORMATTED, LOG_AND_STDOUT, cdelimit, message
use slam_strings, only: toString
use slam_math, only: UNDEFINED, identity_matrix, rad2deg, redang, eps9, mag, outerproduct, uvec1, uvec2, uvec3
use slam_reduction_class, only: Reduction_type
use slam_time, only: time_t, gd2jd
Expand Down Expand Up @@ -402,7 +403,7 @@ subroutine initGravityPotential( this, &

ich = openFile(fileName, SEQUENTIAL, IN_FORMATTED)

call message(' - Reading '//trim(modelName(this%nmodel))//' geopotential coefficients...', LOG_AND_STDOUT)
call message(' - Reading '//trim(modelName(this%nmodel))//' geopotential coefficients with degree '//toString(this%degree)//'...', LOG_AND_STDOUT)

!====================================================================================================
!
Expand Down
55 changes: 42 additions & 13 deletions src/maneuvers.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ module maneuvers
use slam_math, only: mag, eps9
use slam_reduction_class, only: Reduction_type
use slam_time, only: time_t, tokenizeDate, sec_per_day
use slam_strings, only: toString

implicit none

Expand Down Expand Up @@ -465,6 +466,8 @@ subroutine init_maneuvers_file( &
return
end if

call message(' - Initializing maneuvers from file...', LOG_AND_STDOUT)

!** check input
if(date_start%mjd > date_end%mjd) then
call setNeptuneError(E_SGA_INPUT, FATAL)
Expand Down Expand Up @@ -543,6 +546,8 @@ subroutine init_maneuvers_api( this, &
call checkIn(csubid)
end if

call message(' - Initializing maneuvers from API...', LOG_AND_STDOUT)

! Please not that we allow re-initialization

!=====================================================================
Expand Down Expand Up @@ -599,6 +604,8 @@ subroutine init_maneuver_sequence(this, number_phases)

allocate(this%mnv_sequence(size(this%manv))) ! enough elements to account for transitions

call message(" - Initializing maneuver sequences", LOG_AND_STDOUT)

k = 1
do i=1,size(this%manv)

Expand All @@ -609,6 +616,12 @@ subroutine init_maneuver_sequence(this, number_phases)
this%mnv_sequence(k)%mjd_end = this%manv(i)%phase(j)%mjd_end
this%mnv_sequence(k)%trans = NO_TRANSITION

call message(" ... from "//toString(this%mnv_sequence(k)%mjd_start)//" to "//toString(this%mnv_sequence(k)%mjd_end) &
//" with " &
//toString(this%manv(i)%phase(j)%acc(1))//" km/s² x " &
//toString(this%manv(i)%phase(j)%acc(2))//" km/s² x " &
//toString(this%manv(i)%phase(j)%acc(3))//" km/s²", LOG_AND_STDOUT)

k = k + 1 ! +1, as between each maneuver phase, also a transition phase will be added later on (see below)

end do
Expand Down Expand Up @@ -1116,34 +1129,50 @@ end function get_current_index
!! @brief Get the start or end epoch of the next manoeuvre within the mnv_sequence array
!! @author Christopher Kebschull
!!
!! @param[in] mjd current MJD
!! @param[in] mjd current MJD
!! @param[in] using_backwards_propagation Whether propagation is backwards or not. Optional. Defaults to .false.
!!
!! @date <ul>
!! <li> 13.07.2020: initial design</li>
!! </ul>
!!
!-----------------------------------------------------------------------------
real(dp) function get_upcoming_manoeuvre_change_epoch(this, mjd)
real(dp) function get_upcoming_manoeuvre_change_epoch(this, mjd, using_backwards_propagation)

class(Manoeuvres_class) :: this
real(dp), intent(in) :: mjd ! current MJD
logical, intent(in), optional :: using_backwards_propagation

integer :: k ! loop counter
logical :: flag_backwards

flag_backwards = .false.
if (present(using_backwards_propagation)) flag_backwards = using_backwards_propagation

get_upcoming_manoeuvre_change_epoch = 0.d0

! loop through mnv_sequence array to bracket mjd
do k = 1, size(this%mnv_sequence)

if( mjd < this%mnv_sequence(k)%mjd_start) then
get_upcoming_manoeuvre_change_epoch = this%mnv_sequence(k)%mjd_start
return
else if(mjd < this%mnv_sequence(k)%mjd_end) then
get_upcoming_manoeuvre_change_epoch = this%mnv_sequence(k)%mjd_end
return
end if

end do
if (.not. flag_backwards) then ! forward
do k = 1, size(this%mnv_sequence)
if( mjd < this%mnv_sequence(k)%mjd_start) then
get_upcoming_manoeuvre_change_epoch = this%mnv_sequence(k)%mjd_start
return
else if(mjd < this%mnv_sequence(k)%mjd_end) then
get_upcoming_manoeuvre_change_epoch = this%mnv_sequence(k)%mjd_end
return
end if
end do
else ! backward
do k = size(this%mnv_sequence), 1, -1
if( mjd > this%mnv_sequence(k)%mjd_end) then
get_upcoming_manoeuvre_change_epoch = this%mnv_sequence(k)%mjd_end
return
else if(mjd > this%mnv_sequence(k)%mjd_start) then
get_upcoming_manoeuvre_change_epoch = this%mnv_sequence(k)%mjd_start
return
end if
end do
end if

end function

Expand Down
44 changes: 32 additions & 12 deletions src/neptune.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ module libneptune
PAR_MASS, PROP_CONSTANT, PROP_FILE, C_INITIAL_STATE, C_PAR_EARTH_RADIUS
use numint, only: MAX_RESETS
use satellite, only: ORIENT_SPHERE
use maneuvers, only: neptune_maneuver_t => maneuver_t
use slam_strings, only: toString
use slam_time, only: time_t, gd2mjd, date2longstring, operator(>)
use slam_types, only: dp
Expand All @@ -71,6 +72,7 @@ module libneptune
public :: init_neptune ! initialization
public :: propagate ! propagation of state vector and covariance
public :: propagate_set ! propagation of state vector, covariance and state error transition matrix
public :: neptune_maneuver_t ! Maneuver type

contains

Expand Down Expand Up @@ -192,7 +194,6 @@ integer function init_neptune(neptune,initial_state_in)
neptune%satellite_model%surface(1)%reflSpec = neptune%ipar(PAR_CREFL) - 1.d0 ! the specular reflectivity coefficient is saved,
! while the SRP coefficient is used in a cannon-ball model
neptune%satellite_model%surface(1)%reflDiff = 0.d0 ! no diffuse reflectivity in this case

else
cmess = "SRP coefficient not defined."
call setNeptuneError(E_INVALID_SATELLITE, FATAL, (/cmess/))
Expand All @@ -201,7 +202,7 @@ integer function init_neptune(neptune,initial_state_in)

!** Area
if(neptune%isSetArea) then
neptune%satellite_model%surface(1)%area = neptune%ipar(PAR_CROSS_SECTION) ! in km**2
neptune%satellite_model%surface(1)%area = neptune%ipar(PAR_CROSS_SECTION) ! in m**2
else
cmess = "Cross-section area not defined."
call setNeptuneError(E_INVALID_SATELLITE, FATAL, (/cmess/))
Expand Down Expand Up @@ -230,6 +231,14 @@ integer function init_neptune(neptune,initial_state_in)
if(present(initial_state_in)) then
initial_state = initial_state_in
ierr = neptune%setNeptuneVar("INITIAL_STATE", initial_state)
call message(" Using initial state "//date2longstring(initial_state%epoch)// &
" x: "//toString(initial_state%r(1))// &
" y: "//toString(initial_state%r(2))// &
" z: "//toString(initial_state%r(3))// &
" x_dot: "//toString(initial_state%v(1))// &
" y_dot: "//toString(initial_state%v(2))// &
" z_dot: "//toString(initial_state%v(3)), &
LOG_AND_STDOUT)
end if

!** check for start and end epoch being available
Expand Down Expand Up @@ -459,7 +468,7 @@ subroutine propagate( &
type(Neptune_class) ,intent(inout) :: neptune
type(state_t), intent(in) :: state_in
type(covariance_t), intent(in) :: covar_in
type(time_t), dimension(2), intent(in) :: epoch
type(time_t), dimension(:), intent(in) :: epoch
logical, intent(in) :: flag_reset

type(state_t), intent(out) :: state_out
Expand Down Expand Up @@ -584,6 +593,7 @@ subroutine propagate_set( &

restored_request_time = 0.d0
prop_counter = 0.0d0
manoeuvre_change_counter = 0.0d0
propCounterAtReset = 0.0d0
lastPropCounterSuccess = 0.0d0
diff = 0.d0
Expand Down Expand Up @@ -703,8 +713,10 @@ subroutine propagate_set( &
!---------------------------------------------
if(epochs(2)%mjd < epochs(1)%mjd) then
flag_backward = .true.
! call message(" - Set to backward propagation", LOG_AND_STDOUT)
else
flag_backward = .false.
! call message(" - Set to forward propagation", LOG_AND_STDOUT)
end if

!============================================================
Expand Down Expand Up @@ -826,6 +838,7 @@ subroutine propagate_set( &
else
! No previous intermediate call - we can just proceed as planned
request_time = nep_clock%get_next_step(neptune,prop_counter)
! call message(' - Next request time: '//toString(epochs(1)%mjd + request_time/86400.d0)//'('//toString(request_time)//')', LOGFILE)

! Save the current requested time in case we determine that an intermediate call must be done due to an immenant manoeuvre start or end
restored_request_time = request_time
Expand All @@ -841,25 +854,26 @@ subroutine propagate_set( &
intermediate_integrator_call = .false.
if (neptune%derivatives_model%getPertSwitch(PERT_MANEUVERS)) then

!call message(' Maneuver change at '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//')', LOG_AND_STDOUT)
! call message(' Maneuver change at '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//')', LOGFILE)

! Reset when we are starting into the new maneuver or no-maneuver interval
if (abs(prop_counter - manoeuvre_change_counter) < epsilon(1.d0)) then
reset = 1
!force_no_interpolation = .true.
! reset also the count of subroutine calls to the integrator
call neptune%numerical_integrator%resetCountIntegrator()
call message(' - Performing reset of intergator now '//toString(epochs(1)%mjd + prop_counter/86400.d0)//'('//toString(prop_counter)//') for upcoming maneuver '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//')', LOGFILE)
call message(' - Performing reset of integrator now '//toString(epochs(1)%mjd + prop_counter/86400.d0)//'('//toString(prop_counter)//') for upcoming maneuver '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//')', LOGFILE)
call message(' ... state vector going in: '//toString(state_out%r(1))//', '//toString(state_out%r(2))//', '//toString(state_out%r(3))//', '//toString(state_out%v(1))//', '//toString(state_out%v(2))//', '//toString(state_out%v(3)), LOGFILE)
end if

! Get the next manoeuvre change epoch
upcoming_maneuver_epoch_mjd = neptune%manoeuvres_model%get_upcoming_manoeuvre_change_epoch(epochs(1)%mjd + prop_counter/86400.d0)
upcoming_maneuver_epoch_mjd = neptune%manoeuvres_model%get_upcoming_manoeuvre_change_epoch(epochs(1)%mjd + prop_counter/86400.d0, using_backwards_propagation=flag_backward)
if (upcoming_maneuver_epoch_mjd > 0.d0) then
manoeuvre_change_counter = (upcoming_maneuver_epoch_mjd - epochs(1)%mjd) * 86400.d0
!call message(' New maneuver change at: '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//')', LOG_AND_STDOUT)
! call message(' New maneuver change at: '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//')', LOGFILE)
! Check whether the manoeuvre change is more immenent than the step proposed
if (.not. flag_backward .and. (manoeuvre_change_counter < request_time .and. manoeuvre_change_counter > prop_counter) &
.or. flag_backward .and. (manoeuvre_change_counter > request_time .and. manoeuvre_change_counter < prop_counter)) then
if (.not. flag_backward .and. (manoeuvre_change_counter <= request_time .and. manoeuvre_change_counter > prop_counter) &
.or. flag_backward .and. (manoeuvre_change_counter >= request_time .and. manoeuvre_change_counter < prop_counter)) then
if (abs(manoeuvre_change_counter - request_time) > epsilon(1.d0)) then
! The manoeuvre_change_counter does not coincide with a storage and output epoch - surpress the output!
suppressed_output = .true.
Expand Down Expand Up @@ -895,7 +909,7 @@ subroutine propagate_set( &
force_no_interpolation = .true.
! reset also the count of subroutine calls to the integrator
!call neptune%numerical_integrator%resetCountIntegrator()
call message(' - Forcing no interpolation now '//toString(epochs(1)%mjd + prop_counter/86400.d0)//'('//toString(prop_counter)//') for upcoming maneuver due to time step size '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//') int step size: '//toString(neptune%numerical_integrator%get_current_step_size()), LOGFILE)
! call message(' - Enable no interpolation now '//toString(epochs(1)%mjd + prop_counter/86400.d0)//'('//toString(prop_counter)//') for upcoming maneuver due to time step size '//toString(epochs(1)%mjd + manoeuvre_change_counter/86400.d0)//' ('//toString(manoeuvre_change_counter)//') int step size: '//toString(neptune%numerical_integrator%get_current_step_size()), LOGFILE)
end if
end if
end if
Expand Down Expand Up @@ -932,7 +946,11 @@ subroutine propagate_set( &
if(neptune%has_to_write_progress()) then
dtmp = abs(prop_counter - start_epoch_sec)/abs(end_epoch_sec - start_epoch_sec)
call write_progress(neptune%get_progress_file_name(), dtmp, neptune%get_progress_step())
! write(*,*) dtmp, prop_counter
end if

! dtmp = abs(prop_counter - start_epoch_sec)/abs(end_epoch_sec - start_epoch_sec)
! write(*,"(A,F7.2,A,F15.4,A,F15.6)") "progress =",100*dtmp, " // prop_counter =", prop_counter," // mjd =",epochs(1)%mjd + prop_counter/86400.d0

if(hasFailed()) then
call neptune%output%close_open_files(neptune%numerical_integrator)
Expand Down Expand Up @@ -968,6 +986,8 @@ subroutine propagate_set( &
state_out = last_state_out
prop_counter = lastPropCounter
end if
call message(' Resetting to: '//toString(epochs(1)%mjd + prop_counter/86400.d0)//'('//toString(prop_counter), LOGFILE)
call message(' ... restoring state vector: '//toString(state_out%r(1))//', '//toString(state_out%r(2))//', '//toString(state_out%r(3))//', '//toString(state_out%v(1))//', '//toString(state_out%v(2))//', '//toString(state_out%v(3)), LOGFILE)
end if
else if((.not. flag_backward .and. (prop_counter > propCounterAtReset)) .or. &
( flag_backward .and. (prop_counter < propCounterAtReset))) then ! reset 'nresets' flag
Expand All @@ -976,8 +996,8 @@ subroutine propagate_set( &
end if


! diff = abs(prop_counter - request_time)
! write(cmess, '(a, D15.3, a)') 'diff = ', diff, toString(intermediate_integrator_call)
! diff = prop_counter - request_time
! write(cmess, '(a, D15.6, a, a)') 'diff = ', diff, " // intermediate_integrator_call=", toString(intermediate_integrator_call)
! call message(cmess, LOG_AND_STDOUT)

! Exit the integration loop when the desired time is reached
Expand Down
Loading

0 comments on commit 8dbe5f5

Please sign in to comment.