{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce2b568b",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2022-01-17T20:35:05.739Z"
    }
   },
   "outputs": [],
   "source": [
    "# https://gml.noaa.gov/grad/solcalc/solareqns.PDF\n",
    "\n",
    "# https://github.com/SatAgro/suntime\n",
    "\n",
    "# https://stackoverflow.com/questions/19615350/calculate-sunrise-and-sunset-times-for-a-given-gps-coordinate-within-postgresql"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aab13b5",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2022-01-17T20:35:05.742Z"
    }
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import datetime\n",
    "import pytz\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18a435c4",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2022-01-17T20:35:05.744Z"
    }
   },
   "outputs": [],
   "source": [
    "def rise_set_func(lat, lon, day, mo, yr):\n",
    "        \n",
    "# INPUTS\n",
    "#\n",
    "# lat: float32, latitude [degrees] of point of interest\n",
    "# lon: float32, longitude [degrees] of point of interest\n",
    "# day: integer, day of month\n",
    "# mo: integer, month of year\n",
    "# yr: integer, year of interest\n",
    "\n",
    "# OUTPUTS\n",
    "#\n",
    "# sunrise: float32, UTC hour at which sunrise occurs at point of interest\n",
    "# sunset: float32, UTC hour at which sunset occurs at point of interest\n",
    "    \n",
    "    \n",
    "    # set UTC second, minute, hours to 0 for data point\n",
    "    sc, mn, hr = 0, 0, 0\n",
    "            \n",
    "            \n",
    "    # (1) CALCULATE FRACTIONAL YEAR, g [radians]    \n",
    "            \n",
    "    # number of days in a year is 365, unless it is a leap year -> 366\n",
    "    nd = 365\n",
    "    if yr%4 == 0:\n",
    "        nd+=1\n",
    "        \n",
    "    # calculate day of year from datetime\n",
    "    d1 = math.floor( 275*mo / 9)\n",
    "    d2 = math.floor( (mo+9) / 12)\n",
    "    d3 = ( 1 + math.floor( (yr - 4*math.floor(yr/4) + 2) / 3) )\n",
    "        \n",
    "    d = d1 - (d2*d3) + day - 30\n",
    "        \n",
    "    # get g\n",
    "    g = (2*math.pi/nd) * (d - 1 + (hr - 12)/24 )\n",
    "        \n",
    "        \n",
    "    # (2) ESTIMATE EQUATION OF TIME, eqt [minutes]\n",
    "    eqt = 229.18*(0.000075 + 0.001868*math.cos(g) - 0.032077*math.sin(g) \\\n",
    "                    - 0.014615*math.cos(2*g) - 0.040849*math.sin(2*g))\n",
    "        \n",
    "        \n",
    "    # (3) ESTIMATE SOLAR DECLINATION ANGLE, sd [radians]\n",
    "    sd = 0.006918 - 0.399912*math.cos(g) + 0.070257*math.sin(g) \\\n",
    "            - 0.006758*math.cos(2*g) + 0.000907*math.sin(2*g) \\\n",
    "                - 0.002697*math.cos(3*g) + 0.00148*math.sin(3*g)\n",
    "        \n",
    "        \n",
    "    # (4) GET SOLAR HOUR ANGLE, ha, AT SUNRISE/SUNSET [degrees]\n",
    "        \n",
    "    # zenith at sunrise/sunset [degrees]\n",
    "    phi = 90.833\n",
    "    # conversions from degrees to radians\n",
    "    rphi = (math.pi/180)*phi\n",
    "    rlat = (math.pi/180)*lat\n",
    "        \n",
    "    ha = np.arccos(math.cos(rphi)/(math.cos(rlat)*math.cos(sd)) \\\n",
    "                    - math.tan(rlat)*math.tan(sd))\n",
    "    ha = (180/math.pi)*ha\n",
    "\n",
    "    # (6) GET UTC TIME OF SUNRISE AND SUNSET AT POINT OF INTEREST [minutes]\n",
    "    sunrise = (720 - 4*(lon + ha) - eqt)/60\n",
    "    sunset = (720 - 4*(lon - ha) - eqt)/60\n",
    "        \n",
    "    return sunrise, sunset\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95726822",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2022-01-17T20:35:05.746Z"
    }
   },
   "outputs": [],
   "source": [
    "def rise_set(lat, lon, time):\n",
    "        \n",
    "# INPUTS\n",
    "#\n",
    "# lat: float32 xr.DataArray, latitudes of interest\n",
    "# lon: float32 xr.DataArray, longitudes of interest\n",
    "# time: datetime64[ns] xr.DataArray, times of interest \n",
    "        \n",
    "    day = time.dt.day\n",
    "    mo = time.dt.month\n",
    "    yr = time.dt.year\n",
    "        \n",
    "    return xr.apply_ufunc(\n",
    "        rise_set_func,\n",
    "        lat, lon, day, mo, yr,\n",
    "        dask=\"parallelized\",\n",
    "        output_core_dims=[[],[]]\n",
    "    )\n",
    "       "
   ]
  }
 ],
 "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.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}