appsUnsorted.plot_ganymede_tiltedfield

''' A simple code to plot Ganymede B-field at the specified height
'''

import pyana.pep.ganymede_magfield
import matplotlib.pyplot as plt
import numpy as np

if __name__ == "__main__":
    main()

    # Instance the B-field model.
    tf = pyana.pep.ganymede_magfield.TiltedDipole()

    # The longitude and lattitude grid, every 1 degrees.
    lon = np.linspace(0, 2 * np.pi, 360)
    clat = np.linspace(0, np.pi, 180)

    # Ganymede radius
    r = 2631.

    # Altitude of S/C
    h = 200.

    # S/C radius from the center
    rh = r + h

    # Grid on the spacecraft position
    lons, lats = np.meshgrid(lon, clat)

    # x, y and z in Cartesian.
    x = rh * np.sin(lats) * np.cos(lons)
    y = rh * np.sin(lats) * np.sin(lons)
    z = rh * np.cos(lats)

    ### If you want to plot the grided coordinates.
    #plt.figure()
    #plt.pcolor(lons, lats, x)
    #plt.figure()
    #plt.pcolor(lons, lats, y)
    #plt.figure()
    #plt.pcolor(lons, lats, z)

    #Shape of x is (len(clat), len(lon))
    bx = np.zeros_like(x)
    by = np.zeros_like(x)
    bz = np.zeros_like(x)
    bt = np.zeros_like(x)

    # Iterate to calculate the b-field in specific position.
    for ilat in range(len(clat)):
        for ilon in range(len(lon)):
            # Spacecraft position for each grid
            posvec = np.array([x[ilat, ilon], y[ilat, ilon], z[ilat, ilon]])
            # Here, getting the vector
            bvector = tf.get_field(posvec)

            # Save each component.
            bx[ilat, ilon] = bvector[0]
            by[ilat, ilon] = bvector[1]
            bz[ilat, ilon] = bvector[2]

            bt[ilat, ilon] = np.sqrt((bvector ** 2).sum())

    plt.figure()
    plt.contour(lons * 180 / np.pi, 90 - lats * 180 / np.pi, bx)
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.title('Bx [nT]')
    plt.colorbar()

    plt.figure()
    plt.contour(lons * 180 / np.pi, 90 - lats * 180 / np.pi, by)
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.title('By [nT]')
    plt.colorbar()

    plt.figure()
    plt.contour(lons * 180 / np.pi, 90 - lats * 180 / np.pi, bz)
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.title('Bz [nT]')
    plt.colorbar()


    plt.figure()
    plt.contour(lons * 180 / np.pi, 90 - lats * 180 / np.pi, bt)
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.title('Bt [nT]')
    plt.colorbar()