#! /usr/bin/env python """Analyze ACE data for Alfven speed, Elsaesser variables and power spectra Date: December 7, 2008 Author: Bob Wimmer """ from math import pi from numpy import * from numpy.fft import * import numpy import Gnuplot, Gnuplot.funcutils import sys, time from time import time def plot(x,y,desc='plot',line='yes',log=''): #x, y = tsc.fft() data = Gnuplot.Data(x,y,using=(1,2)) #this ensures that t is used as x axis g = Gnuplot.Gnuplot() g('set style data lines') string = 'set title "'+str(desc)+'"' g(string) if line!='yes': string = 'set style data "' + str(line) + '"' g(string) if log =='x': string = 'set logscale x' if log =='y': string = 'set logscale y' if log =='xy': string = 'set logscale xy' if log !='':g(string) g('set ylabel "y-axis [arb. units]"') g('set xlabel "x-axis [arb. units]"') g.plot(data) def plot2(x,y1,y2,title="title",log=''): """Plot the results with Gnuplot.""" data1 = Gnuplot.Data(x,y1,using=(1,2)) #this ensures that t is used as x axis data2 = Gnuplot.Data(x,y2,using=(1,2)) #this ensures that t is used as x axis g = Gnuplot.Gnuplot() if log =='x': string = 'set logscale x' if log =='y': string = 'set logscale y' if log =='xy': string = 'set logscale xy' if log !='':g(string) string = 'set title "'+str(title)+'"' g(string) g('set ylabel "y-axis [arb. units]"') g('set xlabel "x-axis [arb. units]"') g('set style data lines') g.plot(data1,data2) class timeseries(object): """Time series given by x vlaues (e.g. time) and y values (e.g. temperature). Perform various operations on time series wuc as: - treat/replace missing values - smoothing - FFT - and more to come (maybe) Note: requires numpy """ def __init__(self,x,y): if len(x) < 2: print "this is not a time series, but a single point!" if len(y) != len(x): print "len(x) != len(y) -- this needs to be fixed!" self.x = x self.y = y self.beginning = x[0] self.end = x[-1] def max(self): """returns the maximum of x, y""" return self.x.max(), self.y.max() def min(self): """returns the minimum of x, y""" return self.x.min(), self.y.min() def max_arg(self): """returns the maximum of x, y""" arg = self.y.argmax() return arg, self.x[arg], self.y[arg] def min_arg(self): """returns the maximum of x, y""" arg = self.y.argmin() return arg, self.x[arg], self.y[arg] def cut_ind(self, lo,hi): """return subset of time series between lo and hi indices.""" self.xci = self.x[lo:hi] self.yci = self.y[lo:hi] return self.xci, self.yci def cut_x(self, x_lo,x_hi): """return subset of time series between x_lo and x_hi values.""" """find indices of x_lo and x_hi, then call cut_ind""" lo = self.x.searchsorted(x_lo) hi = self.x.searchsorted(x_hi) #print lo, self.x[lo], self.x[lo-1] self.xcx = self.x[lo:hi] self.ycx = self.y[lo:hi] return self.xcx, self.ycx def clean(self): """remove and/or linearly interpolate appropriately missing data values""" dirty = 'true' limit = -9000. inc = 1 if self.y[-1] < limit: self.y[-1] = self.y[-2] if self.y[0] < limit: self.y[0] = self.y[1] while self.y.min() limit: if i < len(self.y)-1: self.y[i] = 0.5*(self.y[i-1] + self.y[i+1]) if i == len(self.y)-1: self.y[i] = self.y[i-1] if i == 0: self.y[i] = self.y[i+1] else: #we have a longer data gap i+=1; inc = 1 while (i+2*inc < len(self.y) and self.y[i+inc]limit and i+inc < len(self.y)): inc /= 2 #print 1, i, inc, i+inc, self.y[i+inc] inc = 1 while (i+2*inc < len(self.y) and self.y[i+inc] freq unit=cycles/day self.fourx = fftfreq(N, d=timestep) self.fourx = fftshift(self.fourx) self.foury = fftshift(self.foury) return self.fourx, self.foury def smooth(self, x,window_len=10,window='hanning'): #print "in smooth method" #print x, window_len, window, len(x) self.s=numpy.r_[2*x[0]-x[window_len:1:-1],x,2*x[-1]-x[-1:-window_len:-1]] #print self.s, len(self.s) #print(len(self.s)) if window == 'flat': #moving average w=ones(window_len,'d') else: w=eval('numpy.'+window+'(window_len)') y=numpy.convolve(w/w.sum(),self.s,mode='same') dy = self.s.copy() dy[window_len-1:-window_len+1] = self.s[window_len-1:-window_len+1] - y[window_len-1:-window_len+1] #print y,dy return y[window_len-1:-window_len+1], dy[window_len-1:-window_len+1] def plot(self,desc='plot',line='yes'): """Plot the results with Gnuplot.""" data = Gnuplot.Data(self.x,self.y,using=(1,2)) #this ensures that t is used as x axis string = 'set title "'+str(desc)+'"' g = Gnuplot.Gnuplot() g(string) g('set ylabel "y-axis [arb. units]"') g('set xlabel "x-axis [arb. units]"') if line=='yes': g('set style data lines') g.plot(data) class elsaesser(timeseries): """Elsaesser variable class zp is Z plus, zm is Z minus, dzp is dZ plus, dzm is dz minus """ def __init__(self,vx,vy,vz,bx,by,bz,vax,vay,vaz,window_length=128,window_type='flat',delta_t=64.): """Initialize Elsaesser variables vx, vy, vz are velocity time series bx, by, bz are B-field time series window length and type are used for smoothing delta t is for normalization of the power spetra, it's 64 seconds for std. ACE data """ vxs, dvxs = vx.smooth(vx.y,window_length,window_type) vys, dvys = vy.smooth(vy.y,window_length,window_type) vzs, dvzs = vz.smooth(vz.y,window_length,window_type) bxs, dbxs = bx.smooth(bx.y,window_length,window_type) bys, dbys = by.smooth(by.y,window_length,window_type) bzs, dbzs = bz.smooth(bz.y,window_length,window_type) #print "now starting work on Alfven speed smoothing..." vaxs, dvaxs = vax.smooth(vax.y,window_length,window_type) vays, dvays = vay.smooth(vay.y,window_length,window_type) vazs, dvazs = vaz.smooth(vaz.y,window_length,window_type) vaxs *= sign(bxs) dvaxs *= sign(bxs) vays *= sign(bxs) dvays *= sign(bxs) vazs *= sign(bxs) dvazs *= sign(bxs) #print "in Elsaesser" #print vx.y, vax.y, vaxs, len(vax.y) #plot2(vx.x,vy.y,vay.y,'vx and vax') #plot2(vx.x,vay.y,vays,'vay and vays') #include sign of bxs for va and vas #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! zpx = vxs + vaxs zpy = vys + vays zpz = vzs + vazs zmx = vxs - vaxs zmy = vys - vays zmz = vzs - vazs dzpx = dvxs + dvaxs dzpy = dvys + dvays dzpz = dvzs + dvazs dzmx = dvxs - dvaxs dzmy = dvys - dvays dzmz = dvzs - dvazs ts_zpx = timeseries(vx.x,zpx) ts_zpy = timeseries(vz.x,zpy) ts_zpz = timeseries(vz.x,zpz) ts_zmx = timeseries(vx.x,zmx) ts_zmy = timeseries(vz.x,zmy) ts_zmz = timeseries(vz.x,zmz) ts_dzpx = timeseries(vx.x,dzpx) ts_dzpy = timeseries(vz.x,dzpy) ts_dzpz = timeseries(vz.x,dzpz) ts_dzmx = timeseries(vx.x,dzmx) ts_dzmy = timeseries(vz.x,dzmy) ts_dzmz = timeseries(vz.x,dzmz) self.zpx = ts_zpx self.zpy = ts_zpy self.zpz = ts_zpz self.zmx = ts_zmx self.zmy = ts_zmy self.zmz = ts_zmz self.dzpx = ts_dzpx self.dzpy = ts_dzpy self.dzpz = ts_dzpz self.dzmx = ts_dzmx self.dzmy = ts_dzmy self.dzmz = ts_dzmz #FFTs of fluctuations in Z+ and Z- (delta Z+, delta -) x,y,z components dzpx_fft_k, dzpx_fft_x = self.dzpx.fft() self.dzpx_fft = timeseries(dzpx_fft_k, dzpx_fft_x) dzpy_fft_k, dzpy_fft_y = self.dzpy.fft() self.dzpy_fft = timeseries(dzpy_fft_k, dzpy_fft_y) dzpz_fft_k, dzpz_fft_z = self.dzpz.fft() self.dzpz_fft = timeseries(dzpz_fft_k, dzpz_fft_z) dzmx_fft_k, dzmx_fft_x = self.dzmx.fft() self.dzmx_fft = timeseries(dzmx_fft_k, dzmx_fft_x) dzmy_fft_k, dzmy_fft_y = self.dzmy.fft() self.dzmy_fft = timeseries(dzmy_fft_k, dzmy_fft_y) dzmz_fft_k, dzmz_fft_z = self.dzmz.fft() self.dzmz_fft = timeseries(dzmz_fft_k, dzmz_fft_z) #Power in the x,y,z components of e+ and e- #First determine scale of x-axis #scale = 86400./(vx.x[len(vx.x)-1] - vx.x[0]) #average x-component of velocity times 64 seconds, the sampling period #scale = -(1./scale - 1./(2.*64.)) scale =1./86400. #This should normalize floating-point day of year to seconds (and thus to Hz in fft x component) pepx = real_if_close(dzpx_fft_x*conjugate(dzpx_fft_x)) nn2 = len(pepx); xx2 = dzpx_fft_k[nn2/2:nn2]*scale; yy2 = 2.*pepx[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.pepx = timeseries(xx2,yy2) #print dzpx_fft_k[nn2/2:nn2] pepy = real_if_close(dzpy_fft_y*conjugate(dzpy_fft_y)) nn2 = len(pepy); xx2 = dzpy_fft_k[nn2/2:nn2]*scale; yy2 = 2.*pepy[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.pepy = timeseries(xx2,yy2) pepz = real_if_close(dzpz_fft_z*conjugate(dzpz_fft_z)) nn2 = len(pepz); xx2 = dzpz_fft_k[nn2/2:nn2]*scale; yy2 = 2.*pepz[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.pepz = timeseries(xx2,yy2) pemx = real_if_close(dzmx_fft_x*conjugate(dzpx_fft_x)) nn2 = len(pemx); xx2 = dzmx_fft_k[nn2/2:nn2]*scale; yy2 = 2.*pemx[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.pemx = timeseries(xx2,yy2) pemy = real_if_close(dzmy_fft_y*conjugate(dzpy_fft_y)) nn2 = len(pemy); xx2 = dzmy_fft_k[nn2/2:nn2]*scale; yy2 = 2.*pemy[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.pemy = timeseries(xx2,yy2) pemz = real_if_close(dzmz_fft_z*conjugate(dzpz_fft_z)) nn2 = len(pemz); xx2 = dzmz_fft_k[nn2/2:nn2]*scale; yy2 = 2.*pemz[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.pemz = timeseries(xx2,yy2) #plot(xx2,yy2,'Power spectral density epx') #calculate power in e+ and e- #print 0.5*self.pepx.y yy2 = 0.5*add(self.pepx.y,self.pepy.y, self.pepz.y) #print yy2, len(yy2),len(self.pepx.y) self.pep = timeseries(self.pepx.x,yy2) pem = -0.5*add(self.pemx.y, self.pemy.y, self.pemz.y) self.pem = timeseries(self.pemx.x,pem) #Now go on to compute coincident and quadrature spectrum, C and Q ciq = conjugate(dzpx_fft_x)*dzmx_fft_x; Cx = real(ciq); Qx = imag(ciq) nn2 = len(Cx); xx2 = dzpx_fft_k[nn2/2:nn2]*scale; yy2 = 2.*Cx[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.Cx = timeseries(xx2,yy2) nn2 = len(Qx); xx2 = dzpx_fft_k[nn2/2:nn2]*scale; yy2 = 2.*Qx[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.Qx = timeseries(xx2,yy2) ciq = conjugate(dzpy_fft_y)*dzmy_fft_y; Cy = real(ciq); Qy = imag(ciq) nn2 = len(Cy); xx2 = dzpy_fft_k[nn2/2:nn2]*scale; yy2 = 2.*Cy[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.Cy = timeseries(xx2,yy2) nn2 = len(Qy); xx2 = dzpy_fft_k[nn2/2:nn2]*scale; yy2 = 2.*Qy[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.Qy = timeseries(xx2,yy2) ciq = conjugate(dzpz_fft_z)*dzmz_fft_z; Cz = real(ciq); Qz = imag(ciq) nn2 = len(Cz); xx2 = dzpz_fft_k[nn2/2:nn2]*scale; yy2 = 2.*Cz[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.Cz = timeseries(xx2,yy2) nn2 = len(Qz); xx2 = dzpz_fft_k[nn2/2:nn2]*scale; yy2 = 2.*Qz[nn2/2:nn2]/nn2 yy2 *= 2*delta_t self.Qz = timeseries(xx2,yy2) #and now the total energy spectrum, cross-helicity spectrum, normalized cross-helicity spectrum, and Elsaesser ratio self.pe_y = 0.5*add(self.pep.y, self.pem.y) self.pe = timeseries(xx2,self.pe_y) self.pec_y = 0.5*add(self.pep.y, -self.pem.y) self.pec = timeseries(xx2,self.pec_y) self.sigma_c_y = self.pec.y/self.pe.y self.sigma_c = timeseries(xx2,self.sigma_c_y) self.r_e_y = self.pem.y/self.pep.y self.r_e = timeseries(xx2,self.r_e_y) #residual energy or symmetric part of the cross-correlation spectrum self.per_y = 0.5*add(self.Cx.y, self.Cy.y, self.Cz.y) self.per = timeseries(xx2,self.per_y) #and the corresponding expression for the imaginary part (symmetric part of the cross-correlation) self.pes_y = 0.5*add(self.Qx.y, self.Qy.y, self.Qz.y) self.pes = timeseries(xx2,self.pes_y) #and their normalized versions self.sigma_r_y = self.per_y/self.pe_y self.sigma_r = timeseries(xx2,self.sigma_r_y) self.sigma_c_y = self.pes_y/self.pe_y self.sigma_c = timeseries(xx2,self.sigma_c_y) #and the Alfven ratio spectrum self.r_a_y = add(1.,self.sigma_r_y)/add(1.,-self.sigma_r_y) self.r_a = timeseries(xx2,self.r_a_y) def main(): # t_start = 101. # t_stop = 151. # t_start = 96. # t_stop = 97. # t_start = 91. # t_stop = 97. t_start = 92.5 t_stop = 95.5 # t_start = 94.5 #outward-propagating waves? # t_stop = 95. t_begin_read = time() infilename = "/home/et085/wimmer/work/ace/2007.dat/2007_sw_cut.dat" #print infilename data = loadtxt(infilename) t_end_read = time() fy = data[:,0] #floating point year fd = data[:,1] #floating point day of year np = data[:,2] #proton density tp = data[:,3] #proton temperature ap = data[:,4] #alpha to proton ratio vp = data[:,5] #proton speed vx = data[:,6] #x-component of proton velocity vy = data[:,7] #y-component vz = data[:,8] #z-component bx = data[:,9] #B_x by = data[:,10] #B_y bz = data[:,11] #B_z bm = data[:,12] #B db = data[:,13] #dB (fluctuation of magnitude) t_end_distrib = time() print "It took "+str(t_end_read - t_begin_read) + " seconds to read the data file." print "It took "+str(t_end_distrib - t_end_read) + " seconds to distribute the data." ts_bm = timeseries(fd,bm) #time series with floating doy as x, B-magnitude as y ts_bm_cx, ts_bm_cy = ts_bm.cut_x(t_start,t_stop) ts_bm_c = timeseries(ts_bm_cx,ts_bm_cy) #print ts_bm_c.min_arg() ts_bm_c.clean() #bm_c_y, bm_c_dy = ts_bm_c.smooth(ts_bm_cy,128,'flat') #ts_bm_cs = timeseries(ts_bm_cx,bm_c_y) #ts_dbm_cs = timeseries(ts_bm_cx,bm_c_dy) #ts_bm_c.plot() #plot2(ts_bm_cs.x,ts_bm_c.y,ts_bm_cs.y,'B_mag, B_mag smoothed') #Now determine density in solar wind ts_np = timeseries(fd,np) ts_np_cx, ts_np_cy = ts_np.cut_x(t_start,t_stop) ts_np_c = timeseries(ts_np_cx, ts_np_cy) ts_np_c.clean() #Now look at fluctuations in Bmag ts_db = timeseries(fd,db) #time series with ts_db_cx, ts_db_cy = ts_db.cut_x(t_start,t_stop) ts_db_c = timeseries(ts_db_cx,ts_db_cy) ts_db_c.clean() #plot2(ts_bm_cs.x,sqrt(ts_dbm_cs.y**2),ts_db_c.y,'my dBmag and dBmag from file') #plot(ts_db_c.y,sqrt(ts_dbm_cs.y**2),'My dB vs. file dB','no') #Now compute Alfven speed V_A = sqrt(B**2/(mu0*rho)) mu_0 = pi*4.e-7 ts_rho_cy = ts_np_c.y*1.67e-21 ts_va_cy = 1.e-12*ts_bm_c.y/sqrt(mu_0*ts_rho_cy) #ts_va = timeseries(fd,va) #ts_va_cx, ts_va_cy = ts_va.cut_x(t_start,t_stop) ts_va_c = timeseries(ts_db_cx, ts_va_cy) ts_va_c.clean() #and now proton speed, vp ts_vp = timeseries(fd,vp) ts_vp_cx, ts_vp_cy = ts_vp.cut_x(t_start,t_stop) ts_vp_c = timeseries(ts_vp_cx, ts_vp_cy) ts_vp_c.clean() ts_vp_c.plot('solar wind speed') #plot2(ts_vp_c.x, ts_vp_c.y, ts_va_c.y,'Vp and Va') #and now Alfvenic Mach number ts_mach_cy = ts_vp_c.y/ts_va_c.y ts_mach_c = timeseries(ts_vp_c.x,ts_mach_cy) ts_mach_c.plot('Alfven Mach number') #Prepare v and b vector time series for Elsaesser analysis ts_vx = timeseries(fd,vx) ts_vx_cx, ts_vx_cy = ts_vx.cut_x(t_start,t_stop) ts_vx_c = timeseries(ts_vx_cx, ts_vx_cy) ts_vx_c.clean() ts_vx_c.y *= -1. #this is just because vx is given positive in the direction towards the Sun (real GSE coordinates) ts_vy = timeseries(fd,vy) ts_vy_cx, ts_vy_cy = ts_vy.cut_x(t_start,t_stop) ts_vy_c = timeseries(ts_vy_cx, ts_vy_cy) ts_vy_c.clean() ts_vz = timeseries(fd,vz) ts_vz_cx, ts_vz_cy = ts_vz.cut_x(t_start,t_stop) ts_vz_c = timeseries(ts_vz_cx, ts_vz_cy) ts_vz_c.clean() ts_bx = timeseries(fd,bx) ts_bx_cx, ts_bx_cy = ts_bx.cut_x(t_start,t_stop) ts_bx_c = timeseries(ts_bx_cx, ts_bx_cy) ts_bx_c.clean() ts_by = timeseries(fd,by) ts_by_cx, ts_by_cy = ts_by.cut_x(t_start,t_stop) ts_by_c = timeseries(ts_by_cx, ts_by_cy) ts_by_c.clean() ts_bz = timeseries(fd,bz) ts_bz_cx, ts_bz_cy = ts_bz.cut_x(t_start,t_stop) ts_bz_c = timeseries(ts_bz_cx, ts_bz_cy) ts_bz_c.clean() #and vector Alfven velocity ts_vax_cy = 1.e-12*ts_bx_c.y/sqrt(mu_0*ts_rho_cy) ts_vax_c = timeseries(ts_db_cx, ts_vax_cy) ts_vax_c.clean() ts_vay_cy = 1.e-12*ts_by_c.y/sqrt(mu_0*ts_rho_cy) ts_vay_c = timeseries(ts_db_cx, ts_vay_cy) ts_vay_c.clean() ts_vaz_cy = 1.e-12*ts_bz_c.y/sqrt(mu_0*ts_rho_cy) ts_vaz_c = timeseries(ts_db_cx, ts_vaz_cy) ts_vaz_c.clean() #print ts_va_c.y, ts_vaz_c.y bx_cs_x, bx_cs_dx = ts_bx_c.smooth(ts_bx_c.y) ts_bx_cs = timeseries(ts_bx_c.x,bx_cs_x) ts_dbx_cs = timeseries(ts_bx_c.x,bx_cs_dx) plot2(ts_bx_cs.x,ts_bx_cs.y,ts_dbx_cs.y,'Bx and dBx smoothed') el = elsaesser(ts_vx_c,ts_vy_c,ts_vz_c,ts_bx_c,ts_by_c,ts_bz_c,ts_vax_c,ts_vay_c,ts_vaz_c,128,'flat') plot2(el.zpx.x,el.zpx.y,ts_vx_c.y,'Elsaesser Z+ x and vx component') plot2(el.zpx.x,el.zmx.y,ts_vx_c.y,'Elsaesser Z- x and vx component') plot2(el.zpx.x,el.zpx.y,el.zmx.y,'Elsaesser Z+ and Z- x component') plot2(el.dzpx.x,el.dzpx.y,el.dzmx.y,'Elsaesser dZ+ and dZ- x component') plot(el.dzpx.x,el.dzpx.y,'Elsaesser dZ+ x component','yes','') el_dzpx_s, el_dzpx_ds = el.dzpx.smooth(el.dzpx.y) plot2(el.dzpx.x,el.dzpx.y,el_dzpx_s,'Elsaesser dZ+ and smoothed') plot2(el.pemx.x,el.pep.y,el.pem.y,'power in e+ and e-','xy') plot(el.r_a.x,el.r_a.y,'Alfven ratio spectrum','yes','x') plot(el.r_e.x, el.r_e.y,'Elsaesser ratio','yes','x') # Now calculate Fourier coefficients raw_input('Please press the return key to continue...\n') main()