{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Defining Aircraft 'Speed'\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Measuring translations and velocities for a ground-based vehicle (your car) is relatively simple as the vehicle is moving with respect to an inertially-fixed reference frame (the Earth), that it is *always* connected to. To determine how long it will take to get somewhere, distance/time/velocity calculations are simple, and the engine performance is chiefly a function of the vehicle speed.\n", "\n", "By comparison, an aircraft is moving with respect to the inertially-fixed reference frame, but it is not connected to it - thus is it difficult to make measurements between the two. Furthermore, aerodynamic performance (how much lift and drag is being produced) is a function of the aircraft's orientation and translation with respect to the incident wind, and not to the Earth. Similarly, the powerplant performance will be a function of the incident wind and not the groundspeed.\n", "\n", "This gives us our first set of 'speeds', and the relationship between them:\n", "\n", "$$\\begin{aligned}\n", " V_{gnd} &= V + V_{wind}\\label{eq:Vgnd1}\n", "\\end{aligned}$$\n", "\n", "where\n", "\n", "$$\\begin{aligned}\n", "&V_{gnd} : \\text{Ground Speed}\\nonumber\\\\\n", " &V : \\text{True Airspeed (TAS)}\\nonumber\\\\\n", " &V_{wind} : \\text{Wind Speed}\\nonumber\\end{aligned}$$\n", "\n", "where we can simply add the two together via vector addition. In the two figures below, the two situations will require the exact same amount of thrust, and the same configuration of control surfaces - but the there is a groundspeed difference of 100kn.\n", "\n", "```{figure} ../Images/HeadWind.png\n", "---\n", "height: 300px\n", "name: HeadWind\n", "---\n", "Aircraft With Headwind\n", "```\n", "\n", "```{figure} ../Images/TailWind.png\n", "---\n", "height: 300px\n", "name: TailWind\n", "---\n", "Aircraft With Tailwind\n", "```\n", "\n", "\n", "This is why LHR-ORD takes typically 8h50m, but ORD-LHR takes 7h30m, yet in each case, the aircraft is flying at typical cruise velocity, requiring the same amount of total thrust.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "## Airspeed Measurement\n", "\n", "Clearly aerospace engineers require a means of determining the true airspeed, in order to determine aircraft powerplant/aerodynamic performance, and in order to calculate groundspeed and thus facilitate speed/time calculations.\n", "\n", "Equation {eq}`TASequation` may be developed from the isentropic flow equations, and allows us to determine **true airspeed** as a function of density and pressure:\n", "\n", "```{math}\n", ":label: TASequation\n", "V =\\sqrt{7\\cdot\\frac{p}{\\rho}\\cdot\\left\\{\\left(\\frac{\\Delta p}{p}+1\\right)^\\frac{\\gamma-1}{\\gamma}-1\\right\\}}\n", "```\n", "\n", "where \n", "\n", "$$\\rho = \\text{ Local air density}\n", "\n", "p = \\text{ Local static pressure}\n", "\n", "\\gamma = \\text{ Specific heat ratio = 1.4 for diatomic gases}\n", "\n", "\\Delta p = \\text{ Dynamic pressure/pressure difference} : p_o-p\n", "\n", "p_o = \\text{ Stagnation/total pressure} p_o = \\tfrac{1}{2}\\rho V^2 + p$$\n", "\n", "Equation {eq}`TASequation` is often simply given in textbooks with no derivation or explanation. You are encouraged to *understand* the derivation, but I won't test you on it.\n", "\n", "```{admonition} Click to show the derivation...\n", ":class: dropdown\n", "This is presented {cite}`Yechout:2014vg`, and other textbooks with no derivation given. It is easily derived from the isentropic flow equations.\n", "\n", "The following should be standard from your compressible flow classes, but I want to ensure you understand the legacy of the equations you're using - so you understand limitations and their development. This is adapted from my own undergraduate notes.\n", "\n", "We have the isentropic relationship showing that the ratio of total to static temperature is a function of the Mach number and the ratio of specific heats:\n", "\n", "$$\\frac{T_{0}}{T}=1+\\frac{\\gamma-1}{2} M^{2}$$ - (x)\n", "\n", "and we have the isentropic relations\n", "\n", "$$\\frac{p_{2}}{p_{1}}=\\left(\\frac{\\rho_{2}}{\\rho_{1}}\\right)^{\\gamma}=\\left(\\frac{T_{2}}{T_{1}}\\right)^{\\frac{\\gamma}{\\gamma-1}}$$\n", "\n", "since from total to stagnation we have an isentropic compression of the flow, we can state\n", "\n", "$$\n", "\\frac{p_{0}}{p}=\\left(\\frac{\\rho_{0}}{\\rho}\\right)^{\\gamma}=\\left(\\frac{T_{0}}{T}\\right)^{\\frac{\\gamma}{\\gamma-1}}$$\n", "\n", "the equations above can be combined to give\n", "\n", "$$\\frac{p_{0}}{p}=\\left(1+\\frac{\\gamma-1}{2} M^{2}\\right)^{\\frac{\\gamma}{\\gamma-1}}\\label{eq:ISEN1}$$\n", "\n", "the fraction may be expressed as:\n", "\n", "$$\\frac{p_{0}}{p} = \\frac{p+\\Delta p}{p}=\\frac{\\frac{p}{p}+\\frac{\\Delta p}{p}}{\\frac{p}{p}}=\\frac{\\Delta p}{p} + 1$$\n", "\n", "where $\\Delta p$ is the \\textsl{impact pressure}, and is the difference between total and static in a pitot-static - for incompressible flow, this would be the same as the dynamic pressure, but for compressible flow it is not - hence\n", "\n", "$$\\frac{\\Delta p}{p} + 1=\\left(1+\\frac{\\gamma-1}{2} M^{2}\\right)^{\\frac{\\gamma}{\\gamma-1}}$$\n", "\n", "$$\\frac{2}{\\gamma-1}\\left[\\left(\\frac{\\Delta p}{p}+1\\right)^{\\frac{\\gamma-1}{\\gamma}}-1\\right]=M^2=\\left(\\frac{V}{\\sqrt{\\gamma\\,R\\,T}}\\right)^2=\\frac{V^2}{{\\gamma\\,R\\,T}}$$\n", " \n", "$$V^2=\\frac{\\gamma\\,R\\,T\\cdot 2}{\\gamma-1} \\left[\\left(\\frac{\\Delta p}{p}+1\\right)^{\\frac{\\gamma-1}{\\gamma}}-1\\right]$$\n", " \n", "from the equation of state, we can express $R\\,T$ as $\\frac{p}{\\rho}$\n", "\n", "$$ V^2=\\frac{\\gamma}{\\gamma-1}\\cdot\\left(\\frac{p}{\\rho}\\right)^2\\cdot 2 \\left[\\left(\\frac{\\Delta p}{p}+1\\right)^{\\frac{\\gamma-1}{\\gamma}}-1\\right]$$\n", "\n", "the *heat capacity ratio*, $\\gamma$, is **exactly** 1.4 for diatomic gases - see footnote - so we can write\n", "\n", "$$V = \\sqrt{7\\cdot\\frac{p}{\\rho}\\left[\\left(\\frac{\\Delta p}{p}+1\\right)^{\\frac{\\gamma-1}{\\gamma}}-1\\right]}$$\n", "\n", "QED\n", " \n", " \n", "Footnote - Wikipedia isn't a valid source in most places, but I don't want to repeat a bit of derivation for an bit of a tangent - https://en.wikipedia.org/wiki/Heat_capacity_ratio\\#Relation_with_degrees_of_freedom\n", "```\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Equation {eq}`TASequation` allows us to determine the velocity with respect to a volume of air, via measurement of three quantities - density $\\rho$, static pressure $p$, and impact pressure $\\Delta p$. 'Impact pressure' is a term that is probably new to you, and refers to the difference between total and static pressure - and, for an incompressible fluid, is the same as dynamic pressure, but not for a compressible fluid.\n", "\n", "We have the relationship $V(p, \\rho, \\Delta p)$, but the first two quantities are relatively difficult to measure, as they are *absolute* quantities. Hence, transducers for measurement of absolute pressure and density would require regular calibration, and would add complexity to a measurement system. The reasons for this are complex, and will be elaborated upon in an instrumentation/measurement course - suffice to say that, in general, all transducers are relative transducers, and thus measurement of any absolute quantity require a datum against which they can be compared, calibrated for changes in non-measured quantities, such as temperature.\n", "\n", "By contrast, $\\Delta p$ is a *relative* quantity. Whilst a transducer for measurement will still require calibration, measurement of $\\Delta p$ is relatively easy - for this purpose, we use a pitot-static device. Originally these were analogue devices, measuring the pressure difference between two sides of a pressure chamber via a mechanical diaphragm, calibrated to provide an accurate reading of $\\Delta p$:\n", "\n", "\n", "```{figure} ../Images/PitotStatic.png\n", "---\n", "height: 300px\n", "name: PitotStatic\n", "---\n", "Aircraft Pitot-Static Measurement Device\n", "```\n", "\n", "\n", "### Calibrated Airspeed \n", "\n", "To avoid having to measure $p$ and $\\rho$, rather than using the pitot-static to measure $\\Delta p$, aircraft Engineers calibrated their devices to provide the correct value of true airspeed at sea-level ISA density and pressure. This provides us with **calibrated airspeed**, $V_c$/CAS, which is only the same as true airspeed at sea-level ISA conditions. $V_c$ can be calculated via Equation {eq}`CASequation`:\n", "\n", "```{math}\n", ":label: CASequation\n", " V_c=\\sqrt{7\\cdot\\frac{p_{SL}}{\\rho_{SL}}\\cdot\\left\\{\\left(\\frac{\\Delta p}{p_{SL}}+1\\right)^\\frac{\\gamma-1}{\\gamma}-1\\right\\}}\n", "```\n", "\n", "### Indicated Airspeed\n", "\n", "In reality, such a pitot-static device will measure the **calibrated airspeed** at a point on the aircraft surface, where the flow has already been disturbed by the aircraft. The device is NOT measuring the freestream velicity - thus there is a *position error*, based on where the pitot-static is on the aircraft. During flight testing, these inaccuracies can be quantified, and a difference between what is *indicated* by the device, **indicated airspeed**/IAS, and calibrated airspeed is determined - $\\Delta V_p$, the position error.\n", "\n", "```{math}\n", ":label: IASequation\n", "V_c = V_I + \\Delta V_p\n", "```\n", "\n", "At this stage - we have three airspeeds; Indicated (what is shown on the measurement device, for ISA SL conditions), Calibrated (what _should_ be shown on the measurement device, for ISA SL conditions), and True Airspeed. We don't have a means to convert Calibrated to True yet, so we adopt a two-step approach.\n", "\n", "### Equivalent Airspeed\n", "\n", "The first correction is for the actual pressure at a given altitude. This gives Equivalent Airspeed, EAS.\n", "\n", "```{math}\n", ":label: EASequation\n", " V_e=\\sqrt{7\\cdot\\frac{p}{\\rho_{SL}}\\cdot\\left\\{\\left(\\frac{\\Delta p}{p}+1\\right)^\\frac{\\gamma-1}{\\gamma}-1\\right\\}}\n", "```\n", "\n", "where we can see that the sea-level pressure has been replaced with the actual pressure. In practice, this correction is applied as a multiplier between CAS and EAS:\n", "\n", "$$V_e = f\\cdot V_c$$\n", "\n", "where\n", "\n", "$$ f = \\frac{V_e}{V_c} $$\n", "$$ = \\frac{\\sqrt{7\\cdot\\frac{p}{\\rho_{SL}}\\cdot\\left\\{\\left(\\frac{\\Delta p}{p}+1\\right)^\\frac{\\gamma-1}{\\gamma}-1\\right\\}}}{\\sqrt{7\\cdot\\frac{p_{SL}}{\\rho_{SL}}\\cdot\\left\\{\\left(\\frac{\\Delta p}{p_{SL}}+1\\right)^\\frac{\\gamma-1}{\\gamma}-1\\right\\}}}$$\n", "$$= f(\\Delta p, p)$$\n", " \n", "Since $f$ is only dependent on the $\\Delta p$ the aircraft speed, and the $p$ the aircraft altitude, it can be calculated and tabulated - see the table below for the pressure correction factor:\n", "\n", "\n", "| h, ft | Calibrated Airspeed, kn | | | | | | | | |\n", "|:------|:------|:------|:------|:------|:------|:------|:------|:------|:------|\n", "| | 100 | 125 | 150 | 175 | 200 | 225 | 250 | 275 | 300 |\n", "| 5000 | 0.999 | 0.999 | 0.999 | 0.998 | 0.998 | 0.997 | 0.997 | 0.996 | 0.995 |\n", "| 10000 | 0.999 | 0.998 | 0.997 | 0.996 | 0.995 | 0.994 | 0.992 | 0.991 | 0.989 |\n", "| 15000 | 0.998 | 0.997 | 0.995 | 0.994 | 0.992 | 0.990 | 0.987 | 0.985 | 0.982 |\n", "| 20000 | 0.997 | 0.995 | 0.993 | 0.990 | 0.987 | 0.984 | 0.981 | 0.977 | 0.973 |\n", "| 25000 | 0.995 | 0.993 | 0.990 | 0.986 | 0.982 | 0.978 | 0.973 | 0.968 | 0.963 |\n", "| 30000 | 0.993 | 0.990 | 0.986 | 0.981 | 0.975 | 0.970 | 0.963 | 0.957 | 0.950 |\n", "| 35000 | 0.991 | 0.986 | 0.981 | 0.974 | 0.967 | 0.959 | 0.951 | 0.943 | 0.934 |\n", "| 40000 | 0.988 | 0.982 | 0.974 | 0.966 | 0.957 | 0.947 | 0.937 | 0.926 | 0.916 |\n", "| 45000 | 0.984 | 0.976 | 0.966 | 0.956 | 0.944 | 0.932 | 0.920 | 0.907 | 0.895 |\n", "| 50000 | 0.979 | 0.969 | 0.957 | 0.944 | 0.930 | 0.915 | 0.901 | 0.886 | 0.871 |\n", "\n", "Intermediate values are found via linear interpolation - you should be comfortable with doing this by hand, or by writing code. An example code is provided below that provides two-dimensional interpolation for any values - you can use this function to find out the pressure correction factor - use the syntax ```f_correction(VC, h)``` (both of these have default values, so you don't need to provide either).\n" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "tags": [ "hide-input", "thebe-init" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.999 0.999 0.999 0.998 0.998 0.997 0.997 0.996 0.995]\n", " [0.999 0.998 0.997 0.996 0.995 0.994 0.992 0.991 0.989]\n", " [0.998 0.997 0.995 0.994 0.992 0.99 0.987 0.985 0.982]\n", " [0.997 0.995 0.993 0.99 0.987 0.984 0.981 0.977 0.973]\n", " [0.995 0.993 0.99 0.986 0.982 0.978 0.973 0.968 0.963]\n", " [0.993 0.99 0.986 0.981 0.975 0.97 0.963 0.957 0.95 ]\n", " [0.991 0.986 0.981 0.974 0.967 0.959 0.951 0.943 0.934]\n", " [0.988 0.982 0.974 0.966 0.957 0.947 0.937 0.926 0.916]\n", " [0.984 0.976 0.966 0.956 0.944 0.932 0.92 0.907 0.895]\n", " [0.979 0.969 0.957 0.944 0.93 0.915 0.901 0.886 0.871]]\n" ] } ], "source": [ "import numpy as np\n", "import scipy.interpolate\n", "from scipy.interpolate import interp2d\n", "presscorr = np.array([[5000,0.999,0.999,0.999,0.998,0.998,0.997,0.997,0.996,0.995],\n", " [10000,0.999,0.998,0.997,0.996,0.995,0.994,0.992,0.991,0.989],\n", " [15000,0.998,0.997,0.995,0.994,0.992,0.99,0.987,0.985,0.982],\n", " [20000,0.997,0.995,0.993,0.99,0.987,0.984,0.981,0.977,0.973],\n", " [25000,0.995,0.993,0.99,0.986,0.982,0.978,0.973,0.968,0.963],\n", " [30000,0.993,0.99,0.986,0.981,0.975,0.97,0.963,0.957,0.95],\n", " [35000,0.991,0.986,0.981,0.974,0.967,0.959,0.951,0.943,0.934],\n", " [40000,0.988,0.982,0.974,0.966,0.957,0.947,0.937,0.926,0.916],\n", " [45000,0.984,0.976,0.966,0.956,0.944,0.932,0.92,0.907,0.895],\n", " [50000,0.979,0.969,0.957,0.944,0.93,0.915,0.901,0.886,0.871]])\n", "\n", "h = presscorr[:, 0]\n", "VC = np.array([100, 125, 150, 175, 200, 225, 250, 275, 300])\n", "table = presscorr[:, 1:]\n", "\n", "\n", "\n", "f_corr = interp2d(VC, h, table, kind='linear')\n", "\n", "def f_correction(VC=125, h=45000, ):\n", " f = f_corr(VC, h)[0]\n", " return f\n" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "For a CAS of 182, at an altitude of 39500 feet, the pressure correction factor is 0.964336\n" ] } ], "source": [ "VC = 182\n", "h = 39500\n", "f = f_correction(VC, h)\n", "print(f\"For a CAS of {VC}, at an altitude of {h} feet, the pressure correction factor is {f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "### True Airspeed\n", "\n", "Finally, **True Airspeed** may be calculated from **Equivalent Airspeed**, by accounting for the actual density at the correct altitude:\n", "\n", "$$V = V_e\\cdot\\sqrt{\\frac{\\rho_{SL}}{\\rho}}$$\n", "\n", "$$V = V_e\\cdot\\sqrt{\\frac{1}{\\sigma}}$$\n", " \n", "$$\\sigma = \\frac{\\rho}{\\rho_{SL}}$$\n", "\n", "Where we may get the density at different altitudes from standard tables of the atmospheric properties. These are available in the original PDF notes. What we'll use more regularly is either ```atmosisa``` in MATLAB or ```ambiance``` in Python.\n", "\n", "If you can't get the following to work - then you might need to install ambiance. Try ```pip install ambiance``` in your terminal window. See here for more instructions: https://pypi.org/project/ambiance/\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At sea level, the pressure is 101325.000Pa, the density is 1.225kg/m^3, the tempertaure is 15.0C, and the viscosity is 1.461e-05kg/s\n", "\n", "At an altitude of 12000m, the pressure is 19399.392Pa,the density is 0.312kg/m^3,the tempertaure is -56.5C,and the viscosity is 4.557e-05kg/s\n" ] } ], "source": [ "from ambiance import Atmosphere\n", "\n", "sealevel = Atmosphere(0)\n", "\n", "h = 12000\n", "\n", "altmospere = Atmosphere(h)\n", "\n", "print(f\"At sea level, the pressure is {sealevel.pressure[0]:6.3f}Pa,\\\n", " the density is {sealevel.density[0]:1.3f}kg/m^3,\\\n", " the tempertaure is {sealevel.temperature[0]-273.15:1.1f}C, \\\n", " and the viscosity is {sealevel.kinematic_viscosity[0]:1.3e}kg/s\")\n", "\n", "print(\"\")\n", "\n", "print(f\"At an altitude of {h}m, the pressure is {altmospere.pressure[0]:6.3f}Pa,\\\n", "the density is {altmospere.density[0]:1.3f}kg/m^3,\\\n", "the tempertaure is {altmospere.temperature[0]-273.15:1.1f}C,\\\n", "and the viscosity is {altmospere.kinematic_viscosity[0]:1.3e}kg/s\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using ```ambiance``` to define the density, we can create a function, ```sigma_density``` (it's a good idea to make function handles descriptive, and there tend to be lots of sigma functions that mean important things in programming languages). This function will return the density ratio for a given altitude provided in metres or feet - bearing in mind than the units for ```ambiance``` are metres." ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At 39500ft, the density ratio is 0.253064\n", "[216.65]\n" ] } ], "source": [ "from ambiance import Atmosphere\n", "\n", "def sigma_density(h=10000, alt_units='m'):\n", " if alt_units == 'm':\n", " factor = 1\n", " elif alt_units == 'km':\n", " factor = 1000\n", " elif alt_units == 'ft':\n", " factor = 0.3048\n", " elif alt_units == 'miles':\n", " factor = 1609.34\n", " \n", " sea_level = Atmosphere(0)\n", " atmospheric_properties = Atmosphere(h*factor)\n", " \n", " sigma = atmospheric_properties.density / sea_level.density\n", " return sigma\n", " \n", "h_in = 39500\n", "\n", "h_unit = 'ft'\n", "\n", "print(f\"At {h_in}{h_unit}, the density ratio is {sigma_density(h_in, h_unit)[0]:1.6f}\")\n", "\n", "atm = Atmosphere(39000*.3048)\n", "print(atm.temperature)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} A word about atmospheres\n", ":class: dropdown\n", "\n", "The ambiance module for Python implments the ICAO (International Civial Aviation Organisation) atmopshere, which is an extension of ISA (International Standard Atmosphere).\n", "\n", "Printed in the old notes is the US Standard Atmosphere, which is another extension of the ISA model.\n", "\n", "Both are equal in the ranges used in this course - but in the printed tables, the values have been truncated. This can lead to differences when interpolating from the table, and when using the routines. See the code dropdown below to understand this further.\n", "```" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "tags": [ "hide-input", "hide-output" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The density at 35,000ft from the table is 7.3800e-04slugs/ft**3 and from ICAO it is 7.3821e-04slugs/ft**3\n", "The density at 40,000ft from the table is 5.8700e-04slugs/ft**3 and from ICAO it is 5.8728e-04slugs/ft**3\n", "\n", "The density at 39,500ft from ICAO is 6.0151e-04slugs/ft**3 and interpolated between 35000 and 40000 ICAO it is 6.0237e-04slugs/ft**3\n", "\n", "Sigma from the table is 0.2533, \n", " from ambiance (at 39000) is is 0.2531, \n", " and from ambiance (interpolated between 35000 and 40000) is is 0.2534\n" ] } ], "source": [ "# Get four different atmospheres\n", "atmo_sl = Atmosphere(0*.3048)\n", "atmo_35kft = Atmosphere(35000*.3048)\n", "atmo_40kft = Atmosphere(40000*.3048)\n", "atmo_395kft = Atmosphere(39500*.3048)\n", "\n", "\n", "# Conversions\n", "kg_to_slugs = 1/14.59390\n", "m_to_feet = 3.28084\n", "\n", "# Convert the densities\n", "dens_imperial_sl = atmo_sl.density * kg_to_slugs / m_to_feet**3\n", "dens_imperial_35kft = atmo_35kft.density * kg_to_slugs / m_to_feet**3\n", "dens_imperial_40kft = atmo_40kft.density * kg_to_slugs / m_to_feet**3\n", "dens_imperial_395kft = atmo_395kft.density * kg_to_slugs / m_to_feet**3\n", "\n", "print(f\"The density at 35,000ft from the table is {7.38*1e-4:1.4e}slugs/ft**3 and from ICAO it is {dens_imperial_35kft[0]:1.4e}slugs/ft**3\")\n", "print(f\"The density at 40,000ft from the table is {5.87*1e-4:1.4e}slugs/ft**3 and from ICAO it is {dens_imperial_40kft[0]:1.4e}slugs/ft**3\")\n", "print(\"\")\n", "\n", "# Interpolate the densities\n", "dens_interpolated_39500_ambiance = dens_imperial_35kft + 4500/5000 * (dens_imperial_40kft - dens_imperial_35kft)\n", "dens_interpolated_39500_table = 7.38e-4 + 4500/5000 * (5.87e-4 - 7.38e-4)\n", "\n", "print(f\"The density at 39,500ft from ICAO is {dens_imperial_395kft[0]:1.4e}slugs/ft**3 and interpolated between 35000 and 40000 ICAO it is {dens_interpolated_39500_ambiance[0]:1.4e}slugs/ft**3\")\n", "\n", "print(\"\")\n", "# Get sigmas\n", "sigma_ambiance_direct = (atmo_395kft.density/atmo_sl.density)[0]\n", "sigma_table = dens_interpolated_39500_table/(23.77*1e-4)\n", "sigma_ambiance = dens_interpolated_39500_ambiance/dens_imperial_sl\n", "\n", "\n", "print(f\"Sigma from the table is {sigma_table:1.4f}, \\n from ambiance (at 39000) is is {sigma_ambiance_direct:1.4f}, \\n and from ambiance (interpolated between 35000 and 40000) is is {sigma_ambiance[0]:1.4f}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inputs can also be ```numpy``` arrays " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At an altitude of 1 mile, the density ratio is 0.8544179176392871\n", "At an altitude of 2 miles, the density ratio is 0.7256946830480345\n", "At an altitude of 3 miles, the density ratio is 0.6124148561281167\n", "At an altitude of 4 miles, the density ratio is 0.5132337621724308\n", "At an altitude of 5 miles, the density ratio is 0.42687647082233393\n", "At an altitude of 6 miles, the density ratio is 0.3521367607463307\n", "At an altitude of 7 miles, the density ratio is 0.28579654351527745\n", "At an altitude of 8 miles, the density ratio is 0.22195323343439954\n", "At an altitude of 9 miles, the density ratio is 0.17239372263976577\n" ] } ], "source": [ "alts = np.arange(1, 10)\n", "sigmas = sigma_density(alts, 'miles')\n", "\n", "for alt, sigma in zip(alts, sigmas):\n", " if alt == 1:\n", " print(f\"At an altitude of {alt} mile, the density ratio is {sigma}\")\n", " else:\n", " print(f\"At an altitude of {alt} miles, the density ratio is {sigma}\")\n", " \n", "# Note that the syntax to write mile/miles for single vs. any other quatity is completely unnecessary but nice" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "#### Approximate density correction\n", "\n", "For altitudes below about 16km, the approximation can be used\n", "\n", "$$\\sigma=\\frac{20-H}{20+H}$$\n", "\n", "where $H$ is the altitude in $km$.\n", "\n", "The error in this approximation is easily shown using the functions created above:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3wVVfr48c+TTkISIIQSAoTeq6EJIkoR2K9gwaVYQF2wY1d0lXXR3Z8N24IFEUVEqg0Li6uCiooQeu8tlBBCEhIgkHJ+f8wELiE9N5nc5Hm/XveVe2fOnXnmzs0895wzc0aMMSillFLu5OV0AEoppSoeTS5KKaXcTpOLUkopt9PkopRSyu00uSillHI7TS5KKaXcTpOLKhdEpI+IxDodR0mJSAMRSRURb6djKSo77sZOx+FpRORdEXnW6TjKG00uHkJE9onIGfsAkP2YUorrK3cHe7HsEZEtTseSF2PMAWNMVWNMpruXLSLPicgnJXh/NRGZISJHRSRFRHaIyJPZ8+2497gn2pIRkVEiEmN/z4+IyGIR6VUO4hojIstdpxlj7jbGPO9UTOWVj9MBqCK51hjzg9NBOKg3UAvwEZEuxphVpbESEfExxmSUxrId9joQBLQCkoHmQFtHI8qFiDwCTADuBpYA54CBwFBgeT5vzW1Zl+zLCrx/yxdjjD484AHsA/rlMe8dYKHL65eAHwEBqgPfAPFAov080qVsDeBD4LA9/0usA9AZIAtItR8Ruaz3L8Ba4CRwEHjOZV4UYIDRwAHgOPB3l/lVgI/sdW4BHgdiC/gMZgCzgc+BKTnmLQP+H7AS68D5FVAjRyzj7O08Ajzq8t7ngIXAJ/a2/A3wB96wyx+2n/vb5Z8EVgA+9ut7gM1AgMu6fFziegH43f4cvwbC7O04CawColxiedP+LE8Cq4Er7OkDsQ6y6fZy1tvTQ4EP7G06ZK/LO4/PbxNwXT6frwGa2s/D7FizY3wBWJ6j7L3ATiAFeB5oAvxhv2c+4GeXzfc7mCOGUHv7bsonzvz2TR8g1t5HR4FZuU2zy/4fsA5IsvdPe5d11Mf6nsUDCcAUrKScBmTaMSbZZT8CXnB571hgF3ACWITL/479ud1tf26JwFRA7HlNgZ+xvr/HgXlOH3dKdMxyOgB9FHJH5Z9cAoEdwBjgCvuLGWnPCwNutMsEAwuAL13e+y0wzz4A+AJX2tP7UPDBvg/QDqt5tT0Qh33w4sJB9n2sRNIBOAu0sue/CPyKldzqYx348lyfHf9JYLC9PcexD172/GVYB9e2WMnxM+CTHLHMsee1sw8a/ez5z2EdtK+zt6UKMAkrgdQCwu2Dz/N2eS/gF/t9zeyDRKcc63JNLruwDryhWIl0B9APq+XgY+BDl+24xd5nPsCjWAfDAJc4P8nxuXwJvGdvVy2s5HpXHp/hdKwkeDvQLJf5rsllrv0IBFpjJbycyWUREAK0sfftj0Bjl+0cXZjvYI4YBgIZ2Z9fHmXy2zd97Pe/hJWEquQxrTNwDOgGeGP9CNpnz/cG1nOhphcA9LKXP8b1c7CnfYSdXICrsb6bne1l/Qf4Jcfn9g1QDWiA9T0caM+bA/wd6/t1fp2e+nA8AH0UckdZX/xUrF9Z2Y+xLvO7Yv1S2g+MzGc5HYFE+3ldrNpJ9VzK9aGA5JLLe94AXrefR9n/SK61pJXACPv5nux/Kvv1uPzWh3XQjcc66Prb23+9y/xlwIsur1tj/dL3domlpcv8l4EP7OfPuR4A7Gm7gcEur68B9rm8jrI/763AUzmm50wurjW2ycBil9fXAuvy2e5EoINLnJ+4zKuNdVCv4jJtJLA0j2VVAZ7GqhGlYyW9QS7zDdavZ297fguXebnVXHq6vF4NPJljO98o6DuYy7ybgaMFfM/y3Df29/YcdkLOZ9o72AnJZdp24EqgR/Z3LZd1jyH/5PIB8LLLvKr2Zxnl8rn1cpk/H5hgP/8YmEYetTpPe2iHvme5zhhTzeXxfvYMY8xKrAO2YH1hARCRQBF5T0T2i8hJrF/c1eyzmeoDJ4wxicUJRkS6ichSEYkXkWSs6n7NHMWOujw/jfXPBhCB9Ws42/4CVjcamG+MyTDGnMVqshido0zO5fnmiCfn/Ig85mXH5xrTReWNMfuApVjJZGoBsce5PD+Ty+vszwQReVREtopIsogkYdUCcn6m2RpibeMREUmyy7+H9Yv+EsaYM8aYfxtjLsOqTcwHFohIjRxFw7GSuOtnkvPzKfR2FfAdzCkBqCki+fUH57tvgHhjTFqO9+Sc1hB4NPtzsz+7+vZy6gP7TfH6ZS6KzRiTirVN9VzK5PU/8QTW/+9KEdksIncUY/3lhiaXCkJE7sP6RX8Y60ua7VGgBdDNGBOC1SkO1pf4IFBDRKrlskhTiNV+itU0Ut8YEwq8ay+3MI5g/RNna5BXQRGJxGpuuMU+0+koMAwYLCKuB96cy0vHaqLIa/5hl9c5t/cw1gEo1/IiMhjrF+6PwCt5xV4UInIFVr/AX7Fqk9Ww2t+zP9OcMR7EqrnUdPnBEWKMaVPQuowxJ4F/YzX7NMoxOx6rGSnSZVp9ii+/72BOf2D1a1yXz/Ly3Tfk/t3N7bP7V44fa4HGmDn2vAZ5JLiC/i8uik1EgrAS+aEC3ocx5qgxZqwxJgK4C3hbRJoW9L7ySpNLBSAizbGaLW4BbgWeEJGO9uxgrF+RSfYv1H9kv88YcwRYjPUlri4iviKS/Y8fB4SJSGg+qw7GqvmkiUhXYFQRwp4PPGWvNxJ4IJ+yt2L1U7TAalLpiHWmUyxWM1C2W0SktYgEYrXLLzQXnxL8rP0rug1Wv8O8fNY5B3hGRMLtBDYRq8Mf+/UHWB3/o4Fr7WRTUsFYB/V4rDPiJmL1aWSLA6JExAvO77/vgckiEiIiXiLSRESuzG3hIvKsiHQRET8RCQAexGpe3O5azv7MPgeesz+vlsBtJdyuXL+DORljkrE+66kicp29fl8RGSQiL9vF8tw3RfA+cLdd+xYRCRKRv4hIMFbz7RHgRXt6gIj0tN8XB0SKiF8ey/0UuF1EOoqIP1YC/9Ou6eZLRG6y/xfAag41WCcPeCRNLp7l6xzXuXxh/7r6BHjJGLPeGLMTq119lv3lfgOrrf04Vifof3Ms81asX/jbsDo4HwIwxmzD+ifeYzcbRHCpe4FJIpKC9Q8+P5cyefknVvPBXqwD5Kx8yo4G3rZ/2Z1/YNWUXJvGZmG1fx/F6hAdn2M5P2P1M/wIvGqM+T6fdb4AxAAbgI3AGnsaWO3iXxljvjPGJAB3AtNFJCz/TS7QEqxkvwPrs0nj4uaoBfbfBBFZYz+/DfDD6kBPxDrrrW4eyzdYZwYex/qF3R/4i910k9P9WE1y2WdczcGqJRVHQd/Bi4M05jXgEeAZrER70I7nS7tIfvumUIwxMVhndU3B+tx2YfWnZCfXa7H6nw5g/YgZbr/1J6yTIo6KyHFyMMb8CDyLdULJEawTOUYUMqwuwJ8ikorVIvCgMWZvUbarPMk+BU4pjyYiy7A6u6fnMi8KK4n5FrMdvdITkZeAOsaYnP1cSuVKay5KqUuISEsRaW83GXXFqp194XRcynPoFfpKqdwEYzWFRWA1l07GujBVqULRZjGllFJup81iSiml3E6bxYCaNWuaqKgop8NQSimPsnr16uPGmPDc5mlyAaKiooiJiXE6DKWU8igikufIGtosppRSyu00uSillHI7TS5KKaXcTvtclFJlLj09ndjYWNLScg5erMqjgIAAIiMj8fX1LfR7NLkopcpcbGwswcHBREVFIVLYgbSVE4wxJCQkEBsbS6NGOQfQzpvHNYuJyAwROSYim/KYLyLylojsEpENItK5rGNUSuUvLS2NsLAwTSweQEQICwsrci3T45IL1qi3A/OZPwjr1rPNsO5u+E4ZxKSUKiJNLJ6jOPvK45KLMeYXrNvL5mUo8LGxrMC6411eQ5CXSFp6Ji//dxv7j+UXjlJKVT4el1wKoR4X3wMjlotvMQqAiIwTkRgRiYmPjy/WipJOp9Pgj2c5Of16sjKzihetUsoR3t7edOzY8fzjxRdfdNuy161bx3fffee25Xmiitihn1v97ZLROY0x07Bu+kR0dHSxRu+sExpA03bdaLfpeX77cio9b8zvZopKqfKkSpUqrFu3rlSWvW7dOmJiYhg82B03KPVMFbHmEsvF9/uO5OL7a7vVZTc8zA6/1rTe8BJHjsSW1mqUUmUgOTmZFi1asH27defnkSNH8v777wNwzz33EB0dTZs2bfjHPy7cqXnVqlVcfvnldOjQga5du5KcnMzEiROZN28eHTt2ZN68/O6mXXFVxJrLIuB+EZkLdAOS7XuNlwrx8ib4prep+klfNn/yEHUeW6AdlUoVwT+/3syWwyfduszWESH849o2+ZY5c+YMHTt2PP/6qaeeYvjw4UyZMoUxY8bw4IMPkpiYyNixYwH417/+RY0aNcjMzKRv375s2LCBli1bMnz4cObNm0eXLl04efIkgYGBTJo0iZiYGKZMmeLW7fIkHpdcRGQO0AeoKSKxwD8AXwBjzLvAd8BgrHtinwZuL+2Y6jbrxPomd9Bu92y++WM9117eseA3KaUclVezWP/+/VmwYAH33Xcf69evPz99/vz5TJs2jYyMDI4cOcKWLVsQEerWrUuXLl0ACAkJKbP4yzuPSy7GmJEFzDfAfWUUznntRjzPXe/34c/v4+jWNo1aIQFlHYJSHqmgGkZZy8rKYuvWrVSpUoUTJ04QGRnJ3r17efXVV1m1ahXVq1dnzJgxpKWlYYzRloo8VMQ+F0d4+VXhqb/24VxGBu/PXYje4VMpz/T666/TqlUr5syZwx133EF6ejonT54kKCiI0NBQ4uLiWLx4MQAtW7bk8OHDrFq1CoCUlBQyMjIIDg4mJSXFyc1wnCYXN2ocXpWZzf/giUPj+fnXZU6Ho5TKR3afS/ZjwoQJ7Nixg+nTpzN58mSuuOIKevfuzQsvvECHDh3o1KkTbdq04Y477qBnz54A+Pn5MW/ePB544AE6dOhA//79SUtL46qrrmLLli2VukNf9Be2dSqyu24WlpmaQMrkThwyNanzyK+EhQS5ZblKVSRbt26lVatWToehiiC3fSYiq40x0bmV15qLm3lXDeN0vxdpw25+mfkPbR5TSlVKmlxKQUTPUewK78vg4x+y9NdfnA5HKaXKnCaXUtLotnc56NuIGT9t4Giy3rNCKVW5aHIpJd7BtfC5axmrM5vx+ML12jymlKpUNLmUoqjwqjwzsDHd9k5l0Y/LnA5HKaXKjCaXUjaqfSij/ZYS9etj7Dvm3iEulFKqvNLkUsokuDYZA1+mg+zi54+fIzNLm8eUKi+++OILRIRt27Y5FsPEiRP54YcfSrycpKQk3n777fOvDx8+zLBhw0q83OLS5FIGqncZweG6/RiR8jELvlvidDhKKducOXPo1asXc+fOdcvyMjIyivyeSZMm0a9fvxKvO2dyiYiIYOHChSVebnFpcikLItQd9TbnfIJoufLvbDmU7HRESlV6qamp/Pbbb3zwwQfnk8uyZcvo3bs3119/Pa1bt+buu+8mK8u6EWDVqlV59NFH6dy5M3379iX7JoN9+vTh6aef5sorr+TNN99k//799O3bl/bt29O3b18OHDgAwNChQ/n4448BeO+997j55psBGDNmzPkkEBUVxdNPP02PHj2Ijo5mzZo1XHPNNTRp0oR33333fNx9+/alc+fOtGvXjq+++gqACRMmsHv3bjp27Mjjjz/Ovn37aNu2LQBpaWncfvvttGvXjk6dOrF06VIAPvroI2644QYGDhxIs2bNeOKJJ9z2+XrcwJWeSoJrY26cwf/7bD+J89ex6P5eBPh6Ox2WUuXDh3+5dFqb66DrWDh3GmbfdOn8jqOg081wKgHm33bxvNu/LXCVX375JQMHDqR58+bUqFGDNWvWALBy5Uq2bNlCw4YNGThwIJ9//jnDhg3j1KlTdO7cmcmTJzNp0iT++c9/nh9SPykpiZ9//hmAa6+9lttuu43Ro0czY8YMxo8fz5dffsm0adPo2bMnjRo1YvLkyaxYsSLXuOrXr88ff/zBww8/zJgxY/jtt99IS0ujTZs23H333QQEBPDFF18QEhLC8ePH6d69O0OGDOHFF19k06ZN50d63rdv3/llTp06FYCNGzeybds2BgwYwI4dOwDrxmZr167F39+fFi1a8MADD1C/fv1L4ioqrbmUoZDWfbl3+LXsiEvl9UW5f7GUUmVjzpw5jBgxAoARI0YwZ84cALp27Urjxo3x9vZm5MiRLF++HAAvLy+GDx8OwC233HJ+OnB+OsAff/zBqFGjALj11lvPl6tduzaTJk3iqquuYvLkydSoUSPXuIYMGQJAu3bt6NatG8HBwYSHhxMQEEBSUhLGGJ5++mnat29Pv379OHToEHFxcflu6/Lly7n11lsBa7DNhg0bnk8uffv2JTQ0lICAAFq3bs3+/fuL8CnmTWsuZezK5uF8ErWEqA3fsqzpEvp0aOp0SEo5L7+ahl9g/vODwgpVU3GVkJDATz/9xKZNmxARMjMzEREGDx58yRD6eQ2p7zo9KCjvMQRdy23cuJGwsDAOH8775rj+/v6Alcyyn2e/zsjIYPbs2cTHx7N69Wp8fX2JiooiLS3/C7Xzu87OdR3e3t7F6jfKjdZcHNB14CjqyglSvnyEYyf16n2lytrChQu57bbb2L9/P/v27ePgwYM0atSI5cuXs3LlSvbu3UtWVhbz5s2jV69egHWfl+y+kU8//fT89Jwuv/zy8304s2fPPl9u5cqVLF68mLVr1/Lqq6+yd+/eYsWenJxMrVq18PX1ZenSpedrGvkN89+7d29mz54NwI4dOzhw4AAtWrQo1voLS5OLA/yiupMcPZ5rzc/MnfkfsvT0ZKXK1Jw5c7j++usvmnbjjTfy6aef0qNHDyZMmEDbtm1p1KjR+XJBQUFs3ryZyy67jJ9++omJEyfmuuy33nqLDz/8kPbt2zNr1izefPNNzp49y9ixY5kxYwYRERFMnjyZO+64o1gjd9x8883ExMQQHR3N7NmzadmyJQBhYWH07NmTtm3b8vjjj1/0nnvvvZfMzEzatWvH8OHD+eijjy6qsZQGHXIf9w65X2iZ6Rx/qw/eSfv4rtdn3Ny/e9muXykHldch95ctW8arr77KN998c8m8qlWrkpqa6kBU5YMOue8pvH0Ju/Uj8PFn8c+/sklPT1ZKVSCaXBwkNZvB+PXsDOzE+LlrOX3OPR1pSqni6dOnT661FqBS11qKQ5OLw6qHBvPaTR24MnEhM+e45yphpTyBNsl7juLsK00u5UDPhoGMr7qUIXsm8s2fm5wOR6lSFxAQQEJCgiYYD2CMISEhgYCAgCK9T69zKQ/8ggi+ZRbB0/sT8N1D7Ir6nKa1Q5yOSqlSExkZSWxs7PkhVFT5FhAQQGRkZJHeo8mlnPCJ7MzJKyfS7+eJTPnwX9z5yL+p4qfDw6iKydfXl0aNGjkdhipFHtcsJiIDRWS7iOwSkQm5zG8gIktFZK2IbBCRwU7EWRwhfcaTEHEVd56ZwSufLy/4DUopVU55VHIREW9gKjAIaA2MFJHWOYo9A8w3xnQCRgBv4ylECLv5Axa1eZ0Z61L5fE2s0xEppVSxeFRyAboCu4wxe4wx54C5wNAcZQyQ3WERCuQ9iE95FBTGjTeOomujGsz64ht2xendK5VSnsfTkks94KDL61h7mqvngFtEJBb4DnggtwWJyDgRiRGRmPLWqejj7cW7V2bymfcEfvjwOc6cy3Q6JKWUKhJPSy65DU+a81zGkcBHxphIYDAwS0Qu2U5jzDRjTLQxJjo8PLwUQi2ZGi16ciKyP3878yEzPp2tp2wqpTyKpyWXWMD1LjaRXNrsdScwH8AY8wcQANQsk+jcSYSat0wnpUo9hu19ls9/XeN0REopVWiellxWAc1EpJGI+GF12C/KUeYA0BdARFphJZfy1e5VWAGhhI6eQzWvM0T+cB9r9nnmZiilKh+PSi7GmAzgfmAJsBXrrLDNIjJJRIbYxR4FxorIemAOMMZ4cJuSV912ZAx+nZX+Pbh39nriU846HZJSShVIh9zHoSH3i2jL4ZPc8M5vXFYvkJljr8DH26N+FyilKiAdcr8CaB0RwrTeabxyZAzTP1/sdDhKKZUvTS4epHfXroT4wjUbH2bxqq1Oh6OUUnnS5OJJQuvhf/OnRHolEPzNOLYfTnQ6IqWUypUmFw/j26gHZ/q9SC/ZwJoPHyLx1DmnQ1JKqUtocvFAIb3GcqzVbfifPcH9s2NIz8xyOiSllLqIJhcPVeumNzBD3ua3PYlM+nqL0+EopdRFNLl4Ki9vboyuz1NdvBi4ehwLf17tdERKKXWeJhcP97ceEUT77Kbpj2NZsf1gwW9QSqkyoMnFw3lHdCDr+mm099pDypy/sf94itMhKaWUJpeKoEr7oST1fJb+rOD3aQ+SkpbudEhKqUpOk0sFUaPfIxxtOoIWaRt4dM5KMvQMMqWUgzS5VBQi1Bk5he3XzOb77UlM+maL3gNGKeUYTS4VibcvI3u24IHLw+kW8zALl/zkdERKqUpKk0sF9HCvWlzht5Muf9zDj6s3Ox2OUqoS0uRSAXnViML/1rnUlURqLrqV9bsPOR2SUqqS0eRSQflHdefsde/TVvaSMusWDhxLdjokpVQlosmlAgvpeB0Jff4fTdnPkzOXkHRaB7lUSpUNTS4VXK0+d3Nw5FJWJwUzbtZq0tIznQ5JKVUJaHKpBLq0aMgrw9ow6ODrfPHeP8nM0lOUlVKlS5NLJTG0fQRX1znD8Pi3mDtzql4Do5QqVZpcKgtvHxqOm8fh4LYM2/ccCz6b63RESqkKTJNLZeIXSL17viIpIIJBGx/m2yXfOR2RUqqC0uRSyUhQGGF3f0eKXzif/7KabzcccTokpVQFpMmlEvKpHkn1R2NIjuzLw/PW8cd2vchSKeVeHpdcRGSgiGwXkV0iMiGPMn8VkS0isllEPi3rGD1BlQB/PhjdhdtDVhHx6dVs3rbN6ZCUUhWIRyUXEfEGpgKDgNbASBFpnaNMM+ApoKcxpg3wUJkH6iFCA30Ze90AakoyAXNvYMeevU6HpJSqIDwquQBdgV3GmD3GmHPAXGBojjJjganGmEQAY8yxMo7Ro9Rs0Z1TN84mkngyP76BPQcPOx2SUqoC8LTkUg9wvVF8rD3NVXOguYj8JiIrRGRgbgsSkXEiEiMiMfHx8aUUrmeo1a4vJ/4ynWbsJ3XGdRw8luh0SEopD+dpyUVymZbzakAfoBnQBxgJTBeRape8yZhpxphoY0x0eHi42wP1NHW7DOXogHdYbjowcsYajiSfcTokpZQH87TkEgvUd3kdCeRsx4kFvjLGpBtj9gLbsZKNKkDk5cPpeecrJJ3J4Kn3FhKfkOB0SEopD+VpyWUV0ExEGomIHzACWJSjzJfAVQAiUhOrmWxPmUbpwTrUr8bMW1rz6qm/c+TtISQmahOZUqroPCq5GGMygPuBJcBWYL4xZrOITBKRIXaxJUCCiGwBlgKPG2P0J3gRXNasPsd7PUebjM3snzqEE0lJToeklPIwogMYQnR0tImJiXE6jHJn25LpNPv9MTb4tqfh/V9To1qo0yEppcoREVltjInObZ5H1VxU2Wp5zd/YdfnLdEjfwHfvPElC6lmnQ1JKeQhNLipfLa4Zx5arP+ClU4MY9f6fHNcEo5QqBE0uqkBtr7yRd8f0JPFEHCvfvIWE43pdqlIqf5pcVKH0bFqTDwf40e/cj5x4+xqOx+mV/EqpvGlyUYXW5oqh7O47nfqZBzn53kCOHNrvdEhKqXJKk4sqkla9b2D/wJnUyTzKuekDObBvp9MhKaXKIU0uqsha9PgLR4bOIc34cc+sGLYeOel0SEqpckaTiyqWJp374n3PryR412LUe7+xedM6p0NSSpUjmlxUsTWtHcKCu3vwsM9C6i8YxMbfFzsdklKqnNDkokqkfo1ABo15imTv6jRbcitrf5jjdEhKqXJAk4sqsfDIpgTf8yMHfBvR7td7+fPzt5wOSSnlME0uyi2qhdel3vjv2ValI+3Wv8D73/2BjlunVOWlyUW5TVBIdVo8/B3vNX6Tf/1ygqe/2ERGRqbTYSmlHKDJRbmVr38VHrptBPf2aQKrP2T16zdy5vRpp8NSSpUxTS7K7USEJwa25MbWIXQ7tZQ9r/cj6Xic02EppcqQJhdVaqJvfo61XV+l2bntnHz7KmJ3b3E6JKVUGck3uYiINpirEuk0eCy7B35CaFYSAbMGsnrbXqdDUkqVgYJqLlImUagKrVWPQaTe/B2f+I9gxKytfLY61umQlFKlrKDkoueSKreo16wjtz/4L7pE1WDuwnn8+NHzZGXp10upisptfS7ahKYKEhroy8w7uvJknRj67nuV5W/eRlpamtNhKaVKgTs79LUJTRXI19uLy+6fxdqGt9M7eRE7JvfnmN54TKkKx53J5Xwbh4hcLiLBbly2qkDE24dOt7/Bhq4v0fLcFs6+cxXrt253OiyllBuV1qnIZ4HbROQhEXmslNahPFz7wXdz5PrPWOPTkZs+2cmnfx5wOiSllJsUKbmISJCIdBGRoPzKGWNWA5lACFC1BPGpCq5hxz70eeQTujepxdtf/sRX057jbHqG02EppUqoqDWXT4GewGQRGVlA2XjgB8CtDeoiMlBEtovILhGZkE+5YSJiRCTanetX7hca6MuHY7rwalQMQw+/zq+TR3LsRJLTYSmlSqCoyWWPMeYNY8zdQNfcCmTXboD/GmN+N8ZMK3GUF5btDUwFBgGtgZEi0jqXcsHAeOBPd61blS5vL6H72LfY2eIe+qV9T/x/+rJh8yanw1JKFVNRk8smEXldRMbkU2YOha/dFFVXYJcxZo8x5hwwFxiaS7nngZcBPc/Vk3h50Wzki8ReM40oc4h68wfx+XeLdeh+pTxQkZKLMeYD4BUgA/ASkU9yKba7oNpNCdQDDrq8jrWnnScinYD6xphv8luQiIwTkRgRiYmPj3dzmKokInsMJ+tvP7Gvakf+/stpxssVbvgAAB0/SURBVM1aTfLpdKfDUkoVQZHPFjPGHDbGfGKMedAYc0suRQpTuymu3K6lcT0F2gt4HXi0oAUZY6YZY6KNMdHh4eFuDFG5Q3Bkazo/tojH/q8zf247wNevjWXjXr0eRilP4faLKAtZuymuWKC+y+tILj5hIBhoCywTkX1Ad2CRdup7JhHhzl6N+GJwBiMzvsLvowF88cMv2kymlAeQ/P5RRSQLeAE4DnxnjNlVVoHlEY8PsAPoCxwCVgGjjDGb8yi/DHjMGBOT33Kjo6NNTEy+RZTDUjZ/j3x2B1mZmcyNmMDwW+8lNNDX6bCUqtREZLUxJtcf7/nWXIwxXsaYicAnQA8RmSQiT4tITxEp8+FejDEZwP3AEmArMN8Ys9mOa0hZx6PKTnCbAQTev5y0kEaMOzKRD157kj/3JDgdllIqD/nWXHJ9g4gvcBXQAys5rQOWGGM89l62WnPxIBlnOfbVRMbv6sTKpGDu7dOUB/s1w9db73unVFnLr+ZS5OSSy8I7AgOwrsSPBb4xxnhUz6smF89z6mwGz321ie4bn2FH9SsZNfpeGoblO3CEUsrN8ksuPiVduDFmHVbtBRGJAK4VkXpACvC9MWZ9SdehVE5B/j68cm1Dkg6f4MakF5j51lqq/t8L3BDdCAdabJVSObi1LcGusXwLxAHXAf8VkbvcuQ6lzqtSnWr3LyW1452Mlm9p8vUwnvnwO46nnnU6MqUqPbckFxHpICL/EJHVWJ3tkcBjQIQx5j13rEOpXPn4U/W618i8aSatfOO4Y/9jXPPaMr7dcMTpyJSq1ErcLCYiXwJNgJnASGPMjhJHpVQRebe5Du+Ijvgd2Ee9X30Y/+kqlq6vztM3dKNGkJ/T4SlV6bijQ785UMcY84t7Qip72qFfsWRkZrH24yeos+8rJvo8yIgbhnFNmzpOh6VUhVPs61wKwxizw5MTi6p4fLy96NL3JmqH+DM981k2f/p3Hp6zisRT55wOTalKQy8OUBVTg2743fcb0vZGHvFdyM1b7+O21xby/eajTkemVKWgyUVVXAGheA2bDje8T8fAeKICzzJu1moemrtWazFKlbISd+grVe61/ys+Lf/CZK8qNFm2i4Rl7zB8Z08eu6EXA7QvRqlSoTUXVTn4BeHn48VD0VX4p/9sFmQ9zBez3+ZBrcUoVSo0uajKpVp9vO7+leA6TXjH7036bn6aG1/7VvtilHIzTS6q8glvgdfffoCrnuFan5XMNM9wz6yVPDh3rd7xUik30T4XVTl5+8CVjyPNB1D3xH4eONKSKT/tZOPuWJ4f3oOeTWs6HaFSHk1rLqpyq9sBnzZDeKhfc5ZdfYD5GQ8wc8YUJn29hbT0TKejU8pjaXJRyhbZujs1wiOY5vc6nVc+xK1vfs3mw8lOh6WUR9LkolS2iE543fUzXP0Mg33XMD31PqZOfZ3X/reDsxlai1GqKDS5KOXK2xd6P47XPb8TVK81HRrX5a0fd/KXt5azen+i09Ep5TFKPHBlRaADV6pcGQMiLN1+jG3zn+PYWS/oMo7HBrYmyF/PhVGqVO9EqVSFZd/R8qrm4VzROAGfnf9lw+rljN98L8OHDqF/69p610ul8qDNYkoVRASfUXNh2AxaBaXy/rknODLnfu75YBl74lOdjk6pckmTi1KFIQJtb8T3wdWYLmO52e9n4g7u4po3fuGl/27j9LkMpyNUqlzRPhe0z0UVQ0ocxwjlpcXbCVv/DusCezJq8NUM6RCBl5c2lanKoVRvFqZUpRRcm1rBAUweVIcnAr9hdvpD7P/sGW6c8hO/7zrudHRKOc7jkouIDBSR7SKyS0Qm5DL/ERHZIiIbRORHEWnoRJyqkgiujc/41fi0vZ4HfT5nSuI9zJwxhTEz/mTb0ZNOR6eUYzwquYiINzAVGAS0BkaKSOscxdYC0caY9sBC4OWyjVJVOsG1kWHT4bZF1K1ZnbcC32fXgVgGvfkrjy1Yz4GE005HqFSZ87RTkbsCu4wxewBEZC4wFNiSXcAYs9Sl/ArgljKNUFVeja/E6+7f8I/fyjchLZjy4w58Vr7D0LV9GNC5Bfdf3ZT6NQKdjlKpMuFRNRegHnDQ5XWsPS0vdwKLc5shIuNEJEZEYuLj490YoqrUvH2gTjuqBfrxTMfTPOk9m+WBj+G7/mP6vvoTEz7bwMETWpNRFZ+nJZfcTsPJ9XQ3EbkFiAZeyW2+MWaaMSbaGBMdHh7uxhCVstXvitz1M0ERrXjB+32WhU5i39qfuOrVZTy+YD0741KcjlCpUuNpySUWqO/yOhI4nLOQiPQD/g4MMcacLaPYlLpU3Q5w+2K48QMifFL4pOZH3NK1Hl9vOEz/13/hbzNXsWrfCfSSAFXReNR1LiLiA+wA+gKHgFXAKGPMZpcynbA68gcaY3YWZrl6nYsqE+dOQeJ+qN2aE0lJ7PzsnzwR24v9ZwLo3KAa43o3oX/r2njrdTLKQ1SY61yMMRnA/cASYCsw3xizWUQmicgQu9grQFVggYisE5FFDoWr1MX8gqC2dXJjjbg/6HbwQ5b5P8zCDqtJTjnF3Z+s5spXlvLuz7tJPHXO4WCVKhmPqrmUFq25KEcc3QT/exZ2/4SpHsXa5g/x8oEWrNibiJ+PF0M6RHBbj4a0j6zmdKRK5Sq/mosmFzS5KIft+gG+fxaqVIcx37I9LpVZK/bx+ZpDnD6XSYf61bi5awMGt69LVR3qX5UjmlwKoMlFOS4rE06fgKrhcPIwLH6S1B6PsfBgCLNW7Gd3/CkC/bwZ3K4uN10WSddGNXS4f+U4vZ+LUuWdl7eVWADiNsOen6m69WvGtLuJ0bc/zZqUaiyIOcg3G46wcHUsDcMCuemySG7oHElEtSrOxq5ULrTmgtZcVDl0+gT89gb8+R5kZUD0HTDoZU6nZ7J441EWrD7Iij0nAOgSVZ1rO0QwqG1dwoP9HQ5cVSbaLFYATS6q3Dp5BH55GTLTYegUa9q50+AXyIGE0yxaf4hF6w+zIy4VL4GeTWtybfsIrmlTh9BAX2djVxWeJpcCaHJR5Z4x1g3LDq+Dj4dCz/HQ9S7wrwrA9qMpfL3+MF9vOMz+hNP4eguXN6lJ/9a16d+6NrVDAhzeAFURaXIpgCYX5TGO74Tvn4Ed/4XAMLh8PHT52/kkY4xh46FkvtlwhCWbj7LfHpG5Q/1qDGhdmwGta9O0VlU9GUC5hSaXAmhyUR4nNgaWvQi7/gfVGsIDa6xBM10YY9h5LJX/bYnj+81HWR+bDEBUWCB9WtTiyubhdGtcg0A/Pa9HFY8mlwJoclEe6+AqOLEbOoywms7WzoI2N5yvybg6mpzG/7bG8cOWOFbsSeBsRhZ+3l50aVSdK5uH07t5OC1qB2utRhWaJpcCaHJRFcKBP2HGAKu5rOeDVnOZX1CuRdPSM1m59wS/7Ijnl53x7IhLBaB2iD89m9SkW+MadG8cRoMagZpsVJ40uRRAk4uqMA6utJrLdv8IgTXh8vuh293gm/+1MEeSz/DrjuP8vDOeP/ckcDzVGtusTkgA3RvXoFvjMLo3DiMqTJONukCTSwE0uagKJzvJxG+H8WvBxw8yMy7pl8mNMYbd8ams2HOCFXsS+HPvCeJTrDtX1Ar2p3OD6nRqUI1ODarTrl4oVfy8S3trVDmlyaUAmlxUhXUm0RqzLOMsvHM5NO1v1WZCIwu9CGMMe4+fYsWeE6zcm8Dag0nnz0Lz8RJa1Q2xk001OtavrrWbSkSTSwE0uagK7/QJWPJ32DAPxAs6DIeeD0PNpsVaXELqWdYdTGLtgSTWHkxk3YEkTp3LBCA4wIfWdUNoWy+UNhEhtIkIpUl4ED7eHnWHD1UImlwKoMlFVRpJB+D3/8Caj63azF2/QN32JV5sZpZh57EU1h5IYvPhZDYdOsm2oydJS88CwN/Hi5Z1Q2hrJ5sWdYJpXrsqwQE6ioAn0+RSAE0uqtJJPQYbF0D3e60r/9d8bDWVNb7Keu0GGZlZ7Dl+6nyy2Xw4mc2HT5KSlnG+TL1qVWheuyrN6wTTonYwzWsH07RWVQJ8tR/HE2hyKYAmF1WpZWXB1C6QsAtqt4Ue90PbG62TANzMGENs4hm2H01he1wKO+NS2B6Xyu5jqZzLtGo5XgINagTSOLwqjWoGEVUziMY1g2hUM4g6IQF46W2gyw1NLgXQ5KIqvYyzVk3m9ykQvxWC68L/vQEtBpbN6jOz2Jdwmh1xKWw/msLOYynsPX6avcdTzzetAQT4ehEVZiWaRjWDaBgWSGT1QCKrVyGiWhV8tV+nTGlyKYAmF6VsxljXyPz+H+j/vNUfk7jfmle9YZmHk5VliEtJY2/8KfYmnLL+HrceB06cJiPrwvHLS6zrcrKTTWQN+2/1KkRWC6RWiL82t7mZJpcCaHJRKh9f3G2dZdZqiNVkFhnttn6ZkkjPzOJochoHE08Tm3jGepzIfn6aIyfTyHl4qx7oS+2QAOqEBlAnJCDX59UDffVU6kLSO1EqpYrv6mehai2I+Qi2fAkRnaDnQ9DmOkfD8vX2on6NQOrXCMx1/rmMC8nnUNIZ4pLTOHoyjbiT1t9Nh06ScOrsJQnIz8eL8Kr+hFX1IyzIj5pV/Qmr6k/NqtnPL/ytEeinp1jnQWsuaM1FqUI5mwLr58LKadbFmAP/bTWjpcZBcB2noyuW9MwsjqWc5WiynXTsv8dTz3E89SwJp85yPOUcCafOkp556bFSBKoH+lE90JdqgX5Uq2L/DfS1n/sSen66L9UD/QgN9CXY36dC1I60WawAmlyUKgJjIP0M+AXC7qUwexi0HmrdvKx+13LRZOZuxhhOpmVYCSc78aSePZ+EEk+fI+l0Okmn00k+k07S6XPnLyrNjbeXEFrFl+AAH4IDfKjq70NwgJV0qp6f5kvVAB9CXOZbf61HkL+P4ycwaLOYUsp9RKzEAhDW1BoYc80s2PQZ1O0I3e6CdjeBd8W5QFLESgahVXxpEl6495zLyCL5TDrJZy4knsTT5+zkk07SmXOkpGWQkpZBaloGsYlnSElLJ/WsNS0zq+Af/r7eQhVfbwL9fAj09ybQz5tAXx+q+FnPs/8G+vnY5S68zp5fI8iP9pHVSvgJXcrjai4iMhB4E/AGphtjXswx3x/4GLgMSACGG2P25bdMrbkoVUJnU61O/z/fg7RkeGijdZ3M2dRc7y2j8meMIS09i5S0dFLOWsknJS2D1LPpnLSTUerZDE6fy+TMOevv6fRMzpzL5HT263MXvz6bkZXrutrWC+GbB64oVpwVpuYiIt7AVKA/EAusEpFFxpgtLsXuBBKNMU1FZATwEjC87KNVqhLxrwpd7oToOyD54IVRmKd2g/Dm1vTmgwo1KrOyakpV7JpFLTctMzPLcCbdSjZn7ORz+lwmvt6l04zpaXu6K7DLGLMHQETmAkMB1+QyFHjOfr4QmCIiYjytiqaUJxKBag2s55nnoPNtsGYmzLvFujCz821w2e0QUtfZOCshby+hqr/Vf1MWPO0cunrAQZfXsfa0XMsYYzKAZCAs54JEZJyIxIhITHx8fCmFq1Ql5hcIfZ6EBzfAiDnW0DI/vwzHNlvzz52CrLw7vZVn87SaS271t5w1ksKUwRgzDZgGVp9LyUNTSuXK2wdaDrYeifshtL41fdmL1nUznUdDx1EQEuFsnMqtPK3mEgvUd3kdCRzOq4yI+AChwIkyiU4plb/qDcHLPuw07AnVGsJPz8PrbWD2TbDtO2fjU27jacllFdBMRBqJiB8wAliUo8wiYLT9fBjwk/a3KFUOtRgIY76BB9ZAr0fg6CbY/MWF+Qm7nYtNlZhHNYsZYzJE5H5gCdapyDOMMZtFZBIQY4xZBHwAzBKRXVg1lhHORayUKlBYE+j7LFz1tHUaM1iJ5t2eUO8y6HSLdQuAgFBn41RF4nHXuZQGvc5FqXLm9AlrqJm1s+DYFvCpYo0C0O85PdOsHKkw17kopSqJwBrQ417ofg8cXgNrP4Htiy9ckHlgBQSGQc1mzsap8qTJRSlVfolYTWP1LoPBr4KXfT+WxU/CkXUQ0Rk6jLCazYJqOhuruoindegrpSorL5cbfY2aBwNegMx0WPwETG4BS//tXGzqEppclFKeJ7gOXP4A3LMc7vkdut9rXaQJkHoMFo2H/b9DVu7jaanSp81iSinPVrsNDHj+wusjG2DjQmvYmdAG0PYGaDfMSj4V8HYA5ZXWXJRSFUuzfvD4Trh+mjVo5u//gXevsG5qBta9aFSp05qLUqri8QuCDsOtx6kEOPD7hbtlzr8NUo5A22FWrSZ7oE3lVlpzUUpVbEFh0OraC6+bDQBvf/jhH/BGO5je37rRmXIrTS5Kqcql61gY+yM8uB76/sNqJjux15p37hSsfB9SjjobYwWgV+ijV+grVellZVqnOu9YAp/+FRBo0B1aDbFqPdXqF7iIyii/K/S15qKUUtnX0DS/Bu79E/o8BWknYclT8EZbiN9hzc9Mdy5GD6Md+kop5apWS+vR50k4vgt2/3hhmJmvH4KjG6D1EGg11DobTeVKay5KKZWXmk2h210Xro+JjAafAPjpBZjaBaZ2hz/fczbGckprLkopVVjRt1uP5EOw7RvYsujCyQBZmdYZaI2vgqgrwMfP2Vgdph36aIe+UqoEjLFqNsd3wnu9If00+IdA037QYjA0H1Bh70WjHfpKKVVaspvMajaDJ/bAyLnWvWf2/Qqf/w1iV1nzTx6BpAPOxVnGtFlMKaXcxbcKtBhkPbKy4FAM1O1gzVs5DZa/BnXaWTWaZtdARCfwqpi/8TW5KKVUafDygvpdL7zufCtUqQ7bv4NfXoGfX4LqjeCBNVbZzAzwrjiH5IqzJUopVZ7VaAw9x1uPUwnWKc6n4i/UXN6/CvyDoVl/q1ZTq5VHj+JcMetjSilVngWFQfu/Qo/7rNdZmdYFnGdPwg/PwTs9rHHPVs90NMyS0JqLUko5zcsbrn7Gepw8DDv/Bzu/t0Z3BkjcD988bNVqmvS1Th4o57UaTS5KKVWehETAZaOtR7aUI5C0H/47wXodWh+aXAW9nyi3455pclFKqfKuQXd4YDUk7oPdS63+mi2L4OqJ1vzNX0DcZqtWExkN3r6OhguaXJRSynNUj7owSkD2SM4AsTGw4m3rLDT/EGjU22pCu2yMY6F6TIe+iNQQkf+JyE77b/VcynQUkT9EZLOIbBCR4U7EqpRSpS47sQBc8y/rAs6/fmzdXfPIBlg7+8L8P6bClq/g9IkyC89jhn8RkZeBE8aYF0VkAlDdGPNkjjLNAWOM2SkiEcBqoJUxJim/ZevwL0qpCsUY68yzgFDrNgGvNoczJwCxLupsfCW0vg7qdS7RavIb/sWTmsWGAn3s5zOBZcBFycUYs8Pl+WEROQaEA/kmF6WUqlBELoxn5u0Lj+2AQ6thz8+wZxn88TYEVCtxcsk3BA+quSQZY6q5vE40xlzSNOYyvytWEmpjjMnKZf44YBxAgwYNLtu/f38pRK2UUuXQ2VTIyoAq1Qoumw+PqbmIyA9AnVxm/b2Iy6kLzAJG55ZYAIwx04BpYDWLFTFUpZTyXP5VS30V5Sq5GGP65TVPROJEpK4x5oidPI7lUS4E+BZ4xhizopRCVUoplQ+POVsMWARkX1U0GvgqZwER8QO+AD42xiwow9iUUkq58KTk8iLQX0R2Av3t14hItIhMt8v8FegNjBGRdfajozPhKqVU5eUxHfqlSU9FVkqpotM7USqllCpTmlyUUkq5nSYXpZRSbqfJRSmllNtphz4gIvFASS7Rrwkcd1M4Tqoo2wG6LeVVRdmWirIdULJtaWiMCc9thiYXNxCRmLzOmPAkFWU7QLelvKoo21JRtgNKb1u0WUwppZTbaXJRSinldppc3GOa0wG4SUXZDtBtKa8qyrZUlO2AUtoW7XNRSinldlpzUUop5XaaXJRSSrmdJpcSEJGBIrJdRHaJyASn4ykJEdknIhvtkaQ9ahRPEZkhIsdEZJPLtBoi8j8R2Wn/zfOupeVJHtvynIgcchnpe7CTMRaGiNQXkaUislVENovIg/Z0j9sv+WyLJ+6XABFZKSLr7W35pz29kYj8ae+XefbtS0q2Lu1zKR4R8QZ2YA3/HwusAkYaY7Y4Glgxicg+INoY43EXholIbyAV6z4+be1pLwMnjDEv2om/ujHmSSfjLIw8tuU5INUY86qTsRWFfUO/usaYNSISDKwGrgPG4GH7JZ9t+Suet18ECDLGpIqIL7AceBB4BPjcGDNXRN4F1htj3inJurTmUnxdgV3GmD3GmHPAXGCowzFVSsaYX4ATOSYPBWbaz2diHQzKvTy2xeMYY44YY9bYz1OArUA9PHC/5LMtHsdYUu2XvvbDAFcDC+3pbtkvmlyKrx5w0OV1LB76hbMZ4HsRWS0i45wOxg1qG2OOgHVwAGo5HE9J3S8iG+xms3LflORKRKKATsCfePh+ybEt4IH7RUS8RWQd1q3i/wfsBpKMMRl2EbccyzS5FJ/kMs2T2xh7GmM6A4OA++zmGVU+vAM0AToCR4DJzoZTeCJSFfgMeMgYc9LpeEoil23xyP1ijMk0xnQEIrFaYFrlVqyk69HkUnyxQH2X15HAYYdiKTFjzGH77zHgC6wvnSeLs9vKs9vMjzkcT7EZY+LsA0IW8D4esm/sNv3PgNnGmM/tyR65X3LbFk/dL9mMMUnAMqA7UE1EfOxZbjmWaXIpvlVAM/ssCz9gBLDI4ZiKRUSC7I5KRCQIGABsyv9d5d4iYLT9fDTwlYOxlEj2wdh2PR6wb+yO4w+ArcaY11xmedx+yWtbPHS/hItINft5FaAfVh/SUmCYXcwt+0XPFisB+9TDNwBvYIYx5l8Oh1QsItIYq7YC4AN86knbIiJzgD5YQ4fHAf8AvgTmAw2AA8BNxphy31Gex7b0wWp6McA+4K7sfovySkR6Ab8CG4Ese/LTWH0VHrVf8tmWkXjefmmP1WHvjVW5mG+MmWQfA+YCNYC1wC3GmLMlWpcmF6WUUu6mzWJKKaXcTpOLUkopt9PkopRSyu00uSillHI7TS5KKaXcTpOLUvkQketFxIhIS5dpUdmjFotIR9fRcEVkSPYI2SJynYi0LsY6l4lIdBHKTxKRfkUoP0ZEphQ1LqWKQpOLUvkbiTVy7Ig85ncEzicXY8wiY8yL9svrgCInl6Iyxkw0xvxQ2utRqig0uSiVB3ssqZ7AneSSXOyRGSYBw+37eQzPrhWIyOXAEOAVe14T1xqJiNS0b3OAiFQRkbn2AIjzgCou6xggIn+IyBoRWWDHlDOOj0RkmP18n4j80y6/0bXGlcc2/sVefk17Oe+Ide+SPSJypT0g41YR+ah4n6KqrDS5KJW364D/GmN2ACdEpLPrTPtWCxOBecaYjsaYeS7zfsca6uRxe97ufNZzD3DaGNMe+BdwGVgJCHgG6GcPKhqDdd+Nghy3y78DPJZXIRG5HpgADHa5j091rOHXHwa+Bl4H2gDtRKRjIdatFKDJRan8jMQaEgP778hSWk9v4BMAY8wGYIM9vTtWs9pv9hDpo4GGhVhe9iCRq4GoPMpcBTwJ/MUYk+gy/WtjDduxEYgzxmy0B2bcnM+ylLqET8FFlKp8RCQM6xd8WxExWGMxGRF5ogSLzeDCD7qAHPNyG4dJgP8ZY4qa1LLHhMok7//xPUBjoDlWjSjne7Ncnme/1uOFKjStuSiVu2FYtxpuaIyJMsbUB/YCvXKUSwGC81hGznn7sJu8uDACLcAvwM0AItIWaG9PXwH0FJGm9rxAEWlevM25xH7gBuBjEWnjpmUqdZ4mF6VyN5ILI0Vn+wwYlWPaUqB1dod+jnlzgcdFZK2INAFeBe4Rkd+xRj3O9g5QVUQ2AE8AKwGMMfFY95yfY89bAeTbQV8UxpjtWEltgR2fUm6joyIrpZRyO625KKWUcjtNLkoppdxOk4tSSim30+SilFLK7TS5KKWUcjtNLkoppdxOk4tSSim3+/+4yPMct3EWqgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def sigma_approx(H):\n", " sigma = (20 - H ) / (20 + H)\n", " return sigma\n", "\n", "Hs = np.linspace(0, 30e3, int(30e3))\n", "\n", "sigma_exact = sigma_density(Hs)\n", "sigma_approximation = sigma_approx(Hs/1e3)\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "plt.figure()\n", "plt.plot(Hs/1e3, sigma_exact, '-', label=\"Exact\")\n", "plt.plot(Hs/1e3, sigma_approximation, '--', label=\"Approximation\")\n", "plt.legend();\n", "plt.title('Exact and Approximate Sigma Corrections')\n", "plt.xlabel('Altitude in km')\n", "plt.ylabel('$\\\\sqrt{\\\\frac{\\\\rho}{\\\\rho_{sl}}}$');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Why use EAS?\n", "\n", "Equivalent airspeed is actually very useful (though pilots tend not to actually use it, annoyingly). We tend to define $q_\\infty$ as $\\tfrac{1}{2}\\rho_{SL}V_e^2$, which means that aerodynamic coefficients remain the same for a given $V_e$ and $\\alpha/\\beta$.\n", "\n", "It is useful to see that the Equivalent Air Speed (EAS) for a given flight condition, $V_e$, is the **speed which if flown at standard sea level density ($\\rho_{sl} \\sim 1.225 kg\\cdot m^{-3}$) would give the same aerodynamic loads, for the same aerodynamic configuration, constant $C_L$**.\n", "\n", "This is useful for simplicity of calculations, but is also useful for ease of flight - stall will always occur at the same angle, for the same EAS, regardless of flight altitude. Structural limits are always defined in EAS for the same reason, as the loads are constant at any altitude for the same EAS." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary of Corrections\n", "\n", "We can remember the order of corrections by the mnemonic \"ICE-T\", with the relative magnitudes of the velocities given by the shape that loosely looks like a square root - \n", "\n", "```{figure} ../Images/IceT.png\n", "---\n", "height: 300px\n", "name: IceT\n", "---\n", "Airspeed Corrections\n", "```\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problems:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below are a range of numerical and theoretical questions. Some of the numerical questions have random numbers that change each time the notes get updated - so you're not going crazy if you thought it was different last time.\n", "\n", "Try and go through the questions on your own before looking at the solutions - questions you'll face in tests/exams will be _similar_ but not exactly the same, so you need to be able to understand the *why* of the solution rather than just following an algorithm." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [ "remove-input", "remove-outputs", "remove-output", "thebe-init" ] }, "outputs": [ { "data": { "application/papermill.record/text/plain": "35000" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "altitude" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "195" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "vc" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "-10" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Vt" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "10" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Vtabs" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'headwind'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "wind" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'subtracted from'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "winddir" } }, "output_type": "display_data" } ], "source": [ "from myst_nb import glue\n", "from random import seed\n", "from random import randint\n", "\n", "# Make altitudes\n", "alts = np.arange(15000, 50000, 5000)\n", "vcs = np.arange(125, 300, 5)\n", "Vts = np.arange(-20, 20, 2)\n", "Vts[Vts == 0] = 3\n", "\n", "# seed(randseed)\n", "\n", "\n", "alt = alts[randint(0, len(alts)-1)]\n", "vc = vcs[randint(0, len(vcs)-1)]\n", "Vt = Vts[randint(0, len(Vts)-1)]\n", "if Vt > 0:\n", " wind = 'tailwind'\n", " winddir = 'added to'\n", "else:\n", " wind = 'headwind'\n", " winddir = 'subtracted from'\n", "\n", "\n", "\n", "# my_variable = \"here is some text!\"\n", "glue(\"altitude\", (alt), display=False);\n", "glue(\"vc\", (vc), display=False);\n", "glue(\"Vt\", (Vt), display=False);\n", "glue(\"Vtabs\", abs(Vt), display=False);\n", "glue(\"wind\", (wind), display=False);\n", "glue(\"winddir\", (winddir), display=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1.1 - Conversion between airspeeds\n", "\n", "An aircraft is flying at an altitude of {glue:text}`altitude`ft, with a calibrated airspeed of {glue:text}`vc`kn, with a {glue:text}`Vtabs`kn {glue:text}`wind`:\n", "\n", "a) How long will it take to cover 100 miles?\n", "\n", "b) How long will it take to cover 200km?\n", "\n", "c) If instead of $V_C$ = {glue:text}`vc`kn, you have $V_I$ = {glue:text}`vc`kn, with a position error of $\\Delta V_P=+2kn$, what do the above answers change to?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try and tackle the problem yourself before you see the solution below. The numbers in this problem will change each time the notes are updated. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "tags": [ "remove-input", "remove-output" ] }, "outputs": [ { "data": { "application/papermill.record/text/plain": "0.9672799999999999" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "press_correction" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "0.3105757644543481" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "dens_correction" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "192.48871999999997" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Ve" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "345.399339813558" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "V" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "335.399339813558" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Vg" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "172.5443268883529" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Vms" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "160934.0" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "dist1" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "200000.0" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "dist2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "932.7110482405745" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "time1" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "1159.1224330975115" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "time2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'15 minutes, 32 seconds'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "text_time1" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'19 minutes, 19 seconds'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "text_time2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "201" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Vc2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "194.30267999999998" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Ve2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "348.6542868382367" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "V2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "358.6542868382367" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Vg2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "184.5077051807132" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "Vms2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "872.234576016084" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "time1_2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "1083.965570999396" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "time2_2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'14 minutes, 32 seconds'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "text_time1_2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'18 minutes, 3 seconds'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "text_time2_2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'1 minutes, 0 seconds'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "dT1" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'1 minutes, 15 seconds'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "dT2" } }, "output_type": "display_data" }, { "data": { "application/papermill.record/text/plain": "'faster'" }, "metadata": { "scrapbook": { "mime_prefix": "application/papermill.record/", "name": "time_text" } }, "output_type": "display_data" } ], "source": [ "import math\n", "\n", "def hms(s):\n", " s = abs(s)\n", " hours = math.floor(s/3600)\n", " minutes = math.floor((s - hours*3600)/60)\n", " seconds = math.floor((s - hours*3600 - minutes*60))\n", " \n", " if hours > 0:\n", " textstr = f\"{hours:1.0f} hours, {minutes:1.0f} minutes, {seconds:1.0f} seconds\"\n", " elif minutes > 0:\n", " textstr = f\"{minutes:1.0f} minutes, {seconds:1.0f} seconds\"\n", " else:\n", " textstr = f\"{seconds:1.0f} seconds\"\n", " \n", " \n", " return textstr\n", "\n", "\n", "### Solution work\n", "press_correction = f_correction(VC=vc, h=alt)\n", "dens_correction = sigma_density(alt*0.3048)[0]\n", "\n", "\n", "\n", "Ve = press_correction * vc\n", "V = Ve * np.sqrt(1/dens_correction)\n", "\n", "\n", "glue(\"press_correction\", press_correction, display=False);\n", "glue(\"dens_correction\", dens_correction, display=False);\n", "glue(\"Ve\", (Ve), display=False);\n", "glue(\"V\", (V), display=False);\n", "\n", "Vtabs = abs(Vt)\n", "\n", "if wind == \"tailwind\":\n", " Vg = V + Vtabs\n", "elif wind == \"headwind\":\n", " Vg = V - Vtabs\n", " \n", "# V to SI\n", "Vms = Vg * 0.514444444\n", " \n", "# Miles to m\n", "dist = 100 * 1609.34\n", "\n", "time1 = dist / Vms\n", "\n", "# 200km \n", "time2 = 200 * 1e3 / Vms\n", "\n", "glue(\"Vg\", Vg, display=False);\n", "glue(\"Vms\", Vms, display=False);\n", "glue(\"dist1\", dist, display=False);\n", "glue(\"dist2\", 200 * 1e3, display=False);\n", "glue(\"time1\", time1, display=False);\n", "glue(\"time2\", time2, display=False);\n", "\n", "glue(\"text_time1\", hms(time1), display=False);\n", "glue(\"text_time2\", hms(time2), display=False);\n", "\n", "\n", "\n", "# Position error\n", "vc = vc + 2\n", "\n", "press_correction = f_correction(VC=vc, h=alt)\n", "dens_correction = sigma_density(alt*0.3048)[0]\n", "\n", "Ve = press_correction * vc\n", "V = Ve * np.sqrt(1/dens_correction)\n", "\n", "\n", "if wind == \"tailwind\":\n", " Vg = V + Vt\n", "elif wind == \"headwind\":\n", " Vg = V - Vt\n", "\n", "# V to SI\n", "Vms = Vg * 0.514444444 \n", " \n", "glue(\"Vc2\", vc, display=False);\n", "glue(\"Ve2\", Ve, display=False);\n", "glue(\"V2\", V, display=False);\n", "glue(\"Vg2\", Vg, display=False);\n", "glue(\"Vms2\", Vms, display=False);\n", " \n", "# Miles to m\n", "dist = 100 * 1609.34\n", "\n", "time1_2 = dist / Vms\n", "\n", "\n", "# 200km \n", "time2_2 = 200 * 1e3 / Vms\n", "\n", "glue(\"time1_2\", time1_2, display=False);\n", "glue(\"time2_2\", time2_2, display=False);\n", "\n", "glue(\"text_time1_2\", hms(time1_2), display=False);\n", "glue(\"text_time2_2\", hms(time2_2), display=False);\n", "\n", "# Time deltas\n", "dT1 = time1_2 - time1\n", "dT2 = time2_2 - time2\n", "if dT1 > 0:\n", " time_text = 'slower'\n", "else:\n", " time_text = 'faster'\n", " \n", "glue(\"dT1\", hms(dT1), display=False);\n", "glue(\"dT2\", hms(dT2), display=False);\n", "\n", "glue(\"time_text\", time_text, display=False);\n", "\n", " \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} Click to show the solution... \n", ":class: dropdown \n", "\n", "The first step is to take Vc to Ve via the pressure correction.\n", "\n", "For the altitude of {glue:text}`altitude`ft, with a calibrated airspeed of {glue:text}`vc`kn, the pressure correction can be interpolated from table using the function created earlier with syntax press_correction = f_correction({glue:text}`altitude`, {glue:text}`vc`).\n", "\n", "This gives a value of {glue:text}`press_correction:1.3f` for the pressure correction, and this determination pf the equivalent airspeed from $V_e = V_c\\cdot f$ so, $V_e=${glue:text}`Ve:1.3f`kn.\n", "\n", "For EAS to TAS, the density ratio is required - $\\sigma=\\frac{\\rho}{\\rho_{sl}}$, and this is found to be $\\sigma=${glue:text}`dens_correction:1.3f`kn. Since this is an altitude _above_ sea level, this has to be a number smaller than one - if it's not, then you've done something wrong.\n", "\n", "True airspeed is obtained from the relationship $V=V_e\\cdot\\sigma^{-\\frac{1}{2}}$ which yields $V$={glue:text}`V:1.3f`kn.\n", "\n", "Note that at no point have the units had to be converted into SI to make the corrections - that's the smart part of the density and pressure corrections; they are both non-dimensional corrections. If you have velocity in knots, you use them and get velocity in knots out.\n", "\n", "Since the wind is a {glue:text}`wind`, this is {glue:text}`winddir` the value above to yield the groundspeed in knots as {glue:text}`Vg:1.3f`kn. This gives the speed of the aircraft relative to the ground, and hence the speed we require for distance/time calculations.\n", "\n", "With the velocity in knots it may be easily converted to m/s by multiplying by 0.5144444 which yields $V$={glue:text}`Vms:1.3f`m/s. \n", "\n", "100 miles is 160,934m, and hence this is covered in {glue:text}`time1:1.1f`s or {glue:text}`text_time1`\n", "\n", "200 km is 200,000, and hence this is covered in {glue:text}`time2:1.1f`s or {glue:text}`text_time2`\n", "\n", "----\n", "\n", "For the last part of the question, this is solved by recalling that $V_c = V_I + \\Delta V_p$ and hence the speeds calculated are now:\n", "\n", "$V_c$ = {glue:text}`Vc2:1.3f`kn\n", "\n", "$V_e$ = {glue:text}`Ve2:1.3f`kn\n", "\n", "$V$ = {glue:text}`V2:1.3f`kn\n", "\n", "$V_g$ = {glue:text}`Vg2:1.3f`kn = {glue:text}`Vms2:1.3f`m/s \n", "\n", "The times taken are found to be {glue:text}`text_time1_2` and {glue:text}`text_time2_2`, for a respective {glue:text}`dT1` and {glue:text}`dT2` {glue:text}`time_text`.\n", "\n", "```\n", "\n", "\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Problem 1.2 - Theory\n", "\n", "a) When converting between different airspeeds, how can you 'sense check' your numbers?\n", "\n", "```{admonition} Click to show the solution... \n", ":class: dropdown \n", "\n", "Look at the square-root symbol for ICE-T {ref}`IceT`. This shows the respective magnitude of the velocities you should expect to find.\n", "\n", "Once you have the true airspeed, you should also compare it to the type of aircraft listed (if one is). The speed of sound is about 340m/s at sea-level - so if you've been given a light commuter aircraft, and you've found a true airspeed of 500m/s then something is wrong. (You probably got the knots to metre/s conversion wrong - [see here](https://www.google.com/search?q=how+do+i+convert+knots+to+metres+per+second&rlz=1C5CHFA_enUS892US892&oq=how+do+i+convert+knots+to+metres+per+second&aqs=chrome..69i57j6.5772j1j7&sourceid=chrome&ie=UTF-8)\n", "\n", "```\n", "\n", "b) What is the significance of \"equivalent airspeed\", and why is it preferable to use when defining limiting speeds for aircraft\n", "\n", "```{admonition} Click to show the solution... \n", ":class: dropdown \n", "\n", "For a given angle of attack, the dimensional lift on a given aircraft will be the same for constant $V_e$ at *any* altitude. This is not the same for any of the other airspeeds.\n", "\n", "Hence - the *stall speed* (which we will calculate soon) is *constant* with altitude when defined in EAS, as is the \"never exceed\" speed (provided $V_{ne}$ has been defined due to aerodymamic structural limitations.\n", "\n", "```\n", "\n", "c) In what cases could you find that your true airspeed is *less* than your calibrated airspeed? \n", "\n", "```{admonition} Click to show the solution... \n", ":class: dropdown \n", "\n", "Flight below sea-level, or if you've made a mistake.\n", "\n", "```\n", "\n", "d) For what conditions is EAS equivalent to TAS?\n", "\n", "\n", "```{admonition} Click to show the solution... \n", ":class: dropdown \n", "\n", "Flight at sea-level.\n", "\n", "```\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }