# This scripts is doing...
# Patrick Kuehl, date
import time
from pylab import *
from datetime import datetime
import matplotlib.pyplot as plt
import numpy as np
import numpy.ma as ma
from matplotlib.colors import LogNorm
from scipy.optimize import leastsq
from matplotlib import dates
import time
import datetime as date
from datetime import datetime
from datetime import timedelta
import os
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning) 
from matplotlib.backends.backend_pdf import PdfPages


thisyear=int(raw_input("Year: "))
thisdoy=int(raw_input("DoY: "))

if thisyear <2000: 
  thisyear2d=thisyear-1900
  prefix='eph'
else: 
  thisyear2d=thisyear-2000
  prefix='epi'

data=np.loadtxt("/data/missions/soho/costep/level2/pha/%s/%s%02d%03d.pl2" %(thisyear,prefix, thisyear2d,thisdoy))


# time
time_bool=int(raw_input("Entire day (0) or limited time period (1): "))
if time_bool==1:
  time= data[:,0]-data[0,0]
  starthour=int(raw_input("Starthour: "))
  startmin=int(raw_input("Startmin: "))
  endhour=int(raw_input("Endhour: "))
  endmin=int(raw_input("Endmin: "))
  starttime=((starthour*60)+startmin)*60000.
  endtime=((endhour*60)+endmin)*60000.
  data=data[time>starttime,:]
  time=time[time>starttime]
  data=data[time<endtime,:]

# coincidence
aseg=data[:,2]
bseg=data[:,3]
coinc=data[:,1]
thiscoinc=raw_input("Designated Coincidence (-1 for all,-2 for stopping, -3 for e-, -4 for p, -5 for he): ")
if thiscoinc=="-1":
  data=data
elif thiscoinc=="-2":
  data=data[coinc<12]
elif thiscoinc=="-3":
  data=data[coinc<4]
elif thiscoinc=="-4":
  data=data[(coinc>=4)&(coinc<8)]
elif thiscoinc=="-5":
  data=data[(coinc>=8)&(coinc<12)]
else:
  data=data[coinc==float(thiscoinc)]


# values to plot
data=data[data[:,5]>0]
thisy=int(raw_input("Hist Data: a=0, b=1, c=2, d=3, e=4, Etot=5, Etot*a=6, Etot/a=7 "))
labels=["Energy in A / MeV","Energy in B / MeV","Energy in C / MeV","Energy in D / MeV","Energy in E / MeV","Total Energy / MeV", "Total Energy * Energy in A / MeV**2","Total Energy / Energy in A"]

if thisy==6:  ydata=data[:,10]*data[:,5]
elif thisy==7:  ydata=data[:,10]/data[:,5]
else: ydata=data[:,5+thisy]


data=ydata

#print data
# plot stuff
bins=np.logspace(np.log10(np.min(data[data>0])),np.log10(np.max(data)),101)
#print bins
hist,binedges=np.histogram(data,bins)
bincenter=np.ones(len(binedges)-1)
for i in range(len(binedges)-1):
  bincenter[i]=0.5*(binedges[i]+binedges[i+1])
plt.figure()
plt.step(bincenter,hist,where='mid')
plt.xlabel(labels[thisy])
plt.ylabel("counts")

plt.yscale("log")
plt.xscale("log")
#plt.ylim([max(min(ydata),1e-8),max(ydata)])
#plt.xlim([max(min(xdata),1e-8),max(xdata)])

plt.ylabel("#")
coinclabel=thiscoinc
if coinclabel=="-1": coinclabel="all"
if coinclabel=="-2": coinclabel="stopping"
if coinclabel=="-3": coinclabel="electron"
if coinclabel=="-4": coinclabel="proton"
if coinclabel=="-5": coinclabel="helium"
plt.title("year %s - doy %s : coinc=%s" %(thisyear,thisdoy,coinclabel) )


#import pickle
#pickle.dump(plt.gca(), file('myplot.pickle', 'w'))
plt.show()