import numpy as np
import matplotlib.pyplot as plt
scipath="/data/missions/soho/costep/level1/sci/"

year=int(raw_input("Year: "))
doy=int(raw_input("DoY: "))

##year,doy=2000,2 #2012,138

if year>=2000:  fname="epi%02d%03d.sci"%(year-2000,doy)
else:  fname="eph%02d%03d.sci"%(year-1900,doy)

data=np.loadtxt(scipath+"%04d/"%year+fname)

yearephin=data[:,0]
dayephin=data[:,1]
millisecephin=data[:,2]
hist=data[:,-40:-8]
histbit=data[:,-8]-64

#print np.where(histbit==0)

hists=[]
sumhists=[]
numhists=[]
for q in range(8):
  hists.append(hist[histbit==q])
  numhists.append(len(hist[histbit==q]))
  sumhists.append(np.sum(hist[histbit==q],axis=0))



plt.figure()  
norm=True
axes=[]
for idx in range(4):
  if idx==0:  axes.append( plt.subplot2grid((4, 1), (idx, 0))   )
  else:  axes.append( plt.subplot2grid((4, 1), (idx, 0), sharex=axes[0])   )
  
  if norm==True:  norm1,norm2=numhists[idx*2],numhists[idx*2+1]
  else: norm1,norm2=1,1
  plt.step(range(32),sumhists[idx*2]/norm1,label="histidx= %i"%(idx*2),lw=2)
  plt.step(range(33,65),sumhists[idx*2+1]/norm2,label="histidx= %i"%(idx*2+1),lw=2)
  if idx==0:  plt.title("%i - %03d"%(year,doy)) 
  plt.legend()
  if norm==True: plt.ylabel("normalized counts")
  if norm==False: plt.ylabel("counts")

plt.show()