In [1]:

    

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
from matplotlib.colors import LogNorm
rc('text', usetex=True)
plt.rc('font', family='serif')

from __future__ import print_function


from IPython.html.widgets import interact


    

In [2]:

    

%%bash 
make
./non_axisymmetric_potential


    

    




make: `non_axisymmetric_potential' is up to date.
 energy (initial): 
  -99.999999999999986     
 ---
 x_max:
  0.52554206776627421     
 x_0: 
   5.2033868095670713E-003   1.0406773619134143E-002   1.5610160428701215E-002   2.0813547238268285E-002   2.6016934047835356E-002   3.1220320857402430E-002   3.6423707666969500E-002   4.1627094476536571E-002   4.6830481286103641E-002   5.2033868095670711E-002   5.7237254905237782E-002   6.2440641714804859E-002   6.7644028524371930E-002   7.2847415333939000E-002   7.8050802143506071E-002   8.3254188953073141E-002   8.8457575762640211E-002   9.3660962572207282E-002   9.8864349381774352E-002  0.10406773619134142       0.10927112300090849       0.11447450981047556       0.11967789662004263       0.12488128342960972       0.13008467023917680       0.13528805704874386       0.14049144385831094       0.14569483066787800       0.15089821747744508       0.15610160428701214       0.16130499109657923       0.16650837790614628       0.17171176471571337       0.17691515152528042       0.18211853833484751       0.18732192514441456       0.19252531195398165       0.19772869876354870       0.20293208557311579       0.20813547238268285       0.21333885919224993       0.21854224600181699       0.22374563281138407       0.22894901962095113       0.23415240643051821       0.23935579324008527       0.24455918004965235       0.24976256685921944       0.25496595366878649       0.26016934047835361       0.26537272728792066       0.27057611409748772       0.27577950090705478       0.28098288771662189       0.28618627452618894       0.29138966133575600       0.29659304814532306       0.30179643495489017       0.30699982176445723       0.31220320857402428       0.31740659538359134       0.32260998219315845       0.32781336900272551       0.33301675581229256       0.33822014262185962       0.34342352943142673       0.34862691624099379       0.35383030305056085       0.35903368986012790       0.36423707666969501       0.36944046347926207       0.37464385028882913       0.37984723709839618       0.38505062390796330       0.39025401071753035       0.39545739752709741       0.40066078433666447       0.40586417114623158       0.41106755795579863       0.41627094476536569       0.42147433157493275       0.42667771838449986       0.43188110519406692       0.43708449200363397       0.44228787881320103       0.44749126562276814       0.45269465243233520       0.45789803924190225       0.46310142605146931       0.46830481286103642       0.47350819967060348       0.47871158648017054       0.48391497328973759       0.48911836009930471       0.49432174690887176       0.49952513371843887       0.50472852052800599       0.50993190733757299       0.51513529414714010       0.52033868095670721     
 
 v_y: 
   3.0381311745351813        2.8007046742851394        2.6519834947348935        2.5411911206052049        2.4518085587611278        2.3762832523136650        2.3105022691120416        2.2519675731064264        2.1990434008927791        2.1505977884519520        2.1058134979351277        2.0640803603800215        2.0249302009598864        1.9879955670101483        1.9529824964597351        1.9196519448074321        1.8878067553566193        1.8572822935381119        1.8279395710333637        1.7996601030679478        1.7723419981018542        1.7458969405339730        1.7202478315127212        1.6953269221559091        1.6710743202940190        1.6474367841102600        1.6243667386627507        1.6018214673714768        1.5797624421759022        1.5581547645719307        1.5369666960322292        1.5161692610269693        1.4957359094270479        1.4756422277944596        1.4558656911623722        1.4363854485375085        1.4171821366338448        1.3982377173534362        1.3795353353296995        1.3610591924874709        1.3427944370877858        1.3247270651405076        1.3068438324051554        1.2891321754754230        1.2715801406682474        1.2541763196234927        1.2369097906728697        1.2197700651625851        1.2027470380180805        1.1858309419247866        1.1690123045690612        1.1522819084406664        1.1356307527441636        1.1190500170028015        1.1025310259659711        1.0860652154508128        1.0696440987606199        1.0532592333275224        1.0369021872245616        1.0205645051824592        1.0042376737286332       0.98791308503953701       0.97158199906104359       0.95523550340383367       0.93886447045952526       0.92245951110586288       0.90601092427221153       0.88950864151525000       0.87294216560319760       0.85630050191728879       0.83957208124132376       0.82274467221028469       0.80580528130916129       0.78874003782880164       0.77153406056369167       0.75417130223169937       0.73663436654499870       0.71890429147650492       0.70096029042198671       0.68277944047456207       0.66433630364269725       0.64560246216538941       0.62654594252083640       0.60713049338333491       0.58731466923778297       0.56705065133152244       0.54628270737169771       0.52494514449236129       0.50295953536864135       0.48023087470714426       0.45664211440560210       0.43204615446088446       0.40625367367693704       0.37901381128200495       0.34998178592812140       0.31866071415668218       0.28428696184621849       0.24557295523199454       0.20000333350556657       0.14106970513308573     
 starting body: 
           1
 starting body: 
           2
 starting body: 
           3
 starting body: 
           4
 starting body: 
           5
 starting body: 
           6
 starting body: 
           7
 starting body: 
           8
 starting body: 
           9
 starting body: 
          10
 starting body: 
          11
 starting body: 
          12
 starting body: 
          13
 starting body: 
          14
 starting body: 
          15
 starting body: 
          16
 starting body: 
          17
 starting body: 
          18
 starting body: 
          19
 starting body: 
          20
 starting body: 
          21
 starting body: 
          22
 starting body: 
          23
 starting body: 
          24
 starting body: 
          25
 starting body: 
          26
 starting body: 
          27
 starting body: 
          28
 starting body: 
          29
 starting body: 
          30
 starting body: 
          31
 starting body: 
          32
 starting body: 
          33
 starting body: 
          34
 starting body: 
          35
 starting body: 
          36
 starting body: 
          37
 starting body: 
          38
 starting body: 
          39
 starting body: 
          40
 starting body: 
          41
 starting body: 
          42
 starting body: 
          43
 starting body: 
          44
 starting body: 
          45
 starting body: 
          46
 starting body: 
          47
 starting body: 
          48
 starting body: 
          49
 starting body: 
          50
 starting body: 
          51
 starting body: 
          52
 starting body: 
          53
 starting body: 
          54
 starting body: 
          55
 starting body: 
          56
 starting body: 
          57
 starting body: 
          58
 starting body: 
          59
 starting body: 
          60
 starting body: 
          61
 starting body: 
          62
 starting body: 
          63
 starting body: 
          64
 starting body: 
          65
 starting body: 
          66
 starting body: 
          67
 starting body: 
          68
 starting body: 
          69
 starting body: 
          70
 starting body: 
          71
 starting body: 
          72
 starting body: 
          73
 starting body: 
          74
 starting body: 
          75
 starting body: 
          76
 starting body: 
          77
 starting body: 
          78
 starting body: 
          79
 starting body: 
          80
 starting body: 
          81
 starting body: 
          82
 starting body: 
          83
 starting body: 
          84
 starting body: 
          85
 starting body: 
          86
 starting body: 
          87
 starting body: 
          88
 starting body: 
          89
 starting body: 
          90
 starting body: 
          91
 starting body: 
          92
 starting body: 
          93
 starting body: 
          94
 starting body: 
          95
 starting body: 
          96
 starting body: 
          97
 starting body: 
          98
 starting body: 
          99
 starting body: 
         100
 
 energy (final): 
  -100.00239792713487     
 overall conservation:
 energy err:
   2.3979271348792965E-005

In [3]:

    

n_bodies = 100

i=1
t_i, x_i, y_i, z_i, v_x_i, v_y_i, v_z_i = np.loadtxt('body_' + str(i).zfill(4) + '.dat', unpack=True, skiprows=1)
n_points = t_i.size

t = np.empty( (n_bodies, n_points))
t[i-1,:]   = t_i
x = np.empty( (n_bodies, n_points))
x[i-1,:]   = x_i
y = np.empty( (n_bodies, n_points))
y[i-1,:]   = y_i
z = np.empty( (n_bodies, n_points))
z[i-1,:]   = z_i
v_x = np.empty( (n_bodies, n_points))
v_x[i-1,:] = v_x_i
v_y = np.empty( (n_bodies, n_points))
v_y[i-1,:] = v_y_i
v_z = np.empty( (n_bodies, n_points))
v_z[i-1,:] = v_z_i 


for i in range(2, n_bodies+1):
    t_i, x_i, y_i, z_i, v_x_i, v_y_i, v_z_i = np.loadtxt('body_' + str(i).zfill(4) + '.dat', unpack=True, skiprows=1)

    t[i-1,:]   = t_i
    x[i-1,:]   = x_i
    y[i-1,:]   = y_i
    z[i-1,:]   = z_i
    v_x[i-1,:] = v_x_i
    v_y[i-1,:] = v_y_i
    v_z[i-1,:] = v_z_i


    

In [4]:

    

plot_step_size = 10
def plot_solution(plot_start=0,body=1):
    print("plot range: t= ",plot_start,  " - ", plot_start + plot_step_size)
    print("x_0 = ", x[body-1, 0])
    t_tmp = t[body-1,:]
    x_tmp = x[body-1,:]
    y_tmp = y[body-1,:]
    
    indices = np.where( ((t_tmp > plot_start) & (t_tmp < plot_start + plot_step_size)))
    plt.plot(x_tmp[indices], y_tmp[indices], linestyle='-')
    plt.xlabel(r"$x$")
    plt.ylabel(r"$y$")
    plt.legend([r"Body " + str(1)], loc="upper right")
    plt.axvline(0, linestyle = '--')
    plt.axhline(0, linestyle = '--')
    
    return

interact(plot_solution, plot_start=(0,plot_step_size*np.floor(t[0,-1] / plot_step_size),plot_step_size), body=(1,n_bodies));


    

    




plot range: t=  0.0  -  10.0
x_0 =  0.0051138






    

In [5]:

    

y_shift = np.zeros(y.shape)
y_shift[:, 0] = 0
y_shift[:, 1:] = y[:, :-1]

surface_points = np.where(((y_shift) < 0) & (y > 0))
plt.hist2d(x[surface_points], v_x[surface_points],  cmap='Greys', bins=n_bodies, norm=LogNorm())
plt.xlabel(r"$x$")
plt.ylabel(r"$\ddot{x}$")
plt.title(r"Surface of Section --- $y=0$, $\dot{y}>0$ --- $E=-1$, $e=0.3$")

figure_filename = "surface_of_section"
plt.savefig(figure_filename + ".eps")
plt.savefig(figure_filename + ".pdf")