Hamiltonian fitting

[1]:

import numpy as np
import matplotlib.pyplot as plt
import galport
import agama

Investigate the 2D bar

Load potential. For this purpose, it is best to take a potential symmetrical with respect to the \(xy\)-plane.

[2]:

pot_gal = agama.Potential(file='data/Pot_tri_CylSpline_t400.ini')
pot_gal_sym = agama.Potential(file='data/Pot_sph_CylSpline_t400.ini')

Omega = np.load('data/omega_p.npy')[400*8]

Function for plotting

[3]:

def make_H_mesh(hamiltionian):
    ntheta = 1001
    nJ = 801
    theta = np.linspace(-np.pi,np.pi,ntheta)
    J = np.linspace(0.0,0.8,nJ)

    X, Y = np.meshgrid(theta, J)

    x_sh = X.reshape(ntheta*nJ)
    y_sh = Y.reshape(ntheta*nJ)
    q = np.zeros((ntheta*nJ, 3))
    q[:,0] = x_sh
    q[:,2] = y_sh

    H_sh = hamiltionian(y_sh, x_sh)
    x_stack = np.vstack((y_sh, x_sh))
    dot = hamiltionian.derivative(y_sh, x_sh)
    H = H_sh.reshape((nJ, ntheta))
    dot_J = dot[0].reshape((nJ, ntheta))
    dot_theta = dot[1].reshape((nJ, ntheta))

    return X, Y, H, dot_J, dot_theta

def plot_phase_portrate(hamiltionian, x_data, Hone):

    X, Y, H, dot_J, dot_theta = make_H_mesh(hamiltionian)

    J_data = np.reshape(x_data[:,:,1], len(x_data[:,:,1])*len(x_data[0,:,1]))
    phi_data = np.reshape(x_data[:,:,3], len(x_data[:,:,3])*len(x_data[0,:,3]))
    H_data = hamiltionian.hamiltonian(J_data, phi_data)
    H_data_tab = np.reshape(H_data, (len(x_data[:,:,3]),len(x_data[0,:,3])))
    H_levels = np.nanmean(H_data_tab, axis=1)

    fig, ax = plt.subplots(figsize=(12, 7))

    levels = np.sort(H_levels)
    CS = ax.contour(X, Y, H, levels=levels, colors='black')
    CS = ax.contour(X, Y, dot_theta, levels=[0], colors='blue')
    CS = ax.contour(X, Y, dot_J, levels=[0], colors='green')

    ax.scatter((x_data[:,:,3]+np.pi)%(2*np.pi)-np.pi, x_data[:,:,1], s=1, linewidth=0)
    ax.set_xlabel("$2\\theta_\\varphi - \\theta_R$", fontsize=60)
    ax.set_ylabel("$J_R$", fontsize=60)

    ticks = [-np.pi, -np.pi/2, 0, np.pi/2, np.pi]
    ticks_text = ['$-\\pi$', '$-\\pi/2$', '$0$', '$\\pi/2$', '$\\pi$']
    ax.tick_params(labelsize=30)
    ax.set_xticks(ticks, ticks_text)
    ax.set_ylim([0,np.nanmax(J_data)*1.05])
    ax.set_xlim([-np.pi,np.pi])

    plt.show()
    plt.close()

[4]:

HJ = -1.8

HF =galport.HFitting(potential=pot_gal, Omega=Omega)
H_vilr = HF.fit(H=HJ, Htype='bar_2d')
plot_phase_portrate(H_vilr, HF.phasecoord, HJ)

10 orbits complete (370.4 orbits/s)

Investigate the vertical structure of the bar

Load potential

[5]:

pot_gal = agama.Potential(file='data/Pot_non_CylSpline_t400.ini')
pot_gal_sym = agama.Potential(file='data/Pot_axi_CylSpline_t400.ini')

Omega = np.load('data/omega_p.npy')[400*8]

[6]:

def make_H_mesh(hamiltionian):
    ntheta = 1001
    nJ = 801
    theta = np.linspace(-np.pi,np.pi,ntheta)
    J = np.linspace(0.0,0.8,nJ)

    X, Y = np.meshgrid(theta, J)

    x_sh = X.reshape(ntheta*nJ)
    y_sh = Y.reshape(ntheta*nJ)
    q = np.zeros((ntheta*nJ, 3))
    q[:,0] = x_sh
    q[:,2] = y_sh

    H_sh = hamiltionian(y_sh, x_sh)
    x_stack = np.vstack((y_sh, x_sh))
    dot = hamiltionian.derivative(y_sh, x_sh)
    H = H_sh.reshape((nJ, ntheta))
    dot_J = dot[0].reshape((nJ, ntheta))
    dot_theta = dot[1].reshape((nJ, ntheta))

    return X, Y, H, dot_J, dot_theta

def plot_phase_portrate(hamiltionian, x_data, Hone):

    X, Y, H, dot_J, dot_theta = make_H_mesh(hamiltionian)

    J_data = np.reshape(x_data[:,:,1], len(x_data[:,:,1])*len(x_data[0,:,1]))
    phi_data = np.reshape(x_data[:,:,3], len(x_data[:,:,3])*len(x_data[0,:,3]))
    H_data = hamiltionian.hamiltonian(J_data, phi_data)
    H_data_tab = np.reshape(H_data, (len(x_data[:,:,3]),len(x_data[0,:,3])))
    H_levels = np.nanmean(H_data_tab, axis=1)

    fig, ax = plt.subplots(figsize=(12, 7))

    levels = np.sort(H_levels)
    CS = ax.contour(X, Y, H, levels=levels, colors='black')
    CS = ax.contour(X, Y, dot_theta, levels=[0], colors='blue')
    CS = ax.contour(X, Y, dot_J, levels=[0], colors='green')

    ax.scatter((x_data[:,:,3]+np.pi)%(2*np.pi)-np.pi, x_data[:,:,1], s=1, linewidth=0)
    ax.set_xlabel("$\\theta_z - \\theta_R$", fontsize=60)
    ax.set_ylabel("$J_z$", fontsize=60)

    ticks = [-np.pi, -np.pi/2, 0, np.pi/2, np.pi]
    ticks_text = ['$-\\pi$', '$-\\pi/2$', '$0$', '$\\pi/2$', '$\\pi$']
    ax.tick_params(labelsize=30)
    ax.set_xticks(ticks, ticks_text)
    ax.set_ylim([0,np.nanmax(J_data)*1.05])
    ax.set_xlim([-np.pi,np.pi])

    plt.show()
    plt.close()

Quick start for the Hamiltonian fitting

[7]:

HJ = -1.8

HF =galport.HFitting(potential=pot_gal, axisym_potential=pot_gal_sym, Omega=Omega)
H_vilr = HF.fit(H=HJ, Htype='buckling')
plot_phase_portrate(H_vilr, HF.phasecoord, HJ)

994 orbits complete (292.5 orbits/s)
994 orbits complete (298.8 orbits/s)
998 orbits complete (275.5 orbits/s)
1000 orbits complete (274.4 orbits/s)
10 orbits complete (44.25 orbits/s)

Set the parameters of the Hamiltonian fitting

\[H = h_0 + \sum_i \left[h_i^c\cos(i\theta) + h_i^s\sin(i\theta)\right]\]

\[h = \sum_{j=0}^{j_{max}} \sqrt{J^j}\]

Norb: number of orbits for fitting

n: list of Fourier modes

Positive integers → cosine terms \(h_i^c\cos(i\theta)\)
Negative integers → sine terms \(h_i^s\sin(i\theta)\)

deg: \(j_{max}\) - maximal degree of \(\sqrt{J}\)

Tint: time of integrating orbit

Nint: number of points on every orbit

[9]:

HJ = -1.8

Norb = 20
n = [0, 1, 2, 3, -1, -2, -3]
# h0, h1c, h2c, h2c, h1s, h2s, h2s
deg = 4
Tint = 200
Nint = 20000

HF =galport.HFitting(potential=pot_gal, axisym_potential=pot_gal_sym, Omega=Omega)
H_vilr = HF.fit(H=HJ, Htype='buckling', Norb=Norb, Tint=Tint, Nint=Nint,
                n=n, deg=deg, weight_dphidt=1.)

plot_phase_portrate(H_vilr, HF.phasecoord, HJ)

1989 orbits complete (299.5 orbits/s)
1989 orbits complete (296.8 orbits/s)
1995 orbits complete (271.5 orbits/s)
2000 orbits complete (270.6 orbits/s)
20 orbits complete (69.69 orbits/s)