{ "cells": [ { "cell_type": "markdown", "id": "621c0dc4-4ff6-4df3-afec-985716d0a324", "metadata": {}, "source": [ "# Generating Visibility Plots with Tudat" ] }, { "cell_type": "markdown", "id": "5f5df435-efe6-4dfe-a8e9-77ac35765033", "metadata": {}, "source": [ "# Objectives\n", "Whether you're looking to impress your friends by showing them when and where the **International Space Station** 🛰️ will pass overhead, or you need to verify if an **interplanetary probe** will be above a certain **elevation** for radio antenna communication 📡, this example has you covered. In this notebook, we'll demonstrate how to use **Tudatpy** to generate **visibility plots (azimuth and elevation)** for a specific target, from a selected location on Earth, across a given time period. In order to validate our results, we will compare them with results from the [JPL Horizons tool](https://ssd.jpl.nasa.gov/horizons/app.html#/). " ] }, { "cell_type": "markdown", "id": "e5ba7f21-37e9-496e-8a8a-b3c4fbcee345", "metadata": {}, "source": [ "## Import Relevant Modules\n", "Let's start by importing all relevant modules and libraries. Notably, we will need the [HorizonsQuery class](https://py.api.tudat.space/en/latest/horizons.html#tudatpy.data.horizons.HorizonsQuery), containing wrapper function to the [JPL Horizons tool](https://ssd.jpl.nasa.gov/horizons/app.html#/). We will also make use of Tudat's interface to the NASA's [spice tool](https://py.api.tudat.space/en/latest/spice.html). The Horizons tool provides access to high-precision ephemerides for solar system objects, offering position and velocity data. The SPICE tool, developed by NASA, enables the computation of spacecraft trajectory, orientation, and other mission-related data using detailed planetary and satellite models." ] }, { "cell_type": "code", "execution_count": 1, "id": "bc367d70", "metadata": {}, "outputs": [], "source": [ "# Load required tudatpy modules\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "from tudatpy.data.horizons import HorizonsQuery\n", "from tudatpy.interface import spice\n", "from tudatpy.astro import time_representation, element_conversion\n", "from tudatpy.math import interpolators\n", "from tudatpy.dynamics import environment, environment_setup\n", "from tudatpy.dynamics import propagation_setup, parameters_setup, simulator\n", "from tudatpy.estimation import observable_models_setup, observable_models, observations_setup, observations, estimation_analysis\n", "from tudatpy.astro.time_representation import DateTime\n", "\n", "from datetime import datetime\n", "import matplotlib.dates as mdates\n", "from itertools import zip_longest" ] }, { "cell_type": "markdown", "id": "7b1d93b4-e714-4b49-9681-fcf3d5eff612", "metadata": {}, "source": [ "## Run the Code and Get Your Results\n", "The two functions we set up for you are `plot_combined_elevation` and `read_state_vector`. \n", "\n", "***If you want to run the code and just get your visibility plot, please feel free to skip directly to the [Get Visibility Plot](#Get-Visibility-Plot) of this notebook.*** There, we show how this example might be turn out to be useful for **Planetary Defense Purposes**.\n", "\n", "- `plot_combined_elevation`: this is the core function of this notebook: if you're keen on understanding **how Tudat truly works**, we highly recommend **diving into its comments**. They are well-written and should help you make sense of the process. \n", "\n", "- `read_state_vector`: this function **reads state vector data from a file**, extracting time and state vector components based on the specified column indices. It's just a utility function you might use to parse a file containing ephemeris of a customized object, for which no ephemeris data is found on the Horizons database. The extracted ephemeris can then be used into Tudat." ] }, { "cell_type": "code", "execution_count": 2, "id": "02bd795d-329c-4d3a-9453-bb8524e8dfd3", "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "def plot_combined_elevation(\n", " target, \n", " station_names, \n", " start_epoch, \n", " end_epoch, \n", " time_step = '10m',\n", " time_format = '%Y-%m-%d %H:%M:%S',\n", " global_frame_origin = 'Earth', \n", " global_frame_orientation = 'J2000', \n", " geodetic_positions = None,\n", " custom_ephemeris = None):\n", "\n", " \"\"\"\n", " Plots the combined elevation of a target celestial body as seen from multiple ground stations.\n", "\n", " Parameters:\n", " -----------\n", " target : str\n", " Name of the celestial body being observed.\n", " station_names : list of str\n", " List of ground station names from which the target's elevation is to be computed.\n", " start_epoch : str\n", " Start epoch of the observation in the specified time format.\n", " end_epoch : str\n", " End epoch of the observation in the specified time format.\n", " time_step : str, optional (default='10m')\n", " Time step for ephemeris queries (e.g., '10m' for 10 minutes).\n", " time_format : str, optional (default='%Y-%m-%d %H:%M:%S')\n", " Format of the input time strings.\n", " global_frame_origin : str, optional (default='Earth')\n", " Origin of the global reference frame.\n", " global_frame_orientation : str, optional (default='J2000')\n", " Orientation of the global reference frame.\n", " geodetic_positions : dict, optional (default=None)\n", " Dictionary containing geodetic positions (altitude, latitude, longitude) for ground stations not present in Tudat's default list.\n", " custom_ephemeris : tudatpy.ephemeris.Ephemeris, optional (default=None)\n", " Custom ephemeris to be used for the target body instead of querying Horizons.\n", "\n", " Returns:\n", " --------\n", " None\n", " This function generates a plot of the elevation angles for the target as seen from the specified ground stations.\n", " \"\"\"\n", " \n", " # Load standard spice kernels. These are needed - for instance - to properly create the Tudat body_settings. \n", " spice.load_standard_kernels()\n", "\n", " # Define the list of bodies that will be created in the system model. In this case, we only need the Earth.\n", " bodies_to_create = [\"Earth\"]\n", "\n", " # Convert start and end epochs to Julian Day\n", " jd_start_epoch = time_representation.DateTime.from_python_datetime(datetime.strptime(start_epoch, time_format)).to_julian_day()\n", " jd_end_epoch = time_representation.DateTime.from_python_datetime(datetime.strptime(end_epoch, time_format)).to_julian_day()\n", " n_day_buffer = 1\n", "\n", " # Add a buffer to the user-defined start and end epochs.\n", " # This ensures that the simulation interpolator operates without errors nor warnings.\n", " # While the buffered epochs will differ from the original user-defined range, \n", " # the code will later filter out any dates outside the requested range. \n", " calendar_date_simulation_start_epoch = time_representation.DateTime.from_julian_day(jd_start_epoch - n_day_buffer).to_python_datetime()\n", " calendar_date_simulation_end_epoch = time_representation.DateTime.from_julian_day(jd_end_epoch + n_day_buffer).to_python_datetime()\n", "\n", "\n", " # Convert the start and end epochs to seconds since the epoch for simulation. \n", " # This conversion is needed, as the 'get_default_body_settings_time_limited' function\n", " # (as well as other Tudat functions) - which we will use later on -\n", " # only accept floats (in seconds from J2000) as arguments. \n", " simulation_seconds_start_epoch = time_representation.julian_day_to_seconds_since_epoch(jd_start_epoch - n_day_buffer)\n", " simulation_seconds_end_epoch = time_representation.julian_day_to_seconds_since_epoch(jd_end_epoch + n_day_buffer)\n", "\n", " # Actual (user-queried) start epoch. Later on, we will filter our results based on this epoch.\n", " actual_seconds_start_epoch = time_representation.julian_day_to_seconds_since_epoch(jd_start_epoch)\n", " actual_seconds_end_epoch = time_representation.julian_day_to_seconds_since_epoch(jd_end_epoch)\n", "\n", " # Create default Earth and target settings.\n", " # First, we use 'get_default_body_settings_time_limited' to retrieve the default settings \n", " #for the given set of input bodies, with a limited valid time interval.\n", " # See [Tudatpy API Reference](https://py.api.tudat.space/en/latest/environment_setup.html#tudatpy.dynamics.environment_setup.get_default_body_settings_time_limited)\n", " # Please note that, although it is not necessary to create time limited settings function, the code will run faster if you use it.\n", " body_settings = environment_setup.get_default_body_settings_time_limited(\n", " bodies_to_create, simulation_seconds_start_epoch, simulation_seconds_end_epoch, global_frame_origin, global_frame_orientation)\n", " body_settings.add_empty_settings(target)\n", "\n", " # Add Earth's shape settings. \n", " # We go for an oblate spherical shape, using a radius of 6378 km and the current accepted flattening value\n", " equatorial_radius = 6378*1e3 # in meters\n", " flattening = 1/298\n", " body_settings.get('Earth').shape_settings = environment_setup.shape.oblate_spherical(\n", " equatorial_radius = equatorial_radius,\n", " flattening = flattening,\n", " )\n", " \n", " # Add Earth's rotational settings via the 'environment_setup.rotation_model.gcrs_to_itrs'.\n", " # This function defines high-accuracy Earth rotation model according to the IERS Conventions 2010. \n", " # The model computes the rotation from ITRS to GCRS.\n", " # See [Tudatpy API Reference](https://py.api.tudat.space/en/latest/rotation_model.html#tudatpy.dynamics.environment_setup.rotation_model.gcrs_to_itrs)\n", " body_settings.get('Earth').rotation_model_settings = environment_setup.rotation_model.gcrs_to_itrs(\n", " environment_setup.rotation_model.iau_2006, global_frame_orientation,\n", " interpolators.interpolator_generation_settings(interpolators.cubic_spline_interpolation(),\n", " simulation_seconds_start_epoch, simulation_seconds_end_epoch, 60),\n", " interpolators.interpolator_generation_settings(interpolators.cubic_spline_interpolation(),\n", " simulation_seconds_start_epoch, simulation_seconds_end_epoch, 60),\n", " interpolators.interpolator_generation_settings(interpolators.cubic_spline_interpolation(),\n", " simulation_seconds_start_epoch, simulation_seconds_end_epoch, 60))\n", "\n", " # As you've already seen in the documentation of this function, users can input custom_ephemeris via the custom_ephemeris flag.\n", " # This is useful for objects that are not present on the Horizons database, \n", " # or for any other custom/fake object for which an ephemeris file is made available by the user.\n", " # If the custom_ephemeris flag is not on, then we create ephemeris by performing an Horizons Query; if the custom_ephemeris flag is on, \n", " # the ephemeris settings are set accordingly as the ones provided by the user.\n", "\n", " if not custom_ephemeris:\n", " # Horizons query to retrieve ephemeris.\n", " query_ephemerides = HorizonsQuery(\n", " query_id=target,\n", " location = '@399',\n", " epoch_start=calendar_date_simulation_start_epoch,\n", " epoch_end=calendar_date_simulation_end_epoch,\n", " epoch_step= time_step\n", " )\n", "\n", " # Allows for ephemeris retrieval\n", " horizons_ephemeris = query_ephemerides.create_ephemeris_tabulated(\n", " frame_origin=global_frame_origin,\n", " frame_orientation= global_frame_orientation\n", " )\n", "\n", " # Set the fetched ephemeris as the target's ephemeris\n", " body_settings.get(target).ephemeris_settings = horizons_ephemeris\n", "\n", " else:\n", " # Use the custom ephemeris provided by the user\n", " body_settings.get(target).ephemeris_settings = custom_ephemeris\n", "\n", " # Tudat comes with the definition of a bunch of relevant antenna locations spread around the world. \n", " # These include - but are not limited to: DSN and EStrack Antennas, EVN/VLBA Antennas. \n", " # The function 'radio_telescope_stations()' of the 'environment_setup.ground_station' class \n", " # allows to add all Tudat radio stations in the Tudat environment. \n", " # See [Tudatpy API Reference](https://py.api.tudat.space/en/latest/ground_station.html#tudatpy.dynamics.environment_setup.ground_station.radio_telescope_stations)\n", " all_tudat_radio_telescopes_settings = environment_setup.ground_station.radio_telescope_stations()\n", " all_tudat_radio_telescopes_names = [telescope.station_name for telescope in all_tudat_radio_telescopes_settings]\n", "\n", " # Initialize the (sub)plots.\n", " fig, axes = plt.subplots(len(station_names), 1, figsize=(17, 15))\n", " fig.tight_layout(pad=5) # Adjust padding between subplots\n", "\n", " # As we already mentioned in this Notebook's Objectives section, users can select their preferred location on Earth \n", " # as observation location. The previous ('radio_telescope_stations') function makes Tudat aware of all Tudat stations locations\n", " # (in cartesian coordinates, on Earth). If a user-requested Antenna is not found among the tudat list of radio telescope stations, \n", " # it is considered as \"new\" and it is added to environemnt via the function 'environment_setup.ground_station.basic_station'.\n", " # See [Tudatpy API Reference](https://py.api.tudat.space/en/latest/ground_station.html#tudatpy.dynamics.environment_setup.ground_station.basic_station)\n", " existing_stations = [station_name for station_name in station_names if station_name in all_tudat_radio_telescopes_names]\n", " new_stations = [station_name for station_name in station_names if station_name not in all_tudat_radio_telescopes_names]\n", "\n", " new_ground_stations_settings = list()\n", " geodetic_stations_dict = dict()\n", " horizons_coordinates_dict = dict()\n", "\n", " # The following code snippet allows to create or reset settings for station geodetic coordinates, and to create the geodetic coordinate dictionary \n", " # required by Horizons (with a given order ['lon', 'lat', 'elevation'] and in given units).\n", " # Three cases are consdiered:\n", " # 1) The considered station is part of the Tudat list of stations, and it is not part of the geodetic_positions user-defined dictionary\n", " # In this case, the geodetic position is retrieved from the Tudat-defined cartesian position via the 'convert_cartesian_to_geodetic_coordinates' function.\n", " # See:[TUDATPY API REFERENCE](https://py.api.tudat.space/en/latest/element_conversion.html#tudatpy.astro.element_conversion.convert_cartesian_to_geodetic_coordinates)\n", " # 2) The considered station is part of the Tudat list of stations, and it is also part of the geodetic_positions user-defined dictionary.\n", " # In this case, the Tudat-default station location is replaced by the cartesian position corresponding to the user-defined geodetic position\n", " # This is done via the 'convert_position_elements' function. \n", " # See: https://py.api.tudat.space/en/latest/element_conversion.html#tudatpy.astro.element_conversion.convert_position_elements\n", " # 3) The considered station is not in the Tudat list of stations. In this case, it is added via the 'environment_setup.ground_station.basic_station' function.\n", " # See: https://py.api.tudat.space/en/latest/ground_station.html#tudatpy.dynamics.environment_setup.ground_station.basic_station\n", " # Please note the deg2rad and rad2deg conversions throughout this snippet, as well as the order for the station location coordinates. \n", " # (Tudat order: elevation, latitude, longitude. Horizons order: longitude, latitude, elevation). \n", " if len(existing_stations) > 0:\n", " for existing_station in existing_stations:\n", " index = next(\n", " (i for i, telescope in enumerate(all_tudat_radio_telescopes_names) if existing_station == telescope),\n", " None\n", " )\n", " if geodetic_positions and existing_station in geodetic_positions.keys(): #Allow reset of ground station position for an already existing Tudat station\n", " print(f\"Stations: {existing_station} already in Tudat's list of stations. Resetting its geodetic position.\")\n", " geodetic_position = [\n", " geodetic_positions[existing_station][0], # in meters\n", " np.deg2rad(geodetic_positions[existing_station][1]), # latitude from radians to degrees\n", " np.deg2rad(geodetic_positions[existing_station][2]), # longitude from radians to degrees\n", " ]\n", "\n", " # Create a temporary bodies object, just to retrieve the shape model.\n", " # The shape model is required to allow for conversion from geodetic to cartesian elements.\n", " # This is required, since all_tudat_radio_telescopes_settings[index].station_position\n", " # must have cartesian coordinates.\n", " temp_bodies = environment_setup.create_system_of_bodies(body_settings)\n", " all_tudat_radio_telescopes_settings[index].station_position = element_conversion.convert_position_elements(\n", " geodetic_position,\n", " element_conversion.geodetic_position_type,\n", " element_conversion.cartesian_position_type,\n", " temp_bodies.get('Earth').shape_model,\n", " 10\n", " )\n", "\n", " horizons_coordinates_dict[existing_station] = {\n", " 'lon': geodetic_positions[existing_station][2], # in degrees\n", " 'lat': geodetic_positions[existing_station][1], # in degrees\n", " 'elevation': geodetic_positions[existing_station][0]/1000, # in km\n", " }\n", "\n", " else:\n", " radians_geodetic_position = element_conversion.convert_cartesian_to_geodetic_coordinates(\n", " all_tudat_radio_telescopes_settings[index].station_position,\n", " equatorial_radius = equatorial_radius,\n", " flattening = flattening,\n", " tolerance = 10\n", " )\n", "\n", " geodetic_position = [\n", " radians_geodetic_position[0], # in meters\n", " np.rad2deg(radians_geodetic_position[1]), # latitude from radians to degrees\n", " np.rad2deg(radians_geodetic_position[2]), # longitude from radians to degrees\n", " ]\n", " geodetic_stations_dict[existing_station] = geodetic_position\n", "\n", " horizons_coordinates_dict[existing_station] = {\n", " 'lon': geodetic_position[2], # in degrees\n", " 'lat': geodetic_position[1], # in degrees\n", " 'elevation': geodetic_position[0]/1000, # in km\n", " }\n", "\n", " if len(new_stations) > 0:\n", " for new_station in new_stations:\n", " if not geodetic_positions:\n", " print(f'Station: {new_station} not found in Tudat list of stations, and no custom geodetic position found in the geodetic_position dictionary.\\n'\n", " f'Please provide it.\\n Aborting...')\n", " exit()\n", " elif geodetic_positions and not geodetic_positions[new_station]:\n", " print(f'Geodetic Position for {new_station} not found in the geodetic_position dictionary. Please provide it.\\n Aborting...')\n", " exit()\n", " if new_station not in geodetic_positions.keys():\n", " print(f'Station Location for station: {new_station} not found in Tudat list of stations, '\n", " f'nor in ground_station_positions. Please provide a valid geodetic location.'\n", " )\n", " exit()\n", "\n", " else:\n", " geodetic_position = geodetic_positions[new_station]\n", " new_ground_stations_settings.append(environment_setup.ground_station.basic_station(\n", " new_station,\n", " [geodetic_position[0], # in meters\n", " np.deg2rad(geodetic_position[1]), # latitude in radians\n", " np.deg2rad(geodetic_position[2]) # longitude in radians\n", " ],\n", " element_conversion.geodetic_position_type\n", " ))\n", "\n", " geodetic_stations_dict[new_station] = geodetic_position\n", "\n", " horizons_coordinates_dict[new_station] = {\n", " 'lon': geodetic_position[2], # in degrees\n", " 'lat': geodetic_position[1], # in degrees\n", " 'elevation': geodetic_position[0]/1000, # in km\n", " }\n", "\n", " # In the end, all ground station settings are put together in a list. \n", " if not new_ground_stations_settings:\n", " body_settings.get('Earth').ground_station_settings = all_tudat_radio_telescopes_settings\n", " else:\n", " body_settings.get('Earth').ground_station_settings = all_tudat_radio_telescopes_settings + new_ground_stations_settings\n", "\n", " # In order to visualize the plot better, we also decide to record the rise and set times (if any) over the user-specified timespan.\n", " # These are defined as the points in time where the target rises above or sets below the local horizon, on a given topocentric location.\n", " rise_set_times_dict = dict()\n", " for idx, station_name in enumerate(station_names):\n", " horizons_coord = horizons_coordinates_dict[station_name]\n", "\n", " # The following query to Horizons and subsequent application of the Tudatpy 'interpolated_station_angles' method \n", " # allow to retrieve inteprolated azimuths and elevations.\n", " # This call to horizons differs from the ones we have done before, in that it does not provide Ephemeris, \n", " # but rather a list of azimuths and elevations (observables). \n", " if not custom_ephemeris:\n", " query_az_el = HorizonsQuery(\n", " query_id=target,\n", " location=horizons_coord,\n", " epoch_start=calendar_date_simulation_start_epoch,\n", " epoch_end=calendar_date_simulation_end_epoch,\n", " epoch_step= time_step\n", " )\n", " \n", " horizons_antenna_az_el = query_az_el.interpolated_station_angles(degrees = True)\n", " start_end_condition = np.logical_and(\n", " horizons_antenna_az_el[:, 0] >= actual_seconds_start_epoch,\n", " horizons_antenna_az_el[:, 0] <= actual_seconds_end_epoch\n", " ) \n", "\n", " horizons_azimuth = horizons_antenna_az_el[:,1][start_end_condition]\n", " horizons_elevation = horizons_antenna_az_el[:,2][start_end_condition]\n", "\n", " # The 0th column of horizons_antenna_az_el represents the queried times (in seconds).\n", " # Earlier, we added a buffer to the user-defined start and end epochs in order for our simulation to run smoothly.\n", " # Notice that we could have added a buffer to the start_epoch only, but in that case we would have gotten some \n", " # extrapolation warnings related to the interpolator. Although these do not compromise the final results, \n", " # it is good to be aware of it.\n", " # This is the moment where we filter out the unwanted data, only keeping the times >= actual_seconds_start_epoch\n", " horizons_seconds_ephemeris = horizons_antenna_az_el[:,0][start_end_condition]\n", "\n", " # Apply conversions for convenience\n", " horizons_calendar_ephemeris = [time_representation.seconds_since_epoch_to_julian_day(horizons_second_ephemeris) for\n", " horizons_second_ephemeris in horizons_seconds_ephemeris]\n", " horizons_datetime_times = [time_representation.DateTime.from_julian_day(horizons_calendar_time).to_python_datetime() for\n", " horizons_calendar_time in horizons_calendar_ephemeris]\n", "\n", " # We would like the times at which the elevation and azimuth is computed in Tudat to be the same as the Horizons ones\n", " # This is done for ease of comparison. \n", " tudat_seconds_ephemeris = horizons_seconds_ephemeris\n", "\n", " else:\n", " # If the custom_ephemeris flag is on, the times have to be retrieved from the provided ephemeris file. \n", " # Moreover, in this case, we cannot make a call to Horizons, hence we will only get the Tudat elevation and azimuth plots.\n", " tudat_seconds_ephemeris = np.linspace(simulation_seconds_start_epoch, simulation_seconds_end_epoch, int(time_step.strip('m'))*10)\n", " tudat_seconds_ephemeris = tudat_seconds_ephemeris[start_end_condition]\n", "\n", " tudat_calendar_ephemeris = [time_representation.seconds_since_epoch_to_julian_day(tudat_second_ephemeris) for tudat_second_ephemeris in tudat_seconds_ephemeris]\n", " tudat_datetime_times = [time_representation.DateTime.from_julian_day(tudat_calendar_time).to_python_datetime() for tudat_calendar_time in tudat_calendar_ephemeris]\n", "\n", "\n", " # We struggled quite a bit to set all the body settings, but now we can finally create our bodies object.\n", " bodies = environment_setup.create_system_of_bodies(body_settings)\n", "\n", " # Notice how - so far - we have not run any simulation yet. \n", " # As we approach doing that, it is time to create the link ends for the simulation. \n", " # The given antenna is set as a receiver through the function: 'body_reference_point_link_end_id' of the observation module. \n", " # (See: https://py.api.tudat.space/en/latest/estimation/observable_models_setup/links.html#tudatpy.estimation.observable_models_setup.links.body_reference_point_link_end_id)\n", " link_ends = {\n", " observable_models_setup.links.receiver: observable_models_setup.links.body_reference_point_link_end_id('Earth', station_name),\n", " }\n", " # Create a single link definition from the link ends\n", " link_definition = observable_models_setup.links.LinkDefinition(link_ends)\n", "\n", " # Finally, the Tudat function 'compute_target_angles_and_range' can be used to simulate the observed azimuth and elevation \n", " # from the given receiver, at the given observation times (notice how we filter tudat_seconds_ephemeris by start_end_condition)\n", " # as we only want to cover the user-defined timespan. \n", " station_id = ('Earth', station_name)\n", " angles_range_dictionary = observations.observations_geometry.compute_target_angles_and_range(\n", " bodies,\n", " station_id,\n", " target,\n", " observation_times=tudat_seconds_ephemeris,\n", " is_station_transmitting=False\n", " )\n", "\n", " # Initialize lists for azimuth and elevation\n", " tudat_azimuth_list = []\n", " tudat_elevation_list = []\n", "\n", " # The following code snippet allows to compute the rise and set times, as seen from each station, \n", " # and to populate the rise_set_dict dictionary. \n", " set_times = []\n", " rise_times = []\n", " keys_list = sorted(angles_range_dictionary.keys())\n", "\n", " initial_elevation = np.rad2deg(angles_range_dictionary[keys_list[0]][0])\n", " flag = initial_elevation <= 0 # True if initially below horizon\n", "\n", " for idx_key, key in enumerate(keys_list): # Ensure iteration follows sorted order\n", " azimuth_deg = np.rad2deg(angles_range_dictionary[key][1]) % 360\n", " elevation_deg = (np.rad2deg(angles_range_dictionary[key][0]) + 90) % 180 - 90 # Normalize elevation\n", " tudat_azimuth_list.append(azimuth_deg)\n", " tudat_elevation_list.append(elevation_deg)\n", "\n", " if elevation_deg > 0: # Object above horizon\n", " if flag: # Transition from below horizon to above\n", " rise_times.append(tudat_datetime_times[idx_key])\n", " flag = False # Update flag\n", " else: # Object below horizon\n", " if not flag: # Transition from above horizon to below\n", " set_times.append(tudat_datetime_times[idx_key])\n", " flag = True # Update flag\n", "\n", " rise_set_times_dict[station_name] = [rise_times, set_times]\n", " rise_times = rise_set_times_dict[station_name][0]\n", " set_times = rise_set_times_dict[station_name][1]\n", "\n", " # Convert to numpy arrays\n", " azimuth_array = np.array(tudat_azimuth_list)\n", " elevation_array = np.array(tudat_elevation_list)\n", "\n", " # The last bit of this Jupyter Notebook deals with visualization of the obtained data. \n", " # Antennas observing at low elevation quickly lose performance. Moreover, if the elvation is negative, the target is not observable.\n", " # That is why we plot horizontal lines at 15 degrees and 0 degrees.\n", " if len(station_names) == 1:\n", " if not custom_ephemeris:\n", " axes.scatter(np.array(horizons_datetime_times)[horizons_elevation > 15], horizons_elevation[horizons_elevation > 15], label=f\"Horizons Elevation\",\n", " s=10, c='pink', alpha=1)\n", " axes.scatter(np.array(horizons_datetime_times)[horizons_elevation > 15], horizons_azimuth[horizons_elevation > 15], label=f\"Horizons Azimuth\",\n", " s=10, c='green', alpha=1)\n", " axes.scatter(np.array(horizons_datetime_times)[horizons_elevation <= 15], horizons_elevation[horizons_elevation <= 15],\n", " s=7, c='pink', alpha=0.3)\n", " axes.scatter(np.array(horizons_datetime_times)[horizons_elevation <= 15], horizons_azimuth[horizons_elevation <= 15],\n", " s=7, c='green', alpha=0.3)\n", " axes.scatter(np.array(horizons_datetime_times)[horizons_elevation <= 0], horizons_elevation[horizons_elevation <= 0],\n", " s=5, c='pink', alpha=0.1)\n", " axes.scatter(np.array(horizons_datetime_times)[horizons_elevation <= 0], horizons_azimuth[horizons_elevation <= 0],\n", " s=5, c='green', alpha=0.1)\n", "\n", " axes.scatter(np.array(tudat_datetime_times)[elevation_array > 15], elevation_array[elevation_array > 15],\n", " s=5, label='Tudat Elevation', c='blue')\n", " axes.scatter(np.array(tudat_datetime_times)[elevation_array > 15], azimuth_array[elevation_array > 15],\n", " s=5, label='Tudat Azimuth', c='red')\n", " axes.scatter(np.array(tudat_datetime_times)[elevation_array <= 15], elevation_array[elevation_array <= 15],\n", " s=3, c='blue', alpha=0.5)\n", " axes.scatter(np.array(tudat_datetime_times)[elevation_array <= 15], azimuth_array[elevation_array <= 15],\n", " s=3, c='red', alpha=0.5)\n", " axes.scatter(np.array(tudat_datetime_times)[elevation_array <= 0], elevation_array[elevation_array <= 0],\n", " s=1, c='blue', alpha=0.1)\n", " axes.scatter(np.array(tudat_datetime_times)[elevation_array <= 0], azimuth_array[elevation_array <= 0],\n", " s=1, c='red', alpha=0.1)\n", "\n", " # Add a horizontal line for elevation = 15\n", " axes.axhline(y=15, color='k', linestyle='--', alpha=0.5, label = '15 deg')\n", " axes.axhline(y=0, color='green', linestyle='--', alpha=0.5, label = '0 deg')\n", "\n", " for i in range(len(set_times)):\n", " if i == 0:\n", " axes.axvline(x = set_times[i], color='grey', linestyle='--', alpha=0.5, label = 'Set times')\n", " else:\n", " axes.axvline(x = set_times[i], color='grey', linestyle='--', alpha=0.5)\n", "\n", " for i in range(len(rise_times)):\n", " if i == 0:\n", " axes.axvline(x = rise_times[i], color='purple', linestyle='--', alpha=0.5, label = 'Rise times')\n", " else:\n", " axes.axvline(x = rise_times[i], color='purple', linestyle='--', alpha=0.5)\n", "\n", " # Add labels, title, and legend\n", " axes.set_xlabel('UTC Time')\n", " axes.set_ylabel('Elevation/Azimuth (Degrees)')\n", " axes.set_title(f'{target} from {station_name}')\n", " axes.legend(loc='upper right')\n", "\n", " # Use zip_longest to pair rise_times and set_times, filling missing values with None\n", " for rise_time, set_time in zip_longest(rise_times, set_times, fillvalue=None):\n", " \n", " if not set_time and rise_time: # If only rise_time exists\n", " print(f'Rise from {station_name}:', rise_time)\n", " elif not rise_time and set_time: # If only set_time exists\n", " print(f'Set from {station_name}:', set_time)\n", " elif not rise_time and not set_time: # If neither exists\n", " print(f'No rise nor set for {station_name} in the selected timespan.')\n", " else: # If both exist\n", " print(f'Rise from {station_name}:', rise_time, f'Set from {station_name}:', set_time)\n", "\n", " # Customize x-axis ticks (datetime formatting and positioning)\n", " axes.xaxis.set_major_formatter(mdates.DateFormatter('%y-%m-%d\\n%H:%M'))\n", " axes.xaxis.set_major_locator(mdates.AutoDateLocator(minticks=30, maxticks=30))\n", " axes.grid()\n", " axes.tick_params(axis='x', labelsize=7) # Adjust 'labelsize' to the desired font size\n", "\n", " else:\n", " # Horizons will be able to give an Elevation plot only when custom_ephemeris == False or None\n", " if not custom_ephemeris: \n", " axes[idx].scatter(np.array(horizons_datetime_times)[horizons_elevation > 15], horizons_elevation[horizons_elevation > 15], label=f\"Horizons Elevation\",\n", " s=20, c='pink', alpha=0.5, marker = 's')\n", " axes[idx].scatter(np.array(horizons_datetime_times)[horizons_elevation > 15], horizons_azimuth[horizons_elevation > 15], label=f\"Horizons Azimuth\",\n", " s=20, c='green', alpha=0.5, marker = 's')\n", " axes[idx].scatter(np.array(horizons_datetime_times)[horizons_elevation <= 15], horizons_elevation[horizons_elevation <= 15],\n", " s=15, c='pink', alpha=0.3, marker = 's')\n", " axes[idx].scatter(np.array(horizons_datetime_times)[horizons_elevation <= 15], horizons_azimuth[horizons_elevation <= 15],\n", " s=15, c='green', alpha=0.3, marker = 's')\n", " axes[idx].scatter(np.array(horizons_datetime_times)[horizons_elevation <= 0], horizons_elevation[horizons_elevation <= 0],\n", " s=10, c='pink', alpha=0.1, marker = 's')\n", " axes[idx].scatter(np.array(horizons_datetime_times)[horizons_elevation <= 0], horizons_azimuth[horizons_elevation <= 0],\n", " s=10, c='green', alpha=0.1, marker = 's')\n", " \n", " axes[idx].scatter(np.array(tudat_datetime_times)[elevation_array > 15], elevation_array[elevation_array > 15],\n", " s=5, label='Tudat Elevation', c='blue')\n", " axes[idx].scatter(np.array(tudat_datetime_times)[elevation_array > 15], azimuth_array[elevation_array > 15],\n", " s=5, label='Tudat Azimuth', c='black')\n", " axes[idx].scatter(np.array(tudat_datetime_times)[elevation_array <= 15], elevation_array[elevation_array <= 15],\n", " s=3, c='blue', alpha=0.5)\n", " axes[idx].scatter(np.array(tudat_datetime_times)[elevation_array <= 15], azimuth_array[elevation_array <= 15],\n", " s=3, c='black', alpha=0.5)\n", " axes[idx].scatter(np.array(tudat_datetime_times)[elevation_array <= 0], elevation_array[elevation_array <= 0],\n", " s=1, c='blue', alpha=0.1)\n", " axes[idx].scatter(np.array(tudat_datetime_times)[elevation_array <= 0], azimuth_array[elevation_array <= 0],\n", " s=1, c='black', alpha=0.1)\n", "\n", " # Add a horizontal line for elevation = 15\n", " axes[idx].axhline(y=15, color='k', linestyle='--', alpha=0.5, label = '15 deg')\n", " axes[idx].axhline(y=0, color='green', linestyle='--', alpha=0.5, label = '0 deg')\n", "\n", " for i in range(len(set_times)):\n", " if i == 0:\n", " axes[idx].axvline(x = set_times[i], color='grey', linestyle='--', alpha=0.5, label = 'Set times')\n", " else:\n", " axes[idx].axvline(x = set_times[i], color='grey', linestyle='--', alpha=0.5)\n", "\n", " for i in range(len(rise_times)):\n", " if i == 0:\n", " axes[idx].axvline(x = rise_times[i], color='purple', linestyle='--', alpha=0.5, label = 'Rise times')\n", " else:\n", " axes[idx].axvline(x = rise_times[i], color='purple', linestyle='--', alpha=0.5)\n", "\n", " # Add labels, title, and legend\n", " axes[idx].set_xlabel('UTC Time')\n", " axes[idx].set_ylabel('Elevation/Azimuth (Degrees)')\n", " axes[idx].set_title(f'{target} Az/El from {station_name}')\n", " axes[idx].legend(loc='upper right')\n", "\n", " # Use zip_longest to pair rise_times and set_times, filling missing values with None\n", " for rise_time, set_time in zip_longest(rise_times, set_times, fillvalue=None):\n", " \n", " if not set_time and rise_time: # If only rise_time exists\n", " print(f'Rise from {station_name}:', rise_time)\n", " elif not rise_time and set_time: # If only set_time exists\n", " print(f'Set from {station_name}:', set_time)\n", " elif not rise_time and not set_time: # If neither exists\n", " print(f'No rise nor set for {station_name} in the selected timespan.')\n", " else: # If both exist\n", " print(f'Rise from {station_name}:', rise_time, f'Set from {station_name}:', set_time)\n", " \n", " # Customize x-axis ticks (datetime formatting and positioning)\n", " axes[idx].xaxis.set_major_formatter(mdates.DateFormatter('%y-%m-%d\\n%H:%M'))\n", " axes[idx].xaxis.set_major_locator(mdates.AutoDateLocator(minticks=10, maxticks=20))\n", " axes[idx].grid()\n", " axes[idx].tick_params(axis='x', labelsize=7) # Adjust 'labelsize' to the desired font size\n", "\n", "\n", " # Show the combined plot\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": 3, "id": "fcc50d5c-b376-421b-aaa1-47c09032ff2c", "metadata": {}, "outputs": [], "source": [ "def read_state_vector(file_path, time_col=0, state_cols=(1, 2, 3, 4, 5, 6)):\n", " \"\"\"\n", " Reads state vector from a text file with an optional selection of time and state columns.\n", "\n", " Parameters:\n", " file_path (str): Path to the file.\n", " time_col (int): Column index for time (default: 0).\n", " state_cols (tuple): Column indices for x, y, z, vx, vy, vz (default: (1, 2, 3, 4, 5, 6)).\n", "\n", " Returns:\n", " times (list): List of time values as strings.\n", " state_vectors (list): List of NumPy arrays containing state vectors.\n", " \"\"\"\n", " # Initialize empty lists to store state vectors and their corresponding timestamps\n", " state_vectors = []\n", " times = []\n", " \n", " # Open the file at the specified path in read mode\n", " with open(file_path, 'r') as f:\n", " # Loop through each line in the file\n", " for line in f:\n", " # Split the line into components (assumes space-separated values)\n", " parts = line.strip().split()\n", " \n", " # Skip lines that don't have enough columns to extract the required data\n", " if len(parts) <= max(state_cols):\n", " continue\n", " \n", " try:\n", " # Extract the timestamp from the specified column (time_col)\n", " time_str = parts[time_col]\n", " \n", " # Extract the state vector by picking values from the specified state_cols and converting them to floats\n", " state_vector = np.array([float(parts[i]) for i in state_cols])\n", " \n", " # Append the extracted timestamp to the list of times\n", " times.append(time_str)\n", " \n", " # Append the extracted state vector to the list of state vectors\n", " state_vectors.append(state_vector)\n", " except ValueError:\n", " # Skip lines that contain invalid (non-numeric) data that can't be converted to floats\n", " continue\n", " \n", " # Return the list of timestamps and state vectors as output\n", " return times, state_vectors" ] }, { "cell_type": "markdown", "id": "f26f3afc-c89a-46ba-853f-8454c218f134", "metadata": {}, "source": [ "## Get Visibility Plot\n", "Imagine we are engaged in planetary defense efforts and aim to determine whether the University of Tasmania's radio antennas (Australia) could conduct radar observations of the Potentially Hazardous Asteroid 2014 HK129 on the date of its closest approach, which occurred on December 20, 2022. We select the following stations (names are taken from [this list of radio telescope stations](https://gitlab.com/gofrito/pysctrack/-/blob/master/cats/glo.sit?ref_type=heads)):\n", "\n", "* Tidbinbilla (64m)\n", "* Katherine (12m)\n", "* Yarragadee (12m)\n", "* Hobart (26m)\n", "\n", "From the visibility plots, we can see that the asteroid is **initially below the horizon for all locations**, and it rises up some hours later.\n", "\n", "*NOTE: Notice how the **custom_ephemeris flag is set to None**. In the case where you had `your_own_ephemeris` file for 2024 HK129 (or for any other object, you could set that flag to the dictionary created by reading the ephemeris file, through the function `read_state_vector` (or a similar home-made function)*. Try it yourself!" ] }, { "cell_type": "code", "execution_count": 4, "id": "6f1ecd96", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Warning, the function calendar_date_to_julian_day is deprecated, and will be removed in a future version. Please use the functionally equivalent DateTime.from_python_datetime(...).to_julian_day() instead\n", "Warning, the function calendar_date_to_julian_day is deprecated, and will be removed in a future version. Please use the functionally equivalent DateTime.from_python_datetime(...).to_julian_day() instead\n", "Warning, the function julian_day_to_calendar_date is deprecated, and will be removed in a future version. Please use the functionally equivalent DateTime.from_julian_day(...).to_python_datetime() instead\n", "Warning, the function julian_day_to_calendar_date is deprecated, and will be removed in a future version. Please use the functionally equivalent DateTime.from_julian_day(...).to_python_datetime() instead\n" ] }, { "ename": "AttributeError", "evalue": "module 'tudatpy.kernel.math.interpolators' has no attribute 'interpolator_generation_settings_float'", "output_type": "error", "traceback": [ "\u001B[31m---------------------------------------------------------------------------\u001B[39m", "\u001B[31mAttributeError\u001B[39m Traceback (most recent call last)", "\u001B[36mCell\u001B[39m\u001B[36m \u001B[39m\u001B[32mIn[4]\u001B[39m\u001B[32m, line 17\u001B[39m\n\u001B[32m 14\u001B[39m start_epoch = \u001B[33m'\u001B[39m\u001B[33m2022-12-19 00:00:00\u001B[39m\u001B[33m'\u001B[39m\n\u001B[32m 15\u001B[39m end_epoch = \u001B[33m'\u001B[39m\u001B[33m2022-12-20 00:00:00\u001B[39m\u001B[33m'\u001B[39m\n\u001B[32m---> \u001B[39m\u001B[32m17\u001B[39m \u001B[43mplot_combined_elevation\u001B[49m\u001B[43m(\u001B[49m\n\u001B[32m 18\u001B[39m \u001B[43m \u001B[49m\u001B[33;43m'\u001B[39;49m\u001B[33;43m2014 HK129\u001B[39;49m\u001B[33;43m'\u001B[39;49m\u001B[43m,\u001B[49m\n\u001B[32m 19\u001B[39m \u001B[43m \u001B[49m\u001B[43mstation_names\u001B[49m\u001B[43m,\u001B[49m\n\u001B[32m 20\u001B[39m \u001B[43m \u001B[49m\u001B[43mstart_epoch\u001B[49m\u001B[43m \u001B[49m\u001B[43m=\u001B[49m\u001B[43m \u001B[49m\u001B[43mstart_epoch\u001B[49m\u001B[43m,\u001B[49m\n\u001B[32m 21\u001B[39m \u001B[43m \u001B[49m\u001B[43mend_epoch\u001B[49m\u001B[43m \u001B[49m\u001B[43m=\u001B[49m\u001B[43m \u001B[49m\u001B[43mend_epoch\u001B[49m\u001B[43m,\u001B[49m\n\u001B[32m 22\u001B[39m \u001B[43m \u001B[49m\u001B[43mtime_step\u001B[49m\u001B[43m=\u001B[49m\u001B[43m \u001B[49m\u001B[33;43m'\u001B[39;49m\u001B[33;43m20m\u001B[39;49m\u001B[33;43m'\u001B[39;49m\u001B[43m,\u001B[49m\n\u001B[32m 23\u001B[39m \u001B[43m \u001B[49m\u001B[43mgeodetic_positions\u001B[49m\u001B[43m=\u001B[49m\u001B[43mgeodetic_positions\u001B[49m\u001B[43m,\u001B[49m\n\u001B[32m 24\u001B[39m \u001B[43m \u001B[49m\u001B[43mcustom_ephemeris\u001B[49m\u001B[43m \u001B[49m\u001B[43m=\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;28;43;01mNone\u001B[39;49;00m\n\u001B[32m 25\u001B[39m \u001B[43m)\u001B[49m\n", "\u001B[36mCell\u001B[39m\u001B[36m \u001B[39m\u001B[32mIn[2]\u001B[39m\u001B[32m, line 98\u001B[39m, in \u001B[36mplot_combined_elevation\u001B[39m\u001B[34m(target, station_names, start_epoch, end_epoch, time_step, time_format, global_frame_origin, global_frame_orientation, geodetic_positions, custom_ephemeris)\u001B[39m\n\u001B[32m 87\u001B[39m body_settings.get(\u001B[33m'\u001B[39m\u001B[33mEarth\u001B[39m\u001B[33m'\u001B[39m).shape_settings = environment_setup.shape.oblate_spherical(\n\u001B[32m 88\u001B[39m equatorial_radius = equatorial_radius,\n\u001B[32m 89\u001B[39m flattening = flattening,\n\u001B[32m 90\u001B[39m )\n\u001B[32m 92\u001B[39m \u001B[38;5;66;03m# Add Earth's rotational settings via the 'environment_setup.rotation_model.gcrs_to_itrs'.\u001B[39;00m\n\u001B[32m 93\u001B[39m \u001B[38;5;66;03m# This function defines high-accuracy Earth rotation model according to the IERS Conventions 2010. \u001B[39;00m\n\u001B[32m 94\u001B[39m \u001B[38;5;66;03m# The model computes the rotation from ITRS to GCRS.\u001B[39;00m\n\u001B[32m 95\u001B[39m \u001B[38;5;66;03m# See [Tudatpy API Reference](https://py.api.tudat.space/en/latest/rotation_model.html#tudatpy.dynamics.environment_setup.rotation_model.gcrs_to_itrs)\u001B[39;00m\n\u001B[32m 96\u001B[39m body_settings.get(\u001B[33m'\u001B[39m\u001B[33mEarth\u001B[39m\u001B[33m'\u001B[39m).rotation_model_settings = environment_setup.rotation_model.gcrs_to_itrs(\n\u001B[32m 97\u001B[39m environment_setup.rotation_model.iau_2006, global_frame_orientation,\n\u001B[32m---> \u001B[39m\u001B[32m98\u001B[39m \u001B[43minterpolators\u001B[49m\u001B[43m.\u001B[49m\u001B[43minterpolator_generation_settings_float\u001B[49m(interpolators.cubic_spline_interpolation(),\n\u001B[32m 99\u001B[39m simulation_seconds_start_epoch, simulation_seconds_end_epoch, \u001B[32m60\u001B[39m),\n\u001B[32m 100\u001B[39m interpolators.interpolator_generation_settings_float(interpolators.cubic_spline_interpolation(),\n\u001B[32m 101\u001B[39m simulation_seconds_start_epoch, simulation_seconds_end_epoch, \u001B[32m60\u001B[39m),\n\u001B[32m 102\u001B[39m interpolators.interpolator_generation_settings_float(interpolators.cubic_spline_interpolation(),\n\u001B[32m 103\u001B[39m simulation_seconds_start_epoch, simulation_seconds_end_epoch, \u001B[32m60\u001B[39m))\n\u001B[32m 105\u001B[39m \u001B[38;5;66;03m# As you've already seen in the documentation of this function, users can input custom_ephemeris via the custom_ephemeris flag.\u001B[39;00m\n\u001B[32m 106\u001B[39m \u001B[38;5;66;03m# This is useful for objects that are not present on the Horizons database, \u001B[39;00m\n\u001B[32m 107\u001B[39m \u001B[38;5;66;03m# or for any other custom/fake object for which an ephemeris file is made available by the user.\u001B[39;00m\n\u001B[32m 108\u001B[39m \u001B[38;5;66;03m# If the custom_ephemeris flag is not on, then we create ephemeris by performing an Horizons Query; if the custom_ephemeris flag is on, \u001B[39;00m\n\u001B[32m 109\u001B[39m \u001B[38;5;66;03m# the ephemeris settings are set accordingly as the ones provided by the user.\u001B[39;00m\n\u001B[32m 111\u001B[39m \u001B[38;5;28;01mif\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m custom_ephemeris:\n\u001B[32m 112\u001B[39m \u001B[38;5;66;03m# Horizons query to retrieve ephemeris.\u001B[39;00m\n", "\u001B[31mAttributeError\u001B[39m: module 'tudatpy.kernel.math.interpolators' has no attribute 'interpolator_generation_settings_float'" ] } ], "source": [ "# Select station names\n", "station_names =['TIDBIN64', 'KATH12M', 'YARRA12M', 'HOBART26']\n", "\n", "# Here, we trick tudat by changing the (already defined in the tudat list of stations) position for KATH12M.\n", "# Hence, we will expect to have the same exact plot for both Yarragadee and Katherine.\n", "# We do this in order to show how an already-existing Tudat station location (such as the one of KATH12M) can be overwritten by the user. \n", "geodetic_positions = {'YARRA12M': [250,-29.0464, 115.3456], 'KATH12M': [250,-29.0464, 115.3456]}\n", "\n", "# Select the global frame and orientation\n", "global_frame_origin = 'Earth'\n", "global_frame_orientation = 'J2000'\n", "\n", "# Retrieving Horizons + Tudat Az/El plot for given start and end epochs\n", "start_epoch = '2022-12-19 00:00:00'\n", "end_epoch = '2022-12-20 00:00:00'\n", "\n", "plot_combined_elevation(\n", " '2014 HK129',\n", " station_names,\n", " start_epoch = start_epoch,\n", " end_epoch = end_epoch,\n", " time_step= '20m',\n", " geodetic_positions=geodetic_positions,\n", " custom_ephemeris = None\n", ")" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "-all", "main_language": "python", "notebook_metadata_filter": "-all" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.13" } }, "nbformat": 4, "nbformat_minor": 5 }