initialize_the_mission Subroutine

    subroutine initialize_the_mission(me,x)

    implicit none

    class(mission_type),intent(inout) :: me
    real(wp),dimension(:),allocatable,intent(out),optional :: x !! initial guess

    logical :: status_ok                        !! status flag
    integer :: i,j                              !! counter
    integer :: n_segs                           !! number of segments
    integer :: iseg                             !! segment counter
    real(wp) :: t_periapsis                     !! time of last periapsis crossing (in days from ref epoch)
    integer :: irev                             !! rev counter
    real(wp),dimension(8) :: t0                 !! [days]
    real(wp),dimension(8) :: tf                 !! [days]
    real(wp),dimension(6,8) :: x0_rotating      !! rotating frame
    integer :: n                                !! number of opt vars
    integer :: m                                !! number of constraints
    integer :: n_nonzero                        !! number of nonzero elements in the jacobian
    real(wp),dimension(:),allocatable :: dpert  !! perturbation step size array
    real(wp),dimension(:),allocatable :: xlow   !! lower bounds array
    real(wp),dimension(:),allocatable :: xhigh  !! upper bounds
    character(len=10) :: seg_name               !! segment name
    integer,dimension(:),allocatable :: irow    !! sparsity pattern nonzero elements row indices
    integer,dimension(:),allocatable :: icol    !! sparsity pattern nonzero elements column indices
    logical :: use_openmp !! if OpenMP is being used
    real(wp) :: et0 !! initial et for splined ephemeris
    real(wp) :: etf !! final et for splined ephemeirs

    use_openmp = .false.
!$  use_openmp = .true.

    ! problem size:
    call me%define_problem_size(n=n,m=m,n_segs=n_segs,n_nonzero=n_nonzero)

    write(*,*) ''
    write(*,'(A,1X,I10)') '    n:                                     ', n
    write(*,'(A,1X,I10)') '    m:                                     ', m
    write(*,'(A,1X,I10)') '    total number of elements in Jacobian:  ', n*m
    write(*,'(A,1X,I10)') '    number of nonzero elements in Jacobian:', n_nonzero
    write(*,'(A,1X,I10)') '    number of segments:                    ', n_segs
    write(*,'(A,1X,I10)') '    number of revs:                        ', me%n_revs
    write(*,*) ''

    ! arrays for the gradient computations:
    allocate(dpert(n))
    allocate(xlow(n))
    allocate(xhigh(n))
    dpert = 1.0e-06_wp

    ! need to construct some bounds so the sparsity can be computed.
    ! [note: these are not used since we are inputing the sparsity below]
    xlow  = -10.0_wp
    xhigh = 10.0_wp

    ! first set up the gradient computation in the base class:
    call me%initialize(n,m, &
                       xlow                       = xlow,&
                       xhigh                      = xhigh,&
                       dpert                      = dpert,&
                       problem_func               = my_func,&
                       info                       = halo_grad_info,&
                       sparsity_mode              = 3,&        ! specified below
                       jacobian_method            = 3,&        ! standard central diff
                       perturb_mode               = 1,&        ! absolute mode
                       partition_sparsity_pattern = .true.,&   ! partition the pattern
                       cache_size                 = 1000 )

    ! generate and set the sparsity pattern for this problem:
    call me%get_sparsity_pattern(irow,icol)
    call me%set_sparsity_pattern(irow,icol)

    ! set up the ephemeris:
    !write(*,*) 'loading ephemeris file: '//trim(me%ephemeris_file)
    if (me%use_splined_ephemeris) then
        ! have to compute et0, etf based on mission range.
        ! use a 2 period buffer
        et0 = me%et_ref - 2.0*me%period*day2sec
        etf = me%et_ref + (me%period*day2sec * me%n_revs) + 2.0*me%period*day2sec
        allocate(jpl_ephemeris_splined :: me%eph)
        associate (eph => me%eph)
            select type(eph)
            class is (jpl_ephemeris_splined)
                call eph%initialize_splinded_ephemeris(filename=me%ephemeris_file,&
                                                       status_ok=status_ok,&
                                                       et0=et0,dt=me%dt_spline_sec,etf=etf)
            end select
        end associate
        if (.not. status_ok) error stop 'error initializing splined ephemeris'
    else
        allocate(jpl_ephemeris :: me%eph)
        call me%eph%initialize(filename=me%ephemeris_file,status_ok=status_ok)
        if (.not. status_ok) error stop 'error initializing ephemeris'
    end if

    if (.not. me%pointmass_central_body) then
        ! set up the force model [main body is moon]:
        allocate(me%grav)
        call me%grav%initialize(me%gravfile,grav_n,grav_m,status_ok)
        if (.not. status_ok) error stop 'error initializing gravity model'
    end if

    ! now, we set up the segment structure for the problem we are solving:
    ! This is for the "forward-backward" method from the paper (see Figure 2b):
    allocate(me%segs(n_segs))

    ! set up the integrators:
    do i = 1, n_segs

        call me%segs(i)%initialize(n_eoms,maxnum,ballistic_derivs,&
                                   [me%rtol],[me%atol],&
                                   report=trajectory_export_func)

        ! make a copy of this global variable in the segment.
        ! [this was set in read_config_file]
        ! [figure out a way to avoid this..]
        me%segs(i)%data%et_ref = me%et_ref

        if (use_openmp) then
            ! make a copy for each segment, so they can run in parallel
            if (.not. me%pointmass_central_body) allocate(me%segs(i)%grav, source = me%grav)  ! maybe not necessary? (is this threadsafe?)
            allocate(me%segs(i)%eph,  source = me%eph)
        else
            ! for serial use, each seg just points to the global ones for the whole mission
            if (.not. me%pointmass_central_body) me%segs(i)%grav => me%grav
            me%segs(i)%eph  => me%eph
        end if
        me%segs(i)%pointmass_central_body = me%pointmass_central_body

        ! name all the segments:
        write(seg_name,'(I10)') i
        me%segs(i)%name = trim(adjustl(seg_name))

    end do

    ! load the initial guess from the patch point file:
    i = 0
    iseg = 0
    t_periapsis = 0.0_wp
    do irev = 1, me%n_revs

        ! these are all unscaled values:

        t0(1) = t_periapsis + me%periapsis%t

        if (me%fix_initial_r .and. allocated(me%initial_r)) then
            x0_rotating(1:3,1) = me%initial_r(1:3) ! use the user-specified r vector
            x0_rotating(4:6,1) = me%periapsis%rv(4:6)
        else
            x0_rotating(:,1) = me%periapsis%rv
        end if

        tf(1) = t0(1) + me%period8

        t0(2) = t_periapsis + me%quarter%t
        x0_rotating(:,2) = me%quarter%rv
        tf(2) = tf(1)

        t0(4) = t_periapsis + me%apoapsis%t
        x0_rotating(:,4) = me%apoapsis%rv
        tf(4) = t0(4) - me%period8

        t0(3) = t0(2)
        x0_rotating(:,3) = x0_rotating(:,2)
        tf(3) = tf(4)

        t0(5) = t0(4)
        x0_rotating(:,5) = x0_rotating(:,4)
        tf(5) = t0(5) + me%period8

        t0(6) = t_periapsis + me%apoapsis%t + me%quarter%t
        x0_rotating(:,6) = [ me%quarter%rv(1),-me%quarter%rv(2),me%quarter%rv(3),&
                            -me%quarter%rv(4),me%quarter%rv(5),-me%quarter%rv(6)]
        tf(6) = tf(5)

        t0(8) = t_periapsis + me%period
        x0_rotating(:,8) = me%periapsis%rv
        tf(8) = t0(8) - me%period8

        t0(7) = t0(6)
        x0_rotating(:,7) = x0_rotating(:,6)
        tf(7) = tf(8)

        do j = 1, 8

            ! put the data for this rev into the segments:
            call me%segs(iseg+j)%set_input(t0(j),tf(j),x0_rotating(:,j))

            ! also set the scales:
            ! the t0 gradually grows larger so use the initial value for each segment.
            me%segs(iseg+j)%data%t0_scale = abs(magnitude(2.0_wp*me%segs(iseg+j)%data%t0))    ! magic number

            ! do the same for the states, but just in case, specify the min values:
            me%segs(iseg+j)%data%x0_rotating_scale = abs(magnitude(2.0_wp*me%segs(iseg+j)%data%x0_rotating, xscale_x0))

            me%segs(iseg+j)%data%xf_rotating_scale = fscale_xf  ! these are all the same for the constraints

            ! scale the dperts
            me%segs(iseg+j)%data%x0_rotating_dpert = me%segs(iseg+j)%data%x0_rotating_dpert / me%segs(iseg+j)%data%x0_rotating_scale
            me%segs(iseg+j)%data%t0_dpert = me%segs(iseg+j)%data%t0_dpert / me%segs(iseg+j)%data%t0_scale

        end do

        ! for next rev:
        i = i + 35
        iseg = iseg + 8
        t_periapsis = t_periapsis + me%period ! add another period

    end do

    allocate(me%xscale(n))
    allocate(me%xname(n))
    allocate(me%dpert_(n))
    allocate(me%fscale(m))

    ! also set the scale factors:
    call me%get_scales_from_segs()

    ! update with the new dperts:
    call me%set_dpert(me%dpert_)

    if (present(x)) then
        ! get the initial guess from the mission:
        allocate(x(n))    ! size the opt var vector
        call me%get_problem_arrays(x=x)
    end if

    end subroutine initialize_the_mission