notebooks/ConstellationGeneration/ConstellationGeneration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "41b99f83",
   "metadata": {},
   "source": [
    "# Constellation Generation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ef7e733",
   "metadata": {},
   "source": [
    "## Quick guide:\n",
    "* Initialize the constellation by executing the first cell and providing metadata in the second cell (1)\n",
    "* Create a shell by providing the 4 shell arguments, and further parameters (extra params) if necessary (2.1).\n",
    "Store the temporary shell data by executing the cell (2.2)\n",
    "* Turn satellites into position and velocity vectors by executing cell (3)\n",
    "* Write the files by executing cell (4) and save them by executing cell (5)\n",
    "\n",
    "The following guide contains more detail\n",
    "\n",
    "## Detailed guide:\n",
    "\n",
    "### 0. Execute the first cell to import libraries and initialize variables\n",
    "\n",
    "### 1. Initialize Constellation and reset planet_list and other lists\n",
    "Provide the constellation name, the time of launch `startTime` as well as\n",
    "the time it takes to fully deploy the constellation `duration` and execute the cell. \n",
    "`startTime` can be specified as a date string with the format \"YYYY-MM-DD\". `duration` is specified as a string that contains a number of days followed by the 'd' character. If the variables are assigned plain integers, the values will be interpreted as simulation iterations.\n",
    "### 2.1 Generate shell = ( altitude, inclination, nPlanes, nSats )\n",
    "The execution of this cell creates one temporary shell (each satellite is described by orbital elements) that is not yet a part of the constellation. It is possible to visualize the generated satellites by executing the cell `Plot current shell` before adding it to the constellation. Append the shell to the constellation by executing cell `2.2`. Each time the code of `2.1` is run, the temporary data is overridden by the new data.\n",
    "\n",
    "Create a shell by specifying the parameters of an orbital shell which are the \n",
    "following:\n",
    "* `altitude` the altitude of each satellite in kilometers, relative to the\n",
    "earths surface\n",
    "* `inclination` the inclination of each orbital plane in degrees\n",
    "* `nPlanes` the number of orbital planes\n",
    "* `nSats` the number of satellites in each plane\n",
    "\n",
    "In some cases it is necessary to specify more parameters: `offsetM`, `argPeriapsis`,\n",
    "`startingW`, and `W_area`. It is best to modify them only if necessary and to use\n",
    "the default values otherwise. Refer to the end of the notebook to see their usage\n",
    "information (additional parameters).\n",
    "\n",
    "### 2.2 Append to elements_list and other\n",
    "Append the temporary set of satellites (elements_list_tmp) to the\n",
    "constellation (elements_list). The same happens to lists with the shells\n",
    "meta information. At this stage elements_list entries are the orbital elements of the satellite.\n",
    "\n",
    "### 3. Create position and velocity vectors\n",
    "Execute the code in order to create position and velocity vectors from the list of orbital elements. They are stored in the 'satellites' variable.\n",
    "\n",
    "### 4. Plot and store results\n",
    "Prepare .csv and .yaml output and simultaneously plot the constellation by\n",
    "executing this cell.\n",
    "\n",
    "### 5. Save Constellation as files in directory\n",
    "Create a directory with the same name that was specified in `1` in the variable\n",
    "`constellation_name` in the `data` directory of the project. It contains three\n",
    "files: Two .csv files with positions and velocities and a .yaml file with the\n",
    "metadata of the constellation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b5c8ba3",
   "metadata": {},
   "source": [
    "## 0. Import and initialize variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a48e2c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Imports\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "# Append main folder\n",
    "import sys\n",
    "sys.path.append(\"../\")\n",
    "\n",
    "import pykep as pk\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import os\n",
    "import yaml\n",
    "\n",
    "starting_t = pk.epoch_from_string('2022-01-01 00:00:00.000')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8754e98",
   "metadata": {},
   "source": [
    "# 1. Initialize Constellation and reset planet_list and other lists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "906f7bda",
   "metadata": {},
   "outputs": [],
   "source": [
    "#---Constellation data----------------------------------------------------------------------------------------------\n",
    "constellation_name = \"Constellation\"\n",
    "startTime = 0\n",
    "duration = 0\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "nShells = 0\n",
    "elements_list = []\n",
    "shell_list = []\n",
    "xparams_list = []\n",
    "print(\"Planet list is now empty\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e78dc60e",
   "metadata": {},
   "source": [
    "# 2.1 Generate shell = ( altitude, inclination, nPlanes, nSats )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fd559cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# creates a temporary shell that has to be appended to planet_list in cell \"Append to planet_list and other\"\n",
    "# to be stored and not get overwritten (this also concerns additional shell data shell_list, xparams_list)\n",
    "\n",
    "#---input-----------------------------------------------------------------------------------------------------------\n",
    "altitude = 700\n",
    "inclination = 40.0\n",
    "nPlanes = 20\n",
    "nSats = 10\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "#---extra-params----------------------------------------------------------------------------------------------------\n",
    "\n",
    "#offsetM = offset for Mean Anomaly added after each plane (relative phasing)\n",
    "# walker constellation: offsetM = F * 360 / (nPlanes * nSats) ; F element {0, ... , nPlanes - 1}\n",
    "offsetM = 0    # default 0\n",
    "\n",
    "#argPeriapsis = argument of periapsis\n",
    "# starting point of satellite placement for each plane\n",
    "# argPeriapsis = pi avoids collisions in planes that intersect at reference plane\n",
    "argPeriapsis = math.pi  # default math.pi \n",
    "\n",
    "#startingW = offset for W that is not accumulating (W = longitude of ascending node)\n",
    "# formula for overlapping shells (same altitude, same inclination):\n",
    "# (360 / G) / 2 ; G = smallest common multiple of the overlapping nPlanes\n",
    "startingW = 0  # default 0\n",
    "\n",
    "# W_area: orbital planes are distributed evenly within range [startingW,startingW + maximumW)\n",
    "W_area = 360     # default 360\n",
    "#-------------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "minW = startingW\n",
    "maxW = W_area + startingW\n",
    "\n",
    "a = altitude * 1000 + 6371000    # in [m], earth radius included\n",
    "e = 0\n",
    "i = inclination * pk.DEG2RAD\n",
    "W = pk.DEG2RAD * minW\n",
    "w = argPeriapsis * pk.DEG2RAD       \n",
    "M = 0             \n",
    "\n",
    "plane_count = 0\n",
    "\n",
    "when = starting_t\n",
    "mu_central_body = pk.MU_EARTH\n",
    "mu_self = 1\n",
    "radius = 1\n",
    "safe_radius = 1\n",
    "\n",
    "pStep = pk.DEG2RAD * W_area / nPlanes  # W goes from startingW to startingW+W_area\n",
    "sStep = 2 * math.pi / nSats              # M goes from 0° to 360°\n",
    "sExtraStep = pk.DEG2RAD*offsetM\n",
    "\n",
    "planet_list_tmp = []\n",
    "elements_list_tmp = []\n",
    "for x in range(nPlanes):\n",
    "    for y in range(nSats):\n",
    "        planet_list_tmp.append(pk.planet.keplerian(when,[a,e,i,W,w,M], mu_central_body, mu_self,radius, safe_radius ,\"sat\"))\n",
    "        elements_list_tmp.append([a,e,i,W,w,M])\n",
    "        M = M + sStep\n",
    "    plane_count = plane_count + 1\n",
    "    W = W + pStep\n",
    "    M = plane_count * sExtraStep     #equals 0 + count*0 = 0 in the usual case\n",
    "    \n",
    "shell = (altitude,inclination,nPlanes,nSats)\n",
    "xparams = (offsetM,argPeriapsis,startingW,W_area)\n",
    "\n",
    "print(\"Added \" + str(len(elements_list_tmp)) + \" planets\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a4e87f4",
   "metadata": {},
   "source": [
    "# 2.2 Append to planet_list and other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "202135a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# stores and appends temporary shell in planet_list (including lists with additional data)\n",
    "nShells += 1\n",
    "elements_list = elements_list + elements_list_tmp\n",
    "shell_list.append(shell)\n",
    "xparams_list.append(xparams)\n",
    "print(\"Added \" + str(len(elements_list_tmp)) + \" to planet_list now totalling \" + str(len(elements_list)) + \" planets.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15e632e6",
   "metadata": {},
   "source": [
    "### Plot current shell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82683ed1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,6),dpi=100)\n",
    "ax = plt.axes(projection='3d')\n",
    "for i in range (nPlanes*nSats):\n",
    "    pk.orbit_plots.plot_planet(planet_list_tmp[i],axes=ax,s=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbd01df1",
   "metadata": {},
   "source": [
    "# 3. Create position and velocity vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3106348d",
   "metadata": {},
   "outputs": [],
   "source": [
    "satellites = []\n",
    "for elements in elements_list:\n",
    "    pos,v = pk.par2ic(elements,pk.MU_EARTH)\n",
    "        \n",
    "    # convert to km and numpy\n",
    "    pos = np.asarray(pos) / 1000.0 \n",
    "    v = np.asarray(v) / 1000.0\n",
    "        \n",
    "    satellites.append((pos,v))\n",
    "        \n",
    "print(\"Successfully propagated \",len(satellites),\" satellites.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c77e85b",
   "metadata": {},
   "source": [
    "# 4. Plot and store results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8aae75ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6,6),dpi=100)\n",
    "ax = plt.axes(projection='3d');\n",
    "# pos , v (.csv)\n",
    "positions = np.array([pos for pos,_ in satellites])\n",
    "velocities = np.array([v for _,v in satellites])\n",
    "ax.scatter(positions[:,0],positions[:,1],positions[:,2],\".\",alpha=0.25)\n",
    "# parameters (.yaml)\n",
    "constellation_cfg = {\n",
    "    \"constellation\": {\n",
    "        \"name\": constellation_name,\n",
    "        \"startTime\": startTime,\n",
    "        \"duration\": duration,\n",
    "        \"nShells\": nShells\n",
    "    }\n",
    "}\n",
    "for i in range(nShells):\n",
    "    constellation_cfg[\"shell\" + str(i+1)] = {\n",
    "        \"altitude\": shell_list[i][0],\n",
    "        \"inclination\": shell_list[i][1],\n",
    "        \"nPlanes\": shell_list[i][2],\n",
    "        \"nSats\": shell_list[i][3]\n",
    "    }\n",
    "for i in range(nShells):\n",
    "    constellation_cfg[\"special\" + str(i+1)] = {\n",
    "        \"offsetM\": xparams_list[i][0],\n",
    "        \"argPeriapsis\": xparams_list[i][1],\n",
    "        \"startingW\": xparams_list[i][2],\n",
    "        \"W_area\": xparams_list[i][3]\n",
    "    }"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91f07380",
   "metadata": {},
   "source": [
    "# 5. Save Constellation as files in directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf6eb7d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# creates directory named as the constellation_name value within the data folder and stores output there\n",
    "\n",
    "os.mkdir(\"../../data/\" + constellation_name);\n",
    "\n",
    "np.savetxt(\"../../data/\" + constellation_name + \"/pos_\"+constellation_name+\".csv\",positions,delimiter=\",\")\n",
    "np.savetxt(\"../../data/\" + constellation_name + \"/v_\"+constellation_name+\".csv\",velocities,delimiter=\",\")\n",
    "\n",
    "with open(\"../../data/\" + constellation_name + \"/shells_\"+constellation_name+\".yaml\", \"w\") as fh:  \n",
    "    yaml.dump(constellation_cfg, fh,sort_keys=False)\n",
    "    fh.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19cfec21",
   "metadata": {},
   "source": [
    "### Additional parameters\n",
    "### offsetM\n",
    "`offsetM` is the step size for an offset that is added to each satellite's mean anomaly. The offset accumulates, meaning it increases by `offsetM` each time a new plane is generated. In a Walker Constellation, \n",
    "which is a typical model for constellations, the parameter is assigned a value of the following form:\n",
    "\n",
    "`offsetM = F * 360 / (nPlanes*nSats)` \n",
    "\n",
    "Note: F is a natural number: F is element of {0, ... , nPlanes - 1}\n",
    "\n",
    "Note: the property of a set phasing difference\n",
    "between neighboring planes gets lost if the shell is added over time (`duration` > 0)\n",
    "\n",
    "\n",
    "\n",
    "### argPeriapsis\n",
    "`argPeriapsis` corresponds to the orbital element \"argument of periapsis\". Choosing an irrational number\n",
    "prevents satellites from two symmetrical planes to have the same position.\n",
    "\n",
    "### startingW\n",
    "`startingW` is the base value for the shell's longitude of ascending node (W). The parameter can be utilized if planes of different shells overlap. This occurs if shells have the same altitude and the same or similar inclination.\n",
    "\n",
    "The shells should follow the formula:\n",
    "\n",
    "`(360/G)*(i/N)`\n",
    "\n",
    "Note: G stands for the smallest common denominator of the nPlanes parameters of\n",
    "each shell\n",
    "\n",
    "Note: N is the number of shells overlapping\n",
    "\n",
    "Note: i is an index in {0,...,N-1}. each shell should have a distinct i\n",
    "\n",
    "### W_area\n",
    "this parameter determines which degree the planes of the shell span. The longitudes of ascending node are distributed in the range \\[startingW , startingW + W_area\\]."
   ]
  }
 ],
 "metadata": {
  "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.10.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}