{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Version 1 (incorrect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "knot_year =\n",
      "\n",
      " Columns 1 through 11:\n",
      "\n",
      "   2001   2002   2003   2004   2005   2006   2007   2008   2009   2010   2011\n",
      "\n",
      " Columns 12 and 13:\n",
      "\n",
      "   2012   2013\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model.knots = [365 730 1095 1461 1826 2191 2556 2922 3287 3652 4017 4383 4748]';\n",
    "\n",
    "frac4=model.knots./1461;\n",
    "complete4= floor(frac4);\n",
    "mod4=frac4-complete4;\n",
    "daysleft4=model.knots - 1461*complete4;\n",
    "\n",
    "for i=1:length(model.knots)\n",
    "if (mod4(i) >= 1096/1461)    % Test for leap year\n",
    "     knot_year(i) = 2000 + 4*complete4(i) + 3 + (daysleft4(i)-1095)/366;  % Leap year\n",
    "else\n",
    "    knot_year(i)= 2000 + 4*complete4(i) + daysleft4(i)/365;   % Non-leap year\n",
    "end\n",
    "end\n",
    "\n",
    "knot_year"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Version 2 (correct)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      " Columns 1 through 11:\n",
      "\n",
      "   2000   2001   2002   2003   2004   2005   2006   2007   2008   2009   2010\n",
      "\n",
      " Columns 12 through 16:\n",
      "\n",
      "   2011   2012   2013   2100   1901\n",
      "\n"
     ]
    }
   ],
   "source": [
    "function decimal_year = mjd2000_to_decimal_year_v2(mjd2000)\n",
    "%\n",
    "% decimal_year = mjd2000_to_decimal_year_v2(mjd2000)\n",
    "%\n",
    "% the conversion is valid for MJD2000 from 36525 to -36159\n",
    "% (i.e., decimal years between 1901.0 and 2100.0)\n",
    "%\n",
    "% Convert MJD2000 to decimal years.\n",
    "%\n",
    "    frac4 = mjd2000 ./ 1461;\n",
    "    complete4 = floor(frac4);\n",
    "    mod4 = frac4 - complete4;\n",
    "    daysleft4 = mjd2000 - 1461*complete4;\n",
    "\n",
    "    mask_leap_years = mod4 <= 366/1461;\n",
    "    mask_regular_years = ~mask_leap_years;\n",
    "    \n",
    "    decimal_year = 2000 + 4*complete4;\n",
    "    decimal_year(mask_leap_years) += daysleft4(mask_leap_years)/366;\n",
    "    decimal_year(mask_regular_years) += 1 + (daysleft4(mask_regular_years) - 366)/365;\n",
    "\n",
    "end % mjd2000_to_decimal_year_2\n",
    "\n",
    "model.knots = [0 366 731 1096 1461 1827 2192 2557 2922 3288 3653 4018 4383 4749 36525 -36159];\n",
    "mjd2000_to_decimal_year_v2(model.knots)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference (VirES implementation compiled as a mex file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "knot_year =\n",
      "\n",
      " Columns 1 through 11:\n",
      "\n",
      "   2000   2001   2002   2003   2004   2005   2006   2007   2008   2009   2010\n",
      "\n",
      " Columns 12 through 16:\n",
      "\n",
      "   2011   2012   2013   2100   1901\n",
      "\n",
      "knot_day =\n",
      "\n",
      " Columns 1 through 10:\n",
      "\n",
      "       0     366     731    1096    1461    1827    2192    2557    2922    3288\n",
      "\n",
      " Columns 11 through 16:\n",
      "\n",
      "    3653    4018    4383    4749   36525  -36159\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model.knots = [0 366 731 1096 1461 1827 2192 2557 2922 3288 3653 4018 4383 4749 36525 -36159];\n",
    "knot_year = mjd2000_to_decimal_year(model.knots)\n",
    "knot_day = decimal_year_to_mjd2000(knot_year)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of the V2 with the reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N =  20000\n",
      "mjd2000_min = -36159\n",
      "mjd2000_max =  36525\n",
      "max_diff = 0\n"
     ]
    }
   ],
   "source": [
    "N = 20000\n",
    "mjd2000_min = decimal_year_to_mjd2000(1901)\n",
    "mjd2000_max = decimal_year_to_mjd2000(2100)\n",
    "\n",
    "mjd2000 = (mjd2000_max - mjd2000_min)*rand(1, N) + mjd2000_min;\n",
    "\n",
    "decimal_year = mjd2000_to_decimal_year_v2(mjd2000);\n",
    "decimal_year_ref = mjd2000_to_decimal_year(mjd2000);\n",
    "\n",
    "# the values should be equal and the deffierence close to zero\n",
    "max_diff = max(abs(decimal_year .- decimal_year_ref))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "4.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}