IPN map plotting script

make_ipn_plot.py

from numpy import sin, cos, arcsin, arctan2, deg2rad, rad2deg
import numpy as np

import matplotlib as mpl
mpl.use('Agg')
mpl.rcParams['xtick.direction'] = 'in'
mpl.rcParams['ytick.direction'] = 'in'
mpl.rcParams['xtick.top'] = True
mpl.rcParams['ytick.right'] = True


import matplotlib.pyplot as pl 
from matplotlib.path import Path
import matplotlib.patches as patches
from matplotlib.ticker import  MultipleLocator 
from mpl_toolkits.axes_grid.anchored_artists import AnchoredText

class Annulus:

    def __init__(self, str_annulus, n_step):
    
        lst_annulus_pars = [float(x) for x in str_annulus.split()]
        lst_annulus_pars.append(n_step)

        #print lst_annulus_pars
    
        self.arr_c, self.arr_p, self.arr_m = self.get_annulus(*lst_annulus_pars)
    
    def cartesian_to_spherical(self, r):
    
        Rabs = np.linalg.norm(r)
        fDec = arcsin(r[2] / Rabs)
        fRA = arctan2(r[1], r[0])
        
        if (fRA < 0): 
            fRA += 2 * np.pi;
    
        return fRA, fDec
    
    def get_small_circle(self, fRA_c, fDec_c, fR, nStep):
    
        # all angles in degrees!
        
        ec = np.zeros(3) 
        e0 = np.zeros(3)
        e1 = np.zeros(3) 
        
        fRA = deg2rad(fRA_c)
        fDec = deg2rad(fDec_c)
        fTheta = deg2rad(fR)
        
        ec[0] = cos(fDec) * cos(fRA) # x
        ec[1] = cos(fDec) * sin(fRA) # y
        ec[2] = sin(fDec) #z
        
        e0[0] = -sin(fDec) * cos(fRA)
        e0[1] = -sin(fDec) * sin(fRA)
        e0[2] = cos(fDec)
    
        e1 = np.cross(ec, e0)
        
        et = np.zeros(3)  
        el = np.zeros(3) 
        e = np.zeros(3)
    
        c = 0
        
        arr_annulus_data = np.zeros((nStep,2))
    
        for i in xrange(nStep):
        
            fFi = 2.0 * np.pi * i / nStep

            s1 = e0 * cos(fFi)
            s2 = e1 * sin(fFi)
    
            el = ec * cos(fTheta)
            et = s1 + s2
            et = et * sin(fTheta)
            e = et + el
            
            fRA_Ann, fDec_Ann = self.cartesian_to_spherical(e)

            #print fRA_Ann, fDec_Ann
            
            arr_annulus_data[i,0] = rad2deg(fRA_Ann)
            arr_annulus_data[i,1] = rad2deg(fDec_Ann)   
            
        return arr_annulus_data

    def get_annulus(self, fRA_c, fDec_c, fR, fR_minus, fR_plus, nStep):
    
        arr_c = self.get_small_circle(fRA_c, fDec_c, fR, nStep)
        arr_p = self.get_small_circle(fRA_c, fDec_c, fR+fR_plus, nStep)
        arr_m = self.get_small_circle(fRA_c, fDec_c, fR+fR_minus, nStep)
    
        return arr_c, arr_p, arr_m

    def write(self, file_name):

        f = open(file_name, 'w')

        for i in xrange(self.arr_c[:,0].size):
            f.write("{:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f} {:8.3f}\n".format(
                self.arr_c[i,0], self.arr_c[i,1], self.arr_p[i,0],  self.arr_p[i,1], self.arr_m[i,0], self.arr_m[i,1]))

        f.close()

class Plot:

    def __init__(self):

        # fonts
        pl.rc('font',family='serif')
        pl.rc('font',serif='Arial')
        pl.rc('font',size=18)
        pl.rc('mathtext',default='regular')

        self.f, self.axes = pl.subplots(edgecolor='w', facecolor='w',figsize=(9, 9))

        self.axes.set_xlabel(r'$\alpha$ (deg)')
        self.axes.set_ylabel(r'$\delta$ (deg)')

    def plot_annulus(self, annulus, color_):

        self.axes.plot(annulus.arr_c[:,0], annulus.arr_c[:,1], color=color_, linestyle='--', linewidth=0.5)
        self.axes.plot(annulus.arr_m[:,0], annulus.arr_m[:,1], color=color_, linestyle='-', linewidth=0.5)
        self.axes.plot(annulus.arr_p[:,0], annulus.arr_p[:,1], color=color_, linestyle='-', linewidth=0.5)

    def plot_polygon(self, lst_verts):

        lst_verts.append((0.0, 0.0))

        codes = [Path.MOVETO,
            Path.LINETO,
            Path.LINETO,
            Path.LINETO,
            Path.CLOSEPOLY,
        ]

        path = Path(lst_verts, codes)
        patch = patches.PathPatch(path, fill=None, edgecolor='k', lw=1)
        self.axes.add_patch(patch)

    def save(self, file_name, zoom=None):

        self.axes.invert_xaxis()

        if zoom:
            self.axes.set_xlim(zoom[0][0], zoom[0][1])
            self.axes.set_ylim(zoom[1][0], zoom[1][1])

        self.f.savefig(file_name)

def main():

    figure_name = 'GRB190501A_IPN_map_.png'

    # Wind (Konus) - INTEGRAL (SPI-ACS) annulus
    str_annulus = '202.0038 -3.5143 69.1117 -0.5121 +0.5121'
    ann1 = Annulus(str_annulus, 2*360)
    ann1.write('KW_INT_annulus.txt')
    
    # Wind (Konus) - Mars-Odyssey (HEND) annulus
    str_annulus = '79.3230 24.1523 70.7404 -0.0651 +0.0650'
    ann2 = Annulus(str_annulus, 2*360)
    ann2.write('KW_HEND_annulus.txt')

    zoom = [(178, 168),(60, 66)]

    pl = Plot()
    pl.plot_annulus(ann1, 'r')
    pl.plot_annulus(ann2, 'b')
    pl.save(figure_name, zoom)

main()