initialize_globals Subroutine

    subroutine initialize_globals(me, et0, dt, etf)

        !! initialize the earth and sun ephemeris.

        class(jpl_ephemeris_splined),intent(inout) :: me
        real(wp),intent(in) :: et0 !! initial ephemeris time [sec]
        real(wp),intent(in) :: dt  !! ephemeris time step [sec]
        real(wp),intent(in) :: etf !! final ephemeris time [sec]

        integer :: i !! counter
        real(wp) :: et !! ephemeris time
        integer :: nx !! number of points
        real(wp),dimension(6) :: rv_earth, rv_sun
        logical :: status_ok
        real(wp),dimension(:),allocatable :: et_vec
        integer :: iflag

        ! write(*,*) 'initialize_globals...'
        ! write(*,*) '  et0 = ', et0
        ! write(*,*) '  dt  = ', dt
        ! write(*,*) '  etf = ', etf

        call me%earth_eph_interface%destroy()
        call me%sun_eph_interface%destroy()
        call me%ssb_eph_interface%destroy()

        if (allocated(earth_eph%tx)) deallocate(earth_eph%tx)
        if (allocated(sun_eph%tx)) deallocate(sun_eph%tx)
        if (allocated(ssb_eph%tx)) deallocate(ssb_eph%tx)
        if (allocated(earth_eph%bcoef)) deallocate(earth_eph%bcoef)
        if (allocated(sun_eph%bcoef)) deallocate(sun_eph%bcoef)
        if (allocated(ssb_eph%bcoef)) deallocate(ssb_eph%bcoef)
        if (allocated(earth_eph%f)) deallocate(earth_eph%f)
        if (allocated(sun_eph%f)) deallocate(sun_eph%f)
        if (allocated(ssb_eph%f)) deallocate(ssb_eph%f)
        if (allocated(et_vec)) deallocate(et_vec)

        ! first, count the number of points and allocate the arrays:
        nx = 0
        et = et0
        do
            if (et>etf) exit
            nx = nx + 1
            et = et + dt
        end do
        allocate(earth_eph%tx(   nx+kx,6))  ! columns are for each state element (rx,ry,rz,vx,vy,vz)
        allocate(sun_eph%tx(     nx+kx,6))
        allocate(ssb_eph%tx(     nx+kx,6))
        allocate(earth_eph%bcoef(nx,6))
        allocate(sun_eph%bcoef(  nx,6))
        allocate(ssb_eph%bcoef(  nx,6))
        allocate(earth_eph%f(    nx,6))
        allocate(sun_eph%f(      nx,6))
        allocate(ssb_eph%f(      nx,6))
        allocate(et_vec(         nx))

        ! function calls from et0 to etf:
        i = 0 ! index in the arrays
        et = et0
        do
            if (et>etf) exit
            i = i + 1
            et = et + dt
            et_vec(i) = et
            ! get_rv_from_jpl_ephemeris(me,et,targ,obs,rv,status_ok)
            call me%jpl_ephemeris%get_rv(et,body_earth,body_moon,earth_eph%f(i,:),status_ok)
            if (.not. status_ok) error stop 'earth eph error'
            call me%jpl_ephemeris%get_rv(et,body_sun,body_moon,sun_eph%f(i,:),status_ok)
            if (.not. status_ok) error stop 'sun eph error'
            call me%jpl_ephemeris%get_rv(et,body_sun,body_moon,ssb_eph%f(i,:),status_ok)
            if (.not. status_ok) error stop 'ssb eph error'
        end do

        ! create the splines (one for each coordinate):
        do i = 1, 6
            call db1ink(et_vec, nx, earth_eph%f(:,i), kx, iknot, earth_eph%tx(:,i), earth_eph%bcoef(:,i), iflag)
            if (iflag/=0) then
                write(*,*) 'db1ink iflag = ', iflag
                error stop 'db1ink error'
            end if
            call db1ink(et_vec, nx, sun_eph%f(:,i), kx, iknot, sun_eph%tx(:,i), sun_eph%bcoef(:,i), iflag)
            if (iflag/=0) then
                write(*,*) 'db1ink iflag = ', iflag
                error stop 'db1ink error'
            end if
            call db1ink(et_vec, nx, ssb_eph%f(:,i), kx, iknot, ssb_eph%tx(:,i), ssb_eph%bcoef(:,i), iflag)
            if (iflag/=0) then
                write(*,*) 'db1ink iflag = ', iflag
                error stop 'db1ink error'
            end if
        end do
        deallocate(earth_eph%f) ! don't need these anymore
        deallocate(sun_eph%f)
        deallocate(ssb_eph%f)

        earth_eph%nx = nx
        sun_eph%nx = nx
        ssb_eph%nx = nx

    end subroutine initialize_globals