textbook example code: plot_hydro_disk_with_bump.pyΒΆ

[source code]

"""
   Visualization for simple N-body integration.
   Reads particle set from file (nbody.hdf5) and prints subsequent frames.
"""
import os
import sys
import numpy

import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot
#from matplotlib import pyplot
from amuse.plot import scatter, xlabel, ylabel
from amuse.lab import *
from amuse.io import store
from optparse import OptionParser
from amuse.units.optparse import OptionParser
from time import sleep
from amuse.ext.orbital_elements import orbital_elements_from_binary

def read_planetary_system(filename): #lim, snapshot_id):
    source = Particles(0)
    gas = Particles(0)
    time = 0 | units.yr
    bodies = read_set_from_file(filename, "amuse")
    i = 0
    for bi in bodies.history:
        i+=1
        if i%2:
            source = bi
        else:
            gas = bi
            break
    print "N=", len(source), len(gas)
    return source, gas

def mu(X = None, Y = 0.25, Z = 0.02, x_ion = 0.1):
    """
    Compute the mean molecular weight in kg (the average weight of particles in a gas)
    X, Y, and Z are the mass fractions of Hydrogen, of Helium, and of metals, respectively.
    x_ion is the ionisation fraction (0 < x_ion < 1), 1 means fully ionised
    """
    if X is None:
        X = 1.0 - Y - Z
    elif abs(X + Y + Z - 1.0) > 1e-6:
        raise Exception("Error in calculating mu: mass fractions do not sum to 1.0")
    return constants.proton_mass / (X*(1.0+x_ion) + Y*(1.0+2.0*x_ion)/4.0 + Z*x_ion/2.0)

def plot_single_image(source, disk):
    lim = 150#|units.AU
    index = 4

    left, width = 0.1, 0.65
    bottom, height = 0.1, 0.65
    bottom_h = left_h = left+width+0.05
    rect_scatter = [left, bottom, width, height]
    rect_histx = [left, bottom_h, width, 0.2]
    rect_histy = [left_h, bottom, 0.2, height]

    fig = pyplot.figure(figsize=(16,14))
    time = disk.get_timestamp()
    xy = pyplot.axes(rect_scatter)
    xz = pyplot.axes(rect_histx)
    yz = pyplot.axes(rect_histy)
    xy.set_xlabel("X [AU]")
    xy.set_ylabel("Y [AU]")
    xz.set_ylabel("Z [AU]")
    yz.set_xlabel("Z [AU]")

    positions = disk.position
    x, y, z = positions.x.value_in(units.AU), positions.y.value_in(units.AU), positions.z.value_in(units.AU)
    vx, vy, vz = disk.vx.value_in(units.kms), disk.vy.value_in(units.kms), disk.vz.value_in(units.kms)
    v = numpy.sqrt(vx**2+vy**2+vz**2)
    v_max = max(v)
    v_min = min(v)

    stars = source
    xs, ys, zs = stars.position.x.value_in(units.AU), stars.position.y.value_in(units.AU), stars.position.z.value_in(units.AU)

#    sizes = 1000
#    sizes = 1000*disk.rho/disk.rho.max()
    alpha = 0.2
    us = disk.u
    xion = disk.xion
    print xion
    u_min, u_max = min(us), max(us)
    xion_min, xion_max = min(xion), max(xion)
    print "u=", u_min, u_max
    Ts = mu() / constants.kB * disk.u
    T_min = max(20|units.K, mu() / constants.kB * u_min)
    T_max = min(1800|units.K, mu() / constants.kB * u_max)
    print "T=", T_min, T_max
    print "X=", xion_min, xion_max

    log_u = numpy.log((us / u_min)) / numpy.log((u_max / u_min))
    clipped_log_u = numpy.minimum(numpy.ones_like(log_u), numpy.maximum(numpy.zeros_like(log_u), log_u))

    log_v = numpy.log((v / v_min)) / numpy.log((v_max / v_min))
    clipped_log_v = numpy.minimum(numpy.ones_like(log_v), numpy.maximum(numpy.zeros_like(log_v), log_v))


    log_T = numpy.log((Ts / T_min)) / numpy.log((T_max / T_min))
    clipped_log_T = numpy.minimum(numpy.ones_like(log_T), numpy.maximum(numpy.zeros_like(log_T), log_T))

    ps = disk.rho
    p_min, p_max = min(ps), max(ps)
    log_p = numpy.log((ps / p_min)) / numpy.log((p_max / p_min))
    clipped_log_p = numpy.minimum(numpy.ones_like(log_p), numpy.maximum(numpy.zeros_like(log_p), log_p))

    blue = clipped_log_T
    red  = clipped_log_p
    green = clipped_log_v
    colors = numpy.transpose(numpy.array([red, green, blue]))

    sizes = 2*2000*disk.h_smooth/disk.h_smooth.max()
    xy.scatter(x, y, sizes, c=colors, edgecolors = "none", alpha = alpha)
    s = 50*stars.mass/stars.mass.min()
    xy.scatter(xs, ys, s, c='b', edgecolors = "none", alpha = 1)

    xy.set_xlim( (-lim, lim) )
    xy.set_ylim( (-lim, lim) )
    sizes = 20*200*disk.h_smooth/disk.h_smooth.max()
    xz.scatter(x, z, sizes, c=colors, edgecolors = "none", alpha = alpha)
    xz.scatter(xs, zs, s, c='b', edgecolors = "none", alpha = 1)

    yz.scatter(z, y, sizes, c=colors, edgecolors = "none", alpha = alpha)
    yz.scatter(zs, ys, s, c='b', edgecolors = "none", alpha = 1)

    xz.set_xlim( xy.get_xlim() )
    yz.set_ylim( xy.get_xlim() )
    yz.set_xlim( (-0.2*lim, 0.2*lim) )
    xz.set_ylim( (-0.2*lim, 0.2*lim) )
#    yz.set_xlim( (-lim, lim) )
#    xz.set_ylim( (-lim, lim) )

    filename = "fig_hydro_disk_with_bump_i{0:04}".format(index)
    fig.savefig(filename)
#    pyplot.show()

def main(filename=None): #, lim=-1|units.AU, image_id=-1):

    output_single_image(filename)
    """
    if filename:
        output_single_image(lim, image_id)
    else:
        output_multiple_images(lim)
    """

def calculate_orbital_elements(star, planet):
    from amuse.ext.orbital_elements import orbital_elements_from_binary

    p = Particles()
    p.add_particle(star)
    p.add_particle(planet)
    M, m, a, e, ta_out, inc, lan_out, aop_out = orbital_elements_from_binary(p, G=constants.G)
    #print "Orbital elements:", M, m, a, e, ta_out, inc, lan_out, aop_out
    return a, e, inc

def output_single_image(filename):
    source, disk = read_planetary_system(filename)
    print disk[0].position.in_(units.AU)
    snapshot_id = 1
    lim = 150|units.AU
    plot_single_image(source, disk)

def output_multiple_images(lim):
    snapshot_id = 0
    filename = "planetary_system_i{0:04}.amuse".format(snapshot_id)
    while os.path.exists(filename):
        bodies = read_set_from_file(filename, "amuse")
        for bi in bodies.history:
            if len(bi)<=20:
                planets = bi.copy()
                time = bi.get_timestamp()
                print "Orbits at t=", time, planets.semimajor_axis.in_(units.AU), planets.eccentricity
            else:
                disk = bi.copy()

                time = bi.get_timestamp()
                print "Snapshot=", snapshot_id, time
                if image_id<0 or image_id == snapshot_id:
                    plot_single_image(planets, disk)
                if image_id == snapshot_id:
                    print "Stop plotting"
                    break
        snapshot_id += 1
        filename = "planetary_system_i{0:04}.amuse".format(snapshot_id)

def new_option_parser():
    result = OptionParser()
    result.add_option("-f", dest="filename", default = "hydro_disk_with_bump_i0004.amuse",
                      help="output filename [%default]")
    return result

if __name__ in ('__main__', '__plot__'):
    o, arguments  = new_option_parser().parse_args()
    main(**o.__dict__)

Keywords: python, amuse, astrophysics, matplotlib, pylab, example, codex (see how-to-search-examples)