pro smooth_wavelet,wv,wid,gaussian=gaussian
dim = [size(/dimen,wv),1,1]
nk = dim[2]
nj = dim[1]
nt = dim[0]
;printdat,dim,'dim'

if keyword_set(gaussian) then begin
  for k=0,nk-1 do $
     for j=0,nj-1 do  if wid[j] gt 1 then begin
       print,k,j,wid[j]
       kernal = exp(-(dgen(wid[j],range=[-2,2]))^2)
       kernal = kernal / total(kernal)
       wv[*,j,k] = convol(wv[*,j,k],kernal,/edge_truncate)
     endif
  return
endif

for k=0,nk-1 do $
   for j=0,nj-1 do $
      if wid[j] gt 1 then wv[*,j,k] = smooth(wv[*,j,k],(wid[j] < (nt-1)),/edge,/nan)
end



pro cross_corr_wavelet,wa,wb,wid,gam,coinc,quad,crat
   p = wa*conj(wb)
   coinc = float(p)
   quad  = imaginary(p)
   smooth_wavelet,coinc,wid
   smooth_wavelet,quad ,wid
   p  = abs(wa)^2
   smooth_wavelet,p,wid
   gam  = abs(wb)^2
   smooth_wavelet,gam,wid
   if arg_present(crat) then  crat = coinc/gam

   p = p*gam
   gam = (coinc^2 + quad^2)/p
   p = sqrt(p)
   coinc = coinc/p
   quad = quad/p
end


function interp_wavelet,wv,time,newtime
dim1 = size(/dimen,wv)
nt = dim1[0]
jv1 = dim1[1]
nk = (n_elements(dim1) eq 3) ? dim1[2] : 1
dim2 = dim1
dim2[0] = n_elements(newtime)
wv2 = make_array(/complex,dimen=dim2,/noz)

for j = 0,jv1-1 do   $
  for k = 0,nk-1 do   $
    wv2[*,j,k] = interpol(wv[*,j,k],time,newtime)

return,wv2
end






function rotate_wavelet2,wv,dir=B,wid=wid,rotmats=rotmats, $
   get_rotmats=get_rotmats

dim = dimen(wv)
nt  = dim[0]
jv  = dim[1]-1
if dim[2] ne 3 then message,'Must have 3 dimensions'

rot = fltarr(dim[0],3,3)
wvp = make_array(dim[0],dim[1],dim[2],/complex,/noz)

message,'Rotating wavelets',/info
time0 = systime(1)

use_rotmat =  keyword_set(rotmats) and (keyword_set(get_rotmats) eq 0)
if arg_present(rotmats) and keyword_set(get_rotmats) then rotmats=fltarr(dim[0],dim[1],3,3)

help,rotmats
help,use_rotmat
help,get_rotmats

for j=0,jv do begin
  if keyword_set(verbose) then print,j
  if use_rotmat then rot=reform(rotmats[*,j,*,*]) $
  else begin
    for k=0,2 do $
       rot[*,k,2] = smooth(B[*,k],wid[j],/edge)
    rot[*,*,2] = rot[*,*,2] / (sqrt(total(rot[*,*,2]^2,2)) # [1,1,1] )
    rot[*,*,1] = crossp3(rot[*,*,2],[1.,0.,0.])
    rot[*,*,1] = rot[*,*,1] / (sqrt(total(rot[*,*,1]^2,2)) # [1,1,1] )
    rot[*,*,0] = crossp3(rot[*,*,1],rot[*,*,2])
    if keyword_set(get_rotmats) then rotmats[*,j,*,*] = rot
  endelse

  wvp[*,j,2] = wv[*,j,0]*rot[*,0,2] + wv[*,j,1]*rot[*,1,2] + wv[*,j,2]*rot[*,2,2]
  wvp[*,j,1] = wv[*,j,0]*rot[*,0,1] + wv[*,j,1]*rot[*,1,1] + wv[*,j,2]*rot[*,2,1]
  wvp[*,j,0] = wv[*,j,0]*rot[*,0,0] + wv[*,j,1]*rot[*,1,0] + wv[*,j,2]*rot[*,2,0]
endfor
print,systime(1)-time0,' Seconds'

return,wvp
end


function rotate_wavelet,wv,B,w0,w1,w2,wid=wid

dim = dimen(wv)
nt  = dim[0]
jv = dim[1]-1
if dim[2] ne 3 then message,'Must have 3 dimensions'
wv = reform(/over,wv,dim[0]*dim[1],dim[2])

if keyword_set(wid)  then begin
message,/info,'Computing Rotation matrices.'
  w2 = fltarr(dim[0],dim[1],dim[2])
  for j=0,jv do  for k=0,2 do  w2[*,j,k] = smooth(B[*,k],wid[j],/edge)

  w2 = reform(w2,/over,dim[0]*dim[1],dim[2])
  w2 = w2/( sqrt(total(w2^2,2)) # [1,1,1] )
  w1 = crossp3(w2,[1.,0.,0.])
  w1 = w1/( sqrt(total(w1^2,2)) # [1,1,1] )
  w0 = crossp3(w1,w2)
endif
message,/info,'Rotating wavelets.'

wvp = make_array(dim[0]*dim[1],dim[2],/complex,/noz)
wvp[*,2] = wv[*,0]*w2[*,0] + wv[*,1]*w2[*,1] + wv[*,2]*w2[*,2]
wvp[*,1] = wv[*,0]*w1[*,0] + wv[*,1]*w1[*,1] + wv[*,2]*w1[*,2]
wvp[*,0] = wv[*,0]*w0[*,0] + wv[*,1]*w0[*,1] + wv[*,2]*w0[*,2]
wvp = reform(/over,wvp,dim[0],dim[1],dim[2])
wv =  reform(/over,wv,dim[0],dim[1],dim[2])
return,wvp

end





function wave_data,varname,frange=frange,param=param,magrat=magrat, $
  trange=tr,data=data,dimennum=dimennum, tplot_prefix=tplot_prefix, $
  verbose=verbose

if size(/type,varname) eq 7  then begin
name=keyword_set(tplot_prefix) ? tplot_prefix : varname

get_data,varname,time,data

;if keyword_set(normname) then normval=data_cut(normname,time)

if n_elements(dimennum) eq 1 then begin
   data=data[*,dimennum]
   name=name+strcompress(/rem,string('(',dimennum,')'))
endif

if keyword_set(rbin) then begin
   n = n_elements(time)
   nrbin = n/rbin
   time = rebin(time[0:nrbin*rbin-1],nrbin)
   data    = rebin(data[0:nrbin*rbin-1,*],nrbin,3)
endif

if not keyword_set(maxpoints) then maxpoints = 32000l ; maxpoints = 2l^16

if n_elements(time) gt maxpoints and not keyword_set(tr) then begin
   Print,'Too many time samples; Select a time range:'
   ctime,tr
endif

if  keyword_set(tr) then trange = tr    ; else get_timespan,trange

if keyword_set(trange) then begin
  if n_elements(trange) eq 1 then begin
     mm = min(abs(time-trange[0]),w)
     rr = [0,maxpoints-1]- maxpoints/2 + w
     rr = 0 > rr < (n_elements(time)-1)
     time=time[rr[0]:rr[1]]
     data=data[rr[0]:rr[1],*]
  endif else begin
     trange=minmax(time_double(trange))
     w = where(time le trange[1] and time ge trange[0],c)
     if c eq 0 then message,'No data in that time range'
     time = time[w]
     data = data[w,*]
  endelse
endif

endif else name = keyword_set(tplot_prefix) ? tplot_prefix : ''

dtime=(time-shift(time,1))[1:*]
dt = average(dtime)
printdat,dt,'dt'
printdat,minmax(dtime),'dt range'
if total(abs(minmax(dtime)/dt-1)) gt .05 then message,'invalid time sampling'


interp_gap,time,data,index=bad_index ,count=nbad_index  ,verbose=verbose   ; remove nans
if 1 then printdat,bad_index,'index'

pad = 2


wv=wavelet2(data,dt,pad=pad,period=period,frange=frange,prange=prange,$
          verbose=verbose,param=param)


data_ptr = {name:name, $
            dt:dt, $
            param: param, $
            period:period, $
            freq:1/period, $
            time: ptr_new(), $
            badtime: ptr_new(), $
            data: ptr_new(), $
            wv: ptr_new()  }

dim = size(/dimen,wv)
nt = dim[0]
jv1 = dim[1]
nk = (n_elements(dim) eq 3) ? dim[2] : 1

badtime = bytarr(nt)
if keyword_set(nbad_index)  then   badtime[bad_index] = 1

if keyword_set(magrat) then begin
message,/info,'Doing magnitude now'
  wvmag= wavelet2(sqrt(total(data^2,2)),dt,pad=pad,period=period,frange=frange,prange=prange,$
          verbose=verbose,param=param)
  str_element,/add,data_ptr,'wvmag',ptr_new(/no_copy,wvmag)
endif

data_ptr.data = ptr_new(/no_copy,data)
data_ptr.time = ptr_new(/no_copy,time)
data_ptr.wv = ptr_new(/no_copy,wv)
data_ptr.badtime = ptr_new(/no_copy,badtime)

return,data_ptr
end





pro wav_tplot,p,avg_period,tplot_prefix=tplot_prefix,rotate_pow=rotate_pow

name = p.name
period = p.period
dt = p.dt
yax = keyword_set(per_axis) ? period : p.freq
ytitle = keyword_set(per_axis) ? 'Seconds' : name+'!cf (Hz)'
mm = minmax(yax)

dim = size(/dimen,*p.wv)
nt = dim[0]
jv1 = dim[1]
nk = (n_elements(dim) eq 3) ? dim[2] : 1

if not keyword_set(avg_period) then avg_period=12.
if 1 then printdat,avg_period,/val,'Num period'

if not keyword_set(wid) then $
   wid = 3 > round(period/2/dt*avg_period)*2+1 < (nt-1)
if 1 then printdat,wid,'width'

if 1 then begin
  mask = *p.badtime # replicate(1.,jv1)
  bad_index = where(*p.badtime,nbad_index)
  if nbad_index ne 0 then  mask[bad_index,*] = 1.
  mask[0,*] = 1.
  mask[nt-1,*] = 1.
  smooth_wavelet,mask,wid  ;,/gaussian
  msk = mask gt .2
  mask = ([1.,!values.f_nan])[msk]
endif else mask=1.

if not keyword_set(normconst) then normconst=1.
;printdat,normconst,'normconst'

if keyword_set(foo) then stop

if keyword_set(normval) then begin
   if n_elements(normval) eq nt then begin
message,/info,'Calculating Normalization'
      normpow = normval # replicate(normconst,jv1)
      smooth_wavelet,normpow,wid
   endif
   if n_elements(normval) eq 1 then normpow=normval
endif else normpow=normconst

if keyword_set(foo) then stop
;printdat,normpow,'Normpow'


tsfx=''

if keyword_set(frac) then begin
  normpow = ( (nk eq 1) ? B^2 : total(B^2,2) ) # replicate(1,jv1)
message,/info,'Fraction Normalization'
  smooth_wavelet,normpow,wid
  normpow = 1/normpow
  tsfx=tsfx+'/<B!u2!n>'
endif

;printdat,normpow,'Normpow'

if keyword_set(kolom) then begin
  normpow = normpow / (replicate(1.,nt) # (kolom*period^(5./3.)) )
  tsfx=tsfx+'/P!dK!n'
  if not keyword_set(zrange) then  zrange = [.1,10]
endif else begin
  if keyword_set(tint_pow) then begin
     normpow = normpow / (replicate(1.,dim[0]) # (period^(tint_pow)) )
     tsfx=tsfx+'*f'
     if tint_pow ne 1 then $
        tsfx=tsfx+'!u'+string(tint_pow,format='(f4.2)')+'!d'
     if not keyword_set(zrange) then  zrange = [.01,1]
  endif
endelse
;printdat,normpow,'Normpow'

for k=0,nk-1 do wv[*,*,k] = wv[*,*,k] * sqrt(normpow)

if keyword_set(wvmag) then wvmag = wvmag * sqrt(normpow)

if nk eq 1 then pow = abs(wv)^2 else pow = total(abs(wv)^2,3)

if not keyword_set(zrange) then  $
    zrange = roundsig(10^(average(alog(pow*mask),/nan)/alog(10)),sigfi=.2) * [.1,10]



ztitle=''
r = keyword_set(resolution) ? resolution : 0
rdtime = reduce_tres(time,r)
polopts = {spec:1,yrange:mm,ylog:1,ystyle:1,no_interp:1,zrange:[-1,1],zlog:0,zstyle:1,ztitle:ztitle,ytitle:ytitle}
powopts = {spec:1,yrange:mm,ylog:1,ystyle:1,no_interp:1,zlog:1,zstyle:1,zrange:zrange,ztitle:ztitle,ytitle:ytitle}

powopts.ztitle='P!dTot!n'+tsfx
store_data,name+'_wv_pow',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts

if nk ne 3 then return     ; scaler quantities...

if keyword_set(magrat) then begin
   powb = abs(wvmag)^2
   ratopts=polopts
   ztitle='P!dmag!n/P!dtot!n'
   pol = powb/pow
   smooth_wavelet,pol,wid    & ztitle = '<'+ztitle+'>'
   ratopts.ztitle=ztitle
   ratopts.zrange=[0,1]
   store_data,name+'_wv_rat_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
endif

if 0 then begin
smooth_wavelet,pow,wid
powopts.ztitle='<P!dTot!n>'+tsfx
store_data,name+'_wvs_pow',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
if keyword_set(magrat) then begin
   smooth_wavelet,powb,wid
   ratopts.ztitle='<P!dmag!n>/<P!dtot!n>'
   pol = powb/pow
   store_data,name+'_wvs_rat_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
endif
endif


if keyword_set(rotate_pow) then begin
;   wv = rotate_wavelet(wv,B,w0,w1,w2,wid= keyword_set(w2) ? 0 :wid)
   wv = rotate_wavelet2(wv,dir=B,rotmat=rotmat,wid=wid,get_rotmat=arg_present(rotmat) and keyword_set(get_rotmat))

   if 1 then message,/info,'Computing power and polarization for TPLOT'

   i=complex(0,1)
   powr = abs(wv[*,*,0]+i*wv[*,*,1])^2
   powl = abs(wv[*,*,0]-i*wv[*,*,1])^2
   powb = abs(wv[*,*,2])^2
;   smooth_wavelet,powb,wid
;   pow = (abs(wv[*,*,0])^2 + abs(wv[*,*,1])^2)
;   smooth_wavelet,pow,wid

   pol = powb/(powb+powl+powr)
   smooth_wavelet,pol,wid
   ratopts=polopts
   ratopts.ztitle='<P!d||!n/P!dtot!n>'
   ratopts.zrange=[0,1]
   store_data,name+'_wv_pol_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts

;   pol = imaginary(wv[*,*,0]*conj(wv[*,*,1]))/pow*2
   pol = (powr-powl)/(powl+powr)
;   polopts.ztitle='!19s!x!dp!n'
;   store_data,name+'_wv_pol2_perp',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=polopts
   smooth_wavelet,pol,wid
   polopts.ztitle='<!19s!x!dp!n>'
   store_data,name+'_wv_pol_perp',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=polopts

   smooth_wavelet,powr,wid
   smooth_wavelet,powl,wid
   smooth_wavelet,powb,wid

if 0 then begin
   pol = powb/(powb+powl+powr)
   ratopts.ztitle='<P!d||!n>/<P!dtot!n>'
   store_data,name+'_wvs_pol_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts

   pol = (powr-powl)/(powr+powl)
   polopts.ztitle='!19s!x!d<p>!n'
   store_data,name+'_wvs_pol_perp',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=polopts
endif

if keyword_set(all) then begin
   i=complex(0,1)
   powopts.ztitle='P!dB!n'+tsfx
   store_data,name+'_wv_pow_par',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=powopts
   powopts.ztitle='P!dperp!n'+tsfx
   store_data,name+'_wv_pow_perp',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
   pow = abs(wv[*,*,0]+i*wv[*,*,1])^2/2
   powopts.ztitle='P!dR!n'+tsfx
   store_data,name+'_wv_pow_r',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
   pow = abs(wv[*,*,0]-i*wv[*,*,1])^2/2
   powopts.ztitle='P!dL!n'+tsfx
   store_data,name+'_wv_pow_l',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
endif

if keyword_set(cross_cor1) then begin
   i=complex(0,1)
   cross_corr_wavelet,wv[*,*,0]+i*wv[*,*,1],wv[*,*,0]-i*wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dl!n'
   store_data,name+'_wv_gam_lin',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dl!n'
   store_data,name+'_wv_coin_lin',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dl!n'
   store_data,name+'_wv_quad_lin',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts

   cross_corr_wavelet,wv[*,*,0],wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dp!n'
   store_data,name+'_wv_gam_cir',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dp!n'
   store_data,name+'_wv_coin_cir',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dp!n'
   store_data,name+'_wv_quad_cir',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts
endif
if keyword_set(cross_cor2) then begin
   cross_corr_wavelet,wv[*,*,2],wv[*,*,0]+i*wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dpr!n'
   store_data,name+'_wv_gam_pr',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dpr!n'
   store_data,name+'_wv_coin_pr',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dpr!n'
   store_data,name+'_wv_quad_pr',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts

   cross_corr_wavelet,wv[*,*,2],wv[*,*,0]-i*wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dpl!n'
   store_data,name+'_wv_gam_pl',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dpl!n'
   store_data,name+'_wv_coin_pl',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dpl!n'
   store_data,name+'_wv_quad_pl',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts

endif


endif
   if keyword_set(stop) then stop


end









pro wav_data,varname, period=period,prange=prange,$
  param=param,avg_period=avg_period, $
  tplot_prefix=tplot_prefix, $
  data=b, wavelet=wv, time=time, mask=msk,$
  magrat=magrat, $
  per_axis=per_axis, $
  kolom=kolom, normval=normval, normconst=normconst,normname=normname, $
  tint_pow = tint_pow, $
  fraction = frac, $
  rotmat=rotmat, get_rotmat=get_rotmat, $
  wid=wid, $
  dimennum=dimennum,rotate_pow=rotate_pow, $
  maxpoints=maxpoints,rbin=rbin,$
  cross1=cross_cor1, cross2=cross_cor2, $
  trange=tr,resolution=resolution,verbose=verbose


if size(/type,varname) eq 7  then begin
name=keyword_set(tplot_prefix) ? tplot_prefix : varname

get_data,varname,time,B

if keyword_set(normname) then normval=data_cut(normname,time)

if n_elements(dimennum) eq 1 then begin
   B=B[*,dimennum]
   name=name+strcompress(/rem,string('(',dimennum,')'))
endif

if keyword_set(rbin) then begin
   n = n_elements(time)
   nrbin = n/rbin
   time = rebin(time[0:nrbin*rbin-1],nrbin)
   B    = rebin(B[0:nrbin*rbin-1,*],nrbin,3)
endif

if not keyword_set(maxpoints) then maxpoints =2l^15; maxpoints =  32000l

if n_elements(time) gt maxpoints and not keyword_set(tr) then begin
   Print,'Too many time samples; Select a time range:'
   ctime,tr
endif

if  keyword_set(tr) then trange = tr    ; else get_timespan,trange

if keyword_set(trange) then begin
  if n_elements(trange) eq 1 then begin
     mm = min(abs(time-trange[0]),w)
     rr = [0,maxpoints-1]- maxpoints/2 + w
     rr = 0 > rr < (n_elements(time)-1)
     time=time[rr[0]:rr[1]]
     B=B[rr[0]:rr[1],*]
     trange= minmax(time)
     printdat,rr,'rr'
  endif else begin
     trange=minmax(time_double(trange))
     w = where(time le trange[1] and time ge trange[0],c)
     if c eq 0 then message,'No data in that time range'
     time = time[w]
     B = B[w,*]
  endelse
endif

;if keyword_set(interp_time) then begin
;  nt = size(/n_elements,interp_time)
;  dim = size(/dimen,b)
;  b2 = replicate(b[1],nt, (n_elements(dim) eq 2) ? dim[1] : 1)
;  for d=0,d2-1 do $
;    b2[*,d] = interp(time,b[*,d],interp_time)
;  b = temporary(b2)
;  help,b
;endif



endif else name = keyword_set(tplot_prefix) ? tplot_prefix : ''


dtime=(time-shift(time,1))[1:*]
dt = average(dtime)
printdat,dt,'dt'
printdat,minmax(dtime),'dt range'
if total(abs(minmax(dtime)/dt-1)) gt .01 then printdat,minmax(dtime),'dt range'


interp_gap,time,B,index=bad_index ,count=nbad_index  ,verbose=verbose   ; remove nans
if 1 then printdat,bad_index,'index'

pad = 2


wv=wavelet2(B,dt,pad=pad,period=period,prange=prange,$
          verbose=verbose,param=param)
dim = size(/dimen,wv)
nt = dim[0]
jv1 = dim[1]
nk = (n_elements(dim) eq 3) ? dim[2] : 1


if keyword_set(magrat) then begin
message,/info,'Doing magnitude now'
  wvmag= wavelet2(sqrt(total(B^2,2)),dt,pad=pad,period=period,prange=prange,$
          verbose=verbose,param=param)
endif


yax = keyword_set(per_axis) ? period : 1/period
ytitle = keyword_set(per_axis) ? 'Seconds' : name+'!cf (Hz)'
mm = minmax(yax)

if not keyword_set(avg_period) then avg_period=12.
if 1 then printdat,avg_period,/val,'Num period'

if not keyword_set(wid) then $
   wid = 3 > round(period/2/dt*avg_period)*2+1 < (nt-1)
if 1 then printdat,wid,'width'

if 1 then begin
  mask = fltarr(nt,jv1)
  if nbad_index ne 0 then  mask[bad_index,*] = 1.
  mask[0,*] = 1.
  mask[nt-1,*] = 1.
  smooth_wavelet,mask,wid  ;,/gaussian
  msk = mask gt .2
  mask = ([1.,!values.f_nan])[msk]
endif else mask=1.

if not keyword_set(normconst) then normconst=1.
;printdat,normconst,'normconst'

if keyword_set(foo) then stop

if keyword_set(normval) then begin
   if n_elements(normval) eq nt then begin
message,/info,'Calculating Normalization'
      normpow = normval # replicate(normconst,jv1)
      smooth_wavelet,normpow,wid
   endif
   if n_elements(normval) eq 1 then normpow=normval
endif else normpow=normconst

if keyword_set(foo) then stop
;printdat,normpow,'Normpow'


tsfx=''

if keyword_set(frac) then begin
  normpow = ( (nk eq 1) ? B^2 : total(B^2,2) ) # replicate(1,jv1)
message,/info,'Fraction Normalization'
  smooth_wavelet,normpow,wid
  normpow = 1/normpow
  tsfx=tsfx+'/<B!u2!n>'
endif

;printdat,normpow,'Normpow'

if keyword_set(kolom) then begin
  normpow = normpow / (replicate(1.,nt) # (kolom*period^(5./3.)) )
  tsfx=tsfx+'/P!dK!n'
  if not keyword_set(zrange) then  zrange = [.1,10]
endif else begin
  if keyword_set(tint_pow) then begin
     normpow = normpow / (replicate(1.,dim[0]) # (period^(tint_pow)) )
     tsfx=tsfx+'*f'
     if tint_pow ne 1 then $
        tsfx=tsfx+'!u'+string(tint_pow,format='(f4.2)')+'!d'
     if not keyword_set(zrange) then  zrange = [.01,1]
  endif
endelse
;printdat,normpow,'Normpow'

for k=0,nk-1 do wv[*,*,k] = wv[*,*,k] * sqrt(normpow)

if keyword_set(wvmag) then wvmag = wvmag * sqrt(normpow)

if nk eq 1 then pow = abs(wv)^2 else pow = total(abs(wv)^2,3)

if not keyword_set(zrange) then  $
    zrange = roundsig(10^(average(alog(pow*mask),/nan)/alog(10)),sigfi=.2) * [.1,10]



ztitle=''
r = keyword_set(resolution) ? resolution : 0
rdtime = reduce_tres(time,r)
polopts = {spec:1,yrange:mm,ylog:1,ystyle:1,no_interp:1,zrange:[-1,1],zlog:0,zstyle:1,ztitle:ztitle,ytitle:ytitle}
powopts = {spec:1,yrange:mm,ylog:1,ystyle:1,no_interp:1,zlog:1,zstyle:1,zrange:zrange,ztitle:ztitle,ytitle:ytitle}

powopts.ztitle='P!dTot!n'+tsfx
store_data,name+'_wv_pow',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts

if nk ne 3 then return     ; scaler quantities...

if keyword_set(magrat) then begin
   powb = abs(wvmag)^2
   ratopts=polopts
   ztitle='P!dmag!n/P!dtot!n'
   pol = powb/pow
   smooth_wavelet,pol,wid    & ztitle = '<'+ztitle+'>'
   ratopts.ztitle=ztitle
   ratopts.zrange=[0,1]
   store_data,name+'_wv_rat_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
endif

if 0 then begin
smooth_wavelet,pow,wid
powopts.ztitle='<P!dTot!n>'+tsfx
store_data,name+'_wvs_pow',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
if keyword_set(magrat) then begin
   smooth_wavelet,powb,wid
   ratopts.ztitle='<P!dmag!n>/<P!dtot!n>'
   pol = powb/pow
   store_data,name+'_wvs_rat_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
endif
endif

   ratopts=polopts
   ratopts.ztitle='<P!d||!n/P!dtot!n>'
   ratopts.zrange=[0,1]

if keyword_set(rotate_pow) then begin
;   wv = rotate_wavelet(wv,B,w0,w1,w2,wid= keyword_set(w2) ? 0 :wid)
   wv = rotate_wavelet2(wv,dir=B,rotmat=rotmat,wid=wid,get_rotmat=arg_present(rotmat) and keyword_set(get_rotmat))

   if 1 then message,/info,'Computing power and polarization for TPLOT'

   i=complex(0,1)
   powr = abs(wv[*,*,0]+i*wv[*,*,1])^2
   powl = abs(wv[*,*,0]-i*wv[*,*,1])^2
   powb = abs(wv[*,*,2])^2
;   smooth_wavelet,powb,wid
;   pow = (abs(wv[*,*,0])^2 + abs(wv[*,*,1])^2)
;   smooth_wavelet,pow,wid

   pol = powb/(powb+powl+powr)
   smooth_wavelet,pol,wid
   store_data,name+'_wv_pol_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts

;   pol = imaginary(wv[*,*,0]*conj(wv[*,*,1]))/pow*2
   pol = (powr-powl)/(powl+powr)
;   polopts.ztitle='!19s!x!dp!n'
;   store_data,name+'_wv_pol2_perp',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=polopts
   smooth_wavelet,pol,wid
   polopts.ztitle='<!19s!x!dp!n>'
   store_data,name+'_wv_pol_perp',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=polopts

   smooth_wavelet,powr,wid
   smooth_wavelet,powl,wid
   smooth_wavelet,powb,wid

if 0 then begin
   pol = powb/(powb+powl+powr)
   ratopts.ztitle='<P!d||!n>/<P!dtot!n>'
   store_data,name+'_wvs_pol_par',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts

   pol = (powr-powl)/(powr+powl)
   polopts.ztitle='!19s!x!d<p>!n'
   store_data,name+'_wvs_pol_perp',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=polopts
endif

if keyword_set(all) then begin
   i=complex(0,1)
   powopts.ztitle='P!dB!n'+tsfx
   store_data,name+'_wv_pow_par',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=powopts
   powopts.ztitle='P!dperp!n'+tsfx
   store_data,name+'_wv_pow_perp',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
   pow = abs(wv[*,*,0]+i*wv[*,*,1])^2/2
   powopts.ztitle='P!dR!n'+tsfx
   store_data,name+'_wv_pow_r',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
   pow = abs(wv[*,*,0]-i*wv[*,*,1])^2/2
   powopts.ztitle='P!dL!n'+tsfx
   store_data,name+'_wv_pow_l',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=powopts
endif


endif

if keyword_set(cross_cor1) then begin
   i=complex(0,1)
   cross_corr_wavelet,wv[*,*,0]+i*wv[*,*,0],wv[*,*,1]+i*wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dl!n'
   store_data,name+'_wv_gam_lin',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dl!n'
   store_data,name+'_wv_coin_lin',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dl!n'
   store_data,name+'_wv_quad_lin',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts

   cross_corr_wavelet,wv[*,*,0],wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dp!n'
   store_data,name+'_wv_gam_cir',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dp!n'
   store_data,name+'_wv_coin_cir',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dp!n'
   store_data,name+'_wv_quad_cir',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts
endif
if keyword_set(cross_cor2) then begin
   cross_corr_wavelet,wv[*,*,2],wv[*,*,0]+i*wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dpr!n'
   store_data,name+'_wv_gam_pr',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dpr!n'
   store_data,name+'_wv_coin_pr',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dpr!n'
   store_data,name+'_wv_quad_pr',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts

   cross_corr_wavelet,wv[*,*,2],wv[*,*,0]-i*wv[*,*,1],wid,pol,pow,powb
   ratopts.ztitle='!19g!x!dpl!n'
   store_data,name+'_wv_gam_pl',data={x:rdtime,y:reduce_tres(pol*mask,r),v:yax},dlim=ratopts
   polopts.ztitle='C!dpl!n'
   store_data,name+'_wv_coin_pl',data={x:rdtime,y:reduce_tres(pow*mask,r),v:yax},dlim=polopts
   polopts.ztitle='Q!dpl!n'
   store_data,name+'_wv_quad_pl',data={x:rdtime,y:reduce_tres(powb*mask,r),v:yax},dlim=polopts

endif

   if keyword_set(stop) then stop


end


