CoLPaper.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Role of the Center of Light for the Uncertainty of photogoniometric Measurements"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## Introduction to the Center of Light\n",
    "The Center of Light is a concept that results from the way Luminous Intensity is defined. \n",
    "The quantity of luminous intensity is used to express the emitted light with a directional dependency- integrating over any spacial distribution of the light source. This is a desired simplification of reality. The Luminous Intensity itself is a model with the assumption, that the luminous object is a point source. The Center of Light is the link between a real luminous object and this model: It is the point that the object is reduced to. The Center of Light is therefore a property of the luminous object, that has to be taken into account when measuring luminous intensity.\n",
    "\n",
    "The standard setup for measuring luminous intensity is a farfield goniometer. It uses a detector revolving around the object or vice versa at a constant distance. A luminous intensity distribution (LID) is constructed from multiple measurements performed for different directions relative to the center of the machine. For every measured direction it samples the illuminance caused by the source at this distance and derives from it the luminous intensity. This measurement setup assumes a point source in the center of the machine.\n",
    "The Model is implicitly used, but the actual object position is up to the user. The Center of Light of the luminous object is generally unknown, this opens up the possibility to introduce an error from the placement of the object. \n",
    "\n",
    "There is an inherent error resulting from the size of the object. This has been explored under the topic of the limiting photometric distance. This Error has the same underlying root as the error from translated Center of Light: Light coming from places other than the assumed point location. The result of is a deviating measurement distance and angle for every measurement direction. The unknown luminance characteristic of the object leads to an uncertainty contribution that cannot be prevented. Instead, the measurement distance is increased to reduce uncertainty to an acceptable margin. There are good tools and rules of thumb in place to estimate this uncertainty. \n",
    "In contrast, the Center of Light is not taken into account at all. Therefore, the aim is to sensitize operators to the problems arising from Center of Light. Most importantly, this factor can be integrated into an uncertainty analysis by expressing the confidence in the position of the source."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Impact for a LID measurement\n",
    "How the object is placed and oriented can be described by a pose in the coordinate system of the goniometer, consisting of two components: rotation and translation. The rotation of the object is normally dictated by the measurement procedure. A desired orientation could be to align the optical axis of the object to a certain direction in goniometer coordinates.\n",
    "The effect of deviating rotation has already been investigated. Changing the rotation of the object in the goniometer in turn results in a rotated LID. As long as the region of interest is still inside the sampled solid angle, there is no loss of information. And a LID with the desired rotation can be calculated after the measurement. Another work examined the possibility to compare LID measurements of the same or similar objects with different rotations. This can be accomplished by aligning the two LIDs via correlation.\n",
    "\n",
    "The effects of a translated source are more involved. A translation brings about a deviation of the measured luminous intensity that is different for every direction. This effect cannot be reversed from the measurement data alone.\n",
    "A LID can be thought of as a three-dimensional surface with the distance between origin and the surface translating to the luminous intensity in that direction. A measurement of a translated source results in a deformation of that surface. \n",
    "In the following example the LID of a lambert-source and its deformed counterpart are displayed. The mesh structure displays the expected sphere shape. The colored surface corresponds to the LID a goniometer measurement would produce for a translated source. The deformation is amplified considerably to make it visible. The surface color shows the actual deformation in %.\n",
    "In each case, the source was translated by the same value, along a single axis. The translation shown here amounts to 0.1% of the measurement distance. For a goniometer with 10m distance to the detector, this would mean a 10mm offset. \n",
    "Except for edge regions, where the LID approaches zero, the resulting error is small. There are multiple factors that can lead to an increase of this error, which are explored in the code sections of this notebook. The following images show simulations of this effect. \n",
    "<table>\n",
    "    <tr>\n",
    "    </tr>\n",
    "    <tr>\n",
    "      <td>\n",
    "      <img src='./figures/graphs_rel_i_r1_e0.001/x_error.PNG'width=400>\n",
    "      </td>\n",
    "      <td>\n",
    "      <img src='./figures/graphs_rel_i_r1_e0.001/y_error.PNG'width=400>\n",
    "      </td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "    </tr>\n",
    "    <tr>\n",
    "      <td>\n",
    "      <img src='./figures/graphs_rel_i_r1_e0.001/z_neg_error.PNG'width=400>\n",
    "      </td>\n",
    "      <td>\n",
    "      <img src='./figures/graphs_rel_i_r1_e0.001/z_pos_error.PNG'width=400>\n",
    "      </td>\n",
    "    </tr>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determining the Center of Light\n",
    "A definition for the Center of Light of an object can be constructed from its spatially distributed emission, more precise the model of light rays:\n",
    "The light emitted by the object consists of rays, that are modelled as lines. Using the point-source model forces the origin of all those rays to a single point, resulting in the error.\n",
    "The Center of Light is the best fitting origin for all rays of the actual object.\n",
    "It is important to note, that only the rays that exit the object are relevant for this approach. If rays originate from a light source within the object and then pass through a system of optics, only the diverted, exiting rays matter for the definition of the Center of Light. This is also the reason that the Center of Light does not need to coincide with the actual light source.\n",
    "\n",
    "It is possible to generate a ray model of an object from a near-field goniometer measurement. The Center of Light can then be solved for. This problem can be described as error-minimization. Center of Light calculation is already an analysis feature, defining it as the point with the closest quadratic distance to all rays. This corresponds to the point, where the rays approximately coincide. The measurement can be modelled as unit rays or rays weighted with an intensity. The in the latter case, the weight has to be considered when estimating the center.\n",
    "Even tough the Center of Light can be measured, this has little practical use. It is infeasible to perform every photogoniometric measurement task with a near field goniometer. For most applications the Center of Light is unknown and has to be estimated.\n",
    "There are special cases, where this estimation is simple. For a source consisting of diffuse emitting surfaces, the Center of Light corresponds to the weighted geometric center of the light emitting surfaces. If the light source is constructed with a rotational symmetry, the Center of Light will lie on the rotation axis. \n",
    "Optics complicate the matter. Any element that either deflects, reflects  or absorbs rays moves the Center of Light away from the actual light source. The Center of Light is not constrained to the dimensions of the object and can lie outside of it. A good intuition for the position of the Center of Light can in some cases be gained by looking into the source.\n",
    "In a setup, where a single light source is visible through the optical system, the source will appear roughly at the location of the Center of Light. This follows from the reversibility of optical paths. By tilting the object and tracking the depth of the apparent source is a good way to locate the plane, in which the Center of lIght lies.\n",
    "In any case the most important thing is to determine an uncertainty for the estimated point. This equivalent with to a confidence statement. This should be made conservatively, lowering the confidence for more complex sources."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Model for the Influence of the Center of Light\n",
    "The model is constructed in the context of a goniometer measurement. It defines a measurement setup based on a detector-based far-field goniometer. The results can also be applied to camera-based goniometers. The aim is to create a model that we can use to evaluate the effects of the center of light or equivalently the position of an object in the goniometer. Every other aspect of the measurement is treated as being ideal. To this purpose, the source is modelled as a point source. As outlined in the introduction, the influence of the object size and the Center of Light are closely linked. Nevertheless the Center of Light is investigated separately. Modeling both the source size and translation at the same time results in a very complex model. The effect of a translated source can also be separated reasonably well. Translating the source in one direction will impact the LID in a similar manner regardless of its size or shape.\n",
    "\n",
    "The LID is generally calculated from the illuminance distribution in a certain distance from the photometric Center: A point light source with the distribution$I_s(\\theta,\\phi)$ placed in the center has an illuminance distribution according to the inverse-square-law.\n",
    "\n",
    "$$E(\\theta,\\phi) = {I_s(\\theta,\\phi) \\over {r_c}^2} \\cdot cos(\\alpha_d)$$\n",
    "\n",
    "With $r_c$ being the measuring distance to the center and $\\alpha_d$ being the incidence angle of the light striking the detector. Since the distance between source and detector is known, and the detector is oriented perpendicular to the source, the luminous intensity distribution can be calculated as follows:\n",
    "\n",
    "$$I_m(\\theta,\\phi) = {E(\\theta,\\phi) \\cdot r_c^2} = I_s$$\n",
    "\n",
    "The goniometer is described by a polar coordinate system. Since the source does no longer coincide with the goniometer origin, it gets its own polar coordinate system.\n",
    "All angle variables are marked according to the coordinate system or origin to which they are measured. $\\theta_c$ and $\\phi_c$ describe a direction originating from the photometric center and therefore correspond with the angles of the C-plane coordinate system. Likewise $\\theta_s$ and $\\phi_s$ describe angles relating to the light source.\n",
    "The measurement $I_m$ is currently the same as the source LID $I_s$. If we now introduce a translation to this source, the measured LID deviates from the LID of the source.\n",
    "The aim of the model is to calculate the measured luminous intensity $I_m(\\theta_c, \\phi_c)$  for a given luminous intensity of the source $I_s(\\theta_s, \\phi_s)$ and a certain position of the source: $\\mathbf{s} = (s_x,s_y,s_z)$.\n",
    "The whole model is first stated and subdivided into parts, that are explained separately.\n",
    "\n",
    "$$I_m(\\theta_c,\\phi_c,\\mathbf{s}) = \\underbrace{I_s(\\theta_s,\\phi_s)}_\\text{direction} \\cdot \\underbrace{r_c^2 \\over r_s^2}_\\text{distance} \\cdot \\underbrace{cos(\\alpha_d)}_\\text{incidence}$$\n",
    "\n",
    "The influence of a translation can be separated in three different components. With the goniometer moved to a certain direction given by $\\theta_c,\\phi_c$ the LID of the source is actually sampled at diverging angles: $\\theta_s, \\phi_s$. This component is designated as the **Direction Error**. Additionally, the distance between source and detector is now different for every direction, introducing a **Distance Error**. Lastly, the incidence angle of the light striking the detector changes, resulting in the **Angle of Incidence Error**.\n",
    "\n",
    "1. **The Distance Error** </br>\n",
    "Because luminous intensity is not measured directly but instead inferred from illuminance, it requires the knowledge of the distance between source and detector. With the source offset from the photometric center, this distance $r_s$ is different for every detector position. When converting to luminous intensity, the constant distance to the photometric center $r_c$ is assumed instead, introducing an error. To describe this error, both source and detector are expressed as points in the cartesian coordinate system of the goniometer. The detector is modelled as revolving around the center and the source standing still.\n",
    "Position of the source: $\\mathbf{s} = (s_x,s_y,s_z)$\n",
    "Position of the detector for a particular measurement direction: $\\mathbf{d(\\theta_c,\\phi_c)} = (d_x,d_y,d_z)$\n",
    "The measurement directions and the measurement distance define the detector position in polar coordinates. They are converted to cartesian coordinates as follows:\n",
    "$$\n",
    "d_x(\\theta_c,\\phi_c) = r_c \\cdot sin(\\theta_c) \\cdot cos(\\phi_c) \\\\\n",
    "d_y(\\theta_c,\\phi_c) = r_c \\cdot sin(\\theta_c) \\cdot sin(\\phi_c) \\\\\n",
    "d_z(\\theta_c,\\phi_c) = r_c \\cdot cos(\\theta_c)\n",
    "$$\n",
    "The actual distance between the source and particular detector position can then be calculated as $r_s(\\theta_c,\\phi_c) = \\vert\\vert \\mathbf{\\overrightarrow{sd}} \\vert\\vert$.\n",
    "Taking into account the actual distance for the illuminance calculation results in the factor $r_c^2 \\over r_s^2$.\n",
    "\n",
    "$$I_{m1} = {I_s \\cdot {r_c^2 \\over {r_s}^2} }$$\n",
    "\n",
    "2. **The Direction Error** </br>\n",
    "The angles $\\theta_s, \\phi_s$ describe the direction at which the detector is pointed at the source. To calculate these angles, the detector positions are required in the coordinate system of the source position. The coordinate transform is done by substracting the cartesian coordinates of source and detectors: $\\mathbf{sd} = \\mathbf{d} - \\mathbf{s} = (sd_x, sd_y, sd_z)$\n",
    "These cartesian coordinates are then converted back into the polar-coordinates as follows:\n",
    "$$\\theta_s(\\theta_c,\\phi_c,\\mathbf{s}) = cos^{-1}({sd_z \\over r_s})$$\n",
    "$$\\phi_s(\\theta_c,\\phi_c,\\mathbf{s}) = tan^{-1}({sd_y \\over sd_x})$$\n",
    "These angles can now be inserted into the model function.\n",
    "$$I_{m2} = {I_s(\\theta_s, \\phi_s) \\cdot {r_c^2 \\over {r_s}^2} }$$\n",
    "Unlike the other two error components, the directional error is not only defined by the geometry of the measurement setup but also depends on the LID of the source that is measured. In the simplest case, when the LID is isotropic, this error does not occur, because the luminous intensity $I_s(\\theta_s,\\phi_s)$ is the same as $I_s(\\theta_c,\\phi_c)$. The magnitude of this error is defined by the difference in luminous intensity between the assumed and the actual measurement direction. It therefore scales with the gradient of the object's LID and is larger for directions with a large gradient in luminous Intensity.\n",
    "\n",
    "3. **The Angle of Incidence Error** </br>\n",
    "Lastly, there is an error because a translation of the source changes the incidence angle with which the light hits the detector. This angle is denoted as $\\alpha_d$ and can be obtained as the angle between two vectors. One vector being $\\mathbf{\\overrightarrow{cd}}$ between photometric center to the detector and the other vector being $\\mathbf{\\overrightarrow{sd}}$. between the source and the detector. It is calculated as: \n",
    "$$cos(\\alpha)={\\langle \\mathbf{\\overrightarrow{cd}},\\mathbf{\\overrightarrow{sd}} \\rangle \\over \\vert\\vert \\mathbf{\\overrightarrow{cd}} \\vert\\vert \\cdot \\vert\\vert \\mathbf{\\overrightarrow{sd}} \\vert\\vert}$$\n",
    "This factor reduces the measured illumninace for an increasing incidence angle. This assumes a cosine characteristic of the detector, which is very reasonable for the typically small incidence angles. With the addition of this factor the final model function becomes:\n",
    "$$I_m(\\theta_c,\\phi_c,\\mathbf{s}) = I_s(\\theta_s,\\phi_s) \\cdot {r_c^2 \\over r_s^2} \\cdot cos(\\alpha_d)$$\n",
    "\n",
    "Since all parts of the error are invariant to scaling, this model can also be used with relative units. For the following investigations all distances are expressed relative to the measurement distance $r_c$. This allows for generalized insights, that are applicable for all goniometer setups.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This model function is now expressed as python code. The mentioned coordinate transforms from cartesian to polar coordinates and back are done by functions defined in the Helperfunctions module.\n",
    "The arguments are:\n",
    "```\n",
    "x_s (float):      x-Coordinate of the source expressed relative to the measurement distance\n",
    "y_s (float):      y-Coordinate of the source expressed relative to the measurement distance\n",
    "z_s (float):      z-Coordinate of the source expressed relative to the measurement distance\n",
    "LID (function):   light intensity distribution: LID(theta, phi) -> I (angles in radians)\n",
    "theta_c(array):   theta angles at which the LID should be probed (radians)\n",
    "phi_c(array):     phi angles at which the LID should be probed (radians)\n",
    "```\n",
    "It returns a two-dimensional array of luminous intensity for every provided ```theta_c,phi_c``` pair"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def goniometer_model(x_s, y_s, z_s, LID, theta_c, phi_c):\n",
    "    import numpy as np\n",
    "    from Helperfunctions import polar2cart, cart2polar\n",
    "    # s: Position of the source as a vector\n",
    "    source_pos = np.array([x_s,y_s,z_s]).reshape(3, 1, 1)\n",
    "    \n",
    "    # d_c: Positions of the detector in cartesian coords relative to the center\n",
    "    detector_pos_c = polar2cart(theta_c,phi_c,1)\n",
    "    # d_s: Positions of the detector in cartesian coords relative to the source\n",
    "    detector_pos_s = detector_pos_c - source_pos\n",
    "\n",
    "    # convert Detector positions back to polar coordinates to get the angles theta_s, phi_s and the actual measurement distance r_s\n",
    "    theta_s, phi_s, r_s = cart2polar(*detector_pos_s)\n",
    "\n",
    "    # sample LID of the source at theta_s and phi_s\n",
    "    I_s = LID(theta_s,phi_s)\n",
    "\n",
    "    # calculate the distance factor\n",
    "    distance_factor = 1 / np.square(r_s)\n",
    "\n",
    "    # calculate the incidence angle factor\n",
    "    dot = np.sum(np.multiply(detector_pos_s,detector_pos_c),axis=0) # dot product of cd, sd\n",
    "    cos_a = dot / r_s\n",
    "    \n",
    "    I_m = I_s * distance_factor * cos_a\n",
    "    # the uncertainpy framework expects an array of timestamps as the first return value.\n",
    "    # Since this Model is evaluated over angles and not time, we return None.\n",
    "    return None, I_m\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This model has the LID of the source as an input parameter. The simulations will be performed using cosine distributions. For a power of one they describe a lambertian source. Higher powers of the cosine function are a good approximation for collimated light sources.\n",
    "\n",
    "We need a function that describes the LID of our source. It takes two arbitrary angles for theta and phi and returns the luminous intensity in that direction. The next cell creates such a function for a selectable power. The intensity is normalized to an arbitrary value of 100 cd."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cosn_LID_generator(n):\n",
    "    def LID(theta,phi):\n",
    "        import numpy as np\n",
    "        intensity = 100 * np.power(np.cos(theta),n) + np.zeros_like(phi)\n",
    "        return intensity.clip(min=0)\n",
    "    return LID"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will use the goniometer model to simulate a measurement of our source for a specific offset location. We first define the angles $\\theta_c$ and $\\phi_c$ that should be measured.\n",
    "The parameter ```theta_max``` is used to limit the measurement region. It is sensible to adjust this value depending on the LID of the source. If the source is collimated to a narrow cone, it is unnecessary to simulate the angular outskirts of the measurement, where the luminous intensity is effectively 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "theta_max=89\n",
    "theta_deg = np.linspace(start=0,stop=theta_max,num=theta_max+1)\n",
    "phi_deg = np.linspace(start=0,stop=360,num=37)[:-1]\n",
    "\n",
    "# convert to radians and add empty dimensions so that operations on phi and theta return a 2d array \n",
    "theta = np.deg2rad(theta_deg)[:,None]\n",
    "phi = np.deg2rad(phi_deg)[None,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a source function for the power ```n``` and evaluate it for the previously defined angles to get a ground truth, that the model output can be compared to later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the source function\n",
    "n = 10\n",
    "LID_source_function = cosn_LID_generator(n)\n",
    "# calculate the discrete LID\n",
    "LID_source = LID_source_function(theta, phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All that is left in order to run our model are inputs for the source position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run model\n",
    "x_offset = 0.001 # relative measure\n",
    "y_offset = -0.002\n",
    "z_offset = 0.003\n",
    "LID_measurement = goniometer_model(x_offset, y_offset, z_offset, LID_source_function, theta, phi)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For realistic offsets the resulting error is rather small. it would not be visible when plotting both the true and measured LID of the source. Therefore we construct a LID that amplifies the error by a certain factor. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional step: amplify the difference between the source and measurement LID to make it more visible\n",
    "amp_factor = 10\n",
    "LID_amplified = LID_source + amp_factor * (LID_measurement - LID_source)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The way the arrays theta and phi are defined results in a representation of the LID in \"half\" C-Planes. For every angle $\\phi$ there is a slice from  $\\theta:$ $0° .. 90°$. For Visualization purposes it is preferrable to have a complete C-Plane: For every angle $\\phi$ the slice ranges from $\\theta:$ $-90° .. 90°$. The C-Planes are constructed by combining opposing half-Planes. For this to work the $\\phi$ angles need to line up. \n",
    "The function ```half2cplanes``` reorders a 2D array, a theta and a phi array to its halfplane equivalent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Helperfunctions import half2cplanes\n",
    "\n",
    "# convert both LIDs to c-plane format\n",
    "LID_source_cplane, theta_cplane, phi_cplane = half2cplanes(LID_source, theta, phi)\n",
    "LID_measurement_cplane = half2cplanes(LID_measurement,theta,phi)[0]\n",
    "LID_amplified_cplane = half2cplanes(LID_amplified, theta, phi)[0]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization\n",
    "\n",
    "The visualization of a LID can be done in different ways. It can be represented as a 3D surface or as Slices of this surface. (C-Planes)\n",
    "A C-Plane plots the luminous intensity over angles. This can be done in a polar or carthesian coordinate system. A polar coordinate system is useful, because it visualizes the shape of the LID. Cartesian axes allow for easier reading of numerical values.    \n",
    "If all the content that is plotted is rotationally symmetrical a polar and Cartesian representation can be combined. The LIDs used here are rotationally symmetrical. For this representation to make sense, the error has to share this symmetry. It is therefore only applicable for situations without a x/y-offset.\n",
    "\n",
    "The repository contains a module called ```HalfplaneVisualization```. It defines a small visualization class, which implements a simplified interface for plotting using pyplot.\n",
    "the next cell contains the creation of a visualization object and a method to show the empty axes. The arguments of the ```show_fig``` method control which axes should be displayed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from HalfplaneVisualization import HalfplaneVisu\n",
    "\n",
    "visu = HalfplaneVisu(theta_max=theta_max,theta_tick_dist=15,r_max = 125,r_min=-100)\n",
    "visu.show_fig(cart=True, polar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot the results from running the model using the plot method. It takes an array of theta angles and an array of luminous intensities. which c-plane is plotted is controlled by the index ```c```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose the c-plane to plot\n",
    "c = 0 \n",
    "\n",
    "visu.plot(theta_cplane, LID_source_cplane[c], label=r'$C_0$-Plane of the source LID', c='lightgreen')\n",
    "visu.plot(theta_cplane, LID_amplified_cplane[c], label=r'$C_0$-Plane of measured LID', c='red')\n",
    "\n",
    "visu.show_fig(cart=True, polar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error analysis\n",
    "The plot now contains these two curves. Since it's hard to reason about the error from this plot, in the next cell, we will plot the actual error. For this purpose one can use the auxiliary axis on the right. To plot something to this axis use the  ```plot_aux ``` method. The axis will automatically be scaled to the plotted content and its zero position aligned to the luminous intensity axis.\n",
    "The method ```aux_set_center``` can be used to change the value on the axis that is aligned.\n",
    "We will plot the error in two possible ways: First as a percentage of the current luminous intensity, which describes the relative error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the relative error between luminous intensity of source and measurement in percent\n",
    "error_rel_prcnt = 100 * (LID_measurement - LID_source) / LID_source\n",
    "error_rel_prcnt_cplane = half2cplanes(error_rel_prcnt,theta,phi)[0]\n",
    "\n",
    "\n",
    "visu.plot_aux(theta_cplane, error_rel_prcnt_cplane[c], label=r'Relative Error in %', c='blue')\n",
    "visu.aux_set_label('%')\n",
    "visu.show_fig(cart=True, polar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next up we instead plot the error in absolute terms. For this we give the error as a percentage of the maximum luminous intensity of the source.\n",
    "The method ```clear_aux``` deletes any previous plots to the auxiliary axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the absolute error as % of the luminous intensity along the optical axis\n",
    "error_abs_prcnt = 100 * (LID_measurement - LID_source) / LID_source[0,0]\n",
    "error_abs_prcnt_cplane = half2cplanes(error_abs_prcnt,theta,phi)[0]\n",
    "\n",
    "visu.clear_aux()\n",
    "visu.plot_aux(theta_cplane, error_abs_prcnt_cplane[c], label=r'Absolute Error as % of $I_0$', c='orange')\n",
    "visu.show_fig(cart=True, polar=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this we have all the parts needs to examine the error resulting from a translated source. To make it easier to explore different szenarios, the visualization is extended to make it interactive. For this purpose we define an update-function, that recalcualtes the model and adjusts the visualization.\n",
    "Input parameters are the source position and the LID function's cosine power that are necessary to rerun the model. Additionally there are multiple inputs to adapt the visualization: \n",
    "* The C-Plane to be shown\n",
    "* The error amplification\n",
    "* The axes to be shown\n",
    "* The way the error should be plotted\n",
    "This function combines the contents of the previous cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "\n",
    "def update(x,y,z,c,n,amp,err,axes):\n",
    "    # % to factor -> /100\n",
    "    x_offset, y_offset, z_offset = np.array((x,y,z)) / 100\n",
    "    # update the source LID\n",
    "    LID_source_function = cosn_LID_generator(n)\n",
    "    LID_source = LID_source_function(theta, phi)\n",
    "    LID_source_cplane, theta_cplane, phi_cplane = half2cplanes(LID_source, theta, phi)\n",
    "    # recalculate model\n",
    "    LID_measurement = goniometer_model(x_offset, y_offset, z_offset, LID_source_function, theta, phi)[1]\n",
    "    LID_amplified = LID_source + amp * (LID_measurement - LID_source)\n",
    "    LID_amplified_cplane = half2cplanes(LID_amplified, theta, phi)[0]\n",
    "    # recalculate error\n",
    "    if 'rel' in err.lower():\n",
    "        error_prcnt = 100 * (LID_measurement - LID_source) / LID_source\n",
    "        error_label = r'Relative Error in %'\n",
    "    else:\n",
    "        error_prcnt = 100 * (LID_measurement - LID_source) / LID_source[0,0]\n",
    "        error_label = r'Absolute Error as % of $I_0$'\n",
    "    error_prcnt_cplane = half2cplanes(error_prcnt,theta,phi)[0]\n",
    "\n",
    "    # plot / update the curves\n",
    "    visu.plot(theta_cplane, LID_source_cplane[c], label=r'$C_0$-Plane of the source LID', c='lightgreen')\n",
    "    visu.plot(theta_cplane, LID_amplified_cplane[c], c='red', label=r'$C_0$-Plane of measured LID')\n",
    "\n",
    "    visu.clear_aux()\n",
    "    visu.plot_aux(theta_cplane, error_prcnt_cplane[c], label=error_label, c='orange')\n",
    "\n",
    "    # display the result in the chosen coordinate system\n",
    "    cart = True if 'cart' in axes.lower() else False\n",
    "    polar = True if 'polar' in axes.lower() else False\n",
    "    visu.show_fig(cart=cart, polar=polar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell defines the UI for parametrization. It uses different elements from the ```ipywidgets``` module and is mostly made up of Sliders, which are linked to the previously defined update function with the ```interactive_output``` call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interactive_output, IntSlider, FloatSlider, SelectionSlider, ToggleButtons, HBox, VBox, Label, Output\n",
    "c_max = np.size(phi_cplane)-1 # maximum index for c-plane-array\n",
    "phi_deg_cplane = np.rint(np.rad2deg(phi_cplane.squeeze())).astype(int)\n",
    "x_slider = FloatSlider(min=-0.5, max=0.5, step=0.05, value=0.0, continuous_update=False)\n",
    "y_slider = FloatSlider(min=-0.5, max=0.5, step=0.05, value=0.0, continuous_update=False)\n",
    "z_slider = FloatSlider(min=-0.5, max=0.5, step=0.05, value=0.1, continuous_update=False)\n",
    "c_slider = IntSlider  (min= 0, max=c_max, continuous_update=False, readout=False )\n",
    "# Add an Output for the C-Plane Angle in Degrees\n",
    "c_angle_output = Label(\"{:03}°\".format(phi_deg_cplane[c_slider.value]))\n",
    "c_slider.observe(lambda val: setattr(c_angle_output, 'value', \"{:03}°\".format(phi_deg_cplane[c_slider.value])), names='value')\n",
    "n_slider = IntSlider  (min= 1, max=120, continuous_update=False)\n",
    "amp_slider = FloatSlider(min=1, max=100, step=1, value=1, continuous_update=False)\n",
    "error_toggle = ToggleButtons(options=['Relative', 'Absolute'],description='Error Representation:')\n",
    "visu_style_toggle = ToggleButtons(options=['Cartesian', 'Polar', 'Cartesian & Polar'],description='Axes Style:')\n",
    "\n",
    "ui = VBox([HBox([Label('X-Offset as % of Measurement Distance'), x_slider]),\n",
    "           HBox([Label('Y-Offset as % of Measurement Distance'), y_slider]),\n",
    "           HBox([Label('Z-Offset as % of Measurement Distance'), z_slider]),\n",
    "           HBox([Label('Source LID Cosine Power'), n_slider]),\n",
    "           HBox([Label('C-Plane Index'), c_slider,c_angle_output]),\n",
    "           HBox([Label('Error Amplification'), amp_slider]),\n",
    "           error_toggle,\n",
    "           visu_style_toggle])\n",
    "\n",
    "out = interactive_output(update, {'x': x_slider, 'y': y_slider, 'z': z_slider, 'c': c_slider, 'n': n_slider,\n",
    "                                  'amp': amp_slider, 'err':error_toggle, 'axes':visu_style_toggle})\n",
    "display(ui, out)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this interface it is possible to explore the parameter space and determine where the error becomes problematic.\n",
    "In general the error increases for bigger offsets or synonymously smaller measurement distances.\n",
    "A higher gradient also increases the error. Especially the relative error takes on large values for the outer regions.\n",
    "What makes matters worse is that the influence from the gradient and offset are often compounding, since it is generally harder to estimate the Center of Light for focusing sources.\n",
    "For some procedures the object is positioned with the exit surface at the origin. That can result in a substantial offset along the optical axis.\n",
    "\n",
    "An interesting prospect of this model is that the error is reversible for a known offset. A goniometer measurement could be correct afterwards, if the translation from origin to source would be known."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uncertainty analysis\n",
    "A Luminous Intensity Distribution is constructed of multiple measurements. Each of these measurements is individually afflicted with an uncertainty. The uncertainty comprises various contributions. They stem from the mechanical measurement setup, the sensor and are also effected by the object to be measured.\n",
    "It is an ongoing effort to model the uncertainty of photogoniometric measurements. Currently most measurements lack an uncertainty evaluation entirely, which complicates a comparison of measurements.\n",
    "The uncertainty contribution of the sensor, either photometer or camera-based, can already be modelled. The same is true for the uncertainty contribution from the manipulator, that is positioning and orienting the source or detector. Since the location and orientation of the object on the goniometer is unknown, the uncertainty of the manipulator pose cannot be propagated to the source and its luminous intensity distribution and further through the sensor system. The Center of Light is therefore a missing link in the uncertainty chain, that prevents a holistic uncertainty evaluation. \n",
    "\n",
    "What is necessary to close this gap is to define a position and its uncertainty. This uncertainty contains both the confidence in the Center of Light and the confidence in the mounting accuracy. It will need to be rules of thumb on how to assign this uncertainty. One such rule could be to take the dimensions of the object as a first estimate for the region in which the Center of Light is assumed to be in.\n",
    "The following code uses the model function for an uncertainty analysis.\n",
    "This is done via the Monte Carlo method. Every uncertain input parameter is assigned a distribution. These distributions are then sampled many times, each time feeding them into the model and calculating the output. The large set of outcomes is then subjected to standard statistical analysis to determine the propagated uncertainty and its dependency on the input parameters.\n",
    "A Monte Carlo simulation is computationally expensive. The resolution of every single run and the number of iterations are therefore limited.\n",
    "To allow a simulation in a few seconds, we only simulate a reduced set of phi-angles. They correspond to the C0 and C90 plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_mc_deg = np.linspace(start=0,stop=360,num=5)[:-1]\n",
    "phi_mc = np.deg2rad(phi_mc_deg)[None,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to create distributions for the uncertain input parameters: the three coordinates of the source location.\n",
    "The uncertainty is described separately for each dimension. We get the distributions from the ```Chaospy``` module with the methods:\n",
    "```Uniform(lower limit, upper limit)``` or ```Normal(mean, std deviation)```\n",
    "depending on the desired distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import chaospy as cp\n",
    "# create the distributions for the location of the source\n",
    "x_dist = cp.Uniform(-0.001, 0.001)\n",
    "y_dist = cp.Uniform(-0.001, 0.001)\n",
    "z_dist = cp.Uniform(-0.001, 0.001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell defines and runs the Monte Carlo simulation. First a model is created from our model function. Then its parameters are defined with either distributions or static values. The model and its parameters are combined to a ```UncertaintyQuantification``` instance, which is subsequently run.\n",
    "All results are stored inside this instance. This also includes the input and output of every single iteration, which can take up a lot of memory. Since we do not need this data, it is removed from the data structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uncertainpy as un\n",
    "import numpy as np\n",
    "from Helperfunctions import polar2cart, cart2polar\n",
    "# run a monte carlo uncertainty quantification\n",
    "if __name__ == '__main__':\n",
    "    model = un.Model(run=goniometer_model)\n",
    "    parameters = {\"x_s\": x_dist, \"y_s\": y_dist,\"z_s\": z_dist,\n",
    "                  \"LID\": LID_source_function, \"theta_c\": theta, \"phi_c\": phi_mc}\n",
    "    UQ = un.UncertaintyQuantification(model=model, parameters=parameters)\n",
    "    UQ.quantify(method='mc', plot=None, save=False, nr_mc_samples=2**10, CPUs=4)\n",
    "    UQ.data['goniometer_model'].pop('evaluations')\n",
    "\n",
    "# save the result with UQ.save(filename, folder)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output data contains arrays of various statistical measures. The mean and standard deviation are extracted and prepared for a visualization.\n",
    "An estimate for the 95% confidence interval is constructed from the standard deviation by subtracting or adding it two times to the mean.\n",
    "All uncertainty values will be expressed as a expanded uncertainty with coverage factor k=2.\n",
    "\n",
    "As before the data is restructured into a C-Plane format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = UQ.data['goniometer_model']\n",
    "mean = np.array(results['mean'])\n",
    "std = np.sqrt(np.array(results['variance']))\n",
    "# enlarge the confidence region to make it visible\n",
    "amp = 10\n",
    "u_band_low = mean - 2 * std * amp\n",
    "u_band_high = mean + 2 * std * amp\n",
    "# restructure data to c-planes\n",
    "mean_cplanes, theta_cplane, phi_mc_cplane   = half2cplanes(mean,theta,phi_mc)\n",
    "u_band_low_cplanes   = half2cplanes(u_band_low,theta,phi_mc)[0]\n",
    "u_band_high_cplanes  = half2cplanes(u_band_high,theta,phi_mc)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a new visualization object for the uncertainty plots. The expanded uncertainty is plotted as a region around the mean using the ```fill``` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visu_mc = HalfplaneVisu(theta_max=theta_max,theta_tick_dist=15,r_max = 125,r_min=-50)\n",
    "c = 0\n",
    "\n",
    "visu_mc.plot(theta_cplane, mean_cplanes[c],label=r'$C_0$-Plane of the mean measured LID', color='blue')\n",
    "visu_mc.fill(theta_cplane,u_band_low_cplanes[c],u_band_high_cplanes[c], label='95% Interval', color='blue')\n",
    "visu_mc.show_fig(cart=True, polar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next two cells, the auxiliary axis is used to plot the expanded uncertainty. Just as the error, it can be given in relative or absolute terms. The relative expanded uncertainty is constructed as a percentage of the mean luminous intensity for a certain direction. The absolute expanded uncertainty is expressed as a percentage of the maximum luminous intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Relative representation of the Uncertainty\n",
    "u2_rel_prcnt = std * 200 / mean\n",
    "u2_rel_prcnt_cplanes  = half2cplanes(u2_rel_prcnt,theta,phi_mc)[0]\n",
    "\n",
    "visu_mc.aux_set_label('%')\n",
    "visu_mc.plot_aux(theta_cplane, u2_rel_prcnt_cplanes[c],label=r'Relative Uncertainty in %', color='orange')\n",
    "visu_mc.show_fig(cart=True, polar=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u2_abs_prcnt = std * 200 / mean[0,0]\n",
    "u2_abs_prcnt_cplanes  = half2cplanes(u2_abs_prcnt,theta,phi_mc)[0]\n",
    "\n",
    "visu_mc.clear_aux()\n",
    "visu_mc.plot_aux(theta_cplane, u2_abs_prcnt_cplanes[c],label=r'Absolute Uncertainty as % of $I_0$', color='orange')\n",
    "visu_mc.show_fig(cart=True, polar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can again easily be made interactive. The update function has similar parameters as before. The parameters x,y,z previously describing an offset now discribe a distribution. For simplicity this distribution is defined by a single value. The coordinates are set to get Normal distributions that are centered around zero. The parameters x,y,z describe the spread of this distribution expressed as two standard deviations, to be consistent with the expanded uncertainty. The value is again described as a percentage of the measurement distance. \n",
    "\n",
    "In case the offset or the LID-parameter is changed, the simulation has to run again, which is slow. Other changes are possible without a resimulation. Because of this, the update function can be called in two capacities: if all necessary parameters are provided, it runs the simulation, updates the plot and displays it. Otherwise it it skips the simulation and uses the last available data.\n",
    "The interface is constructed similarly as for the error analysis above. There is an additional button to start the simulation. This button gets assigned a callback function that calls the update_mc function with all parameters when the button is pressed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_mc(c, amp, u, style, x=None, y=None, z=None, n=None):\n",
    "    if not None in (x,y,z,n): # all necessary simulation parameters were supplied\n",
    "        # percent to factor -> /100\n",
    "        # 95% to std -> /2\n",
    "        x_dist = cp.Normal(0, x/200)\n",
    "        y_dist = cp.Normal(0, y/200)\n",
    "        z_dist = cp.Normal(0, z/200)\n",
    "        # update the source LID\n",
    "        LID_source_function = cosn_LID_generator(n)\n",
    "        parameters = {\"x_s\": x_dist, \"y_s\": y_dist,\"z_s\": z_dist,\n",
    "                      \"LID\": LID_source_function, \"theta_c\": theta, \"phi_c\": phi_mc}\n",
    "        UQ.parameters = parameters\n",
    "        UQ.quantify(method='mc', plot=None, save=False, nr_mc_samples=2**10, CPUs=4)\n",
    "        UQ.data['goniometer_model'].pop('evaluations')\n",
    "\n",
    "    results = UQ.data['goniometer_model']\n",
    "    \n",
    "    mean = np.array(results['mean'])\n",
    "    mean_cplanes, theta_cplane, phi_mc_cplane   = half2cplanes(mean,theta,phi_mc)\n",
    "    visu_mc.plot(theta_cplane, mean_cplanes[c],label=r'$C_0$-Plane of the mean measured LID', color='blue')\n",
    "    \n",
    "    # calculate confidence region, amplification for better visibility\n",
    "    std = np.sqrt(np.array(results['variance']))\n",
    "    u_band_low = mean - 2 * std * amp\n",
    "    u_band_high = mean + 2 * std * amp\n",
    "    u_band_low_cplanes   = half2cplanes(u_band_low,theta,phi_mc)[0]\n",
    "    u_band_high_cplanes  = half2cplanes(u_band_high,theta,phi_mc)[0]\n",
    "    visu_mc.fill(theta_cplane,u_band_low_cplanes[c],u_band_high_cplanes[c], label='95% Interval', color='blue')\n",
    "    \n",
    "    # calculate uncertainty plots\n",
    "    if 'rel' in u.lower():\n",
    "        u2_prcnt = std * 200 / mean\n",
    "        u2_label = r'Relative extended Uncertainty in %'\n",
    "    else:\n",
    "        u2_prcnt = std * 200 / mean[0,0]\n",
    "        u2_label = r'Absolute extended Uncertainty as % of $I_0$'\n",
    "    u2_prcnt_cplanes  = half2cplanes(u2_prcnt,theta,phi_mc)[0]\n",
    "\n",
    "    visu_mc.clear_aux()\n",
    "    visu_mc.plot_aux(theta_cplane, u2_prcnt_cplanes[c],label=u2_label, color='orange')\n",
    "\n",
    "    # display the result in the chosen coordinate system\n",
    "    cart = True if 'cart' in style.lower() else False\n",
    "    polar = True if 'polar' in style.lower() else False\n",
    "    visu_mc.show_fig(cart=cart, polar=polar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import Button, Output, fixed\n",
    "c_max = np.size(phi_mc_cplane)-1 # maximum index for c-plane-array\n",
    "phi_mc_deg_cplane = np.rint(np.rad2deg(phi_mc_cplane.squeeze())).astype(int)\n",
    "x_slider = FloatSlider(min=0, max=2, step=0.05, value=0.2, continuous_update=False)\n",
    "y_slider = FloatSlider(min=0, max=2, step=0.05, value=0.2, continuous_update=False)\n",
    "z_slider = FloatSlider(min=0, max=2, step=0.05, value=0.2, continuous_update=False)\n",
    "c_slider = IntSlider  (min= 0, max=c_max, continuous_update=False, readout=False)\n",
    "c_angle_output = Label(\"{:03}°\".format(phi_mc_deg_cplane[c_slider.value]))\n",
    "c_slider.observe(lambda val: setattr(c_angle_output, 'value', \"{:03}°\".format(phi_mc_deg_cplane[c_slider.value])), names='value')\n",
    "n_slider = IntSlider  (min= 1, max=120, value=n, continuous_update=False)\n",
    "amp_Slider = FloatSlider(min=1, max=10, step=0.1, value=amp, continuous_update=False)\n",
    "u_toggle = ToggleButtons(options=['Relative', 'Absolute'],description='Uncertainty Representation:')\n",
    "visu_style_toggle = ToggleButtons(options=['Cartesian', 'Polar', 'Cartesian & Polar'],description='Axes Style:',)\n",
    "run_button = Button(description=\"Start Simulation\")\n",
    "\n",
    "ui = VBox([HBox([Label('X-Direction one-sided expanded Uncertainty as % of Measurement Distance'), x_slider]),\n",
    "           HBox([Label('Y-Direction one-sided expanded Uncertainty as % of Measurement Distance'), y_slider]),\n",
    "           HBox([Label('Z-Direction one-sided expanded Uncertainty as % of Measurement Distance'), z_slider]),\n",
    "           HBox([Label('Source LID Cosine Power'), n_slider]),\n",
    "           run_button,\n",
    "           HBox([Label('C-Plane Index'), c_slider, c_angle_output]),\n",
    "           HBox([Label('Uncertainty Amplification'), amp_Slider]),\n",
    "           u_toggle,\n",
    "           visu_style_toggle])\n",
    "\n",
    "out = interactive_output(update_mc, {'c': c_slider, 'amp': amp_Slider, 'u':u_toggle, 'style': visu_style_toggle}) \n",
    "def run_simulation(b):\n",
    "    out.clear_output()\n",
    "    with out:\n",
    "        update_mc(c=c_slider.value, amp=amp_Slider.value, u=u_toggle.value, style=visu_style_toggle.value,\n",
    "                  x = x_slider.value, y=y_slider.value, z=z_slider.value, n=n_slider.value )\n",
    "run_button.on_click(run_simulation)\n",
    "display(ui,out)"
   ]
  }
 ],
 "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}