#import "@preview/diatypst:0.2.0": *

#show: slides.with(
  title: "N-Body project ",
  subtitle: "Computational Astrophysics, HS24",
  date: "04.02.2024",
  authors: ("Rémy Moll"),
  toc: false,
  // layout: "large",
  // ratio: 16/9,
)

#show footnote.entry: set text(size: 0.6em)
#set footnote.entry(gap: 3pt)
#set align(horizon)


#import "helpers.typ"


// Setup of code location
#let t1 = json("../task1.ipynb")
#let t2 = json("../task2-particle-mesh.ipynb")


// Finally - The real content
= N-body forces and analytical solutions

== Overview - the system
Get a feel for the particles and their distribution
#columns(2)[
  #helpers.image_cell(t1, "plot_particle_distribution")
  // Note: for visibility the outer particles are not shown.
  #colbreak()
  The system at hand is characterized by:
  - $N ~ 10^4$ stars
  - a _spherical_ distribution

  $==>$ treat the system as a *globular cluster*
  #footnote[Unit handling [#link(<task1:function_apply_units>)[code]]]
]

== Density
Compare the computed density
#footnote[Density sampling [#link(<task1:function_density_distribution>)[code]]]
with the analytical model provided by the _Hernquist_ model:

#grid(
  columns: (3fr, 4fr),
  inset: 0.5em,
  block[
    $
      rho(r) = M/(2 pi)  a / (r dot (r + a)^3) 
    $
    where we infer $a$ from the half-mass radius:
    $
      r_"hm" = (1 + sqrt(2)) dot a
    $
  ],
  block[
    #helpers.image_cell(t1, "plot_density_distribution")
  ]
)

// Having more bins means to have shells that are nearly empty
// => the error is large, NBINS = 30 is a good compromise



== Force computation
#grid(
  columns: (3fr, 2fr),
  inset: 0.5em,
  block[
    #helpers.image_cell(t1, "plot_force_radial") 
  ],
  block[
    Discussion:
    - the analytical
      #footnote[Analytical force [#link(<task1:function_analytical_forces>)[code]]]
      method replicates the behavior accurately
    - at small softenings the $N^2$
      #footnote[$N^2$ force [#link(<task1:function_n2_forces>)[code]]]
      method has noisy artifacts 
    - a $1 dot epsilon$
      #footnote[$epsilon$ computation [#link(<task1:function_interparticle_distance>)[code]]]
      softening is a good compromise between accuracy and stability
  ]
)



== Relaxation
We express system relaxation in terms of the dynamical time of the system.
$
  t_"relax" = overbrace(N / (8 log N), n_"relax") dot t_"crossing"
$
where the crossing time of the system can be estimated through the half-mass velocity $t_"crossing" = v(r_"hm")/r_"hm"$.

We find a relaxation of $approx 30 "Myr"$ ([#link(<task1:compute_relaxation_time>)[code]])


#grid(
  columns: (1fr, 1fr),
  inset: 0.5em,
  block[
    #image("relaxation.png")
  ],
  block[
    - Each star-star interaction contributes $delta v approx (2 G m )/b$
    - Shifting by $epsilon$ *dampens* each contribution
    - $=>$ relaxation time increases
  ]
)


= Particle Mesh

== Overview - the system
#page(
  columns: 2
)[
  #helpers.image_cell(t2, "plot_particle_distribution")

  $==>$ use $M_"sys" approx 10^4 M_"sol" + M_"BH"$
]




== Force computation
#helpers.code_reference_cell(t2, "function_mesh_force")

#helpers.image_cell(t2, "plot_force_radial")

#grid(
  columns: (2fr, 1fr),
  inset: 0.5em,
  block[
    #helpers.image_cell(t2, "plot_force_radial_single")
  ],
  block[
    - using the (established) baseline of $N^2$
      #footnote[$N^2$ force [#link(<task1:function_n2_forces>)[code]]]
      with $1 dot epsilon$
      #footnote[$epsilon$ computation [#link(<task1:function_interparticle_distance>)[code]]]
      softening
    - small grids
      #footnote[Mesh force [#link(<task2:function_mesh_force>)[code]]]
      are stable but inaccurate at the center
    - very large grids have issues with overdiscretization

    $==> 75 times 75 times 75$ as a good compromise
  ]
)



#helpers.image_cell(t2, "plot_force_computation_time")


== Time integration
*Integration step*
#helpers.code_reference_cell(t2, "function_runge_kutta")

*Timesteps*
Chosen such that displacement is small (compared to the inter-particle distance) [#link(<task2:integration_timestep>)[code]]:
$
  op(d)t = 10^(-3) dot S / v_"part"
$

*Full integration*

[#link(<task2:function_time_integration>)[code]]


#pagebreak()
== First results
#helpers.image_cell(t2, "plot_system_evolution")


== Varying the softening

#helpers.image_cell(t2, "plot_second_system_evolution")


== Stability [#link("../task2_nsquare_integration.gif")[1 epsilon]]
#page(
  columns: 2
)[
  #helpers.image_cell(t2, "plot_integration_stability")
]


== Particle mesh solver
#helpers.image_cell(t2, "plot_pm_solver_integration")

#page(
  columns: 2
)[
  #helpers.image_cell(t2, "plot_pm_solver_stability")
]


= Appendix - Code <appendix>

== Code
#helpers.code_reference_cell(t1, "function_apply_units")
<task1:function_apply_units>

#pagebreak(weak: true)

#helpers.code_reference_cell(t1, "function_density_distribution")
<task1:function_density_distribution>

#pagebreak(weak: true)

#helpers.code_reference_cell(t1, "function_analytical_forces")
<task1:function_analytical_forces>

#pagebreak(weak: true)

#helpers.code_reference_cell(t1, "function_n2_forces")
<task1:function_n2_forces>

#pagebreak(weak: true)

#helpers.code_reference_cell(t1, "function_interparticle_distance")
<task1:function_interparticle_distance>

#pagebreak(weak: true)

#helpers.code_cell(t1, "compute_relaxation_time")
<task1:compute_relaxation_time>

#pagebreak(weak: true)

#helpers.code_reference_cell(t2, "function_mesh_force")
<task2:function_mesh_force>

#pagebreak(weak: true)

#helpers.code_cell(t2, "integration_timestep")
<task2:integration_timestep>

#pagebreak(weak: true)

#helpers.code_cell(t2, "function_time_integration")
<task2:function_time_integration>




#context {
  counter(page).update(locate(<appendix>).page())
}