import numpy as np
import matplotlib.pyplot as plt

### browser to http://server.sepserver.eu/
### go to "greens functions and download data

####################################### protons #################################################

# load downloaded data (proton)

#Center of the pitch-angle bins: 11.25, 33.75, 56.25, 78.75, 101.25, 123.75, 146.25, 168.75 deg 
#E1 (1-3 MeV), E2 (3-6 MeV), E3 (6-11 MeV)
t0,t1,i1,i2,i3,i4,i5,i6,i7,i8=np.loadtxt("data_sepserver_proton.dat",skiprows=53,unpack=True)

datalen=1200
channels=3

parallel=[]
perp=[]
for q in range(channels):
  parallel.append(i1[q*datalen:((q+1)*datalen)])
  perp.append(i8[q*datalen:((q+1)*datalen)])
time=t0[:datalen]

# make a nice plot
plt.figure()
labels=["1-3 MeV","3-6 MeV","6-11 MeV"]
colors=["k","r","b"]

for q in range(len(parallel)):
  plt.plot(time,parallel[q],"-",label=labels[q]+" parallel",color=colors[q])
  plt.plot(time,perp[q],"--",label=labels[q]+" parallel",color=colors[q])
plt.yscale("log")
plt.xlabel("time / minutes")
plt.ylabel("proton intensity / (s sr cm$^2$ MeV)$^{-1}$")
plt.legend(loc=4)
plt.savefig("plot_sepserver_proton.pdf")

# save onset time
f=open("onset_time_proton.dat","w")
f.write("#E_low E_high onset_parallel onset_perp\n")
elows=[1,3,6]
ehighs=[3,6,11]
for q in range(channels):
  # np.argmax gives the index of the first occurence above threshold
  onset_parallel=time[np.argmax(parallel[q]>0)]
  onset_perp=time[np.argmax(perp[q]>0)]
  f.write("%2.2f %2.2f %i %i \n"%(elows[q],ehighs[q],onset_parallel,onset_perp))
f.close()

##################################################################################################

####################################### electrons #################################################

# load downloaded data (electron)

#Center of the pitch-angle bins: 11.25, 33.75, 56.25, 78.75, 101.25, 123.75, 146.25, 168.75 deg
#E1 (24-30 keV), E2 (30-50 keV), E3 (50-82 keV), E4 (82-135 keV), E5 (135-230 keV), E6 (230-400 keV) 
t0,t1,i1,i2,i3,i4,i5,i6,i7,i8=np.loadtxt("data_sepserver_electron.dat",skiprows=53,unpack=True)

datalen=600
channels=6

parallel=[]
perp=[]
for q in range(channels):
  parallel.append(i1[q*datalen:((q+1)*datalen)])
  perp.append(i8[q*datalen:((q+1)*datalen)])
time=t0[:datalen]

# make a nice plot
plt.figure()
labels=["24-30 keV","30-50 keV","50-82 keV", "82-135 keV", "135-230 keV", "230-400 keV"]
colors=["k","r","b","g","c","y"]
for q in range(len(parallel)):
  plt.plot(time,parallel[q],"-",label=labels[q]+" parallel",color=colors[q])
  plt.plot(time,perp[q],"--",label=labels[q]+" parallel",color=colors[q])
plt.yscale("log")
plt.xlabel("time / minutes")
plt.ylabel("electron intensity / (s sr cm$^2$ MeV)$^{-1}$")
plt.legend(loc=4,ncol=2)
plt.savefig("plot_sepserver_electron.pdf")

# save onset time
f=open("onset_time_electron.dat","w")
f.write("#E_low E_high onset_parallel onset_perp\n")
elows=[24,30,50,82,135,230]
ehighs=[30,50,82,135,230,400]
for q in range(channels):
  # np.argmax gives the index of the first occurence above threshold
  onset_parallel=time[np.argmax(parallel[q]>0)]
  onset_perp=time[np.argmax(perp[q]>0)]
  f.write("%2.2f %2.2f %i %i \n"%(elows[q]/1e3,ehighs[q]/1e3,onset_parallel,onset_perp))
f.close()

##################################################################################################