{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a5306a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "using OrdinaryDiffEq\n",
    "using PyPlot\n",
    "using CSV\n",
    "using DataFrames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a391c2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "\n",
    "The system of two dimensionless differential equations\n",
    "describing the radial distribution of the density, rho(r),\n",
    "and mass, m(r), inside a white dwarf star\n",
    "\n",
    "     m' =   rho r^2\n",
    "   rho' = - m rho /(gamma(\\rho) r^2)\n",
    "\n",
    "where gamma(rho) = rho^(2/3)/(3 sqrt(1 + rho^(2/3))\n",
    "\"\"\"\n",
    "function dwarf_eqs!(dudr, u, p, r)\n",
    "    m = u[1]\n",
    "    rho = u[2]\n",
    "    rho_c = p\n",
    "    if rho >= 0.0\n",
    "        w = rho^(2/3)\n",
    "        gamma = w/(3 * sqrt(1. + w))\n",
    "        dudr[1] = rho * r * r\n",
    "        if (r > 1.e-6)\n",
    "            dudr[2] = -m * rho/(gamma * r * r)\n",
    "        else\n",
    "            dudr[2] = -rho_c/3 * r * rho/gamma\n",
    "        end\n",
    "    else\n",
    "        dudr[1] = 0.0\n",
    "        dudr[2] = 0.0\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "526b70b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_c1 = 0.08      #  min density in the star center\n",
    "rho_c2 = 100000.0  #  max density in the star center\n",
    "np = 100\n",
    "\n",
    "factor = exp(log(rho_c2/rho_c1)/ (np - 1))\n",
    "rho_c = rho_c1\n",
    "radi = zeros(np)\n",
    "mass = zeros(np)\n",
    "rspan = (0.0, 10.0)\n",
    "for i = 1:np\n",
    "    u0 = [0.0, rho_c]\n",
    "    prob = ODEProblem(dwarf_eqs!, u0, rspan, rho_c)\n",
    "    condition(u, t, integrator) = u[2]  # stop integration when rho = 0.\n",
    "    affect!(integrator) = terminate!(integrator)\n",
    "    cb = ContinuousCallback(condition, affect!)\n",
    "    sol = solve(prob, Tsit5(), callback=cb)\n",
    "    radi[i] = sol.t[end]\n",
    "    mass[i] = sol(radi[i])[1]\n",
    "    global rho_c *= factor\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d80fa8d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download white dwarfs observational data from \\\"The VizieR database of astronomical catalogues\n",
    "url = \"http://vizier.u-strasbg.fr/viz-bin/asu-tsv?-source=J/A+A/420/507/tablea1&-out=Mass&-out=Rad\"\n",
    "catalog = download(url)\n",
    "df = CSV.read(catalog, DataFrame, comment=\"#\", skipto=4, \n",
    "              types=[Float64, Float64], delim='\\t', ignorerepeated=true, \n",
    "              silencewarnings=true)\n",
    "dropmissing!(df);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d605175",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(mass*0.71, radi*0.006, label=\"theory\")\n",
    "plot(df.Mass, df.Rad, marker=\"o\", linestyle=\"none\", \n",
    "     fillstyle=\"none\", markersize=4, label=\"observations\")\n",
    "grid(true)\n",
    "xlabel(L\"$M/M_{\\odot}$\")\n",
    "ylabel(L\"$R/R_{\\odot}$\")\n",
    "legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab1fae9a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}