export_rp_ra_json_file Subroutine

    subroutine export_rp_ra_json_file(me,et,rmag,rdot,filename)

    use root_module
    use bspline_module

    class(mission_type),intent(inout) :: me
    real(wp),dimension(:),intent(in) :: et
    real(wp),dimension(:),intent(in) :: rmag
    real(wp),dimension(:),intent(in) :: rdot
    character(len=*),intent(in) :: filename !! file name [without extension]

    type,extends(chandrupatla_solver) :: my_solver
        !! the solver to use for the root finding
    end type my_solver

    real(wp) :: et0, etf, et_root, rdot_root, initial_et, final_et, rdot0, rdotf, rmag_root
    integer :: i, n
    real(wp),dimension(:),allocatable :: et_for_rp, et_for_ra, rp_vec, ra_vec
    real(wp) :: x, f
    integer :: iflag
    type(my_solver) :: solver
    type(bspline_1d) :: rdot_spline, rmag_spline
    type(json_file) :: json
    logical :: finished !! to indicate the arrays are done
    integer :: n_rp, n_ra, n_rp_et, n_ra_et !! counters for array sizes

    integer,parameter :: kx = 4  !! spline order (cubic bspline)
    integer,parameter :: idx = 0 !! interpolate value only
    real(wp),parameter :: et_step = 3600.0_wp * 24.0_wp !! et step size [should be an input] (sec)
    integer,parameter :: chunk_size = 1000 !! chunk size for expanding the arrays

    write(*,'(A)') ' * Generating ra/rp file.'

    ! error check to make sure et is strictly increasing.
    do i = 2, size(et)
        if (et(i)<=et(i-1)) then
            write(*,*) 'error in et point ', i
            write(*,*) 'et(i)  ', et(i)
            write(*,*) 'et(i-1)', et(i-1)
            error stop
        end if
    end do
    ! first spline rmag and rdot as a functions of et:
    call rdot_spline%initialize(et,rdot,kx,iflag,.false.); call check_spline('rdot')
    call rmag_spline%initialize(et,rmag,kx,iflag,.false.); call check_spline('rmag')

    allocate(et_for_rp(0), et_for_ra(0), rp_vec(0), ra_vec(0)) ! arrays to hold results

    ! step through the trajectory and find all the rdot roots (periapsis and apoapsis):
    call solver%initialize(rdot_func)
    n = size(et)
    initial_et = et(1)
    final_et = et(n)
    et0 = initial_et
    n_rp = 0; n_ra = 0; n_rp_et = 0; n_ra_et = 0
    do
        etf = min(final_et, et0+et_step)
        ! if there is a root on this interval (change of sign of rdot):
        rdot0 = rdot_func(solver, et0)
        rdotf = rdot_func(solver, etf)
        finished = etf==et(n) ! the last step
        if (rdot0*rdotf<=zero) then  ! there is a root on this interval
            ! find the root:
            call solver%solve(et0,etf,et_root,rdot_root,iflag)
            rmag_root = rmag_func(et_root)  ! get corresponding rmag at the root
            if (rdotf>rdot0) then
                ! increasing - periapsis
                ! et_for_rp = [et_for_rp, et_root]
                ! rp_vec    = [rp_vec, rmag_root]
                call add_to(et_for_rp, et_root, n_rp_et, chunk_size, finished)
                call add_to(rp_vec, rmag_root, n_rp, chunk_size, finished)
            else
                ! decreasing - apoapsis
                !et_for_ra = [et_for_ra, et_root]
                !ra_vec    = [ra_vec, rmag_root]
                call add_to(et_for_ra, et_root, n_ra_et, chunk_size, finished)
                call add_to(ra_vec, rmag_root, n_ra, chunk_size, finished)
            end if
        end if
        if (finished) exit
        et0 = etf ! set up for next step
    end do

    ! output results to file:
    call json%initialize(compress_vectors=.true.)
    call json%add('et_for_rp', et_for_rp)
    call json%add('rp_vec',    rp_vec)
    call json%add('et_for_ra', et_for_ra)
    call json%add('ra_vec',    ra_vec)
    call json%print(trim(filename)//'_rp_ra.json')
    call json%destroy()

    contains

        pure subroutine add_to(vec, val, n, chunk_size, finished)
            !! resize an array in chunks

            real(wp), dimension(:), allocatable, intent(inout) :: vec
                !! the vector to add to
            real(wp), intent(in) :: val
                !! the value to add
            integer, intent(inout) :: n
                !! counter for last element added to vec.
                !! must be initialized to size(vec)
                !! (or 0 if not allocated) before first call
            integer, intent(in) :: chunk_size
                !! allocate vec in blocks of this size (>0)
            logical, intent(in) :: finished
                !! set to true to return vec
                !! as its correct size (n)

            real(wp), dimension(:), allocatable :: tmp

            if (allocated(vec)) then
                if (n == size(vec)) then
                    ! have to add another chunk:
                    allocate (tmp(size(vec) + chunk_size))
                    tmp(1:size(vec)) = vec
                    call move_alloc(tmp, vec)
                end if
                n = n + 1
            else
                ! the first element:
                allocate (vec(chunk_size))
                n = 1
            end if

            vec(n) = val

            if (finished) then
                ! set vec to actual size (n):
                if (allocated(tmp)) deallocate (tmp)
                allocate (tmp(n))
                tmp = vec(1:n)
                call move_alloc(tmp, vec)
            end if

        end subroutine add_to

        subroutine check_spline(case) !! error checking for splines
            character(len=*),intent(in) :: case !! case name
            if (iflag/=0) then
                select case (iflag)
                case(2); error stop case//' : bspline error: iknot out of range.'
                case(3); error stop case//' : bspline error: nx out of range.'
                case(4); error stop case//' : bspline error: kx out of range.'
                case(5); error stop case//' : bspline error: x not strictly increasing.'
                end select
            end if
        end subroutine check_spline

        function rdot_func(me,x)  !! compute rdot from spline
            class(root_solver),intent(inout) :: me
            real(wp),intent(in) :: x
            real(wp) :: rdot_func
            integer :: iflag
            call rdot_spline%evaluate(x,idx,rdot_func,iflag)
            if (iflag/=0) error stop 'error evaluating bspline for rdot.'
        end function rdot_func
        function rmag_func(x)  !! compute rmag from spline
            real(wp),intent(in) :: x
            real(wp) :: rmag_func
            integer :: iflag
            call rmag_spline%evaluate(x,idx,rmag_func,iflag)
            if (iflag/=0) error stop 'error evaluating bspline for rmag.'
        end function rmag_func

    end subroutine export_rp_ra_json_file