# 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 import matplotlib.dates as mdates 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/" #"/data/etph/kuehl/EPHIN_LEVEL3_FINAL/lvl3_output/" 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: 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]),0,int(minute[date])) for date in range(len(year))] longshort=raw_input("Longterm (>1 year): 'l' ; or Shortterm (<1 year): 's' ") if longshort=="l": mission=True if longshort=="s": mission=False if longshort not in ["s","l"]: print "Wrong Input: can only accept 's' or 'l'!" sys.exit() if longshort=='s': tres=int(raw_input("Timeresolution in minutes (1, 5, 10, 30, 60 or 1440): ")) if tres not in [1, 5, 10, 30, 60, 1440]: print "Wrong Input: can only accept '1', '5', '10', '30', '60' or '1440' minutes!" sys.exit() year=int(raw_input("Year: ")) mysyear,myeyear=year,year mysdoy=raw_input("Start Doy ('mm-dd' allowed as well): ") myedoy=raw_input("Start End ('mm-dd' allowed as well): ") if longshort=='l': tres=int(raw_input("Timeresolution in minutes (60 or 1440): ")) if tres not in [60, 1440]: print "Wrong Input: can only accept '60' or '1440' minutes for longterm plot. Try shortterm for higher time resolution!" sys.exit() mysyear=int(raw_input("Start Year: ")) myeyear=int(raw_input("End Year: ")) mysdoy,myedoy=1,367 try: mysdoy=int(mysdoy) except: mon=int(mysdoy.split("-")[0]) day=int(mysdoy.split("-")[1]) mydate=datetime(int(mysyear),int(mon),int(day)) mysdoy=int(mydate.strftime("%j")) try: myedoy=int(myedoy) except: mon=int(myedoy.split("-")[0]) day=int(myedoy.split("-")[1]) mydate=datetime(int(myeyear),int(mon),int(day)) myedoy=int(mydate.strftime("%j")) 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 = load_level3_l3i(mysyear,1,tres,annual=True,mission=mission,unpack=True) """mask=[(doy>=mysdoy)*(doy<=myedoy)] if longshort=="s": 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=year[mask], month[mask], day[mask], doy[mask], hour[mask], minute[mask], status[mask], accumtime[mask], int_p4[mask], int_p8[mask], int_p25[mask], int_p41[mask], sys_p4[mask], sys_p8[mask], sys_p25[mask], sys_p41[mask], stat_p4[mask], stat_p8[mask], stat_p25[mask], stat_p41[mask], int_h4[mask], int_h8[mask], int_h25[mask], int_h41[mask], sys_h4[mask], sys_h8[mask], sys_h25[mask], sys_h41[mask], stat_h4[mask], stat_h8[mask], stat_h25[mask], stat_h41[mask]""" ringoff=np.zeros(len(status)) for q in range(len(status)): binaries='{0:08b}'.format(int(status[q])) if int(binaries[-2]): ringoff[q]=1 datearray=[datetime(int(year[date]),int(month[date]),int(day[date]),int(hour[date]),int(minute[date])) for date in range(len(year))] # plot variables log=True lfs=10 myncol=4 mycs=0 myhtp=0.1 myds="default" #"steps" myfmt="." #"-" rast=True if True: mycolors=["deepskyblue","mediumslateblue","steelblue","blue","orange","tomato","indianred","red","lightgreen","lawngreen","forestgreen","g","k"] mylabels=["0.25 - 0.70 MeV","0.67 - 3.00 MeV","2.64 - 6.18 MeV","4.80 - 10.40 MeV","4.3 - 7.8 MeV","7.8 - 25 MeV","25 - 40.9 MeV","40.9 - 53 MeV","4.3 - 7.8 MeV/nuc","7.8 - 25 MeV/nuc","25 - 40.9 MeV/nuc","40.9 - 53 MeV/nuc","e$^-$: > 8.70 MeV\np: >53 MeV\nHe: > 53 MeV/nuc"] # first panel (protons) ax1 = plt.subplot2grid((2, 1), (0, 0)) ax1.yaxis.tick_left() plt.plot_date(datearray,int_p4,myfmt,color=mycolors[4],label=mylabels[4],drawstyle=myds,rasterized=rast) plt.plot_date(datearray,int_p8,myfmt,color=mycolors[5],label=mylabels[5],drawstyle=myds,rasterized=rast) plt.plot_date(datearray,int_p25,myfmt,color=mycolors[6],label=mylabels[6],drawstyle=myds,rasterized=rast) plt.plot_date(datearray,int_p41,myfmt,color=mycolors[7],label=mylabels[7],drawstyle=myds,rasterized=rast) plt.ylabel("protons / (cm$^2$ s sr MeV)$^{-1}$") if log==True: plt.yscale("log") plt.legend(fontsize=lfs,ncol=myncol,columnspacing=mycs,handletextpad=myhtp) plt.setp(ax1.get_xticklabels(), visible=False) # do not show tick labels # second panel (helium) ax2 = plt.subplot2grid((2, 1), (1, 0), sharex=ax1) ax2.yaxis.tick_left() #ax2.yaxis.set_label_position("right") #ax2.yaxis.tick_right() plt.plot_date(datearray,int_h4,myfmt,color=mycolors[8],label=mylabels[8],drawstyle=myds,rasterized=rast) plt.plot_date(datearray,int_h8,myfmt,color=mycolors[9],label=mylabels[9],drawstyle=myds,rasterized=rast) plt.plot_date(datearray,int_h25,myfmt,color=mycolors[10],label=mylabels[10],drawstyle=myds,rasterized=rast) plt.plot_date(datearray,int_h41,myfmt,color=mycolors[11],label=mylabels[11],drawstyle=myds,rasterized=rast) plt.ylabel("helium / (cm$^2$ s sr MeV/nuc)$^{-1}$") if log==True: plt.yscale("log") plt.legend(fontsize=lfs,ncol=myncol,columnspacing=mycs,handletextpad=myhtp) #plt.setp(ax2.get_xticklabels(), visible=False) # do not show tick labels myFmt = mdates.DateFormatter('%H:%M \n %d %b \n %Y ') #%H:%M:%S \n #%y %b #ax3.get_xaxis().set_major_locator(dates.HourLocator(byhour=range(0,24,3))) #ax3.get_xaxis().set_minor_locator(dates.HourLocator(byhour=range(0,24,1))) ax2.xaxis.set_major_formatter(myFmt) # fill ring offs for tax in [ax1,ax2]: tax.fill_between(datearray, tax.get_ylim()[0], tax.get_ylim()[1], where=ringoff>0, facecolor='red', interpolate=False,alpha=0.3) plt.xlim(datetime(mysyear,1,1)+timedelta(mysdoy-1),datetime(myeyear,1,1)+timedelta(myedoy-1)) plt.show()