import sys,os
import datetime
from irfpy.swim.ace import *
from irfpy.cy1orb import *
import numpy.ma as ma
import matplotlib
# matplotlib.use('Agg')
from pylab import *
a=AceMoon()
def plot(t0, t1, fig=None):
if fig==None:
fig=figure(figsize=(8,12))
val=a.getList(t0,t1)
print(val.size())
t=[]
np=[]
vx=[]
vy=[]
vz=[]
bx=[]
by=[]
bz=[]
tp=[]
he=[]
for data in val:
prm=data.getData()
_t=prm.getJulday().getDatetime()
_np=prm.getNp()
_vx=prm.getV_gse_x()
_vy=prm.getV_gse_y()
_vz=prm.getV_gse_z()
_bx=prm.getB_gse_x()
_by=prm.getB_gse_y()
_bz=prm.getB_gse_z()
_tp=prm.getTp()
_he=prm.getAlpha_ratio()
t.append(_t)
np.append(_np)
vx.append(_vx)
vy.append(_vy)
vz.append(_vz)
bx.append(_bx)
by.append(_by)
bz.append(_bz)
tp.append(_tp)
he.append(_he)
avg = a.getAverage(t0,t1)
print(avg)
sn = fig.add_subplot(5,1,2)
sn.plot(ma.masked_where(np<0, t), ma.masked_where(np<0, np))
if avg[0]!=float('nan'):
print(avg[0], avg[1])
eb=sn.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[0],avg[0],avg[0],avg[0]], [avg[1],avg[1],avg[1],avg[1]], ecolor='b', color='b', linestyle='--')
sn.text(t0+3*(t1-t0)//4, 0.05, 'N: %.2f+-%.2f'%(avg[0], avg[1]), color='b', rotation=10)
sn.set_yscale('log')
sn.set_ylim(0.05, 30)
sn.set_ylabel('Np [/cc]')
# Plot of B field
sb = fig.add_subplot(5,1,1, sharex=sn)
sb.plot(t, bx, 'b-', t, by, 'g-', t, bz, 'r-')
if avg[12]!=float('nan'):
print(avg[12], avg[13])
eb=sb.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[12],avg[12],avg[12],avg[12]], [avg[13],avg[13],avg[13],avg[13]], ecolor='b', color='b', linestyle='--')
sb.text(t0+(t1-t0)//4, -20, 'Bx: %.1f+-%.1f'%(avg[12], avg[13]), color='b', rotation=10)
if avg[14]!=float('nan'):
print(avg[14], avg[15])
eb=sb.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[14],avg[14],avg[14],avg[14]], [avg[15],avg[15],avg[15],avg[15]], ecolor='g', color='g', linestyle='--')
sb.text(t0+(t1-t0)//2, -20, 'By: %.1f+-%.1f'%(avg[14], avg[15]), color='g', rotation=10)
if avg[16]!=float('nan'):
print(avg[16], avg[17])
eb=sb.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[16],avg[16],avg[16],avg[16]], [avg[17],avg[17],avg[17],avg[17]], ecolor='r', color='r', linestyle='--')
sb.text(t0+3*(t1-t0)//4, -20, 'Bz: %.1f+-%.1f'%(avg[16], avg[17]), color='r', rotation=10)
sb.set_ylim(-20,20)
sb.set_ylabel('B [nT]')
# Plot of velocities
sv = fig.add_subplot(5,1,3, sharex=sn)
sv.plot(t, vx, 'b-', t, vy, 'g-', t, vz, 'r-')
if avg[2]!=float('nan'):
print(avg[2], avg[3])
eb=sv.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[2],avg[2],avg[2],avg[2]], [avg[3],avg[3],avg[3],avg[3]], ecolor='b', color='b', linestyle='--')
sv.text(t0+(t1-t0)//4, -600, 'Vx: %.1f+-%.1f'%(avg[2], avg[3]), color='b', rotation=10)
if avg[4]!=float('nan'):
print(avg[4], avg[5])
eb=sv.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[4],avg[4],avg[4],avg[4]], [avg[5],avg[5],avg[5],avg[5]], ecolor='g', color='g', linestyle='--')
sv.text(t0+(t1-t0)//2, -600, 'Vy: %.1f+-%.1f'%(avg[4], avg[5]), color='g', rotation=10)
if avg[6]!=float('nan'):
print(avg[6], avg[7])
eb=sv.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[6],avg[6],avg[6],avg[6]], [avg[7],avg[7],avg[7],avg[7]], ecolor='r', color='r', linestyle='--')
sv.text(t0+3*(t1-t0)//4, -600, 'Vz: %.1f+-%.1f'%(avg[6], avg[7]), color='r', rotation=10)
sv.set_ylim(-600, 200)
sv.set_ylabel('V [km/s]')
# Plot of temperature
st = fig.add_subplot(5,1,5, sharex=sn)
st.plot(ma.masked_where(tp<0, t), ma.masked_where(tp<0, tp))
if avg[8]!=float('nan'):
print(avg[8], avg[9])
eb=st.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[8],avg[8],avg[8],avg[8]], [avg[9],avg[9],avg[9],avg[9]], ecolor='b', color='b', linestyle='--')
st.text(t0+(t1-t0)//2, 10000, 'T: %.4e+-%.4e'%(avg[8], avg[9]), color='b', rotation=10)
st.set_yscale('log')
st.set_ylim(10000, 500000)
st.set_ylabel('T [K]')
# Plot of H-He ratio
sr = fig.add_subplot(5,1,4, sharex=sn)
if max(he)>0:
sr.plot(ma.masked_where(he<0, t), ma.masked_where(he<0, he))
if avg[10]!=float('nan'):
print(avg[10], avg[11])
eb=sr.errorbar([t0,t0+(t1-t0)//4, t0+3*(t1-t0)//4,t1], [avg[10],avg[10],avg[10],avg[10]], [avg[11],avg[11],avg[11],avg[11]], ecolor='b', color='b', linestyle='--')
sr.text(t0+2*(t1-t0)//4, 1e-5, 'He/H: %.4f+-%.4f'%(avg[10], avg[11]), color='b', rotation=10)
else:
sr.plot( (t[0], t[-1]), (2e10,2e10) ) # Plot unrealistic value
print('Since Alpha ratio is not available, no plot is made. %f'%max(he))
sr.set_yscale('log')
sr.set_ylim(1e-5, 1e-0)
sr.set_ylabel('N(He)/N(H)')
# Formatting the x-axis
formatter=matplotlib.dates.DateFormatter('%FT%T')
sr.xaxis.set_major_formatter(formatter)
locator=matplotlib.dates.AutoDateLocator()
sr.xaxis.set_major_locator(locator)
fig.autofmt_xdate()
return (fig, (sn, sv, st, sr))
if __name__=='__main__':
orb=Cy1OrbitNr.Cy1OrbitNr()
orb.setFromDefaultUri()
for nr in range(645,3160):
t0=datetime.datetime.utcfromtimestamp(orb.getStartTime(nr))
t1=datetime.datetime.utcfromtimestamp(orb.getStopTime(nr))
print(nr, t0, t1)
fig, ss = plot(t0, t1)
fig.savefig('ace_%04d.png'%nr)