apps120911_iotorus.app01_iotorus_viz

''' Visualize the io torus in 2D (slice)
'''

import numpy as np
import matplotlib.pyplot as plt

import logging

from pyana.pep.iotorus import SimpleNeutralCloudModel, TwoPeakModel

def torus_viz_xy(torus, ax, z=0):
    x = np.linspace(-10, 10, 200)
    y = np.linspace(-10, 10, 200)

    xx, yy = np.meshgrid(x, y)
    density = np.zeros_like(xx)

    for ix, x0 in enumerate(x):
        logging.debug('Exec x=%f' % x0)
        for iy, y0 in enumerate(y):
            density[iy, ix] = torus.get_density(x0, y0, z)

    cf = ax.contour(xx, yy, density, levels=[1, 2, 3, 5, 10, 25, 50])
    plt.clabel(cf)

def torus_viz_x_colz(torus, ax, y=0):
    x = np.linspace(-10, 10, 200)
    z = np.linspace(-4, 4, 200)

    density = np.zeros_like(x)

    for ix, x0 in enumerate(x):
        for iz, z0 in enumerate(z):
            density[ix] += (torus.get_density(x0, y, z0) * (0.04 * 7e9))

    ax.plot(x, density)
    ax.set_yscale('log')
    ax.set_ylim(1e7, 1e14)

def torus_viz_xz(torus, ax, y=0):
    x = np.linspace(-10, 10, 200)
    z = np.linspace(-10, 10, 200)

    xx, zz = np.meshgrid(x, z)
    density = np.zeros_like(xx)

    for ix, x0 in enumerate(x):
        for iz, z0 in enumerate(z):
            density[iz, ix] = torus.get_density(x0, y, z0)

    cf = ax.contour(xx, zz, density, levels=[1, 2, 3, 5, 10, 25, 50])
    ax.set_xlim(5.4, 6.4)
    ax.set_ylim(-0.5, 0.5)
    plt.clabel(cf)


def torus_viz_xy_col(torus, ax, z=0):
    x = np.linspace(-7, 7, 70)
    y = np.linspace(-7, 7, 70)
    z = np.linspace(-2, 2, 100)
    logging.getLogger().setLevel(logging.DEBUG)

    xx, yy = np.meshgrid(x, y)
    cdensity = np.zeros_like(xx)

    for ix, x0 in enumerate(x):
        logging.debug('Exec x=%f' % x0)
        for iy, y0 in enumerate(y):
            for z0 in z:
                cdensity[iy, ix] += (torus.get_density(x0, y0, z0) * (0.04 * 7e9))
        logging.debug('   Max = %g' % cdensity[:, ix].max())

    cf = plt.contour(xx, yy, np.log10(cdensity), levels=(9, 10, 11, 12, 13))
    plt.clabel(cf)


def main():
    '''Main script'''

#    torus = SimpleNeutralCloudModel()
    torus = TwoPeakModel()

    fig = plt.figure()
    ax1 = fig.add_subplot(221)
    ax2 = fig.add_subplot(222)
    ax3 = fig.add_subplot(223)
    ax4 = fig.add_subplot(224)

    torus_viz_xy(torus, ax1)
    torus_viz_xz(torus, ax2)
    torus_viz_x_colz(torus, ax3)
    torus_viz_xy_col(torus, ax4)
#    torus_viz_r_col(torus, ax4)


    plt.savefig('app01_iotorus_viz.png')    
    plt.savefig('app01_iotorus_viz.eps')    

if __name__ == "__main__":
    main()