Gadget2 and HDF5

Kemarin saya iseng menginstal Gadget 2A code for cosmological simulations of structure formation.  Cara instalasi dapat dilihat pada website

Setelah berhasil menginstal, lalu dicoba examplenya.. (btw saya instal di server Perseus, Prodi Astronomi, karena dependensinya sudah terinstal semua, jadi mudah, haha). Example yang sudah saya coba adalah tumbukan galaksi. Boleh dibilang problem SPH (Smoothed Particle Hydrodynamics) klasik, dengan dua jenis partikel, perhitungan gravitasinya menggunakan metode Tree. Untuk tujuan pembelajaran saya sekaligus keluarkan outputnya dalam format hdf5, lalu cari tahu cara membaca format file tersebut. HDF5 adalah format file yang katanya paling efektif untuk data besar dan kompleks; banyak bidang sudah mengimplemetasikan format file ini antara lain di meteorologi dan astronomi. :D

Hasil output hdf5 coba saya baca dengan script python sederhana.. ingin coba multiprocess di python namun tidak sempat, ya sudahlah kapan2 saja. Setelah mengulik cara pakainya akhirnya berhasil, haha… berikut script sederhana untuk membaca dan plot tersebut:

#!/usr/bin/python -tt

import h5py
import numpy as numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plotter(ifile):
    f = h5py.File(ifile+".hdf5", 'r')
    halo = f['/PartType1']
    disk = f['/PartType2']

    fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3)
    ax1.scatter(halo['Coordinates'][:,0], halo['Coordinates'][:,1], color='red', s=0.01)
    ax1.scatter(disk['Coordinates'][:,0], disk['Coordinates'][:,1], color='blue', s=0.5)

    ax2.scatter(halo['Coordinates'][:,0], halo['Coordinates'][:,2], color='red', s=0.01)
    ax2.scatter(disk['Coordinates'][:,0], disk['Coordinates'][:,2], color='blue', s=0.5)

    ax2.scatter(halo['Coordinates'][:,1], halo['Coordinates'][:,2], color='red', s=0.01)
    ax2.scatter(disk['Coordinates'][:,1], disk['Coordinates'][:,2], color='blue', s=0.5)

    f.close() # don't forget to close it

if __name__ == "__main__":
    n_snap = 402
    for i in range(n_snap):
        if i < 10:
            j = "00"+str(i)
        elif i>=10 and i<100:
            j = "0"+str(i)
            j = str(i)
        ifile = "snapshot_"+j # without .hdf5

Setelah digabungkan “screenshot” berupa image menjadi movie, maka hasilnya sebagai berikut

ffmpeg -f image2 -r 15 -i 'snapshot_%03d.png' col.mp4

warna biru dan merah mewakili partikel disk dan halo. Lumayan menggambarkan yang di dapat di alam seperti ini (galaksi antennae.. mirip bentuk cintah.. <3 :v ):
Credit & Copyright: Star Shadows Remote Observatory and PROMPT/CTIO
(Jack Harvey, Steve Mazlin, Rick Gilbert, and Daniel Verschatse)

berikutnya ingin saya coba yang simulasi kosmologi.. :D

The Precession of Mercury’s Perihelion (simulation)

One of the classical test of general relativity is the perihelion precession of Mercury (also gravitational lensing & gravitational redshift).

Image source:

Under Newtonian physics, a two-body system consisting of a lone object orbiting a spherical mass would trace out an ellipse with the spherical mass at a focus. The point of closest approach, called the periapsis (or, as the central body in our Solar System is the sun, perihelion), is fixed. A number of effects in our solar system cause the perihelia of planets to precess (rotate) around the sun. The principal cause is the presence of other planets which perturb each other’s orbit. Another (much less significant) effect is solar oblateness.

In 1859, the French mathematician and astronomer Urbain Le Verrier reported that the slow precession of Mercury’s orbit around the Sun could not be completely explained by “Newtonian mechanics” and “perturbations by the known planets“.

Observation: The perihelion precession of Mercury is 574.10±0.65 arcseconds per century relative to the inertial ICFR (International Celestial Reference Frame: barycenter of solar system).

Classic calculation: Newtonian mechanics, taking into account all the effects from the other planets, predicts a precession 531.63±0.69 arcseconds per century.

In the early 20th century, Albert Einstein’s General Theory of Relativity provided the explanation for the observed precession. The effect is very small: the Mercurian relativistic perihelion advance excess is just 42.98 arcseconds per century. Similar, but much smaller, effects operate for other planets: 8.62 arcseconds per century for Venus, 3.84 for Earth, 1.35 for Mars, and 10.05 for 1566 Icarus.

source: wikipedia hahahaha

Sooo…. Lets try this out! :) with numerical simulation..

First, of course the theoretical background, you can find the derivation of the General relativity correction in literature. :D ahaha
Ok, we can calculate that in gravitational system consist of many particles, the Post-Newtonian equation of motion is [ujug-ujug.. mak bedunduk]:
then, become [klothaa..ak]:
(the derivation is a lot more complicated, LOL)
Source: Benitez and Gallardo, Celest Mech Dyn Astr 2008
I really want to calculate the gravity using above equation, but if we just want to know Mercury’s precession, that equation is overkill.. [halah]. (maybe next time, e.g. for binary Neutron star, globular cluster, or S-star around supermassive black hole :) )

Assume that only the central body give effect of general relativity, the general relativity correction given by this simple equation:

WOKE, lets do our simulation.. using Rebound integrator. Oh wait, there is no GR-correction in this integrator yet, so we can add it,
simply give “additional_force” like this:

void relativity_correction_force(){
    for (int i=1; i<N; i++){ // for all planet 
        struct particle* com = &(particles[i]);
        // relative to the Sun
        struct particle sun = particles[0];
        double dx = com->x - sun.x;
        double dy = com->y - sun.y;
        double dz = com->z - sun.z;
        double dvx = com->vx - sun.vx;
        double dvy = com->vy - sun.vy;
        double dvz = com->vz - sun.vz;

        double r = sqrt(dx*dx + dy*dy + dz*dz);
        double _r = 1./r; // reciprocal
        // Benitez and Gallardo 2008
        // alpha = k^2/(r^3*c^2); beta  = 4k^2/r - v.v; gamma = 4*(v.r)
        // delta_a = alpha*(beta*r + gamma*v)
        double alpha = G*_r*_r*_r*_c2;
        double beta  = 4.*G*_r - (dvx*dvx + dvy*dvy + dvz*dvz);
        double gamma = 4.*(dvx*dx + dvy*dy + dvz*dz);
        // add acceleration to comet
        com->ax += alpha * (beta*dx + gamma*dvx);
        com->ay += alpha * (beta*dy + gamma*dvy);
        com->az += alpha * (beta*dz + gamma*dvz);

or similar with this code. :)

using this program and some initial data of planet from HORIZON JPL NASA we can make simulation of “Solar System” for example in 10000 years. We want to know how much the different of the Mercury’s orbit with and without GR correction, hopefully I can get the value given above. ahahaha

Package Rebound
Integrator: IAS15 (similar with RADAU15)
time step: 0.1 days
integration time: 10000 years
including all planets, moon, and 3 massive asteroid.. LOL (why these asteroids here? because I’m too lazy to delete it)

computing time (OpenMP 8 threads)
with GR-correction 1011.007421 s
w/o (Classic) 948.855602 s

First of all lets remind you (if u already know) about Keplerian coordinate, hehe, its easier to look how the orbit change using this coordinate. :D

Mercury’s orbit change after 10000 years:





Different between integration w/o and with GR-correction is clearly visible in argument of perihelion parameter.
assume that the other parameters are exactly the same (not change due to GR-correction) we can simply calculate the precession of Mercury’s perihelion using argument of perihelion. After 10000 years the different is 1.22397536º or 44.063 arcseconds/century. Uh.. not exactly as calculated before (by someone).. hahaha. I think 1 arcseconds is tolerable for this simulation :P
This different perhaps due to:
– that last calculation only using argument of perihelion
– need a smaller “integrator_epsilon” (error)
– the precession is not linear! (look in paper above)
– didn’t include correction from Oblateness of the Sun (quadrupole moment of gravity)


Neville’s algorithm (polynomial interpolation and extrapolation)

Salah satu jenis interpolasi (dan ekstrapolasi) yang penting adalah interpolasi polinom Lagrange. Penurunannya dapat dilihat disini dan disitu (halah).

Interpolasi polinom tidak sering digunakan terutama pada kumpulan data yang besar, orang lebih sering menggunakan Cubic Spline untuk “memuluskan” data.. wkwkwk. Namun interpolasi polinom ini banyak digunakan dalam teknik integrasi juga penyelesaian ODE (ya iyalah).. salah satunya metode terkenal algoritma Bulirsch-Stoer-Gragg (ekstrapolasi ke h = 0) :D

Contoh program sederhana dalam C++

/* Neville's algorithm
 * Copyleft (2014). Ridlo W. Wibowo

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
using namespace std;

double x[100], y[100];
int N = -1;

double neville(double input);
void plot();

int main(){
    ifstream read("data.txt");
    while (!read.eof()){
        N += 1;
        read >> x[N] >> y[N];
        //cout << x[N] << " " << y[N] << endl; 
    cout << "Jumlah data : " << N << endl;
    /* interpolate input "x" */
    double input[10] = {1.1, 1.6, 2.0, 2.3, 2.8, 3.5, 4.2, 4.5, 4.8, 5.2}; // last data is outside region
    /* y-interpolate and extrapolate */        
    ofstream outNeville("data.out");
    for (int i=0; i<10; i++){
        outNeville << input[i] << " " << neville(input[i]) << endl;
    return 0;

double neville(double input){
    double Q[100][100];    
    // initiate first column    
    for (int i=0; i<N; i++){
        Q[i][0] = y[i];
    // iteration
    for (int i=1; i<N; i++){
        for (int j=1; j<i; j++){
            Q[i][j] = ((input-x[i-j])*Q[i][j-1] - (input-x[i])*Q[i-1][j-1])/(x[i] - x[i-j]);
            //cout << Q[i][j] << " ";
        //cout << endl;

    return Q[N-1][N-2]; // output last iteration

void plot(){
    ofstream ploter("");
    ploter << "# gnuplot input file\n";
    ploter << "set xlabel \"x\"\n";
    ploter << "set ylabel \"f(x)\"\n";
    ploter << "plot \"data.txt\" t \"data input\", \"data.out\" t \"interpolation + last data as extrapolation\" \n";
    system("gnuplot -persist <");
    system("rm -f");


0.9 1.3
1.3 1.5
1.9 1.85
2.1 2.1
2.6 2.6
3.0 2.7
3.9 2.4
4.4 2.15
4.7 2.05
5.0 2.1



Tangential of radial vector in velocity direction

How to compute the tangential of radial vector in velocity direction in 3-dimension space?

assume there is a particle move in 3D space relative to the origin, if radial vector from origin and velocity vector of particle are known, how to compute tangential of radial vector in velocity direction?

\vec{t} = (\vec{r} \times \vec{v}) \times \vec{r}

is it correct?