swim_momfit1d

A script program to calculate SWIM moment by 1-D maxwellian fitting.

''' A script program to calculate SWIM moment by 1-D maxwellian fitting.
'''

import matplotlib
#matplotlib.use('Agg')   # To prevent using GUI.

from pylab import *
import os,sys,logging
#logging.basicConfig(level=logging.DEBUG)
from irfpy.swim import *
from irfpy.util import maxwell
import datetime
import dateutil
import numpy

if __name__ == '__main__':

    if len(sys.argv) < 2:
        print('USAGE:   %s   time'%sys.argv[0], file=sys.stderr)
        sys.exit(0)


    t0=dateutil.parser.parse(sys.argv[1])   ### This is the equator crossing.

    # Obtain VDF
    data=numpy.array(swim_start.getvdf(t0).getData())
    idx=data.argmax()     # argmax considers 16x16 matrix as 256 array.  To know which deflection is the maximum, you need to take a residual of 16
    idx_def = idx%16
    print(idx_def)   # It may be 12

    vdf = data[:,idx_def]

    print(vdf)

    # Energy and velocity

    e=SwimEnergy.EnergyTable.table()
    v=SwimEnergy.eV2kms(e.getTable())
    print(v)

    plot(v, vdf, label='Observed flux')

    n,vs,vt = swim_start.simple_moment(v, vdf)

    print(n, vs, vt)

    # Fitting by maxwellian

    fitfunc = maxwell.mkslice(n, vs*1e3, vt*1e3)
    v0 = arange(100,800,10)
    plot(v0, fitfunc(v0*1e3), label='Maxwell Fitting')

    legend()

    text(600, 1e-10, 'n=%.3f [/cc]\nvsw=%.1f [km/s]\nvth=%.1f [km/s]'%(n*1e-6, vs, vt))
    title('SWIM D=%d @ %s'%(idx_def, t0))
    xlabel('Velocity [km/s]')
    ylabel('VDF [s^3/m^6]')
    yscale('log')
    ylim(1e-13,1e-8)

    plt.show()