plot_trajectory Subroutine

public subroutine plot_trajectory(me, filename, export_trajectory, draw_trajectory, only_first_rev)

Plot the trajectory using matplotlib and/or generate a text file (for MKSPK)

Note

It is assumed that all the data is present in the segments needed to propagate.

Type Bound

mission_type

Arguments

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

plot file name [without extension]
logical, intent(in), optional :: export_trajectory

if true, a text trajectory file is also produced (that can be read by MKSPK) [default is False]
logical, intent(in), optional :: draw_trajectory

if true [Default], a matplotlib plot is generated
logical, intent(in), optional :: only_first_rev

to only do the first rev [Default is False]

Calls

proc~~plot_trajectory~~CallsGraph proc~plot_trajectory mission_type%plot_trajectory add_3d_plot add_3d_plot proc~plot_trajectory->add_3d_plot add_sphere add_sphere proc~plot_trajectory->add_sphere initialize initialize proc~plot_trajectory->initialize proc~destroy_trajectory trajectory%destroy_trajectory proc~plot_trajectory->proc~destroy_trajectory proc~get_case_name mission_type%get_case_name proc~plot_trajectory->proc~get_case_name proc~propagate_segment segment%propagate_segment proc~plot_trajectory->proc~propagate_segment savefig savefig proc~plot_trajectory->savefig first_call first_call proc~propagate_segment->first_call icrf_frame icrf_frame proc~propagate_segment->icrf_frame integrate integrate proc~propagate_segment->integrate proc~set_segment_outputs segment%set_segment_outputs proc~propagate_segment->proc~set_segment_outputs transform transform proc~propagate_segment->transform two_body_rotating_frame two_body_rotating_frame proc~propagate_segment->two_body_rotating_frame

Called by

proc~~plot_trajectory~~CalledByGraph proc~plot_trajectory mission_type%plot_trajectory proc~halo_solver_main halo_solver_main proc~halo_solver_main->proc~plot_trajectory

Source Code

    subroutine plot_trajectory(me,filename,export_trajectory,draw_trajectory,only_first_rev)

    implicit none

    class(mission_type),intent(inout) :: me
    character(len=*),intent(in) :: filename !! plot file name [without extension]
    logical,intent(in),optional :: export_trajectory    !! if true, a text trajectory
                                                        !! file is also produced
                                                        !! (that can be read by MKSPK)
                                                        !! [default is False]
    logical,intent(in),optional :: draw_trajectory      !! if true [Default], a matplotlib
                                                        !! plot is generated
    logical,intent(in),optional :: only_first_rev       !! to only do the first rev [Default is False]

    type(pyplot) :: plt  !! for generating the plots
    integer :: iseg  !! segment number counter
    integer :: istat !! pyplot status code
    character(len=10) :: iseg_str !! string version of segment number
    logical :: export  !! if the txt file is to be produced
    logical :: plot    !! if the py file is to be produced
    integer :: iunit !! file unit for the txt file
    integer :: i !! counter for txt file write
    integer :: istart,iend,istep,iendprev  !! index counters
    integer :: nsegs_to_plot !! number of segments to plot

    integer,dimension(2),parameter :: figsize = [20,20] !! figure size for plotting
    logical,parameter :: plot_rotating = .true.  !! if true, the rotating state is plotted.
                                                 !! if false, the inertial state is plotted.

    ! optional arguments:
    if (present(export_trajectory)) then
        export = export_trajectory
    else
        export = .false.
    end if
    if (present(draw_trajectory)) then
        plot = draw_trajectory
    else
        plot = .true.
    end if
    nsegs_to_plot = size(me%segs) ! default export all the segments
    if (present(only_first_rev)) then
        if (only_first_rev) nsegs_to_plot = 8 ! only the first rev (8 segments)
    end if

    if (.not. plot .and. .not. export) return ! nothing to do

    if (plot) then
        ! initialize the plot:
        call plt%initialize(grid=.true.,&
                            xlabel='\n\n\nx (km)',&
                            ylabel='\n\n\ny (km)',&
                            zlabel='\n\n\nz (km)',&
                            figsize         =figsize,&
                            font_size       = 20,&
                            axes_labelsize  = 25,&
                            xtick_labelsize = 20,&
                            ytick_labelsize = 20,&
                            ztick_labelsize = 20,&
                            title='Halo Trajectory: '//trim(filename),&
                            legend=.false.,&
                            axis_equal=.true.,&
                            mplot3d=.true.)
    end if

    if (export) then
        open(newunit=iunit,file=trim(filename)//'_'//me%get_case_name()//&
                                '.txt',status='REPLACE',iostat=istat)
        if (istat/=0) error stop 'error opening trajectory file.'
    end if

    do iseg = 1, size(me%segs)
        ! destroy all trajectories:
        call me%segs(iseg)%traj_inertial%destroy()
        call me%segs(iseg)%traj_rotating%destroy()
        call me%segs(iseg)%traj_se_rotating%destroy()
    end do

    do iseg = 1, nsegs_to_plot

        write(iseg_str,'(I10)') iseg

        ! generate the trajectory for this segment:
        call me%segs(iseg)%propagate(mode=2)    ! [export points]

        if (plot) then
            ! export to plot:
            if (plot_rotating) then
                call plt%add_3d_plot(x=me%segs(iseg)%traj_rotating%x,&
                                     y=me%segs(iseg)%traj_rotating%y,&
                                     z=me%segs(iseg)%traj_rotating%z,&
                                     label='seg'//trim(adjustl(iseg_str)),&
                                     linestyle='r-',linewidth=1,istat=istat)
            else
                call plt%add_3d_plot(x=me%segs(iseg)%traj_inertial%x,&
                                     y=me%segs(iseg)%traj_inertial%y,&
                                     z=me%segs(iseg)%traj_inertial%z,&
                                     label='seg'//trim(adjustl(iseg_str)),&
                                     linestyle='r-',linewidth=1,istat=istat)
            end if
        end if

        if (export) then

            ! write the trajectory data:     ! Note: THIS HAS TO BE THE INERTIAL FRAME FOR THE SPK FILE

            if (me%segs(iseg)%traj_inertial%et(2) > me%segs(iseg)%traj_inertial%et(1)) then
                ! forward propagated
                istart = 1
                iend   = size(me%segs(iseg)%traj_inertial%et)
                istep  = 1
            else
                ! backward propagated
                istart = size(me%segs(iseg)%traj_inertial%et)
                iend   = 1
                istep  = -1
            end if

            if (iseg>1) then
                ! if the last one written is the same as the first one here, don't write it
                ! [MKSPK does not allow duplicate time tags]
                if (me%segs(iseg-1)%traj_inertial%et(iendprev) >= me%segs(iseg)%traj_inertial%et(istart)) then
                    !write(*,*) 'duplicate point: ', me%segs(iseg-1)%traj_inertial%et(iendprev)
                    istart = istart + istep
                end if
            end if

            do i=istart,iend,istep
                ! inertial trajectory for SPK:
                write(iunit,'(*(E30.16E3,A,1X))',iostat=istat) &
                        me%segs(iseg)%traj_inertial%et(i),';',&
                        me%segs(iseg)%traj_inertial%x(i) ,';',&
                        me%segs(iseg)%traj_inertial%y(i) ,';',&
                        me%segs(iseg)%traj_inertial%z(i) ,';',&
                        me%segs(iseg)%traj_inertial%vx(i),';',&
                        me%segs(iseg)%traj_inertial%vy(i),';',&
                        me%segs(iseg)%traj_inertial%vz(i),''
            end do
            iendprev = iend
        end if

    end do

    if (plot) then
        ! add the moon as a sphere:
        call plt%add_sphere(r=rad_moon,xc=zero,yc=zero,zc=zero,istat=istat,color='Grey')

        ! save the plot:
        ! call plt%showfig(pyfile=trim(filename)//'.py',istat=istat)
        call plt%savefig(trim(filename)//'_'//me%get_case_name()//'.png',&
                        pyfile=trim(filename)//'_'//me%get_case_name()//'.py',dpi='300',&
                        istat=istat)

        ! cleanup:
        call plt%destroy()
    end if

    if (export) close(iunit,iostat=istat)

    end subroutine plot_trajectory