{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Retrieval\n",
    "=========\n",
    "\n",
    "Retrievals are done using a series of interchangable \"blocks\" which carry\n",
    "out the many functions required for a DOAS-style retrieval. The hierarchy of\n",
    "classes is shown below.\n",
    "\n",
    "| :ref:`ret`\n",
    "|     :ref:`ret_scd`\n",
    "|         :ref:`ret_scd_cross_section`\n",
    "|             :ref:`ret_scd_instrument`\n",
    "|         :ref:`ret_scd_ring`\n",
    "|             :ref:`ret_scd_solar`\n",
    "|                 :ref:`ret_scd_instrument`\n",
    "|     :ref:`ret_amf`\n",
    "|         :ref:`ret_amf_atmosphere`\n",
    "|         :ref:`ret_amf_box`\n",
    "\n",
    "The following example shows how to simulate radiances using SASKTRAN, package the radiances in a Radiance object that skdoas.retrieval can use, and then perform a total column NO2 retrieval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "from scipy.integrate import trapz\n",
    "\n",
    "import sasktran as sk\n",
    "import skdoas.retrieval as skr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### SIMULATE RADIANCES ###\n",
    "mjd = 54000.6\n",
    "lats, lons = [52.1], [-106.7]  # Saskatoon\n",
    "\n",
    "# create lines of sight\n",
    "observer = sk.Geodetic()\n",
    "observer.from_lat_lon_alt(0., -100., 35786e3)  # approximate TEMPO location\n",
    "geometry = sk.NadirGeometry()\n",
    "geometry.from_lat_lon(mjd=mjd, observer=observer, lats=lats, lons=lons)\n",
    "geometry.reference_point = lats + lons + [0.0, mjd]\n",
    "\n",
    "# create atmosphere\n",
    "atmosphere = sk.Atmosphere()\n",
    "atmosphere.atmospheric_state = sk.MSIS90()\n",
    "atmosphere['air'] = sk.Species(sk.Rayleigh(), sk.MSIS90())\n",
    "atmosphere['ozone'] = sk.Species(sk.O3DBM(), sk.Labow())\n",
    "atmosphere['no2'] = sk.Species(sk.NO2Vandaele1998(), sk.Pratmo())\n",
    "\n",
    "wl = np.arange(405., 465., 0.2)\n",
    "engine = sk.EngineHR(geometry, atmosphere, wl)\n",
    "rad = engine.calculate_radiance('xarray')\n",
    "\n",
    "measurement_l1 = skr.Radiance_SASKTRAN(sk_result_xarray=rad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "setting up retrieval\n",
      "performing retrieval\n",
      "\n",
      "retrieval results:\n",
      "2.2799e+15 | 2.2922e+15 |  -0.54%\n"
     ]
    }
   ],
   "source": [
    "### PERFORM RETRIEVAL ON SIMULATED RADIANCES ###\n",
    "\n",
    "print('setting up retrieval')\n",
    "# create SCD fitting object\n",
    "instrument = skr.DeltaFunction()\n",
    "cross_sections = skr.CrossSection_Default(instrument=instrument, wavelength=np.arange(405., 465.1, 0.1))\n",
    "cross_sections.add_cross_section(sk.NO2Vandaele1998(), 220., 1e5)\n",
    "cross_sections.add_cross_section(sk.O3DBM(), 220., 1e5)\n",
    "cross_sections.add_cross_section(sk.HITRANChemical('H2O'), 220., 1e5)\n",
    "cross_sections.release_weights()\n",
    "solar_spectrum = skr.SolarSpectrum_Default(instrument=instrument, wavelength=np.arange(405., 465.1, 0.1))\n",
    "ring_spectrum = skr.RingSpectrum_Default(solar_spectrum=solar_spectrum, wavelength=np.arange(405., 465.1, 0.1))\n",
    "scd_object = skr.SlantColumn_OMNO2(cross_sections, ring_spectrum)\n",
    "\n",
    "# create AMF object\n",
    "amf_atmosphere = skr.Atmosphere_Default(atmosphere=atmosphere)\n",
    "amf_table = skr.BoxAirMassFactor_440nm50kmClearTable()\n",
    "amf_object = skr.AirMassFactor_OMNO2(amf_table, amf_atmosphere)\n",
    "\n",
    "retrieval = skr.Retrieval_OMNO2_TotalColumnOnly(scd_object, amf_object)\n",
    "\n",
    "print('performing retrieval')\n",
    "vcdr, vcdr_error = retrieval.retrieve(measurement_l1)\n",
    "\n",
    "alts = amf_table.altitude_edges()\n",
    "print('\\nretrieval results:')\n",
    "for i in range(measurement_l1.num_lines_of_sight()):\n",
    "    profile = sk.Pratmo().get_parameter('SKCLIMATOLOGY_NO2_CM3', lats[i], lons[i], alts, mjd)\n",
    "    vcd = trapz(profile, 1e2 * alts)\n",
    "    print(f'{vcdr[i]:10.4e} | {vcd:10.4e} | {1e2*(vcdr[i]/vcd - 1):6.2f}%')\n"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}