# 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

########################################################################
# load level3 l3i
########################################################################

def load_level3_l3i(year,doy,tres,annual=False,mission=False,unpack=True):
  # tres in minutes, e.g. 1, 5, 10, 30, 60, 1440
  # mission: entire mission (for 60 and 1440 min only), annual: entire year
  path="/data/missions/soho/costep/level3/l3i/"
  if mission==True:
    data=np.loadtxt(path+"%imin/entire_misson_%imin.l3i"%(tres,tres))
  elif annual==True:
    data=np.loadtxt(path+"%imin/%i.l3i"%(tres,year))
  else:
    data=np.loadtxt(path+"%imin/%i/%i_%03d.l3i"%(tres,year,year,doy))
  if unpack==True:
    print  len( np.hsplit(data, 1)[0] )
    year, month, day, doy, hour, minute, status, accumtime, int_p4, int_p8, int_p25, int_p41,  sys_p4, sys_p8, sys_p25, sys_p41,  stat_p4, stat_p8, stat_p25, stat_p41,   int_h4, int_h8, int_h25, int_h41,  sys_h4, sys_h8, sys_h25, sys_h41,  stat_h4, stat_h8, stat_h25, stat_h41 = np.hsplit(data, 1)[0].T
    return year, month, day, doy, hour, minute, status, accumtime, int_p4, int_p8, int_p25, int_p41,  sys_p4, sys_p8, sys_p25, sys_p41,  stat_p4, stat_p8, stat_p25, stat_p41,   int_h4, int_h8, int_h25, int_h41,  sys_h4, sys_h8, sys_h25, sys_h41,  stat_h4, stat_h8, stat_h25, stat_h41
  else:
    return data
  #datearray=[datetime(int(year[date]),int(month[date]),int(day[date]),int(hour[date]),int(minute[date])) for date in range(len(year))]


def year_doy_to_name(year,doy):
  if year<2000: 
    pre="eph"
    year=year-1900
  else: 
    pre="epi"
    year=year-2000
  name="%s%02d%03d" %(pre,year,doy)
  return name

########################################################################
# load level1 single counter
########################################################################

def load_soho_ephin_level1(start_doy,end_doy,year,filename):
  startyear,endyear=year,year
  mean=4  #
  init=1
  for thisyear in range(startyear,endyear+1): 
    sdoy=start_doy
    edoy=end_doy
    for thisdoy in range(sdoy,edoy):
      print "Loading DoY %s of %s..." %(thisdoy,thisyear)
      if thisyear <2000: 
        thisyear2d=thisyear-1900
        prefix='eph'
      else: 
        thisyear2d=thisyear-2000
        prefix='epi'
      try:
        ephindata=np.loadtxt("/data/missions/soho/costep/level1/sci/%s/%s%02d%03d.sci" %(thisyear,prefix, thisyear2d,thisdoy))
        yearephin=ephindata[:,0]
        dayephin=ephindata[:,1]
        millisecephin=ephindata[:,2]
        g=ephindata[:,5]
        a=ephindata[:,6]+ephindata[:,7]+ephindata[:,8]+ephindata[:,9]+ephindata[:,10]+ephindata[:,11]
        b=ephindata[:,12]+ephindata[:,13]+ephindata[:,14]+ephindata[:,15]+ephindata[:,16]+ephindata[:,17]
        c=ephindata[:,18]
        d=ephindata[:,19]
        e=ephindata[:,20]
        f=ephindata[:,21]
        thisoutput=np.array(([yearephin,dayephin,millisecephin,g,a,b,c,d,e,f])).T
        if init==0:
          output=np.vstack((output,thisoutput))
        if init==1:
          output=thisoutput
          init=0
      except: dummy=1
  #np.savetxt("data/%s" %filename, output, fmt='%i')
  return output


########################################################################
# load level2 fluxes
########################################################################
def year_doy_to_name(year,doy):
  if year<2000: 
    pre="eph"
    year=year-1900
  else: 
    pre="epi"
    year=year-2000
  name="%s%02d%03d" %(pre,year,doy)
  return name

myyear,mydoy=2012,20
name=year_doy_to_name(myyear,mydoy)
data=np.loadtxt("/data/missions/soho/costep/level2/rl2/%4d/%s.rl2" %(myyear,name))

year=data[:,0]
doy=data[:,1]
msec=data[:,2]

datearray=[(datetime(int(year[date]),1,1,0)+timedelta(int(doy[date])-1,int(msec[date]/1000.))) for date in range(len(year))]

e150=data[:,6]
e300=data[:,7]
e1300=data[:,8]
#e3000=data[:,9]
p4=data[:,10]
p8=data[:,11]
p25=data[:,12]
#p41=data[:,13]
he4=data[:,14]
he8=data[:,15]
he25=data[:,16]
#he41=data[:,17]
  


########################################################################
# load level2 phas
########################################################################

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))
aseg=data[:,2]
bseg=data[:,3]
coinc=data[:,1]
etot=data[:,10]
ea=data[:,5]
eb=data[:,6]
ec=data[:,7]
ed=data[:,8]
ee=data[:,9]

########################################################################
# load level1 sci coinc counter
########################################################################
def load_level1_sci(year,doy):
  if year <2000: 
    thisyear2d=year-1900
    prefix='eph'
  else: 
    thisyear2d=year-2000
    prefix='epi'
  year,doy,msdoy,e1,e2,e3,e4,p1_1,p1_2,p1_3,p2_1,p2_2,p2_3,p3_1,p3_2,p3_3,p4_1,p4_2,p4_3, h1_1,h1_2,h1_3,h1_4,h2_1,h2_2,h2_3,h2_4,h3_1,h3_2,h3_3,h3_4,h4_1,h4_2,h4_3,h4_4, total_int_counts=np.loadtxt("/data/missions/soho/costep/level1/sci/%s/%s%02d%03i.sci"%(year,prefix,thisyear2d,doy),usecols=(0,1,2,36,37,38,39, 22,23,24, 25,26,27, 41,42,43, 44,45,46, 28,29,30,31, 32,33,34,35, 47,48,49,50, 51,52,53,54, 40),unpack=True)
  p1=p1_1+p1_2+p1_3
  p2=p2_1+p2_2+p2_3
  p3=p3_1+p3_2+p3_3
  p4=p4_1+p4_2+p4_3
  h1=h1_1+h1_2+h1_3+h1_4
  h2=h2_1+h2_2+h2_3+h2_4
  h3=h3_1+h3_2+h3_3+h3_4
  h4=h4_1+h4_2+h4_3+h4_4
  lvl1_counts=[year,doy,msdoy, e1,e2,e3,e4, p1,p2,p3,p4, h1,h2,h3,h4, total_int_counts]
  return lvl1_counts