apps160222_jdc.scenario4ΒΆ

""" JUICE senario 4. Moon flyby.

I chose first Europa flyby as scenario 4.

The time is  2030-10-05T10:57:19 with closest approach of 1964 km (from center).
"""
import datetime
import irfpy.pep.pep_attitude as patt
import irfpy.juice.jspice as js
#js.furnsh('./160222.mk')
#js.init()
#js.define_juice()

import matplotlib.patches
import matplotlib.dates as mpldates

import numpy as np

import os

js.spice.furnsh('scenario.mk')

n_jp = js.spice.bodn2c('JUPITER')
n_sun = js.spice.bodn2c('SUN')
n_earth = js.spice.bodn2c('EARTH')
n_ganymede = js.spice.bodn2c('GANYMEDE')
n_io = js.spice.bodn2c('IO')
n_europa = js.spice.bodn2c('EUROPA')
n_callisto = js.spice.bodn2c('CALLISTO')

rj = 71492.
re = 1560.8

import irfpy.juice.attitude0 as jatt

import matplotlib.pyplot as plt

def length(v):
    return np.sqrt(np.inner(v, v))

def unit_vector(v):
    return v / np.sqrt(np.inner(v, v))

def longitude(v):
    return np.rad2deg(np.arctan2(v[1], v[0]))

def latitude(v):
    vv = unit_vector(v)
    return np.rad2deg(np.arcsin(vv[2]))

def get_state(n_target, et, frame, n_origin, unit=rj):
    posvel, lt = js.spice.spkez(n_target, et, frame, 'LT+S', n_origin)
    return posvel[:3] / unit, posvel[3:]


def main():
    fp = open('scenario4.txt', 'w')
    print('#datetime mpldates et vx_jse vy_jse vz_jse vx_sc vy_sc vz_sc', file=fp)
    print('# Velocity is measured in the spacecraft frame', file=fp)

    pathname = 'scenario4'
    try:
        os.mkdir(pathname)
    except:
        pass

    ###  2030-10-05T10:57:19 plus minus 12 hr

    t0 = datetime.datetime(2030, 10, 4, 23)
    t1 = datetime.datetime(2030, 10, 5, 23)
    dt = datetime.timedelta(minutes=10)     # Time difference to be confirmed.


    t = t0 - dt
#    t = t1 - dt

    while t <= t1:
        t = t + dt

        et = js.spice.str2et(t.strftime('%FT%T'))

        ### Ok, now many object positions can be calculated.
        ### The variable names are "target"_"frame"_"origin"

        ### Here SC centered JSE positions
#        jupiter_jse_juice, v_jupiter_jse_juice = get_state(n_jp, et, "JSE", -907)
        sun_jse_juice, v_sun_jse_juice = get_state(n_sun, et, "JSE", -907)
#        earth_jse_juice, __vel = get_state(n_earth, et, "JSE", -907)
#        io_jse_juice, __vel = get_state(n_io, et, "JSE", -907)
#        europa_jse_juice, __vel = get_state(n_europa, et, "JSE", -907)
#        ganymede_jse_juice, __vel = get_state(n_ganymede, et, "JSE", -907)
#        callisto_jse_juice, __vel = get_state(n_callisto, et, "JSE", -907)

        ### Here JSE positions
        juice_jse_jupiter, v_juice_jse_jupiter = get_state(-907, et, "JSE", n_jp)
        juice_jse_sun, v_juice_jse_sun = get_state(-907, et, "JSE", n_sun)
        juice_jse_earth, v_juice_jse_earth = get_state(-907, et, "JSE", n_earth)
        juice_jse_europa, v_juice_jse_europa = get_state(-907, et, "JSE", n_europa)
        io_jse_jupiter, __vel = get_state(n_io, et, "JSE", n_jp)
        europa_jse_jupiter, __vel = get_state(n_europa, et, "JSE", n_jp)
        ganymede_jse_jupiter, __vel = get_state(n_ganymede, et, "JSE", n_jp)
        callisto_jse_jupiter, __vel = get_state(n_callisto, et, "JSE", n_jp)
        sun_jse_jupiter, __vel = get_state(n_sun, et, "JSE", n_jp)
        sun_jse_jupiter_unitvector = unit_vector(sun_jse_jupiter)
        earth_jse_jupiter, __vel = get_state(n_earth, et, "JSE", n_jp)
        earth_jse_jupiter_unitvector = unit_vector(earth_jse_jupiter)

        ### Here System3 positions (many velocity not used)
        juice_sys3_jupiter, v_juice_sys3_jupiter = get_state(-907, et, "IAU_JUPITER", n_jp)
        io_sys3_jupiter, __vel = get_state(n_io, et, "IAU_JUPITER", n_jp)
        europa_sys3_jupiter, __vel = get_state(n_europa, et, "IAU_JUPITER", n_jp)
        ganymede_sys3_jupiter, __vel = get_state(n_ganymede, et, "IAU_JUPITER", n_jp)
        callisto_sys3_jupiter, __vel = get_state(n_callisto, et, "IAU_JUPITER", n_jp)
        sun_sys3_jupiter, __vel = get_state(n_sun, et, "IAU_JUPITER", n_jp)
        sun_sys3_jupiter_unitvector = unit_vector(sun_sys3_jupiter)
        earth_sys3_jupiter, __vel = get_state(n_earth, et, "IAU_JUPITER", n_jp)
        earth_sys3_jupiter_unitvector = unit_vector(earth_sys3_jupiter)

        ### Here Europa frame positions (many velocity not used)
        juice_europa_europa, v_juice_europa_europa = get_state(-907, et, "IAU_EUROPA", n_europa, unit=re)
        io_europa_europa, __vel = get_state(n_io, et, "IAU_EUROPA", n_europa, unit=re)
        europa_europa_europa, __vel = get_state(n_europa, et, "IAU_EUROPA", n_europa, unit=re)
        ganymede_europa_europa, __vel = get_state(n_ganymede, et, "IAU_EUROPA", n_europa, unit=re)
        callisto_europa_europa, __vel = get_state(n_callisto, et, "IAU_EUROPA", n_europa, unit=re)
        sun_europa_europa, __vel = get_state(n_sun, et, "IAU_EUROPA", n_europa, unit=re)
        sun_europa_europa_unitvector = unit_vector(sun_europa_europa)
        earth_europa_europa, __vel = get_state(n_earth, et, "IAU_EUROPA", n_europa, unit=re)
        earth_europa_europa_unitvector = unit_vector(earth_europa_europa)

        ### Conversion matrixes.
        cmat_eclip_jse = js.spice.pxform('ECLIPJ2000', 'JSE', et)   # Ecliptic -> JSE
        cmat_sys3_jse = js.spice.pxform('IAU_JUPITER', 'JSE', et)        ### System III to JSE frame conversion matrix
        cmat_sys3_europa = js.spice.pxform('IAU_JUPITER', 'IAU_EUROPA', et)   # System 3 to Europa fixed.
        cmat_jse_europa = js.spice.pxform('JSE', 'IAU_EUROPA', et)   # JSE 2 Europa

        ### Earth's orbital plane.
        eclip_norm_vec = cmat_eclip_jse.dot(np.array([0, 0, 1]))   # Almost along z by definition.

        ### Spacecraft frame definition
        juice_frame = jatt.get_juice_farflyby(juice_jse_europa, sun_jse_juice)

        ### Objects in Spacecraft frame
#        jupiter_sc_juice = juice_frame.convert_to_nsc(jupiter_jse_juice)
#        sun_sc_juice = juice_frame.convert_to_nsc(sun_jse_juice)
#        earth_sc_juice = juice_frame.convert_to_nsc(earth_jse_juice)
#        io_sc_juice = juice_frame.convert_to_nsc(io_jse_juice)
#        europa_sc_juice = juice_frame.convert_to_nsc(europa_jse_juice)
#        ganymede_sc_juice = juice_frame.convert_to_nsc(ganymede_jse_juice)
#        callisto_sc_juice = juice_frame.convert_to_nsc(callisto_jse_juice)

        fig = plt.figure(figsize=(12, 12))
        axes={}

        ### Orbit plot (xy in JSE)
        ax = fig.add_axes([0.1, 0.4, 0.35, 0.35])
        axes['orbit_jse_xy'] = ax
        # Spacecraft position in JSE
        ax.plot([juice_jse_jupiter[0]], [juice_jse_jupiter[1]], 'ko')
        # SC->JSE conversion, plot
        armlen = 3
        vec_xaxis_jse = juice_frame.convert_to_ref([1, 0 ,0])
        vec_yaxis_jse = juice_frame.convert_to_ref([0, 1 ,0])
        vec_zaxis_jse = juice_frame.convert_to_ref([0, 0 ,1])

        ax.plot(juice_jse_jupiter[0] + [0, vec_xaxis_jse[0] * armlen],
                juice_jse_jupiter[1] + [0, vec_xaxis_jse[1] * armlen],
                'r')
        ax.text(juice_jse_jupiter[0] + vec_xaxis_jse[0] * armlen,
                juice_jse_jupiter[1] + vec_xaxis_jse[1] * armlen,
                'x', color='r')
        ax.plot(juice_jse_jupiter[0] + [0, vec_yaxis_jse[0] * armlen],
                juice_jse_jupiter[1] + [0, vec_yaxis_jse[1] * armlen],
                'g')
        ax.text(juice_jse_jupiter[0] + vec_yaxis_jse[0] * armlen,
                juice_jse_jupiter[1] + vec_yaxis_jse[1] * armlen,
                'y', color='g')
        ax.plot(juice_jse_jupiter[0] + [0, vec_zaxis_jse[0] * armlen],
                juice_jse_jupiter[1] + [0, vec_zaxis_jse[1] * armlen],
                'b')
        ax.text(juice_jse_jupiter[0] + vec_zaxis_jse[0] * armlen,
                juice_jse_jupiter[1] + vec_zaxis_jse[1] * armlen,
                'z', color='b')

        # Decoration
        ax.plot([0, 0], [-100, 100], 'k:')
        ax.plot([-100, 100], [0, 0], 'k:')
        ax.plot([-100, 100], [100, -100], 'k:')
        ax.plot([-100, 100], [-100, 100], 'k:')
        x = np.cos(np.linspace(0, np.pi * 2, 100))
        y = np.sin(np.linspace(0, np.pi * 2, 100))
        for r in np.arange(5, 100, 5): ax.plot(r * x, r * y, 'k:')
        ax.plot(x, y, 'k')
        ax.set_xlim(-30, 30)
        ax.set_ylim(-30, 30)
        ax.set_xlabel('X [JSE]')
        ax.set_ylabel('Y [JSE]')
        ax.set_title('{:%FT%T}'.format(t))
        ax.set_aspect('equal')

        ax.text(io_jse_jupiter[0], io_jse_jupiter[1], 'I', color='g', ha='center', va='center')
        ax.text(europa_jse_jupiter[0], europa_jse_jupiter[1], 'E', color='g', ha='center', va='center')
        ax.text(ganymede_jse_jupiter[0], ganymede_jse_jupiter[1], 'G', color='g', ha='center', va='center')
        ax.text(callisto_jse_jupiter[0], callisto_jse_jupiter[1], 'C', color='g', ha='center', va='center')

        # Sun,earth direction
        ax.plot([0, sun_jse_jupiter_unitvector[0] * 10],
                [0, sun_jse_jupiter_unitvector[1] * 10],
                'r-')
        ax.text(sun_jse_jupiter_unitvector[0] * 10, sun_jse_jupiter_unitvector[1] * 10,
                'S', color='r')
        ax.plot([0, earth_jse_jupiter_unitvector[0] * 10],
                [0, earth_jse_jupiter_unitvector[1] * 10],
                'b-')
        ax.text(earth_jse_jupiter_unitvector[0] * 10, earth_jse_jupiter_unitvector[1]* 10,
                'E', color='b')


        ### Orbit plot (xz in JSE)
        ax = fig.add_axes([0.1, 0.78, 0.35, 0.20])
        axes['orbit_jse_xz'] = ax
#        ax = patt.axis_configuration_xz(juice_frame, axis=ax)
        # Spacecraft position in JSE
        ax.plot([juice_jse_jupiter[0]], [juice_jse_jupiter[2]], 'ko')
        # SC->JSE conversion, plot
        ax.plot(juice_jse_jupiter[0] + [0, vec_xaxis_jse[0] * armlen],
                juice_jse_jupiter[2] + [0, vec_xaxis_jse[2] * armlen],
                'r')
        ax.text(juice_jse_jupiter[0] + vec_xaxis_jse[0] * armlen,
                juice_jse_jupiter[2] + vec_xaxis_jse[2] * armlen,
                'x', color='r')
        ax.plot(juice_jse_jupiter[0] + [0, vec_yaxis_jse[0] * armlen],
                juice_jse_jupiter[2] + [0, vec_yaxis_jse[2] * armlen],
                'g')
        ax.text(juice_jse_jupiter[0] + vec_yaxis_jse[0] * armlen,
                juice_jse_jupiter[2] + vec_yaxis_jse[2] * armlen,
                'y', color='g')
        ax.plot(juice_jse_jupiter[0] + [0, vec_zaxis_jse[0] * armlen],
                juice_jse_jupiter[2] + [0, vec_zaxis_jse[2] * armlen],
                'b')
        ax.text(juice_jse_jupiter[0] + vec_zaxis_jse[0] * armlen,
                juice_jse_jupiter[2] + vec_zaxis_jse[2] * armlen,
                'z', color='b')

        # Decoration
        ax.plot([0, 0], [-1.1, 1.1], 'k-')
        ax.plot([-100, 100], [0, 0], 'k:')
        x = np.cos(np.linspace(0, np.pi * 2, 100))
        y = np.sin(np.linspace(0, np.pi * 2, 100))
        ax.plot(x, y, 'k')
        ax.set_xlim(-30, 30)
        ax.set_ylim(-10, 10)
        ax.set_xlabel('X [JSE]')
        ax.set_ylabel('Z [JSE]')
        ax.set_title('{:%FT%T}'.format(t))
        ax.set_aspect('equal')

        ax.text(io_jse_jupiter[0], io_jse_jupiter[2], 'I', color='g', ha='center', va='center')
        ax.text(europa_jse_jupiter[0], europa_jse_jupiter[2], 'E', color='g', ha='center', va='center')
        ax.text(ganymede_jse_jupiter[0], ganymede_jse_jupiter[2], 'G', color='g', ha='center', va='center')
        ax.text(callisto_jse_jupiter[0], callisto_jse_jupiter[2], 'C', color='g', ha='center', va='center')

        ### Sun direction
        ax.plot([0, sun_jse_jupiter_unitvector[0] * 10],
                [0, sun_jse_jupiter_unitvector[2] * 10],
                'r-')
        ax.text(sun_jse_jupiter_unitvector[0] * 10, sun_jse_jupiter_unitvector[2] * 10,
                'S', color='r')
        ax.plot([0, earth_jse_jupiter_unitvector[0] * 10],
                [0, earth_jse_jupiter_unitvector[2] * 10],
                'b-')
        ax.text(earth_jse_jupiter_unitvector[0] * 10, earth_jse_jupiter_unitvector[2] * 10,
                'E', color='b')


        ### S/C at Europa frame

        ### Orbit plot (xy)
        ax = fig.add_axes([0.5, 0.4, 0.35, 0.35])
        axes['orbit_moon_xy'] = ax
        ax.plot([juice_europa_europa[0]], [juice_europa_europa[1]], 'ko')
        ax.plot([0, 0], [-100, 100], 'k:')
        ax.plot([-100, 100], [0, 0], 'k:')
        ax.plot([-100, 100], [100, -100], 'k:')
        ax.plot([-100, 100], [-100, 100], 'k:')
        x = np.cos(np.linspace(0, np.pi * 2, 100))
        y = np.sin(np.linspace(0, np.pi * 2, 100))
        for r in np.arange(5, 100, 10): ax.plot(r * x, r * y, 'k:')
        ax.plot(x, y, 'k')
        ax.set_xlim(-100, 100)
        ax.set_ylim(-100, 100)
        ax.set_xlabel('X [EUROPA]')
        ax.set_ylabel('Y [EUROPA]')
#        ax.set_xscale('symlog', linthreshx=30, linscalex=2)
#        ax.set_yscale('symlog', linthreshy=30, linscaley=2)
#        import matplotlib.patches as mplp
#        bbox = mplp.Rectangle((-30, -30), 60, 60, facecolor='none', edgecolor='k', linestyle='-')
#        ax.add_patch(bbox)

        if np.abs(juice_europa_europa[0]) < 100 and np.abs(juice_europa_europa[1]) < 100:
            armlen = 15
            vec_xaxis_europa = cmat_jse_europa.dot(vec_xaxis_jse)
            vec_yaxis_europa = cmat_jse_europa.dot(vec_yaxis_jse)
            vec_zaxis_europa = cmat_jse_europa.dot(vec_zaxis_jse)

            ax.plot(juice_europa_europa[0] + [0, vec_xaxis_europa[0] * armlen],
                juice_europa_europa[1] + [0, vec_xaxis_europa[1] * armlen],
                'r')
            ax.text(juice_europa_europa[0] + vec_xaxis_europa[0] * armlen,
                juice_europa_europa[1] + vec_xaxis_europa[1] * armlen,
                'x', color='r')
            ax.plot(juice_europa_europa[0] + [0, vec_yaxis_europa[0] * armlen],
                juice_europa_europa[1] + [0, vec_yaxis_europa[1] * armlen],
                'g')
            ax.text(juice_europa_europa[0] + vec_yaxis_europa[0] * armlen,
                juice_europa_europa[1] + vec_yaxis_europa[1] * armlen,
                'y', color='g')
            ax.plot(juice_europa_europa[0] + [0, vec_zaxis_europa[0] * armlen],
                juice_europa_europa[1] + [0, vec_zaxis_europa[1] * armlen],
                'b')
            ax.text(juice_europa_europa[0] + vec_zaxis_europa[0] * armlen,
                juice_europa_europa[1] + vec_zaxis_europa[1] * armlen,
                'z', color='b')

        ax.set_title('{:%FT%T}'.format(t))

#        ax.text(io_europa_europa[0], io_europa_europa[1], 'I', color='g')
#        ax.text(europa_europa_europa[0], europa_europa_europa[1], 'E', color='g')
#        ax.text(ganymede_europa_europa[0], ganymede_europa_europa[1], 'G', color='g')
#        ax.text(callisto_europa_europa[0], callisto_europa_europa[1], 'C', color='g')

        # Sun,earth direction
        ax.plot([0, sun_europa_europa_unitvector[0] * 10],
                [0, sun_europa_europa_unitvector[1] * 10],
                'r-')
        ax.text(sun_europa_europa_unitvector[0] * 10, sun_europa_europa_unitvector[1] * 10,
                'S', color='r')
        ax.plot([0, earth_europa_europa_unitvector[0] * 10],
                [0, earth_europa_europa_unitvector[1] * 10],
                'b-')
        ax.text(earth_europa_europa_unitvector[0] * 10, earth_europa_europa_unitvector[1] * 10,
                'E', color='b')


        ### Plot of SC in x-z plane.
        ax = fig.add_axes([0.5, 0.78, 0.35, 0.20])
        axes['orbit_moon_xz'] = ax
        ax.plot([juice_europa_europa[0]], [juice_europa_europa[2]], 'ko')
        ax.plot([0, 0], [-1.1, 1.1], 'k-')
        ax.plot([-100, 100], [0, 0], 'k:')
        x = np.cos(np.linspace(0, np.pi * 2, 100))
        y = np.sin(np.linspace(0, np.pi * 2, 100))
        ax.plot(x, y, 'k')
        ax.set_xlim(-100, 100)
        ax.set_ylim(-100 / 3., 100 / 3.)
        ax.set_xlabel('X [EUROPA]')
        ax.set_ylabel('Z [EUROPA]')
        ax.set_aspect('equal')

        if np.abs(juice_europa_europa[0]) < 100 and np.abs(juice_europa_europa[2]) < 33.33:
            armlen = 15
            vec_xaxis_europa = cmat_jse_europa.dot(vec_xaxis_jse)
            vec_yaxis_europa = cmat_jse_europa.dot(vec_yaxis_jse)
            vec_zaxis_europa = cmat_jse_europa.dot(vec_zaxis_jse)

            ax.plot(juice_europa_europa[0] + [0, vec_xaxis_europa[0] * armlen],
                juice_europa_europa[2] + [0, vec_xaxis_europa[2] * armlen],
                'r')
            ax.text(juice_europa_europa[0] + vec_xaxis_europa[0] * armlen,
                juice_europa_europa[2] + vec_xaxis_europa[2] * armlen,
                'x', color='r')
            ax.plot(juice_europa_europa[0] + [0, vec_yaxis_europa[0] * armlen],
                juice_europa_europa[2] + [0, vec_yaxis_europa[2] * armlen],
                'g')
            ax.text(juice_europa_europa[0] + vec_yaxis_europa[0] * armlen,
                juice_europa_europa[2] + vec_yaxis_europa[2] * armlen,
                'y', color='g')
            ax.plot(juice_europa_europa[0] + [0, vec_zaxis_europa[0] * armlen],
                juice_europa_europa[2] + [0, vec_zaxis_europa[2] * armlen],
                'b')
            ax.text(juice_europa_europa[0] + vec_zaxis_europa[0] * armlen,
                juice_europa_europa[2] + vec_zaxis_europa[2] * armlen,
                'z', color='b')


        ### Other objects
#        ax.text(io_europa_europa[0], io_europa_europa[2], 'I', color='g')
#        ax.text(europa_europa_europa[0], europa_europa_europa[2], 'E', color='g')
#        ax.text(ganymede_europa_europa[0], ganymede_europa_europa[2], 'G', color='g')
#        ax.text(callisto_europa_europa[0], callisto_europa_europa[2], 'C', color='g')

        # Sun,earth direction
        ax.plot([0, sun_europa_europa_unitvector[0] * 10],
                [0, sun_europa_europa_unitvector[2] * 10],
                'r-')
        ax.text(sun_europa_europa_unitvector[0] * 10, sun_europa_europa_unitvector[2] * 10,
                'S', color='r')
        ax.plot([0, earth_europa_europa_unitvector[0] * 10],
                [0, earth_europa_europa_unitvector[2] * 10],
                'b-')
        ax.text(earth_europa_europa_unitvector[0] * 10, earth_europa_europa_unitvector[2] * 10,
                'E', color='b')


        ### Calculate the corotation flow.
        # Assume the corotation flow is axisymmetric flow wrt the rotation axis (System 3)
        # The flow vector is perpendicular to the position vector in x-y plane, without z component.
        speed_solid_corotation_sys3 = length(juice_sys3_jupiter) * rj * 2 * np.pi / 36000.   # 10 hour per rotation assumed
        vel_solid_corotation_sys3 = [-juice_sys3_jupiter[1], juice_sys3_jupiter[0], 0]
        vel_solid_corotation_sys3_unitvector = unit_vector(vel_solid_corotation_sys3)
        vel_solid_corotation_sys3 = vel_solid_corotation_sys3_unitvector * speed_solid_corotation_sys3
        print('Corotation/1st order', vel_solid_corotation_sys3, speed_solid_corotation_sys3)

        ### Geometric conversion from SystemIII to JSE. (Velocity measured at spacecraft -> geometric/position conversion,
        ### not velocity conversion)
        vel_solid_corotation_jse = cmat_sys3_jse.dot(vel_solid_corotation_sys3)

        ### Plasma flow at the spacecraft. Expressed in JSE frame.
        vel_solid_corotation_sys3 = vel_solid_corotation_jse - v_juice_jse_jupiter

        ### Plasma flow at the spacecraft. Convert it to SC frame.
        vel_solid_corotation_sc = juice_frame.convert_to_nsc(vel_solid_corotation_sys3)

        print('{0:%FT%T} {4:.5f} {1:.3f} {2[0]:.2f} {2[1]:.2f} {2[2]:.2f} {3[0]:.2f} {3[1]:.2f} {3[2]:.2f}'.format(
                t, et, vel_solid_corotation_sys3, vel_solid_corotation_sc, mpldates.date2num(t)), file=fp)
        print('{0:%FT%T} {4:.5f} {1:.3f} {2[0]:.2f} {2[1]:.2f} {2[2]:.2f} {3[0]:.2f} {3[1]:.2f} {3[2]:.2f}'.format(
                t, et, vel_solid_corotation_sys3, vel_solid_corotation_sc, mpldates.date2num(t)))



        ### Viewing direction axis
        ax = fig.add_axes([0.1, 0.08, 0.35, 0.20])

        ## Boundary of the viewing direction
        import fovplot

        fovplot.frame_plot_limit(ax=ax, alpha=0.1)
        ax.text(0, 0, 'Tilt=0', color='r')
        fovplot.frame_plot_limit(ax=ax, tilt=15, alpha=0.1)
        ax.text(0, 15, 'Tilt=15', color='r')
        fovplot.frame_plot_limit(ax=ax, tilt=30, alpha=0.1)
        ax.text(0, 30, 'Tilt=30', color='r')
        fovplot.frame_plot_limit(ax=ax, tilt=45, alpha=0.1)
        ax.text(0, 45, 'Tilt=45', color='r')

        ## Direction of the corotation flow in SC
        looking_vco_long = longitude(-vel_solid_corotation_sc)
        looking_vco_lat = latitude(-vel_solid_corotation_sc)

        ax.plot([looking_vco_long], [looking_vco_lat], 'bo')
        ax.text(looking_vco_long, looking_vco_lat, 'Corot', color='b')

        ## Galileans / Earth / Sun
        for num, nam, clr in zip((n_jp, n_sun, n_earth, n_io, n_europa, n_ganymede, n_callisto),
                                 ('Jupiter', 'Sun', 'Earth', 'Io', 'Europa', 'Ganymede', 'Callisto'),
                                 ('g', 'r', 'b', 'g', 'g', 'g', 'g')):
            looking_object_jse, __vel = get_state(num, et, 'JSE', -907)
            looking_object_sc = juice_frame.convert_to_nsc(looking_object_jse)
            looking_object_long = longitude(looking_object_sc)
            looking_object_lat = latitude(looking_object_sc)
            if np.abs(looking_object_lat) > 89.5:
                looking_object_long = 0.

            ax.plot([looking_object_long], [looking_object_lat], clr + 'o')
            ax.text(looking_object_long, looking_object_lat, nam, color=clr)

        ## Europa horizon
        distance = length(juice_europa_europa)    # Europa-S/C in re.
        angle = np.rad2deg(np.arcsin(1 / distance))   # Europa angle
        ax.axhline(90 - angle, color='g')
        print(angle)

        plt.gcf().savefig(os.path.join(pathname, 'scenario4_{:d}.png'.format(int(et))))

        plt.close('all')

        fp.flush()

if __name__ == '__main__':
    main()