generate_eclipse_data Subroutine

public subroutine generate_eclipse_data(me, fileprefix, filetype)

Propagate all the segments with a dense time step and export the Eclipsing data w.r.t. Earth. In this file, any point <0 is in eclipse.

Based on plot_trajectory

TODO: use the data from export_trajectory_json_file, so that should be run first

Type Bound

mission_type

Arguments

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

file prefix for the csv file (case name will be added)
integer, intent(in), optional :: filetype

type of output file: 1=csv [default] or 2=json

Calls

proc~~generate_eclipse_data~~CallsGraph proc~generate_eclipse_data mission_type%generate_eclipse_data add add proc~generate_eclipse_data->add icrf_frame icrf_frame proc~generate_eclipse_data->icrf_frame initialize initialize proc~generate_eclipse_data->initialize proc~destroy_trajectory trajectory%destroy_trajectory proc~generate_eclipse_data->proc~destroy_trajectory proc~get_case_name mission_type%get_case_name proc~generate_eclipse_data->proc~get_case_name proc~get_sun_fraction get_sun_fraction proc~generate_eclipse_data->proc~get_sun_fraction proc~propagate_segment segment%propagate_segment proc~generate_eclipse_data->proc~propagate_segment transform transform proc~generate_eclipse_data->transform two_body_rotating_frame two_body_rotating_frame proc~generate_eclipse_data->two_body_rotating_frame proc~apparent_position apparent_position proc~get_sun_fraction->proc~apparent_position proc~cubic_shadow_model cubic_shadow_model proc~get_sun_fraction->proc~cubic_shadow_model proc~from_j2000moon_to_j2000ssb from_j2000moon_to_j2000ssb proc~get_sun_fraction->proc~from_j2000moon_to_j2000ssb proc~propagate_segment->icrf_frame proc~propagate_segment->transform proc~propagate_segment->two_body_rotating_frame first_call first_call proc~propagate_segment->first_call integrate integrate proc~propagate_segment->integrate proc~set_segment_outputs segment%set_segment_outputs proc~propagate_segment->proc~set_segment_outputs axis_angle_rotation axis_angle_rotation proc~apparent_position->axis_angle_rotation cross cross proc~apparent_position->cross proc~get_pos get_pos proc~apparent_position->proc~get_pos unit unit proc~apparent_position->unit acosd acosd proc~cubic_shadow_model->acosd asind asind proc~cubic_shadow_model->asind proc~cubic_shadow_model->unit wrap_angle wrap_angle proc~cubic_shadow_model->wrap_angle proc~from_j2000moon_to_j2000ssb->icrf_frame proc~from_j2000moon_to_j2000ssb->transform get_r get_r proc~get_pos->get_r get_rv get_rv proc~get_pos->get_rv

Called by

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

Source Code

    subroutine generate_eclipse_data(me,fileprefix,filetype)

        class(mission_type),intent(inout) :: me
        character(len=*),intent(in) :: fileprefix !! file prefix for the csv file (case name will be added)
        integer,intent(in),optional :: filetype !! type of output file: 1=csv [default] or 2=json

        integer :: i,iunit,iseg,istat,istart,istep,iend
        real(wp) :: phi !! sunfrac value
        real(wp) :: p
        integer :: ifiletype
        type(json_file) :: json
        character(len=10) :: iseg_str
        character(len=:),allocatable :: full_filename
        real(wp),dimension(:),allocatable :: phi_vec !! vector of `phi` values from a segment (for JSON file)
        real(wp),dimension(:),allocatable :: et_vec !! vector of `et` values from a segment (for JSON file)
        real(wp),dimension(:),allocatable :: phi_vec_total, et_vec_total !! cumulative for all segments

        real(wp),dimension(:),allocatable :: x_se, y_se, z_se, vx_se, vy_se, vz_se, et_se, phi_se
        type(icrf_frame) :: inertial
        type(two_body_rotating_frame) :: se_rotating
        real(wp),dimension(6) :: x_se_rotating
        logical :: status_ok

        integer,parameter :: FILETYPE_CSV = 1
        integer,parameter :: FILETYPE_JSON = 2
        character(len=*),dimension(*),parameter :: file_ext = ['.csv ', &
                                                               '.json']

        if (present(filetype)) then
            ifiletype = filetype
        else
            ifiletype = FILETYPE_JSON
        end if
        full_filename = trim(fileprefix)//'_'//me%get_case_name()//&
                        trim(file_ext(ifiletype))

        write(*,*) 'Generating eclipse file: '//full_filename

        select case (ifiletype)
        case (FILETYPE_CSV)
            open(newunit=iunit,file=full_filename,status='REPLACE',iostat=istat)
            if (istat/=0) error stop 'error opening eclipse csv file.'
            write(iunit,'(A30,A,1X,A30)',iostat=istat) 'ET (sec)', ',', 'PHI'
        case (FILETYPE_JSON)
            call json%initialize(compress_vectors=.true.)
        case default
            error stop 'invalid filetype for eclipse file'
        end select

        allocate(x_se(0));allocate(y_se(0));allocate(z_se(0))
        allocate(vx_se(0));allocate(vy_se(0));allocate(vz_se(0))
        allocate(et_se(0))
        allocate(phi_se(0))
        allocate(phi_vec_total(0))
        allocate(et_vec_total(0))
        do iseg = 1, size(me%segs)

            if (ifiletype==FILETYPE_JSON) then
                if (allocated(phi_vec)) deallocate(phi_vec); allocate(phi_vec(0))
                if (allocated(et_vec))  deallocate(et_vec);  allocate(et_vec(0))
            end if
            write(iseg_str,'(I10)') iseg    ! segment number as string
            iseg_str = adjustl(iseg_str)

            ! destroy all trajectories first:
            call me%segs(iseg)%traj_inertial%destroy()
            call me%segs(iseg)%traj_rotating%destroy()
            call me%segs(iseg)%traj_se_rotating%destroy()

            ! generate the trajectory for this segment:
            ! [export points at a fixed time step]
            call me%segs(iseg)%propagate(mode=2,tstep=me%eclipse_dt_step)

            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

            do i=istart,iend,istep
                associate (et => me%segs(iseg)%traj_inertial%et(i),&
                            rv_moon => [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) ])

                    p = get_sun_fraction(me%eph, et, rv_moon, me%r_eclipse_bubble)
                    phi = min(0.0_wp, p) ! phi >0 mean the spacecraft is in sunlight

                    select case (ifiletype)
                    case (FILETYPE_CSV)
                        write(iunit,'(*(E30.16E3,A,1X))',iostat=istat) et, ',', phi
                    case (FILETYPE_JSON)
                        ! for this we accumulate the data for the whole
                        ! segment and write all at once at the end
                        phi_vec = [phi_vec, phi] ! TODO: could use expand routines to make this more efficient
                        phi_se = [phi_se, p] ! save the actual value
                        et_se = [et_se, et]

                        ! also convert to sun-earth rotating
                        se_rotating = two_body_rotating_frame(primary_body=body_sun,&
                                                              secondary_body=body_earth,&
                                                              center=center_at_secondary_body,&
                                                              et=et)
                        ! from inertial to rotating:
                        inertial = icrf_frame(b=body_moon)
                        call inertial%transform(rv_moon,se_rotating,et,me%segs(iseg)%eph,x_se_rotating,status_ok)
                        if (.not. status_ok) error stop 'transformation error in propagate_segment'
                        x_se = [x_se, x_se_rotating(1)]
                        y_se = [y_se, x_se_rotating(2)]
                        z_se = [z_se, x_se_rotating(3)]

                    end select

                end associate
            end do
            if (ifiletype==FILETYPE_JSON) then
                ! accumulate only the points where phi < 0 (in shadow)
                et_vec  = pack(me%segs(iseg)%traj_inertial%et, mask=phi_vec<0.0_wp)
                phi_vec = pack(phi_vec, mask=phi_vec<0.0_wp)
                if (size(et_vec)>0) then
                    et_vec_total = [et_vec_total, et_vec]
                    phi_vec_total = [phi_vec_total, phi_vec]
                end if
            end if

        end do

        ! save file:
        select case (ifiletype)
        case (FILETYPE_CSV)
            close(iunit,iostat=istat)
        case (FILETYPE_JSON)
            call json%add('et',  et_vec_total)     ! TODO: could remove any duplicate time entries (e.g., at the segment interfaces)
            call json%add('phi', phi_vec_total)

            call json%add('phi_se', phi_se) ! also the full S-E trajectory data for plotting
            call json%add('et_se', et_se)
            call json%add('x_se_rotating', x_se)
            call json%add('y_se_rotating', y_se)
            call json%add('z_se_rotating', z_se)
            call json%print(full_filename)
        end select
        call json%destroy()

    end subroutine generate_eclipse_data