diff --git a/sst/mtm1m3/M1M3_max_forces_histogram.ipynb b/sst/mtm1m3/M1M3_max_forces_histogram.ipynb new file mode 100644 index 0000000..ac2c0cd --- /dev/null +++ b/sst/mtm1m3/M1M3_max_forces_histogram.ipynb @@ -0,0 +1,378 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "9627519e-ef49-41c5-9085-1a7373e7b3c0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "dayObs = 20250514 # format YYYYMMDD, e.g. \"20250514\"\n", + "bin_number = 100" + ] + }, + { + "cell_type": "markdown", + "id": "7c235e45-d932-432a-b473-f72774f97a37", + "metadata": {}, + "source": [ + "# SITCOM-2089 Evaluate M1M3 hard points individual forces (histogram)\n", + "\n", + "This notebook makes a histogram of maximum forces and spread of all slews passing certain conditions given a day_obs.\n", + "\n", + "### Prepare Notebook" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3d229f3-648a-4ee3-ac70-8e22bd90e53e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from astropy.time import Time\n", + "import matplotlib.dates as mdates\n", + "\n", + "from lsst.summit.utils.tmaUtils import TMAEventMaker, TMAState\n", + "from lsst.summit.utils.efdUtils import getEfdData, makeEfdClient\n", + "from lsst_efd_client import EfdClient\n", + "\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\")\n", + "\n", + "from lsst.ts.xml.enums.MTM1M3 import DetailedStates" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fa49a3c-6c77-4d44-939f-b713b1e5ca76", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# create a client to retrieve datasets in the EFD database\n", + "client = EfdClient(\"usdf_efd\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd6cf13c-81be-44a0-9177-6f0e1a5faed9", + "metadata": {}, + "outputs": [], + "source": [ + "OPERATIONAL_LIMIT = 450 # newtons,TBC with Petr\n", + "SAFE_LIMIT = 900\n", + "BREAKAWAY_LIMIT = 3200" + ] + }, + { + "cell_type": "markdown", + "id": "58085fea-42ff-43f0-89be-f3a2a1101b42", + "metadata": { + "tags": [] + }, + "source": [ + "### Define functions" + ] + }, + { + "cell_type": "markdown", + "id": "c741d971-c0c5-48eb-b99d-553650ccc406", + "metadata": { + "tags": [] + }, + "source": [ + "#### get previous logged detailedState event from a given time stamp\n", + "\n", + "The status of M1M3 is not persistent in the EFD. In order to get at a given time what is the current status, use this function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f13c2bfb-3490-4685-a18b-39633d90ee2a", + "metadata": {}, + "outputs": [], + "source": [ + "def get_previous_logged_detailedState(df_state, timestamp):\n", + " \"\"\"\n", + " Get logged detailedState from M1M3 immediately before arbitrary time\n", + " Args:\n", + " df_state (pandas dataframe): pandas dataframe obtained from time series of\n", + " \"lsst.sal.MTM1M3.logevent_detailedState\" covering a wide time frame which includes\n", + " the time stamp\n", + " timestamp (pandas timestamp): a timestamp where we want to probe the current status of M1M3\n", + " Returns:\n", + " prev_state: human readable status of current M1M3 status\n", + " \"\"\"\n", + " df_state_names = df_state[\"detailedState\"].map(lambda x: DetailedStates(x).name)\n", + " previous_index = df_state.index.asof(timestamp)\n", + " try:\n", + " prev = df_state.index.get_loc(previous_index)\n", + " except KeyError:\n", + " return \"KeyError\"\n", + " return df_state_names[prev]" + ] + }, + { + "cell_type": "markdown", + "id": "72fdf187-0dbb-477c-a8dc-6749a9f8c9f8", + "metadata": { + "tags": [] + }, + "source": [ + "### Define relevant settings" + ] + }, + { + "cell_type": "markdown", + "id": "5a56db2c-4172-46db-9155-3ff63fcbc5b9", + "metadata": {}, + "source": [ + "#### Observation day" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "beda404b-dea9-49ba-9544-be6cbf635a63", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "## Insert here the day_obs of interest\n", + "day_obs = dayObs" + ] + }, + { + "cell_type": "markdown", + "id": "decee784-6779-4188-84d7-86727cb7cdb8", + "metadata": { + "tags": [] + }, + "source": [ + "### Load data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b66bab89-b2ef-42ae-a73e-ab34e72ea514", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Select data from a given date\n", + "eventMaker = TMAEventMaker()\n", + "events = eventMaker.getEvents(day_obs)\n", + "\n", + "# Get lists of slew and track events\n", + "slews = [e for e in events if e.type == TMAState.SLEWING]\n", + "tracks = [e for e in events if e.type == TMAState.TRACKING]\n", + "print(f\"There are {len(events)} events\")\n", + "print(f\"Found {len(slews)} slews and {len(tracks)} tracks\")" + ] + }, + { + "cell_type": "markdown", + "id": "d3a7738e-6b42-4a06-8fc9-e8f5db9951c0", + "metadata": {}, + "source": [ + "### Select slews passing certain criteria" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05266ebb-b00b-44a5-a26f-d246ebc11843", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "slews_selected = []\n", + "hp_max_hist = np.array([])\n", + "hp_spread_max_hist = np.array([])\n", + "min_ele_range = -1 # minimum elevation change in slew, degrees\n", + "min_azi_range = -1 # minimum azimuth change in slew, degrees\n", + "hp_threshold = OPERATIONAL_LIMIT\n", + "velocity_threshold = 0.5 # deg/s\n", + "\n", + "df_state = getEfdData(\n", + " client,\n", + " \"lsst.sal.MTM1M3.logevent_detailedState\",\n", + " begin=Time(slews[0].begin, format=\"isot\", scale=\"utc\"),\n", + " end=Time(slews[-1].end, format=\"isot\", scale=\"utc\"),\n", + ") # get an array for all state changes from first slew of day_obs to final one\n", + "\n", + "for i, slew in enumerate(slews):\n", + " if (\n", + " slew.seqNum == 0\n", + " ): # skip first one to avoid problems looking for a previous detailedState outside the df_state range\n", + " continue\n", + " \n", + " df_azi = getEfdData(client, \"lsst.sal.MTMount.azimuth\", event=slew)\n", + " df_ele = getEfdData(client, \"lsst.sal.MTMount.elevation\", event=slew)\n", + " df_hp = getEfdData(client, \"lsst.sal.MTM1M3.hardpointActuatorData\", event=slew)\n", + " timestamp = pd.Timestamp(\n", + " Time(slew.begin, format=\"iso\", scale=\"utc\").value, tz=\"utc\"\n", + " )\n", + " begin_state = get_previous_logged_detailedState(df_state, timestamp)\n", + " \n", + " if begin_state == \"KeyError\":\n", + " continue\n", + " \n", + " timestamp = pd.Timestamp(Time(slew.end, format=\"iso\", scale=\"utc\").value, tz=\"utc\")\n", + " end_state = get_previous_logged_detailedState(df_state, timestamp)\n", + " \n", + " if len(df_hp) > 0: \n", + " hp_max_individual = np.array(\n", + " [\n", + " np.max(abs(df_hp[\"measuredForce0\"].values)),\n", + " np.max(abs(df_hp[\"measuredForce1\"].values)),\n", + " np.max(abs(df_hp[\"measuredForce2\"].values)),\n", + " np.max(abs(df_hp[\"measuredForce3\"].values)),\n", + " np.max(abs(df_hp[\"measuredForce4\"].values)),\n", + " np.max(abs(df_hp[\"measuredForce5\"].values)),\n", + " ]\n", + " )\n", + " hp_max = np.max(hp_max_individual)\n", + " hp_spread_individual = np.array(\n", + " [\n", + " np.max(df_hp[\"measuredForce0\"].values)-np.min(df_hp[\"measuredForce0\"].values),\n", + " np.max(df_hp[\"measuredForce1\"].values)-np.min(df_hp[\"measuredForce1\"].values),\n", + " np.max(df_hp[\"measuredForce2\"].values)-np.min(df_hp[\"measuredForce2\"].values),\n", + " np.max(df_hp[\"measuredForce3\"].values)-np.min(df_hp[\"measuredForce3\"].values),\n", + " np.max(df_hp[\"measuredForce4\"].values)-np.min(df_hp[\"measuredForce4\"].values),\n", + " np.max(df_hp[\"measuredForce5\"].values)-np.min(df_hp[\"measuredForce5\"].values),\n", + " ]\n", + " )\n", + " hp_spread_max = np.max(hp_spread_individual)\n", + " \n", + " # introduce conditions to select specific slews\n", + "\n", + " # condition on minimum azimuth and minimum elevation\n", + " slew_delta_azi = df_azi[\"demandPosition\"].max() - df_azi[\"demandPosition\"].min()\n", + " slew_delta_ele = df_ele[\"demandPosition\"].max() - df_ele[\"demandPosition\"].min()\n", + " slew_ele_condition = slew_delta_ele > min_ele_range\n", + " slew_azi_condition = slew_delta_azi > min_azi_range\n", + "\n", + " # ensure that the HPs are active (mirror is 'raised')\n", + " state_condition = ( \n", + " (begin_state == \"ACTIVE\") or (begin_state == \"ACTIVEENGINEERING\")\n", + " ) and ((end_state == \"ACTIVE\") or (end_state == \"ACTIVEENGINEERING\"))\n", + "\n", + " # condition on maximum azimuth velocity\n", + " velocity_condition = (df_azi[\"actualVelocity\"] < velocity_threshold).all()\n", + "\n", + " # check all conditions to select slews and fill histogram\n", + " if (\n", + " slew_ele_condition\n", + " and slew_azi_condition\n", + " and state_condition\n", + " and velocity_condition\n", + " ):\n", + " hp_max_hist = np.append(hp_max_hist,hp_max) \n", + " hp_spread_max_hist = np.append(hp_spread_max_hist,hp_spread_max)\n" + ] + }, + { + "cell_type": "markdown", + "id": "1fe2c77d-8883-4ee1-af48-f62b6b3a9f9d", + "metadata": {}, + "source": [ + "### Plot histograms" + ] + }, + { + "cell_type": "markdown", + "id": "194e2712-3e93-49c8-8f0f-136b4d4ff5ab", + "metadata": {}, + "source": [ + "#### Plot distribution of maximum forces" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53704b04-cd37-42eb-aa03-970f4e3510a5", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot a distribution of maximum forces\n", + "plt.ylabel(\"Number of slews\")\n", + "plt.xlabel(\"Maximum recorded force in any hardpoint per slew (N)\")\n", + "plt.title(f\"Maximum force on hardpoints for {day_obs}\")\n", + "plt.axvline(x=OPERATIONAL_LIMIT, color='orange', linestyle='--', linewidth=2, label='Operational limit')\n", + "plt.axvline(x=SAFE_LIMIT, color='red', linestyle='--', linewidth=2, label='Safe limit')\n", + "plt.axvline(x=BREAKAWAY_LIMIT, color='black', linestyle='--', linewidth=2, label='Breakaway limit')\n", + "plt.legend()\n", + "h = plt.hist(hp_max_hist, bins=bin_number)" + ] + }, + { + "cell_type": "markdown", + "id": "15170c22-f455-4873-a6f1-fdeaa06a37ab", + "metadata": {}, + "source": [ + "#### Plot distribution of difference of maximum and minimum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80e39543-a4ba-436f-b638-68047fb62de0", + "metadata": {}, + "outputs": [], + "source": [ + "plt.ylabel(\"Number of slews\")\n", + "plt.xlabel(\"Maximum spread of forces in any hardpoint per slew (N)\")\n", + "plt.title(f\"Maximum spread of forces on hardpoints for {day_obs}\")\n", + "plt.legend()\n", + "h = plt.hist(hp_spread_max_hist, bins=bin_number)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8db084a-b066-47c9-a834-097f8bd881b7", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "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.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/sst/mtm1m3/M1M3_max_forces_histogram.yaml b/sst/mtm1m3/M1M3_max_forces_histogram.yaml new file mode 100644 index 0000000..09c1d9c --- /dev/null +++ b/sst/mtm1m3/M1M3_max_forces_histogram.yaml @@ -0,0 +1,18 @@ +title: MTM1M3 Max Forces Histogram +description: A notebook to make max forces histogram to check M1M3 status every morning +authors: + - name: HyeYun Park + slack: hyeyun +tags: + - mtm1m3 + - force actuators + - following error +parameters: + dayObs: + type: integer + description: Observation day in format YYYYMMDD + default: 20250514 + bin_number: + type: integer + description: Bin number for the histogram,, e.g. 100 + default: 100