3_Lets Talk Astrophysics JL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IgTPJ8p5hW3-"
   },
   "source": [
    "# Let's Talk Astrophysics Workshop\n",
    "\n",
    "Оригинал: https://github.com/LTSQueens/Let-s-Talk-Astrophysics  \n",
    "Полная версия урока на русском питоне: Let's Talk Astrophysics Workshop: Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WokkH0j5_LWr"
   },
   "source": [
    "### 2.3 Observing Rotation Curves: CO Lines\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SD5SzLXv6i2b"
   },
   "source": [
    "Помимо GAIA, еще один телескоп собирает данные, которые можно использовать для построения кривых вращения галактики. Это 45-метровый радиотелескоп Нобеяма, наблюдающий за радиоволнами, известными как «миллиметровые волны». В частности, он может наблюдать то, что известно как молекулярные линии.\n",
    "\n",
    "Атомы могут образовывать молекулы, наиболее распространенными атомами во Вселенной являются C, H, O и He, поэтому логично предположить, что наиболее распространенные молекулы, которые могут быть сформированы, состоят из этих атомов. Для изучения Млечного Пути мы сосредоточимся на молекуле монооксида углерода (СО)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "joQ55Gc7onyH"
   },
   "source": [
    "![](images/img10.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qOuRR3Xv1ohs"
   },
   "source": [
    "Молекулярные связи на самом деле не такие прочные, излучение во Вселенной может легко их разорвать, поэтому молекулы, подобные СО, можно найти только там, где они защищены от этого излучения, в облаках газа и пыли.\n",
    "\n",
    "Ядра углерода и кислорода могут вращаться друг вокруг друга с целыми скоростями, точно так же, как электроны могут вращаться вокруг ядра атома в целочисленных энергетических состояниях, которые мы называем орбиталями или оболочками. Внимательно посмотрите на диаграммы ниже."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WT7hEQgaHSMQ"
   },
   "source": [
    "![](images/img11.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dI7aU2v_WJ4Y"
   },
   "source": [
    "![](images/img12.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oEPiJD1cSYGC"
   },
   "source": [
    "Это целочисленное изменение вращения излучает целочисленную единицу энергии, в данном случае свет, который мы называем фотоном с длиной волны, находящейся в диапазоне миллиметров. Вот почему радиотелескоп Нобеяма ищет длины волн в миллиметровом диапазоне!\n",
    "\n",
    "Этот фотон от линий молекулярного вращения, как и фотон от CO, можно использовать для определения физических свойств (таких как плотность и температура) межзвездного газа в галактике Млечный Путь.\n",
    "\n",
    "Эти миллиметровые фотоны межзвездных молекул являются хорошим индикатором плотных, толстых областей межзвездной среды. Эти плотные области обычно являются местами образования звезд."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OjJqGv8VTZDm"
   },
   "source": [
    "![](images/img13.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NXx5kHLaHffT"
   },
   "source": [
    "Это карта СО нашей галактики, цвета обозначают плотность молекулярных облаков. Более теплые цвета означают более высокую плотность, а более холодные тона указывают на менее плотные зоны.\n",
    "\n",
    "Вы заметили что-то странное на этой карте CO Млечного Пути? Похоже, что слева угарного газа больше, чем справа. Забавный факт: это, скорее всего, потому, что левая сторона нашего Млечного Пути взаимодействует с галактикой Андромеды!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "K1SKNlI6Lyn5"
   },
   "source": [
    "**Дополнительно:** Хотя молекулярный водород $H_2$ на сегодняшний день является доминирующей молекулой в молекулярных облаках, мы используем молекулярные линии $CO$.\n",
    " \n",
    "Одна из причин того, насколько сложно или легко обнаружить молекулу, заключается в том, что сила спектральных линий от молекул связана с тем, насколько «асимметрична» молекула.\n",
    "\n",
    "Какую молекулу труднее обнаружить, $H_2$ или $CO$? Почему?\n",
    "$H_2$ труднее обнаружить, поскольку это полностью симметричная молекула, поэтому ее спектральные линии слабее по сравнению с $CO$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "W436U5ra_ZXL"
   },
   "source": [
    "## Section 3: Rotation Curves and Dark Matter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bRvacamT_fDm"
   },
   "source": [
    "### 3.1 Plotting Rotation Curves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TrIxEYJzxoXH"
   },
   "source": [
    "1. На этом семинаре мы рассмотрим свидетельства существования темной материи, рассмотрев кривые вращения галактик, в частности, кривую вращения галактики Млечный Путь. Для этого мы рассмотрим фактические данные с радиотелескопа. Итак, давайте продолжим и загрузим данные кривой вращения."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DelimitedFiles, Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# r (kpc)   δr (kpc)  V (km s−1) δV (km s−1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dat = readdlm(\"data/COdata.txt\", skipstart = 3)\n",
    "(Rs, dRs, Vs, dVs) = [dat[:,i] for i in 1:4 ];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kIslegu5yPJY"
   },
   "source": [
    "2. Все физические величины имеют единицы измерения, как и данные кривой вращения.\n",
    "Радиус измеряется в килопарсеках, а скорость измеряется в $\\rm км\\,с^{-1}$. Давайте посмотрим на массив скоростей - как видите, он имеет 68 элементов в диапазоне от 100 км/с до 250 км/с.\n",
    "В файле данных вы также увидите, что измерения скорости и радиуса имеют связанные с ними ошибки.\n",
    "Ошибки являются неотъемлемой частью проведения экспериментов и говорят нам о том, насколько точны наши измерения.\n",
    "Как правило, ошибка зависит от характера используемого оборудования, но иногда может также заключаться в нашем непонимании физических процессов."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vI3Yw56hy2Ql"
   },
   "source": [
    "Мы можем попытаться построить кривую вращения и проанализировать ее. Мы увидим, что скорость очень быстро возрастает, затем уменьшается и выравнивается. Давайте попробуем построить кривую вращения Млечного Пути и посмотрим, как она выглядит."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot( Rs, Vs, yerr = dVs, lab = \"CO data\", line = (3, :gray), m = (3,:gray, 0.5),\n",
    "                xlab = \"R, kpc\", ylab = \"V, km/s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "savefig(\"MilkyWayCOdata.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "G5Pf2Pjg1QRm"
   },
   "source": [
    "Как видите, кривая вращения выходит на ~1000 кпк, но понять это может быть сложно. Кроме того, размер галактики Млечный Путь составляет всего 200 кпк, а звезды простираются примерно на 60 кпк. (по википедии радиус 16 кпк) Итак, давайте попробуем приблизить внутренние части кривой вращения и посмотреть, как она выглядит."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot( Rs, Vs, yerr = dVs, lab = \"CO data\", line = (3, :gray), m = (3,:gray, 0.5),\n",
    "                xlab = \"R, kpc\", ylab = \"V, km/s\", xlim = (0, 35), ylim = (0,250) )\n",
    "vline!( [8.2], line = (2, :dash, :yellow), lab = \"Position of the Sun\" )\n",
    "plot!(Rs, 1 ./ Rs.^2, line = (2, :dash, :red), lab = \"Keplerian\" )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "savefig(\"MilkyWayCO.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ihWstPhR2LoD"
   },
   "source": [
    "Мы видим, что кривая вращения Млечного Пути достигает очень высоких скоростей в центре, а затем падает и выходит на стационар при значении $220\\,\\text{км}\\, \\text{s}^{-1}$. Горизонтальная часть кривой вращения говорит нам о том, что существует какая-то невидимая материя.\n",
    "Красная пунктирная линия показывает наше ожидание кривой вращения, если бы она следовала закону Кеплера.\n",
    "\n",
    "Золотая пунктирная линия показывает нам расстояние от Солнца до центра Галактики. Как видите, мы находимся в плоской части кривой вращения, живя внутри диска Млечного Пути.\n",
    "Далее мы увидим доказательства этого, разбивая кривую вращения Млечного Пути на компоненты, а именно на Балдж, Диск и гало Темной Материи."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8tZWheKJ_pq-"
   },
   "source": [
    "### 3.2 Темная материя"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_0SjoD2L66Vx"
   },
   "source": [
    "Когда люди начали строить кривые вращения, стало ясно, что масса в спиральных галактиках не соответствует кеплеровским кривым вращения, ожидаемым от видимой массы, такой как звезды и газ.\n",
    "\n",
    "Как видно из данных наблюдений, кривые вращения плоские, насколько можно измерить радиус. Это не то, что ожидается, учитывая только видимую массу. Помните наше уравнение Кеплера, связывающее скорость, расстояние и массу?\n",
    "\n",
    "$v_{circ}(R)=\\sqrt{\\frac{GM(R)}{R}}$\n",
    "\n",
    "Вот почему кривые вращения интерпретируются как убедительное свидетельство «отсутствующей» массы в галактиках: если вы добавляете больше массы по мере увеличения расстояния, то вы можете поддерживать постоянную скорость! Мы называем эту массу, которую не видим, «темной материей».\n",
    "\n",
    "Чтобы доказать это, мы сначала разделим нашу галактику Млечный Путь на ее компоненты: балдж, диск и гало темной материи.\n",
    "\n",
    "Image Credit: ESA (European Space Agency)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vyP9zhdF6-zC"
   },
   "source": [
    "![](images/img14.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e6p_Jqk97VzK"
   },
   "source": [
    "На картинке выше мы можем видеть круговое распределение звезд в центре, называемое галактическим балджем. Именно здесь живет сверхмассивная черная дыра Стрелец A*. Тонкое «дискообразное» распределение звезд, выходящих из центра, называется «Галактический диск». И галактическая балдж, и диск заключены в круглое гало, заполненное темной материей.\n",
    "\n",
    "Прежде чем мы сосредоточимся на темной материи, давайте посмотрим, сможем ли мы объяснить кривую вращения Млечного пути только светящейся материей. Для начала создадим массив радиуса в килопарсеках."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "c_K2odulDEL2"
   },
   "outputs": [],
   "source": [
    "R = range(0.1, 35, length=256);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "t5oupaIf8rNq"
   },
   "source": [
    "Отлично, мы готовы рассчитать скорость звезд в Балдже и на диске. На этом этапе нам нужно отметить, как добавить скорости. В отличие от обычного сложения скорости складываются в квадратуре. Так,\n",
    "$V_{light}=\\sqrt{V_{Bulge}^{2}+V_{Disk}^2}$. Это очень похоже на теорему Пифагора. Посмотрите [здесь](https://colab.research.google.com/drive/15IsaQ0M00UGkzJRKTHyRKW35y22oW_GG#scrollTo=t5oupaIf8rNq&line=2&uniqifier=1).\n",
    "Будьте осторожны, расчет кривой вращения для балджа займет некоторое время.\n",
    "Требуется некоторое время, чтобы вычислить массу с помощью интегрирования (причудливое слово для суммирования). Так что не пугайтесь."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using QuadGK, Bessels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://github.com/JuliaMath/Bessels.jl  \n",
    "https://github.com/JuliaMath/QuadGK.jl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const G = 4.302e-6 # kpc M_sol^-1 (km/s)^2 #6.67408e-11 # m^3 kg^-1 s^-2\n",
    "\n",
    "rho_integrand(x, kappa, Rb, r) = \n",
    "        exp( kappa * (1.0 - (x / Rb)^0.25) ) / ( x^0.75 * sqrt(x^2 - r^2) )\n",
    "\n",
    "rho(r, sigma_be, kappa, Rb) = \n",
    "        kappa * sigma_be / (4pi * Rb^0.25) * \n",
    "            quadgk(x-> rho_integrand(x, kappa, Rb, r), r, Inf, rtol=1e-7)[1]\n",
    "\n",
    "Mb_integrand(r, sigma_be, kappa, Rb) = \n",
    "     rho(r, sigma_be, kappa, Rb) * r^2\n",
    "\n",
    "Mb(R, sigma_be, kappa, Rb) = \n",
    "        4pi * quadgk(r-> Mb_integrand(r, sigma_be, kappa, Rb), 0.0, R, rtol=1e-7)[1]\n",
    "\n",
    "Vb(R, sigma_be, kappa, Rb) = \n",
    "     sqrt(G * Mb(R, sigma_be, kappa, Rb) / R)\n",
    "\n",
    "Vd(R, Rd, sigma_0) = begin\n",
    "    y = 0.5R / Rd\n",
    "    sqrt(4pi*G*sigma_0*Rd*y^2 * (besseli0(y) * besselk0(y) - besseli1(y) * besselk1(y)) )\n",
    "end\n",
    "\n",
    "Vh(R, Vinf, Rh) = Vinf * sqrt(1.0 - (Rh / R) * atan(R / Rh));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Формулы:  \n",
    "http://www.astro.yale.edu/vdbosch/astro610_lecture18.pdf  \n",
    "https://ned.ipac.caltech.edu/level5/Sept16/Sofue/Sofue4.html  \n",
    "https://academic.oup.com/pasj/article/61/2/227/2898209?login=true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Qf542s1O8zsJ"
   },
   "outputs": [],
   "source": [
    "# Here are some constants that you will need \n",
    "sigma_be = 2.7e9#3.2e9\n",
    "kappa = 7.6695\n",
    "Rb = 0.5#0.5\n",
    "sigma_0 = 8.6e8#8.44e8\n",
    "Rd = 3.5#3.5\n",
    "Vinf = 220#200\n",
    "Rh = 16#12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Qf542s1O8zsJ"
   },
   "outputs": [],
   "source": [
    "Vbul = [ Vb(r, sigma_be, kappa, Rb) for r in R] # вложенный интеграл с бесконечным перделом\n",
    "Vdisk = [ Vd(r, Rd, sigma_0) for r  in R];  # считает долго... еще по чаю?\n",
    "Vtot = sqrt.(Vbul.^2 + Vdisk.^2);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    " plot( R, Vtot, lab = \"Vtot\", line = (3, :red) )\n",
    "plot!( Rs, Vs, yerr = dVs, lab = \"obs\", line = (3, :gray), m = (3,:gray, 0.5) )\n",
    "plot!( R, Vdisk, lab = \"disk\", line = (2, :green) )\n",
    "plot!( R, Vbul,  lab = \"bulge\", line = (2, :orange), \n",
    "    xlab = \"R, kpc\", ylab = \"V, km/s\",\n",
    "    xlim = (0, 35), ylim = (0,250))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "savefig(\"MilkyWayRC.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Y4RW6l2CD10d"
   },
   "source": [
    "Как мы видим, только сумма скоростей балджа и диска действительно хорошо соответствует наблюдаемой кривой вращения при малых R.\n",
    "Однако при больших R полная скорость значительно ниже наблюдаемой скорости.\n",
    "Это говорит нам о том, что нам не хватает какой-то дополнительной массы, которая оказывается ореолом темной материи.\n",
    "\n",
    "Далее, к балджу и скорости диска мы добавим профиль скорости гало."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RsVv9n6d_APg"
   },
   "outputs": [],
   "source": [
    "Vinf = 220\n",
    "Vdm = [ Vh(r, Vinf, Rh ) for r in R ]\n",
    "\n",
    "# The halo velocity is calculated. \n",
    "# Try to add this to the bulge and disk velocity in the same way we did above:\n",
    "Vtotal = sqrt.( Vbul.^2 + Vdisk.^2 + Vdm.^2 );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " plot( R, Vtotal, lab = \"Vtot\", line = (3, :red) )\n",
    "plot!( Rs, Vs, lab = \"obs\", line = (3, :gray), m = (3,:gray, 0.5) )\n",
    "plot!( R, Vdm, lab = \"DM\", line = (2, :black) )\n",
    "#plot!( Rs, Vgas,  lab = \"gas\", line = (2, :blue) )\n",
    "plot!( R, Vdisk, lab = \"disk\", line = (2, :green) )\n",
    "plot!( R, Vbul,  lab = \"bulge\", line = (2, :orange), \n",
    "    xlab = \"R, kpc\", ylab = \"V, km/s\",\n",
    "    xlim = (0, 35), ylim = (0,250))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "savefig(\"MilkyWayRC_DM.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RBq34kFVFyN9"
   },
   "source": [
    "Теперь, когда вы закончили, вы можете видеть, что внешняя полная скорость лучше совпадает с наблюдаемой кривой вращения.\n",
    "\n",
    "**Дополнительно:** Если вы закончили раньше, вернитесь назад и попробуйте изменить значение некоторых констант, показанных выше. В частности, изменить ```Vinf```, ```sigma_0``` и ```sigma_be```.\n",
    "\n",
    "**ANS** Все константы, представленны на рис. Эти константы наблюдаются и рассчитываются для Млечного Пути, и их изменение может привести к расхождению между наблюдаемой и расчетной кривой вращения.\n",
    "В любом случае увеличение```Sigma_0```, ```Sigma_be``` или ```V_inf``` увеличит скорость на кривой вращения (фиолетовая сплошная линия) и наоборот.\n",
    "\n",
    "**Дополнительно:** Мы видели, что квадрат скорости равен массе, так что же это означает, когда мы увеличиваем или уменьшаем скорость? Что это означает для конкретных компонентов Млечного Пути и количества вещества в разных компонентах?\n",
    "\n",
    "**ANS** Наблюдение за тем, что квадрат скорости равен массе, означает, что увеличение скорости связано с добавлением массы (которое изменяет гравитацию, ощущаемую массой) в разных частях кривой вращения. Таким образом, увеличение, например, ``V_inf`` означает увеличение количества темной материи."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "VFE8WsogcjE-"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "colab": {
   "collapsed_sections": [],
   "name": "SOLUTIONS: Let’s Talk Astrophysics Workshop.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Julia 1.7.2",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}