function morlet2,t,param
if not keyword_set(param) then param = 1d
i = dcomplex(0,1)
mm = !dpi^.25 * exp(-t^2/2/param^2) * exp(-i*2*!dpi*t)
return,mm
end


n = 2l^12
t = dgen(n,range=[-6,6.])
dt = total(minmax(t)*[-1,1])/(n-1)
 
mm = morlet2(t,param)
plot,t,mm,xtitle='Time domain'
oplot,t,imaginary(mm),/color,line=2

fmm = fft(double(mm))
fmm = fft(conj(mm))
fmm = fft2(conj(mm),dt,f)
fmm = fft2(mm,dt,f)
fmm = fft2(mm,dt,f,/shif)

plot,f,fmm,xrange = [-20,20] ,/xsty
plot,f,abs(fmm),xrange = [-20,20] ,/xsty
oplot,f,imaginary(fmm),/color




plot,fft(mm)


help,mm

plot,x,mm


!p.multi=[0,2,2]
!x.margin=[6,2]
t = dgen(n,range=[-6,6.])
sigma = 1.d
s = 1.
mm = morlet2(t/s,sigma)
plot,t,mm
oplot,t,imaginary(mm),col=2,line=2

ws=s * dgen(n,range=[-2,2]*2*!pi)
sigmaw = 1/sigma/s
fmm = !dpi^.25 *sigma*s*exp(-(ws-2*!pi)^2/2/sigmaw^2)
plot,ws/2/!pi,fmm

sigma = 2.d
s = 1.
mm = morlet2(t/s,sigma)
plot,t,mm
oplot,t,imaginary(mm),col=2,line=2

ws=s * dgen(n,range=[-2,2]*2*!pi)
sigmaw = 1/sigma/s
fmm = !dpi^.25 *sigma*s*exp(-(ws-2*!pi)^2/2/sigmaw^2)
plot,ws/2/!pi,fmm



popen,/port,'morlet',ysize=7.
pclose,/xv






.r chirp


popen,/land,xsize=9,ysize=7,'chirp'

!p.multi=[0,1,2]
!x.margin=10
plot,t,x1,ytitle='Amplitude',xtitle="Time"
ylim,lim,.001,100,1
options,lim,no_interp=1,ytitle='Frequency',xtitle="Time",ztitle='Power'
;box,lim
scale = replicate(1,n_elements(ind)) # period
specplot,t[ind],1/period,abs(w)^2/scale,lim=lim
oplot,t[ind],1/(2*!pi/omega(t[ind]))
oplot,t[ind],1/coi
pclose,/xv

x0=x1
w0=w
x1=x1+shift(x1,20000)/2
w = wavelet(x1[ind],dt*skip,pad=2,period=period,/ver,coi=coi,param=param)




popen,/land,'wvlt_meth'

!p.multi = [0,2,3]
plot,t,x1,ytitle='Amplitude',xtitle="Time"

fx = fft2(x1,dt,f)  & f =f/2/!dpi
;plot,f,sqrt(abs(fx)^2),/ylog,/xlog,xrange=[1,n]/dt
plot,f,fx,/xlog,ytitle='Amplitude',xrange=minmax(/pos,f),xtitle='Frequency',yrange=[-1,1]*max(abs(fx))

s = 4.
ws = f*2 *!dpi*s 
fmm = !dpi^.25 *sigma*s*exp(-(ws-2*!pi)^2/2/sigmaw^2)
!p.multi[0] =3
plot,f,fmm,/xlog,xrange=minmax(/pos,f)

!p.multi[0] =1
plot,f,fmm*fx,/xlog,xrange=minmax(/pos,f),yrange=[-1,1]*max(abs(fmm*fx))

!p.multi[0] =2
wft=fft(/inv,fmm*fx)
plot,t,float(wft)
;oplot,t,imaginary(wft),col=2

pclose,/xv




popen,/land,'wvlt_meth2'

!p.multi=[0,1,3]
xrange=[100,200]

plot,t,x1,ytitle='Amplitude',xtitle="Time",xrange=xrange
j=69
print,period[j]
plot,t,w[*,j],xrange=xrange
oplot,t,imaginary(w[*,j]),col=2
plot,t,abs(w[*,j])^2,xrange=xrange

pclose,/xv




alfven_pow,/verb,save=wvsave
s = wvsave
ctime,t
ind = where(*s.time lt t[1] and *s.time gt t[0])
wv = (*s.wv)[ind,*,*]
wb = (*s.wb)[ind,*,*]
period = *s.period
ts = (*s.time)[ind]

i = complex(0,1)
 vl = wv[*,*,0] - i * wv[*,*,1]
 vr = wv[*,*,0] + i * wv[*,*,1]
 bl = wb[*,*,0] - i * wb[*,*,1]
 br = wb[*,*,0] + i * wb[*,*,1]

k = 0
j = 15
xlim,lim,-2,2
ylim,lim,-2,2

print,period(j)
g = conj(wv)*wb/abs(wb)^2
box,lim
oplot,float(g[*,j,k]),imaginary(g[*,j,k]),ps=3

u = wb[*,j,0]
v = wb[*,j,1]

.r
gaa = replicate(i,jv1)
para = replicate(par,jv1)
for j=0,jv1-1 do begin
b = br[*,j]
v = vr[*,j]

;g = conj(b)*v/abs(b)/abs(v)
;g = conj(b)*v
g = conj(b)*v/abs(b)^2
if plot then  box,lim
if plot then  oplot,g,g/i,ps=3
ga = average(g,std=sga)
;ga = average(g *abs(b) *abs(v))/average(abs(v))/average(abs(b))
;ga = average(g *abs(b)^2 *abs(v)^2)/average(abs(v)^2)/average(abs(b)^2)
plots,ga,ga/i,ps=4,col=4
gaa[j]=ga
if plot then  plots,ga,ga/i,ps=4,col=6

if plot then  box,lim
if plot then  oplot,b,v,ps=3
if plot then  oplot,dgen(),dgen()
fit,float(b),float(v),nam='a0 a1',func='polycurve',param=par,over=plot,p_sig=sig
para[j]=par
print,ga,abs(ga)

endfor
end

plot,period,[-1,1],/xlog,/nodata
oplot,period,para.a1 
oplot,period,gaa,/col

plot,[-2,2],[-2,2],/nod,/iso
oplot,gaa,gaa/i


ga = average(g *abs(b) *abs(v))/average(abs(v))/average(abs(b))
plots,ga,ga/i,ps=4,col=4

ga = average(g *abs(b)^2 *abs(v)^2)/average(abs(v)^2)/average(abs(b)^2)
plots,ga,ga/i,ps=4,col=3

.run
!p.multi=[0,1,3]
j = 20
for j=0,jv-1 do begin
!p.multi[0] = 0
b = br[*,j]
plot,b,yrange =[-1,1]*max(abs(b))
oplot,b/i,/col
b = br[*,j]*exp(-i*(ts-ts[0])/period[j]/1.03*2*!dpi)
plot,b,yrange =[-1,1]*max(abs(b))
oplot,b/i,/col
oplot,abs(b),col=3
dbdt = deriv(ts,b)
plot,dbdt/b,yrange=max(abs(dbdt/b))*[-1,1]
oplot,(dbdt/b)/i,/col
stop ;wait,.1
endfor
end

