{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "569a05dc-d2f8-4d59-89a4-80e79befda51",
   "metadata": {},
   "outputs": [],
   "source": [
    "using CSV\n",
    "using DataFrames\n",
    "using CairoMakie\n",
    "using Distributions\n",
    "using StatsBase\n",
    "using Printf\n",
    "using Dierckx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "307cd583-c97d-48c7-bdd9-2011218d67f0",
   "metadata": {},
   "source": [
    "# Symulacja spowalniania neutronów\n",
    "\n",
    "Korzystając z wyprowadzonych wzorów spróbujemy policzyć liczbę potrzebnych zderzeń, aby neutron został spowolniony z energii $E_0$ do energii $E_f$. Skorzystamy z następujących wzorów:\n",
    "1. Kąt rozproszenia w układzie CM jest izotropowy, zatem $\\theta_C$ pochodzi z rozkładu $$(\\theta_C) = 1/2\\sin\\theta_C$$\n",
    "2. Energia po rozproszeniu zależy od kąta w układzie CM jak\n",
    "$$ E_L' = \\frac{(1+\\alpha)+(1-alpha)\\cos(\\theta_C)}{2}E_L $$\n",
    "gdzie\n",
    "$$\\alpha = \\left(\\frac{A-1}{A+1}\\right)^2$$\n",
    "3. Początkowa energia pochodzi z rozkładu Maxwella o temperaturze $kT = 1.29$ MeV"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2894e3e8-6b2f-41af-a0dc-37fbc0b1a5e9",
   "metadata": {},
   "source": [
    "## Proste rozpraszanie\n",
    "\n",
    "Dla pojedynczego neutronu, rozpraszanie na deuterze"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3b4a5ed8-16da-4421-9e4b-1bc89e575552",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n"
     ]
    }
   ],
   "source": [
    "A = 2\n",
    "α = ((A-1)/(A+1))^2\n",
    "\n",
    "E0 = 1.29e6\n",
    "El = E0\n",
    "Ef = 25e-3 \n",
    "n = 0\n",
    "while El > Ef\n",
    "    n += 1\n",
    "    θc = acos((1 - 2 * rand()))  # Metoda odwrotnej dystrybuanty\n",
    "    El = ((1 + α) + (1 - α) * cos(θc)) / 2 * El\n",
    "end\n",
    "println(n)    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0823d89-ae9e-4965-8edd-defe244a5f12",
   "metadata": {},
   "source": [
    "Zbierzemy to w funkcję i przeprowadzimy losowanie historii wielu neutronów. Rozkład [Maxwella w funkcji energii](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution#Distribution_for_the_energy) ma postać rozkładu [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9e7c4dc7-8f11-4ecc-9a40-b827b8f71bee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "simple_slowdown (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function simple_slowdown(N; kT=1.29e6, A=2, Ef=25e-3)\n",
    "    \n",
    "    ns = Int64[]\n",
    "    maxwell = Gamma(3/2, kT)\n",
    "    α = ((A-1)/(A+1))^2\n",
    "    for i in 1:N\n",
    "        E0 = rand(maxwell)\n",
    "        El = E0\n",
    "        n = 0\n",
    "        while El > Ef\n",
    "            n += 1\n",
    "            θc = acos((1 - 2 * rand()))\n",
    "            El = ((1 + α) + (1 - α) * cos(θc)) / 2 * El\n",
    "        end\n",
    "        push!(ns, n)\n",
    "    end\n",
    "    @printf(\"%.3f\\n\", mean(ns))\n",
    "    h = fit(Histogram, ns, 0:maximum(ns))\n",
    "    fig = Figure(size=(800, 500))\n",
    "    ax = Axis(fig[1, 1]; xlabel=L\"N_\\text{hits}\", xlabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, limits=(0, nothing, 0, nothing))\n",
    "    stairs!(ax, h.edges[1][1:end-1], h.weights ./ sum(h.weights))\n",
    "    fig\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de8fbe82-21e1-429c-aca9-5c5ec90f3dca",
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_slowdown(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f754185-63f4-409d-b68f-361a61b79f72",
   "metadata": {},
   "source": [
    "## Rozpraszanie z wychwytem\n",
    "\n",
    "W bardziej zaawansowanej wersji powinniśmy wziąć pod uwagę, że neutrony mogą podczas spowalniania ulec wychwytowi. Aby to uwzględnić potrzebne są nam przekroje czynne na różne procesy - w szczególności na rozpraszanie elastyczne i na wychwyt. Weźmiemy je z bazy Nuclear Data Section IAEA:\n",
    "\n",
    "https://www-nds.iaea.org/exfor/endf.htm\n",
    "\n",
    "Do wczytania danych użyjemy biblioteki DataFrames (odpowiednik Pandas pythona)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86dcf594-74c3-45a2-b142-be9cd0c56c61",
   "metadata": {},
   "outputs": [],
   "source": [
    "h1_el = DataFrame(CSV.File(\"h1_n_el.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92bc2ef9-fe63-4f83-b3f2-6c62610ee5fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = Figure(size=(1000, 600))\n",
    " ax = Axis(fig[1, 1]; xlabel=L\"$E$ (eV)\", ylabel=L\"\\sigma_{EL}\", xlabelsize=24, ylabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, xscale=log10, yscale=log10, limits=(1e-5, 2e7, 1e-6, 1e4))\n",
    "scatter!(ax, h1_el.E, h1_el.cs, label=L\"$^{1}H$ (EL)\")\n",
    "current_figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9199162-13ca-40c1-a295-7ff16f301c08",
   "metadata": {},
   "source": [
    "Aby móc liczyć przekroje czynne dla dowolnych energii, musimy użyć interpolacji."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51246422-675e-405b-b7fe-3af503f8bea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_h1_el = Spline1D(h1_el.E, h1_el.cs, k=1, bc=\"extrapolate\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f40a9692-e890-4b49-a480-f79c02f89f19",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ei = 10 .^ (-6:0.01:7.5)\n",
    "lines!(ax, Ei, s_h1_el(Ei), color=:blue)\n",
    "current_figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c4a6192-479b-40c3-958d-d16d26466827",
   "metadata": {},
   "source": [
    "Teraz wczytamy przekroje na wychwyt $(n, \\gamma)$ i również stworzymy funkcję interpolującą"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "613a3798-ec5b-4b4d-b165-f200992f88aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "h1_g = DataFrame(CSV.File(\"h1_n_g.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]));\n",
    "scatter!(ax, h1_g.E, h1_g.cs, marker=:utriangle, color=:tomato, label=L\"$^{1}H$ $(n,\\gamma)$\")\n",
    "s_h1_g = Spline1D(h1_g.E, h1_g.cs, k=1, bc=\"extrapolate\")\n",
    "lines!(ax, Ei, s_h1_g(Ei), color=:red)\n",
    "axislegend(position=:lb)\n",
    "current_figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "708781ff-936f-4d4e-aa39-66bd9d7842c1",
   "metadata": {},
   "source": [
    "Te same operacje musimy zrobić dla pozostałych danych (C-12, O-16, U-235), przekroje elastyczne i na wychwyt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae59eab3-c8c6-412c-b084-236bec99bdc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "c12_el = DataFrame(CSV.File(\"c12_n_el.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]))\n",
    "s_c12_el = Spline1D(c12_el.E, c12_el.cs, k=1, bc=\"extrapolate\")\n",
    "c12_g = DataFrame(CSV.File(\"c12_n_g.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]))\n",
    "s_c12_g = Spline1D(c12_g.E, c12_g.cs, k=1, bc=\"extrapolate\")\n",
    "ax2 = Axis(fig[1, 2]; xlabel=L\"$E$ (eV)\", ylabel=L\"\\sigma_{EL}\", xlabelsize=24, ylabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, xscale=log10, yscale=log10, limits=(1e-5, 2e7, 1e-6, 1e4))\n",
    "scatter!(ax2, c12_el.E, c12_el.cs, label=L\"$^{12}C$ (EL)\")\n",
    "scatter!(ax2, c12_g.E, c12_g.cs, marker=:utriangle, color=:tomato, label=L\"$^{12}C$ $(n,\\gamma)$\")\n",
    "lines!(ax2, Ei, s_c12_el(Ei), color=:blue)\n",
    "lines!(ax2, Ei, s_c12_g(Ei), color=:red)\n",
    "axislegend(position=:lb)\n",
    "\n",
    "o16_el = DataFrame(CSV.File(\"o16_n_el.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]))\n",
    "s_o16_el = Spline1D(o16_el.E, o16_el.cs, k=1, bc=\"extrapolate\")\n",
    "o16_g = DataFrame(CSV.File(\"o16_n_g.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]))\n",
    "s_o16_g = Spline1D(o16_g.E, o16_g.cs, k=1, bc=\"extrapolate\")\n",
    "\n",
    "ax3 = Axis(fig[2, 1]; xlabel=L\"$E$ (eV)\", ylabel=L\"\\sigma_{EL}\", xlabelsize=24, ylabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, xscale=log10, yscale=log10, limits=(1e-5, 2e7, 1e-6, 1e4))\n",
    "scatter!(ax3, o16_el.E, o16_el.cs, label=L\"$^{16}O$ (EL)\")\n",
    "scatter!(ax3, o16_g.E, o16_g.cs, marker=:utriangle, color=:tomato, label=L\"$^{16}O$ $(n,\\gamma)$\")\n",
    "lines!(ax3, Ei, s_o16_el(Ei), color=:blue)\n",
    "lines!(ax3, Ei, s_o16_g(Ei), color=:red)\n",
    "axislegend(positon=:lb)\n",
    "\n",
    "# W rejonie rezonansów są powtarzające się energie, które trzeba wyrzucić przed interpolacją!\n",
    "uniqueidx(v) = unique(i->v[i], eachindex(v))\n",
    "fe_el = DataFrame(CSV.File(\"fe0_n_el.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]))\n",
    "#fe_el = fe_el[uniqueidx(fe_el.E), :]\n",
    "s_fe_el = Spline1D(fe_el.E, fe_el.cs, k=1, bc=\"extrapolate\")\n",
    "fe_g = DataFrame(CSV.File(\"fe0_n_g.txt\", comment=\"#\", delim=\" \", ignorerepeated=true, header=[\"E\", \"cs\", \"interpolation\"]))\n",
    "#fe_g = fe_g[uniqueidx(fe_g.E), :]\n",
    "s_fe_g = Spline1D(fe_g.E, fe_g.cs, k=1, bc=\"extrapolate\")\n",
    "\n",
    "ax4 = Axis(fig[2, 2]; xlabel=L\"$E$ (eV)\", ylabel=L\"\\sigma_{EL}\", xlabelsize=24,\n",
    "          ylabelsize=24, xticklabelsize=22, yticklabelsize=22, xscale=log10,\n",
    "          yscale=log10, limits=(1e-5, 2e7, 1e-6, 1e4))\n",
    "scatter!(ax4, fe_el.E, fe_el.cs, label=\"Fe (EL)\")\n",
    "scatter!(ax4, fe_g.E, fe_g.cs, marker=:utriangle, \n",
    "         color=:tomato, label=L\"Fe $(n,\\gamma)$\")\n",
    "lines!(ax4, Ei, s_fe_el(Ei), color=:blue)\n",
    "lines!(ax4, Ei, s_fe_g(Ei), color=:red)\n",
    "axislegend(ax4, position=:lb)\n",
    "\n",
    "current_figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce41c2ae-0d5b-4fa0-a939-7cda83b17009",
   "metadata": {},
   "source": [
    "Mając dane na temat przekrojów czynnych możemy wrócić do obliczeń. Dla każdego kroku będziemy teraz liczyli całkowity przekrój czynny (suma elastycznego i wychwytu, pomniejsze procesy zaniedbujemy), a następnie losowali zachodzący proces. \n",
    "\n",
    "```\n",
    "0 ...                          1 \n",
    "\n",
    "|-----------------------|------|\n",
    "   el                     (n,g)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cfe394c-dbc8-45c2-83ee-8baf86bd1fb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "function slowdown(N, s_el, s_g, A; kT=1.29e6, Ef=25e-3)\n",
    "    ns = Int64[]\n",
    "    nc = Int64[]\n",
    "    Eg = Float64[]\n",
    "    maxwell = Gamma(3/2, kT)\n",
    "    α = ((A-1)/(A+1))^2\n",
    "    for i in 1:N\n",
    "        E0 = rand(maxwell)\n",
    "        El = E0\n",
    "        n = 0\n",
    "        capture = false\n",
    "        while El > Ef\n",
    "            n += 1\n",
    "            cs_tot = s_el(El) + s_g(El)\n",
    "            r = rand()\n",
    "            if r < s_el(El) / cs_tot\n",
    "                θc = acos((1 - 2 * rand()))\n",
    "                El = ((1 + α) + (1 - α) * cos(θc)) / 2 * El\n",
    "            else\n",
    "                capture = true\n",
    "                break\n",
    "            end\n",
    "        end\n",
    "        if capture\n",
    "            push!(nc, n)\n",
    "            push!(Eg, El)\n",
    "        else\n",
    "            push!(ns, n)\n",
    "        end\n",
    "    end\n",
    "   \n",
    "    fig = Figure(size=(1000, 500))\n",
    "    ax = Axis(fig[1, 1]; xlabel=L\"N_\\text{hits}\", xlabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, limits=(0, nothing, 0, nothing))\n",
    "    if length(ns) > 0\n",
    "        @printf(\"n = %.3f\\n\", mean(ns))\n",
    "        h = fit(Histogram, ns, 0:maximum(ns))\n",
    "        stairs!(ax, h.edges[1][1:end-1], h.weights ./ N)\n",
    "    end\n",
    "    if length(nc) > 0\n",
    "        @printf(\"c = %.1f%%\\n\", length(nc)/N * 100)\n",
    "        h = fit(Histogram, nc, 0:maximum(nc))\n",
    "        stairs!(ax, h.edges[1][1:end-1], h.weights ./ N)\n",
    "    end\n",
    "    \n",
    "    ax2 = Axis(fig[1, 2]; xlabel=L\"E_\\text{cap}\", xlabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, xscale=log10,\n",
    "                         limits=(1e-6, 1e7, 0, nothing)) \n",
    "    h = fit(Histogram, Eg, 10 .^ (-6.0:1.0:7.0))\n",
    "    stairs!(ax2, h.edges[1][1:end-1], h.weights ./ N)\n",
    "\n",
    "    fig\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa19df6a-30ec-4763-a6bf-eb9b8085dcf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "slowdown(10000, s_c12_el, s_c12_g, 12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a7f411f-f328-4141-8deb-e319d694e36e",
   "metadata": {},
   "outputs": [],
   "source": [
    "slowdown(10000, s_fe_el, s_fe_g, 56)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdd1ea20-6f4a-4e67-87dc-e22b29255580",
   "metadata": {},
   "source": [
    "Dla wody sytuacja jest trochę bardziej skomplikowana, bo mamy dwa rodzaje atomów. Losowanie przekroju czynnego będzie teraz wyglądało następująco:\n",
    "\n",
    "```\n",
    "0 ...                          1 \n",
    "\n",
    "|---------|------|-------|-----|\n",
    "\n",
    "  h1_el     h1_g   o16_el o16_g\n",
    "```\n",
    "\n",
    "Uwaga - przekrój czynny dla molekuł nie musi być sumą przekrojów czynnych poszczególnych składników. Dla powolnych neutronów ich długość fali jest porównywalna\n",
    "z rozmiarem całej cząsteczki i mogą występować anomalie. Np. dla wody przekrój czynny na rozpraszanie dla małych energii jest większa niż suma 2H + O."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "583685ca-8f57-4e90-9750-9620a2ebaa08",
   "metadata": {},
   "outputs": [],
   "source": [
    "function slowdown_water(N; kT=1.29e6, Ef=25e-3)\n",
    "    \n",
    "    ns = Int64[]\n",
    "    nc = Int64[]\n",
    "    Eg = Float64[]\n",
    "    maxwell = Gamma(3/2, kT)\n",
    "    α = ((A-1)/(A+1))^2\n",
    "    for i in 1:N\n",
    "        E0 = rand(maxwell)\n",
    "        El = E0\n",
    "        n = 0\n",
    "        capture = false\n",
    "        while El > Ef\n",
    "            n += 1\n",
    "            cs_tot = s_h1_el(El) + s_h1_g(El) + s_o16_el(El) + s_o16_g(El)\n",
    "            r = rand()\n",
    "            if r < (s_h1_el(El) + s_h1_g(El)) / cs_tot\n",
    "                A = 1\n",
    "                if r < s_h1_el(El) / cs_tot\n",
    "                    α = ((A-1)/(A+1))^2\n",
    "                    θc = acos((1 - 2 * rand()))\n",
    "                    El = ((1 + α) + (1 - α) * cos(θc)) / 2 * El\n",
    "                else\n",
    "                    capture = true\n",
    "                    break\n",
    "                end\n",
    "            else\n",
    "                A = 16\n",
    "                if r < (s_h1_el(El) + s_h1_g(El) + s_o16_el(El)) / cs_tot\n",
    "                    α = ((A-1)/(A+1))^2\n",
    "                    θc = acos((1 - 2 * rand()))\n",
    "                    El = ((1 + α) + (1 - α) * cos(θc)) / 2 * El\n",
    "                else\n",
    "                    capture = true\n",
    "                    break\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "        if capture\n",
    "            push!(nc, n)\n",
    "            push!(Eg, El)\n",
    "        else\n",
    "            push!(ns, n)\n",
    "        end\n",
    "    end\n",
    "    @printf(\"n = %.3f\\n\", mean(ns))\n",
    "    @printf(\"c = %.1f%%\\n\", length(nc)/N * 100)\n",
    "\n",
    "    \n",
    "    fig = Figure(size=(1000, 500))\n",
    "    ax = Axis(fig[1, 1]; xlabel=L\"N_\\text{hits}\", xlabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, limits=(0, nothing, 0, nothing))\n",
    "    h = fit(Histogram, ns, 0:maximum(ns))\n",
    "    stairs!(ax, h.edges[1][1:end-1], h.weights ./ N)\n",
    "    h = fit(Histogram, nc, 0:maximum(ns))\n",
    "    stairs!(ax, h.edges[1][1:end-1], h.weights ./ N)\n",
    "    \n",
    "    ax2 = Axis(fig[1, 2]; xlabel=L\"E_\\text{cap}\", xlabelsize=24, xticklabelsize=22, \n",
    "                         yticklabelsize=22, xscale=log10, \n",
    "                         limits=(1e-6, 1e7, 0, nothing)) \n",
    "    h = fit(Histogram, Eg, 10 .^ (-6.0:1.0:7.0))\n",
    "    stairs!(ax2, h.edges[1][1:end-1], h.weights ./ N)\n",
    "    \n",
    "    fig\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6daf0561-fc89-4380-9da7-64c4c3e347bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "slowdown_water(10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "441ab362-2110-416c-be03-57ea0f10089c",
   "metadata": {},
   "source": [
    "## Makroskopowy przekrój czynny\n",
    "\n",
    "Opisuje prawdopodobieństwo na jednostkę odległości, że neutron wejdzie w reakcję $i$ z materiałem\n",
    "$$\\Sigma_i = N\\sigma_i = \\frac{\\rho N_A}{A}\\sigma_i$$\n",
    "\n",
    "Dla mieszaniny\n",
    "$$\\Sigma_i = \\sum_k N_k\\Sigma_{ik}$$\n",
    "\n",
    "gdzie $N$ to liczba atomów na cm$^3$.\n",
    "\n",
    "Policzymy makroskopowy przekrój czynny na wychwyt dla neutronów termicznych dla [węglika boru](https://en.wikipedia.org/wiki/Boron_carbide) (używanego m.in. jako prętów kontrolnych). Podany w tabeli $\\sigma_\\gamma$ jest w rzeczywistości dla $^{10}$B reakcją $(n, \\alpha)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "887f6fb9-4f05-4909-b0f6-8eaadbefa7b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "ρ = 2.50\n",
    "NA = 6.022e23\n",
    "MB = 10.8\n",
    "MC = 12\n",
    "MB4C = 4 * 10.8 + 12\n",
    "@show MB4C\n",
    "NB4C = ρ * NA / MB4C\n",
    "@show NB4C\n",
    "NB = 4 * NB4C\n",
    "@show NB\n",
    "NB10 = 0.199 * NB\n",
    "NB11 = 0.801 * NB\n",
    "@show NB10\n",
    "@show NB11\n",
    "NC = NB4C\n",
    "@show NC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b8167be-af75-435f-8aa5-3a80ab225ca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ΣB4C = 3868e-24 * NB10 + 0.005e-24 * NB11 + 0.004e-24 * NC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e45f4c6d-a1bf-403e-a35c-5497c315e31a",
   "metadata": {},
   "outputs": [],
   "source": [
    "λB4C = 1 / ΣB4C"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64dbb413-e417-4efe-9675-948abb147ed9",
   "metadata": {},
   "source": [
    "## Wzbogacenie\n",
    "\n",
    "Rozpowszechnienie pierwiastków jest podawane w procentach cząsteczkowych $\\gamma_i$ (jak wyżej dla B-10 i B-11) i łatwo użyć to w obliczeniach makroskopowego przekroju czynnego.\n",
    "\n",
    "Wzbogacenie uranu podawane jest jednak w procentach wagowych ($w_i$). Ponieważ U-235 i U-238 mają inne masy, przeliczenie wymaga uwzględnienie tego faktu.\n",
    "\n",
    "W 1 gramie materiału liczba atomów izotopu $i$\n",
    "$$n_i = w_i\\frac{N_A}{M_i}$$\n",
    "przeliczając na wzbogacenie cząsteczkowe\n",
    "$$\\gamma_i = \\frac{n_i}{\\sum_in_i} = \\frac{\\frac{w_i}{M_i}}{\\sum_i\\frac{w_i}{M_i}}$$\n",
    "a następnie na liczbę atomów\n",
    "$$N_i = \\gamma_i\\frac{\\rho N_A}{M}$$\n",
    "gdzie $\\rho$ i $M$ są gęstością i masą molową w stanie niewzbogaconym (bo liczba cząstek w materiale wzbogaconym jest taka sama jak w niewzbogaconym!)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5ac4d06-02ab-4030-b8a9-2145ee2c51d0",
   "metadata": {},
   "source": [
    "Funckja przeliczająca wzbogacenie na rozpowszechnienie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f533a4e-7481-4f4f-9ac8-91e0f676e09f",
   "metadata": {},
   "outputs": [],
   "source": [
    "function abundance(w235)\n",
    "    w238 = 1 - w235\n",
    "    M235 = 235.04\n",
    "    M238 = 238.05\n",
    "    w235/M235 / (w235/M235 + w238/M238) \n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b834ab4b-9539-4136-bbc3-0bf899106e9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = Figure(size=(800, 500))\n",
    "ax = Axis(fig[1, 1]; xlabel=L\"$w$ (%)\", ylabel=L\"$\\gamma - w$ (%)\", \n",
    "                     xlabelsize=24, ylabelsize=24, xticklabelsize=22, yticklabelsize=22)\n",
    "w = collect(0.0:0.01:1.0)\n",
    "lines!(ax, w .* 100, (abundance.(w) .- w) .* 100)\n",
    "current_figure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e34fa3e-4393-4e0c-91c8-e1ce3d00e7e4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.11.6",
   "language": "julia",
   "name": "julia-1.11"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}