cena_espec_3orbit

Plot for Figure for three consective energy spectra.

#!/usr/bin/env python

''' Plot for Figure for three consective energy spectra.

'''
import matplotlib
matplotlib.use('Agg')  # If  you want to plot without showing window, enable it.

import getopt

#### Main routine starts.
import sys
import math
import datetime

import matplotlib.pyplot as plt
import numpy as np
from numpy import ma
import scipy.optimize

from irfpy.util import utc
from irfpy.util import julday
from irfpy.util import bipower
import irfpy.cy1orb.Cy1Orbit
from irfpy.cena import energy
import irfpy.cena.cena_mass2 as cenam
import irfpy.cena.cena_flux

def main(start=2019, ch=3, dt=600):
    ''' Just plot the energy spectra of given channel.
    '''
    cy1orb = irfpy.cy1orb.Cy1Orbit.Cy1Orbit()

    orbs = start + np.arange(3)
    # Get the time.
    perlist = [julday.Julday(cy1orb.getDayEqOfOrb(orb)[0]).getDatetime() for orb in orbs]
    print(perlist)    

    deltat = datetime.timedelta(seconds=dt)
    t0list = [per - deltat / 2 for per in perlist]
    t1list = [per + deltat / 2 for per in perlist]

    # Energy level
    etbl= energy.getEnergyE16()

    # One count level
    cnt2flx = irfpy.cena.cena_flux.Count2Flux()
    ocl = np.array([cnt2flx.getFlux(1, ie, plotch) for ie in range(16)])

    # Plotting frame
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.plot(etbl, ocl)

    # Read the data
    symbols = ['k--+', 'k-+', 'k-.+']

    for i in range(3):
        orb = orbs[i]
        t0 = t0list[i]
        t1 = t1list[i]

        print(orb, t0, t1)

        # This is the data of one orbit.
        datalist = []
        tlist = cenam.getobstime(timerange=[t0, t1])

        for t in tlist:
            try:
#                d = cenam.getHdefluxE16(t)   # 16x7 JdObject
                c = cenam.getdataE16(t)      # 16x7x128 JdObject
            except RuntimeError as e:
                print(('Skipped non-supported data %s' % utc.convert(t, datetime.datetime)), file=sys.stderr)
                continue
#            datalist.append(d.getData())
            datalist.append(c.getData())

        if len(datalist) == 0:
            raise RuntimeError('Failed to load the data in orbit %d.  Abort.' % orb)

        datalist = np.ma.array(datalist)  # N x 16x7x128 data.
        datalen = datalist.shape[0]
        datalist = datalist[:, :, :, 0:64].sum(3) # Sum over light mass data. Nx16x7 data
        datalist = datalist[:, :, plotch]   # Only select the plotting channel.  Nx16 array
        datalist = datalist.sum(0)    # Average over the time.  16 array.

        # Errorbar from statistics.

        # Cnt -> Flux
        countmax = (datalist + np.sqrt(datalist)) / datalen
        countmin = (datalist - np.sqrt(datalist)) / datalen
        countave = (datalist) / datalen

#        print countmax
#        print countave
#        print countmin

        fluxave = ma.array([cnt2flx.getFlux(countave[ie], ie, plotch) for ie in range(16)])
        fluxmax = ma.array([cnt2flx.getFlux(countmax[ie], ie, plotch) for ie in range(16)]) - fluxave
        fluxmin = fluxave - ma.array([cnt2flx.getFlux(countmin[ie], ie, plotch) for ie in range(16)]) - 1e-5

        print('ave', fluxave)
        print('min', fluxmin)

        # Get the SV table.
        if datalist.mask[2] == False:
            sv = 1
        elif datalist.mask[-2] == False:
            sv = 3
        else:
            sv = 2

        ax.errorbar(etbl, fluxave, np.array([fluxmin, fluxmax]), fmt=symbols[sv-1], label='#%04d (SV=%d; acc=%d)' % (orb, sv, datalen),
                        elinewidth=0.5)
        
        if sv == 2:   # Fitting.
            enesv2 = etbl.compress(~fluxave.mask)
            flxsv2 = fluxave.compressed()
            fitprm, scode = bipower.fitting(enesv2, flxsv2)
            print(fitprm)
            k0, r0, k1, r1 = fitprm
            fitfnc = bipower.mkfunc(k0, r0, k1, r1)
            enesmooth = np.logspace(0, 4)
            fitflx = [fitfnc(e) for e in enesmooth]
            ax.plot(enesmooth, fitflx, 'k:', label='Fit: (%.1f %.2f) (%.1f %.2f)' % (math.log10(k0), r0, math.log10(k1), r1))
            

    ax.plot(etbl, ocl/datalen, 'b:', label='1 Lvl (%d acc.)' % datalen)
    # Shaping.
    ax.set_ylim(1e-1, 1e7)
    ax.set_ylabel('Flux (H) [#/cm2 sr s eV]')
    ax.set_yscale('log')

    ax.set_xlim(5, 5000)
    ax.set_xlabel('Energy [eV]')
    ax.set_xscale('log')

    ax.legend(loc='lower left')

    fig.savefig('three_conseq_%04d.eps' % start)

    return (fluxave, fluxmax, fluxmin)


def help():
    print("USAGE:  %s  -s START_ORB -t TIME_SEC -c CHANNEL" % sys.argv[0], file=sys.stderr)

if __name__ == "__main__":

    try:
        opts, args = getopt.getopt(sys.argv[1:], "s:c:t:")
    except:
        help()
        raise

    start = 2019
    plotch = 3
    dt = 30 * 60    # in seconds

    for opt, optarg in opts:
        if "-s" == opt:
            start = int(optarg)
        if "-c" == opt:
            plotch = int(optarg)
        if "-t" == opt:
            dt = float(optarg)
        if "-h" == opt:
            help()
            

    (flx, flxupp, flxdwn) = main(start=start, ch=plotch, dt=dt)