apps120811_flux_recalc.mhd_lvalue

''' Calculate L-value, dump to a file, and map it on the surface.

The L-value is calculated following the B-field.
The :func:`get_lvalue` returns the L-value, which is defined here
the fartherest point of the magnetic field line.

The data is saved into two files.

* mhd_lvalue.dat
* ganymede_lvalue.dat.gz

The data is also plotted into a figure.

.. image:: ../../../src/scripts/apps120811_flux_recalc/mhd_lvalue.png

.. autofunction:: get_lvalue
'''

import os
import sys
import logging
logging.basicConfig()
logger = logging.getLogger('mhd_lvalue.py')
import datetime
import math
import gzip
import pickle

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

from pyana.util import streamline
from pyana.pep import mhddata

def get_lvalue(lonr, latr, mf):
    ''' Return the L-value.

    L-value is defined as the "furtherest point along the B-field".

    :param lonr: Longitude in radians.
    :param latr: Latitude in radians.
    :param mf: Magnetic field object.
    :type mf: :class:`pyana.pep.mhddata.MagneticField1201`.
    :returns: L value.
    '''

    ### Bfield function
    def bfield(pos):
        ''' Return the bfield direction in x-z plane
        '''
        bx, by, bz = mf.interpolate3d(pos[0], pos[1], pos[2])
        babs = np.sqrt(bx * bx + by * by + bz * bz)
        return np.array((bx / babs, by / babs, bz / babs))

    l2value = -1

    pos = np.array([np.cos(latr) * np.cos(lonr), np.cos(latr) * np.sin(lonr), np.sin(latr)])

    ### Start tracing.
    tracer = streamline.SimpleTracing(pos, bfield, 0.02)
    for t in range(1000):
        try:
            tracer.step_forward()
            x, y, z = tracer.x, tracer.y, tracer.z
            l2 = x * x + y * y + z * z
            if l2 > l2value:
                l2value = l2
            if l2 < 1:
                break
            
        except IndexError:
            logger.warn('Index error. Stop calculation')
            break

    tracer2 = streamline.SimpleTracing(pos, bfield, -0.02)
    for t in range(1000):
        try:
            tracer2.step_forward()
            x, y, z = tracer2.x, tracer2.y, tracer2.z
            l2 = x * x + y * y + z * z
            if l2 > l2value:
                l2value = l2
            if l2 < 1:
                break
        except IndexError:
            logger.warn('Index error. Stop calculation')
            break

    return np.sqrt(l2value)


def mkdata(filename='mhd_lvalue.dat'):
    '''Main script'''

    mf = mhddata.MagneticField1201()
    mf.set_interpolate_strategy('project')

    ## Start point

    dlat = 1
    dlon = 2

    lats = np.arange(-90, 91, dlat)
    lons = np.arange(0, 361, dlon)
#    lats = np.arange(-40, 41, 10)
#    lons = np.arange(0, 361, 90)

    lvalues = np.zeros([len(lons), len(lats)]) - 1.  # Initialize by -1.

    fp = open(filename, 'w')

    for ilon in range(len(lons)):
        lonr = np.deg2rad(lons[ilon])
        for ilat in range(len(lats)):
            latr = np.deg2rad(lats[ilat])

            lvalues[ilon, ilat] = get_lvalue(lonr, latr, mf)
            print(np.rad2deg(lonr), np.rad2deg(latr), '%.3f' % lvalues[ilon, ilat], file=fp)
            print(np.rad2deg(lonr), np.rad2deg(latr), lvalues[ilon, ilat])

        fp.flush()
    fp.close()

def plotdata(filename='mhd_lvalue.dat'):
    dat = np.loadtxt(filename)   # shape of (N, 3)
    lons = np.unique(dat[:, 0]); lons.sort()
    lats = np.unique(dat[:, 1]); lats.sort()

    lvalues = dat[:, 2].reshape((len(lons), len(lats)))

    lvalues = lvalues.clip(1.1, 8)

    plt.figure()
    cntr = plt.contour(lons, lats, lvalues.T, levels=[1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 3., ])
    plt.clabel(cntr, inline=1, fontsize=10)
#    plt.colorbar(cntr)
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    plt.title('L-value')
    plt.gca().yaxis.set_ticks(np.arange(-90, 91, 10))
    plt.grid()
    
    plt.savefig('mhd_lvalue.png')

def pickledata(filename='mhd_lvalue.dat', outfile='ganymede_lvalue.dat.gz'):
    dat = np.loadtxt(filename)
    lons = np.unique(dat[:, 0]); lons.sort()
    lats = np.unique(dat[:, 1]); lats.sort()
    lvalues = dat[:, 2].reshape((len(lons), len(lats)))
    lvalues = np.ma.masked_outside(lvalues, 0, 5)
    lvalues = lvalues.filled(-1)
    f = gzip.open(outfile, 'w')
    pickle.dump(lons, f)
    pickle.dump(lats, f)
    pickle.dump(lvalues, f)
    f.close()
    

if __name__ == "__main__":
    
    if not os.path.exists('mhd_lvalue.dat'):
        mkdata()
    plotdata()
    pickledata()