read_config_file Subroutine

public subroutine read_config_file(me, filename)

Read the config file that defines the problem to be solved and set all the global variables.

Notes

  • "period" is the normalized period ("jc", the jacobi constant, is also allowed)
  • The states in the patch point file are assumed to be "S" family

Type Bound

my_solver_type

Arguments

Type IntentOptional Attributes Name
class(my_solver_type), intent(inout) :: me
character(len=*), intent(in) :: filename

the JSON config file to read

Calls

proc~~read_config_file~~CallsGraph proc~read_config_file my_solver_type%read_config_file evaluate evaluate proc~read_config_file->evaluate initialize initialize proc~read_config_file->initialize none~get config_file%get proc~read_config_file->none~get proc~generate_patch_points mission_type%generate_patch_points proc~read_config_file->proc~generate_patch_points proc~read_epoch read_epoch proc~read_config_file->proc~read_epoch proc~update_epoch mission_type%update_epoch proc~read_config_file->proc~update_epoch proc~get_char config_file%get_char none~get->proc~get_char proc~get_int config_file%get_int none~get->proc~get_int proc~get_logical config_file%get_logical none~get->proc~get_logical proc~get_real config_file%get_real none~get->proc~get_real proc~get_real_vec config_file%get_real_vec none~get->proc~get_real_vec proc~generate_patch_points->initialize compute_crtpb_parameter compute_crtpb_parameter proc~generate_patch_points->compute_crtpb_parameter first_call first_call proc~generate_patch_points->first_call integrate integrate proc~generate_patch_points->integrate proc~read_epoch->none~get jd_to_et jd_to_et proc~read_epoch->jd_to_et julian_date julian_date proc~read_epoch->julian_date et_to_jd et_to_jd proc~update_epoch->et_to_jd proc~update_epoch->jd_to_et proc~update_epoch->julian_date julian_date_to_calendar_date julian_date_to_calendar_date proc~update_epoch->julian_date_to_calendar_date get get proc~get_char->get proc~get_int->get proc~get_logical->get proc~get_real->get proc~get_real_vec->get

Called by

proc~~read_config_file~~CalledByGraph proc~read_config_file my_solver_type%read_config_file proc~initialize_the_solver my_solver_type%initialize_the_solver proc~initialize_the_solver->proc~read_config_file proc~halo_solver_main halo_solver_main proc~halo_solver_main->proc~initialize_the_solver

Source Code

    subroutine read_config_file(me, filename)

    implicit none

    class(my_solver_type),intent(inout) :: me
    character(len=*),intent(in) :: filename  !! the JSON config file to read

    type(config_file) :: f
    logical :: found, found_jc, found_period, found_initial_r
    real(wp):: jc !! jacobii constant from the file
    real(wp):: p !! period from file (normalized)
    type(patch_point),dimension(3) :: pp !! patch points
    logical :: use_json_file  !! if true, use the JSON file.
                              !! if false, use the CSV file.
    real(wp) :: lstar, tstar
    real(wp),dimension(:),allocatable :: jcvec, normalized_period
    real(wp),dimension(:),allocatable :: x0, z0, ydot0
    real(wp) :: pp_period, pp_x0, pp_z0, pp_ydot0

    write(*,*) '* Loading config file: '//trim(filename)
    call f%open(filename)

    write(*,*) '* Reading config file'

    call f%get('jc',                        jc, found_jc)          ! one of these two must be present
    call f%get('period',                    p, found_period)       !
    call f%get('solve',                             me%mission%solve,                     found)
    call f%get('fix_initial_time',                  me%mission%fix_initial_time,          found)
    call f%get('fix_initial_r',                     me%mission%fix_initial_r,             found)
    call f%get('fix_ry_at_end_of_rev',              me%mission%fix_ry_at_end_of_rev,      found)
    call f%get('fix_final_ry_and_vx',               me%mission%fix_final_ry_and_vx,       found)
    call f%get('generate_plots',                    me%mission%generate_plots,            found)
    call f%get('generate_trajectory_files',         me%mission%generate_trajectory_files, found)
    call f%get('generate_json_trajectory_file',     me%mission%generate_json_trajectory_file,     found)
    call f%get('generate_guess_and_solution_files', me%mission%generate_guess_and_solution_files, found)
    call f%get('generate_kernel',                   me%mission%generate_kernel,           found)
    call f%get('generate_defect_file',              me%mission%generate_defect_file,      found)
    call f%get('generate_eclipse_files',            me%mission%generate_eclipse_files,    found)
    call f%get('run_pyvista_script',                me%mission%run_pyvista_script,        found)

    call f%get('r_eclipse_bubble',    me%mission%r_eclipse_bubble, found)
    call f%get('eclipse_dt_step',     me%mission%eclipse_dt_step,  found)
    call f%get('eclipse_filetype',    me%mission%eclipse_filetype, found)

    call f%get('use_splined_ephemeris',  me%mission%use_splined_ephemeris,  found)
    call f%get('dt_spline_sec',          me%mission%dt_spline_sec,          found)

    call f%get('pointmass_central_body',  me%mission%pointmass_central_body,  found)

    call f%get('initial_guess_from_file', me%mission%initial_guess_from_file, found)

    call f%get('solver_mode', me%mission%solver_mode, found)
    if (.not. found) me%mission%solver_mode = 1

    ! required inputs:
    call f%get('N_or_S',   me%mission%N_or_S   )
    call f%get('L1_or_L2', me%mission%L1_or_L2 )

    ! reach the epoch. can specify either calendar date or et:
    call read_epoch(f,  me%mission%epoch_mode, &
                        me%mission%year, &
                        me%mission%month, &
                        me%mission%day, &
                        me%mission%hour, &
                        me%mission%minute, &
                        me%mission%sec, &
                        me%mission%et_ref)

    ! tolerances:
    call f%get('rtol',          me%mission%rtol          , found )
    call f%get('atol',          me%mission%atol          , found )
    call f%get('nlesolver_tol', me%mission%nlesolver_tol , found )

    ! if fix_initial_r is true, can specify the r to use
    ! [otherwise, the initial guess is used from the patch points]
    call f%get('initial_r', me%mission%initial_r , found_initial_r )

    call f%get('n_revs',          me%mission%n_revs           )
    call f%get('ephemeris_file',  me%mission%ephemeris_file   )
    call f%get('gravfile',        me%mission%gravfile         )
    call f%get('patch_point_file',me%mission%patch_point_file )
    call f%get('patch_point_file_is_periapsis',me%mission%patch_point_file_is_periapsis, found )
    call f%close() ! cleanup

    use_json_file = index(me%mission%patch_point_file, '.json') > 0
    if (.not. use_json_file) error stop 'error: patch point file must be a JSON file.'
    if (use_json_file .and. (found_jc .eqv. found_period)) &
        error stop 'error: must either specify period or jc in the config file to use JSON patch point file.'
        ! csv file no longer supported

    ! set the epoch:
    call me%mission%update_epoch()

    ! read the patch point file and load it:

    ! read the JSON file:
    write(*,*) '* Reading JSON patch point file: '//trim(me%mission%patch_point_file)
    call f%open(me%mission%patch_point_file)

    ! all inputs are required:
    write(*,*) '* getting data... '
    call f%get('lstar', lstar            )
    call f%get('tstar', tstar            )
    call f%get('jc',    jcvec            )
    call f%get('period',normalized_period)
    call f%get('x0',    x0               )
    call f%get('z0',    z0               )
    call f%get('ydot0', ydot0            )
    call f%close()

    ! write(*,'(A, *(F6.2,",",1X))') 'period range in database (days): ', &
    !       normalized_period * tstar * sec2day

    ! which is the independant variable:
    if (found_jc) then
        pp_period = interpolate_point(jcvec, normalized_period, jc)
        pp_x0     = interpolate_point(jcvec, x0, jc)
        pp_z0     = interpolate_point(jcvec, z0, jc)
        pp_ydot0  = interpolate_point(jcvec, ydot0, jc)
    else
        pp_period = p
        pp_x0     = interpolate_point(normalized_period, x0, p)
        pp_z0     = interpolate_point(normalized_period, z0, p)
        pp_ydot0  = interpolate_point(normalized_period, ydot0, p)
    end if

    call me%mission%generate_patch_points(lstar, tstar, pp_period, &
                                          pp_x0, pp_z0, pp_ydot0, pp)
    me%mission%periapsis = pp(1)
    me%mission%quarter   = pp(2)
    me%mission%apoapsis  = pp(3)

    ! write(*,'(a15,1x,F10.6,1x,a)') 'periapsis t: ' , me%mission%periapsis%t, 'days'
    ! write(*,'(a15,1x,F10.6,1x,a)') 'quarter t:   ' , me%mission%quarter%t,   'days'
    ! write(*,'(a15,1x,F10.6,1x,a)') 'apoapsis t:  ' , me%mission%apoapsis%t,  'days'

    ! compute some time variables:
    me%mission%period  = me%mission%apoapsis%t * 2.0_wp  ! Halo period [days]
    me%mission%period8 = me%mission%period / 8.0_wp      ! 1/8 of Halo period [days]

    write(*,*) ''
    write(*,'(A30,1x,f16.6,a,i4,a,i2,a,i2,a,i0.2,a,i0.2,a,f9.6,a)') '   Reference ephemeris time: ', &
                                                                          me%mission%et_ref, &
                                                                    ' (', me%mission%year,'-',&
                                                                          me%mission%month,'-',&
                                                                          me%mission%day,' ',&
                                                                          me%mission%hour,':',&
                                                                          me%mission%minute,':',&
                                                                          me%mission%sec, ')'
    write(*,'(a30,1x,f30.17)') '   Orbit period (days): ', me%mission%period

    contains
!*****************************************************************************************

    !**********************************************
    !>
    !  Iterpolate if necessary to compute `fvec(x)`.
    !  We will allow data sets that are not strictly increasing/decreasing
    !  (`jc` isn't, but `period` is for L2 file),
    !  so, first the input `x` is checked to see if that is a point in `xvec`.
    !  If not, the bspline interpolation is used.

        function interpolate_point(xvec, fvec, x) result(y)

        implicit none

        real(wp),dimension(:),intent(in) :: xvec
        real(wp),dimension(:),intent(in) :: fvec
        real(wp),intent(in) :: x
        real(wp) :: y !! `fvec(x)`

        type(bspline_1d) :: spline
        integer :: i, iflag, irow

        integer,parameter :: kx = 4  !! spline order (cubic bspline)
        integer,parameter :: idx = 0 !! interpolate value only

        ! so, first just look for the exact point:
        irow = 0
        do i=1,size(xvec)
            if (xvec(i)==x) then
                irow = i
                exit
            end if
        end do

        if (irow/=0) then
            ! if found, then use it:
            y = fvec(irow)
        else
            call spline%initialize(xvec,fvec,kx,iflag, extrap=.true.)
            if (iflag/=0) then
                select case (iflag)
                case(2); error stop 'bspline error: iknot out of range.'
                case(3); error stop 'bspline error: nx out of range.'
                case(4); error stop 'bspline error: kx out of range.'
                case(5); error stop 'bspline error: x not strictly increasing.'
                end select
            end if
            call spline%evaluate(x,idx,y,iflag)
            if (iflag/=0) error stop 'error evaluating bspline.'
        end if

        end function interpolate_point
    !**********************************************

    end subroutine read_config_file