function random_field_1D, kmin, kmax, Npoints, VRMS=vrms, SEED=seed
if not keyword_set(seed) then seed=1
if not keyword_set(vrms) then vrms=1.
; Define the number of points and the interval:  
T = 1./FLOAT(Npoints)  
; Midpoint+1 is the most negative frequency subscript:  
N21 = Npoints/2 + 1  
; Wave numbers
K = FINDGEN(Npoints)/(T*Npoints)
K[N21] = N21 -Npoints + FINDGEN(N21-2)  
;select range of wave numbers 
ind= where(  (K ge kmin) and ((K) le kmax))
Nmodes = N_ELEMENTS(ind)
print, 'setting up ', Nmodes, ' modes.'
; flat power spectrum. I.e. all modes from gaussian with same amplitude
amplitude = RANDOMN(seed, Nmodes ,/NORMAL)
phase     = RANDOMN(seed, Nmodes, /UNIFORM)*2.*!Pi
f = complexarr(Npoints)
f[ind] = complex( amplitude*cos(phase), amplitude*sin(phase) )
; random from reverse fourier transform
v = real_part(FFT(f, -1))
; normalize to correct rms value
v = vrms * v/sqrt(mean(v^2))
return, v
end

function random_field_2D, kmin, kmax, Npoints, VRMS=vrms, SEED=seed
if not keyword_set(seed) then seed=1
if not keyword_set(vrms) then vrms=1.
; Define the number of points and the interval:  
T = 1./FLOAT(Npoints)  
; Midpoint+1 is the most negative frequency subscript:  
N21 = Npoints/2 + 1  
; Wave numbers
oK = FINDGEN(Npoints)/(T*Npoints)
oK[N21] = N21 -Npoints + FINDGEN(N21-2)
K = FLTARR(Npoints, Npoints)
for i=0, Npoints-1 do $
   for j=0, Npoints-1 do $
      K[i,j] = sqrt(oK[i]^2+oK[j]^2)
   
;select range of wave numbers 
ind= where( (K ge kmin) and (K le kmax))
Nmodes = N_ELEMENTS(ind)
print, 'setting up ', Nmodes, ' modes.'
;  rising (lineaerly with k) power spectrum. 
; I.e. all modes have same amplitude from Gaussian 
; so for the larger k with more modes there is more power 
amplitude = RANDOMN(seed, Nmodes ,/NORMAL)
phase     = RANDOMN(seed, Nmodes, /UNIFORM)*2.*!Pi
f = complexarr(Npoints, Npoints)
f[ind] = complex( amplitude*cos(phase), amplitude*sin(phase) )
; random from reverse fourier transform
v = real_part(FFT(f, -1))
; normalize to correct rms value
v = vrms * v/sqrt(mean(v^2))
return, v
end

function random_field_3D, kmin, kmax, Npoints, VRMS=vrms, SEED=seed
if not keyword_set(seed) then seed=1
if not keyword_set(vrms) then vrms=1.
; Define the number of points and the interval:  
T = 1./FLOAT(Npoints)  
; Midpoint+1 is the most negative frequency subscript:  
N21 = Npoints/2 + 1  
; Wave numbers
oK = FINDGEN(Npoints)/(T*Npoints)
oK[N21] = N21 -Npoints + FINDGEN(N21-2)
K = FLTARR(Npoints, Npoints, Npoints)
for i=0, Npoints-1 do $
   for j=0, Npoints-1 do $
      for kk=0, Npoints-1 do $
      K[i,j,kk] = sqrt(oK[i]^2+oK[j]^2+oK[kk]^2)
; k^2 rising power spectrum. I.e. all modes have same amplitude from
; Gaussian 
; So for the large k with more modes there is more power.
;select range of wave numbers 
ind= where( (K ge kmin) and (K le kmax))
; delete this 3D array
K = 0.
Nmodes = N_ELEMENTS(ind)
print, 'setting up ', Nmodes, ' modes.'

amplitude = RANDOMN(seed, Nmodes ,/NORMAL)
phase     = RANDOMN(seed, Nmodes, /UNIFORM)*2.*!Pi
f = complexarr(Npoints, Npoints, Npoints)
f[ind] = complex( amplitude*cos(phase), amplitude*sin(phase) )
; random from reverse fourier transform
v = REAL_PART(FFT(f, -1))
; normalize to correct rms value
v /= sqrt(mean(v^2))
v *= vrms 

return, v
end

function PowerSpectrum_1D, array, obin=obin
  npoints = N_elements(array)
  ps = abs(FFT(array,1))^2/npoints^2
  if ARG_PRESENT(obin) then begin
     T = 1./FLOAT(Npoints)  
; Midpoint+1 is the most negative frequency subscript:  
     N21 = Npoints/2 + 1  
; Wave numbers
     OBIN = FINDGEN(N21)/(T*Npoints)
;     OBIN[N21] = N21 -Npoints + FINDGEN(N21-2)  
  endif
; power spectrum of real variable has same poer in positive and
; negative indeces. Hence multiply by two. 
; consequentl fullfills discrete parcevals theorem 
; sum |hk|^2 = 1/N  sum  |Hk|^2   (Numerical recipes 12.1)
; not for some reason in IDL the inverse transform = the transform 
; in Numerical recipes.... so they're switched and we have to
; use inverse above. 
  return, 2.*ps[0:N21-1]
end

function PowerSpectrum_2D, array, obin=obin
; powerspectrum of real 2D data
  npoints = (SIZE(array))[1]  ; assume its square ... 
  N = N_elements(array)
  N21 = Npoints/2 + 1  
  T = 1./FLOAT(Npoints)  
  oK = FINDGEN(Npoints)/(T*Npoints)
  oK[N21] = N21 -Npoints + FINDGEN(N21-2)
  K = FLTARR(Npoints, Npoints)
  for i=0, Npoints-1 do $
;   for j=0, Npoints-1 do $
;      K[i,j] = sqrt(oK[i]^2+oK[j]^2)
     K[i,*] = sqrt(oK^2+oK[i]^2)
  hist = histogram(K, reverse_indices=r, max=sqrt(2.)*(N21-1))
  psf = abs((FFT(array, 1)))^2/DOUBLE(n)^2
  print, 'total:', total(psf)
  ps = DBLARR(N_elements(hist)-2)
  obin = ps
  for i=0,N_elements(hist)-3 do begin
     if hist[i] gt 0 then begin
        ii = [r[r[i]:r[i+1]-1]]
; done such that int ps dk = |array^2| 
        ps[i] = total(psf[ii])
        obin[i] = mean(K[ii])
     endif
  endfor

  return, ps
end


function PowerSpectrum_3D, array, obin=obin
; powerspectrum of real 3D data
; look at http://arxiv.org/abs/0704.3851
; Kritsuk et al 2007
; Their equations 16 thorugh 19 show what we do here. 
; We do not have the factor 1/2 which is specific for energy however.
  npoints = (SIZE(array))[1]  ; assume its cubical ... 
  N = N_elements(array)
  N21 = Npoints/2 + 1  
  T = 1./FLOAT(Npoints)  
  oK = FINDGEN(Npoints)/(T*Npoints)
  oK[N21] = N21 -Npoints + FINDGEN(N21-2)
  K = FLTARR(Npoints, Npoints, Npoints)
  for i=0, Npoints-1 do $
   for j=0, Npoints-1 do $
     K[i,j,*] = sqrt(oK^2+oK[i]^2+oK[j]^2)
  hist = histogram(K, reverse_indices=r, max=sqrt(3.)*(N21-1))
; powerspectrum 
  psf = abs((FFT(array, 1)))^2/DOUBLE(n)^2
; binned such that Sum 
  ps = DBLARR(N_elements(hist))
  obin = ps

  for i=0L,N_elements(hist)-1 do begin
     if hist[i] gt 0 then begin
        ii = [r[r[i]:r[i+1]-1]]
        obin[i] = mean(K[ii])
        ps[i] = total(psf[ii], /double)
     endif
  endfor
  obin[0] = (obin[0]+obin[1])/2 ; bin center
  return, ps
end


;; ; code to test the power-spectrum routines uncomment and run with
;; ; .r powerspectrum.pro on idl command line

;; Np = 128
;; a = random_field_1d( 4, 16, Np, seed=3)
;; ps = PowerSpectrum_1D(a, obin=obin)
;; plot, obin, ps, /ylog
;; oplot, obin, ps, psym=6

;; b = random_field_2d( 0, 3, Np, seed=0)+10
;; bb = abs(fft(b,1))
;; ps2 = PowerSpectrum_2D(b, obin=obin2)
;; plot, obin2, ps2 > 1e-7, /ylog
;; oplot, obin2, ps2 , psym=4
;; print, total(ps2*2.*!DPI*obin2)+ps2[0], mean(b^2)
;; ;tvscl, alog10(shift((bb > 1e-5), np/2,np/2))

;; ;breaki

;; c = random_field_3d( 10, 180, Np, seed=0)
;; cc = abs(fft(c,1))
;; ps3 = PowerSpectrum_3D(c, obin=obin3)
;; plot, obin3, ps3 > 1e-12, /ylog
;; oplot, obin3, ps3 >1e-12 , psym=4
;; print, total(ps3*4.*!DPI*obin3^2)+ps[0], mean(c^2)

;; end

;; Aake Nordlunds power-spectrum functions:

FUNCTION aver,a,axis,array=array,sum=sum
;+
; FUNCTION aver,a
;  Calculate the average.  Avoid accumulation roundoff by
;  using the /double key word.
;-
  sz=size(a)
  if keyword_set(axis) then begin
    if max(axis) gt sz[0] then begin                                            ; guard against last index having been dropped
      w = where(axis le sz[0], nw)                                              ; ignore axes which are not there
      if nw gt 0 then axis=axis(w) else return, a                               ; if no axis left we are done
    end
  end
  ;if keyword_set(cm) then a=arg(3:s(1)-5,*,*)
  s=size(a)

  if n_elements(axis) le 0 then begin
    ;a=reform(a)
    n=n_elements(a)
    if keyword_set(sum) then av=total(a,/double) else av=total(a,/double)/n
  end else if n_elements(axis) eq 1 then begin
    if keyword_set(sum) then begin
      av=total(a,axis,/double)
    end else begin
      av=total(a,axis,/double)/product(s[axis])
    end
  end else if n_elements(axis) eq 2 then begin
    s=size(a)
    axs = axis(sort(axis))
    if keyword_set(sum) then $
      return,total(total(a,axs(1)),axs(0)) $
    else begin
      av = total(a,axs[1])/product(sz[axis])
      sz = size(av)
      if axs[0] le sz[0] then av = total(av,axs[0])
      return, av
    end
  end else begin
    print,'cannot do more than 2 axes yet'
    return,0
  end
  if s(s(0)+1) ne 5 then av=float(av)

  if keyword_set(array) then begin
    if keyword_set(axis) then begin
      if s(0) eq 2 then begin
        if axis eq 1 then av=rebin(reform(av,1,s(2)),s(1),s(2))
        if axis eq 2 then av=rebin(reform(av,s(1),1),s(1),s(2))
      end else if s(0) eq 3 then begin
        if axis eq 1 then av=rebin(reform(av,1,s(2),s(3)),s(1),s(2),s(3))
        if axis eq 2 then av=rebin(reform(av,s(1),1,s(3)),s(1),s(2),s(3))
        if axis eq 3 then av=rebin(reform(av,s(1),s(2),1),s(1),s(2),s(3))
      end
    end else begin
      av=rebin(reform([av],1,1,1),s(1),s(2),s(3))
    end
  endif

  return,av
END


; $Id: power3d.pro,v 1.23 2007/11/26 10:28:14 aake Exp $
;+------------------------------------------------------------------------------
 PRO power3d,a,b,c,oplt=oplt,title=title,xtitle=xtitle,ytitle=ytitle,$
  size=size,spectrum=spectrum,wavenumbers=wns,fit=fit,afit=afit,$
  average=average,corners=corners,debug=debug,kindex=kindex,$
  compensate=compensate,offset=o,noplot=noplot,potsdam=potsdam,_extra=extra
;
;  Calculate 3d power spectrum.  The resulting spectrum satisfies
;  total(a^2)=wavenumber(1)*total(spectrum).
;
;  Only handles the common case where all sides of the periodic
;  domain are equal (=size).
;
; OPTIONS:
;
;    SPECTRUM : return value (OUT)
; WAVENUMBERS : K-values (OUT)
;  COMPENSATE : multiply spectrum with K^COMPENSATE
;
;       TITLE : title of plot
;      XTITLE : x-title of plot
;      YTITLE : y-title of plot
;      OPLT=2 : dashed line overplot (see PSYM in the IDL docs)
;        SIZE : physical size of the region
;
;    /AVERAGE : average of power in shell x volume of the shell (smoother results)
;    /CORNERS : include corners of k-space
;     /KINDEX : place k-values on integers
;      OFFSET : tuning parameter
;     /NOPLOT : don't plot
;    /POTSDAM : Potsdam choices
;
;  FIT=[3,20] : fit the spectrum btw. k=3 and k=20
;        AFIT : 1-st degree polynomial coeffs of fit
;
; Other options are passed to PLOT / OPLOT
;
;-------------------------------------------------------------------------------
if n_elements(title ) eq 0 then title =''
if n_elements(xtitle) eq 0 then xtitle='Wavenumber'
if n_elements(ytitle) eq 0 then ytitle='Power'
if n_elements(o     ) eq 0 then o     =0.18

if keyword_set(kindex) then o=0.

if keyword_set(potsdam) then begin
  kindex=1
  o=0.5
  plo=0.
end else begin
  plo=o
end
print,o,plo

s=size(a) & nx=s(1) & ny=s(2) & nz=s(3)
if n_elements(size) eq 0 then size=2.*!pi
k0=2.*!pi/size

if keyword_set(debug) then print,'doing fft ...'
t0=systime(1)
if s(0) eq 3 then begin
  ta=fft(a,-1)                                          ; complex 3D FFT
  ta=shift(ta,nx/2,ny/2,nz/2)                           ; shift to corner
end else if s(0) eq 4 then begin
  ta=complexarr(s(1),s(2),s(3),s(4))
  for i=0,s(4)-1 do ta(*,*,*,i)=shift(fft(a(*,*,*,i),-1),nx/2,ny/2,nz/2)
end
parseval=aver(a^2)
t1=systime(1)
if keyword_set(debug) then print,t1-t0,' sec'

t0=systime(1)
x1=indgen(nx)-nx/2
y1=indgen(ny)-ny/2
z1=indgen(nz)-nz/2
x1=reform(x1,nx,1,1)
y1=reform(y1,1,ny,1)
z1=reform(z1,1,1,nz)
rad=rebin(x1,nx,ny,nz)^2
rad=temporary(rad)+rebin(y1,nx,ny,nz)^2                 ; conserve memory
rad=temporary(rad)+rebin(z1,nx,ny,nz)^2
rad=temporary(sqrt(rad))                                ; wave vector radius
t1=systime(1)
if keyword_set(debug) then print,t1-t0,' sec'

if keyword_set(corners) then $
  imax=fix(max(rad)) $                                  ; include corners
else $
  imax=min([nx,ny,nz]/2)                                ; cut corners
spectrum=fltarr(imax)
ka=fltarr(imax)
if n_elements(average) ne 1 then average=spectrum

if keyword_set(debug) then print,'looping over shells'
t0=systime(1)
for i=0L,imax-1 do begin                                ; cut the corners?
  ii=where((rad ge i+o-.5) and (rad lt i+o+.5),nw)      ; indices into shell
  if keyword_set(debug) and i le 10 then print,i,nw
  if s(0) eq 3 then begin
    tmp=ta(ii) & tmp=abs(temporary(tmp))^2              ; conserve memory
  end else begin
    tmp=abs(ta(*,*,*,0))^2
    for k=1,s(4)-1 do tmp=temporary(tmp)+abs(ta(*,*,*,k))^2
  end
  ka(i)=aver(rad(ii))                                   ; average wave number
  if keyword_set(kindex) then ka(i)=i+plo
  if n_elements(average) ne 1 then begin                ; average array or not?
    spectrum(i)=total(tmp)                              ; sum over shell
    average(i)=aver(tmp*rad(ii)^2)*4.*!pi               ; average*volume
  end else begin
    spectrum(i)=aver(tmp*rad(ii)^2)*4.*!pi              ; average*volume
  end
endfor
t1=systime(1)
if keyword_set(debug) then print,t1-t0,' sec'

wns=ka*k0                                               ; normalize wavenumbers
spectrum=spectrum/k0                                    ; normalize to k0

if keyword_set(corners) then begin
  spectrum=spectrum*(parseval/sum(spectrum))            ; make sure Parseval is OK
end

if keyword_set(compensate) then spectrum(1:*)=spectrum(1:*)*wns(1:*)^compensate

if n_elements(c) gt 0 then begin
  power3d,c,aver=aver,compensate=compensate,kindex=kindex,sp=sp1,/noplot,offset=o,potsdam=potsdam
  spectrum=spectrum+sp1
end
if n_elements(b) gt 0 then begin
  power3d,b,aver=aver,compensate=compensate,kindex=kindex,sp=sp1,/noplot,offset=o,potsdam=potsdam
  spectrum=spectrum+sp1
end

if keyword_set(fit) then begin
  if n_elements(fit) eq 2 then begin
    i=fit
  end else begin
   i=[1,10]
  end
  kfit=alog(wns(i(0):i(1)))
  pfit=alog(spectrum(i(0):i(1)))
  afit=poly_fit(kfit,pfit,1,yfit)
  if keyword_set(compensate) then afit[1]=afit[1]-compensate
  print,afit[1]
end

if keyword_set(noplot) then return

if n_elements(oplt) eq 0 then begin
  plot_oo,wns(1:*),spectrum(1:*),$              ; skip DC in log-log
    title=title,xtitle=xtitle,ytitle=ytitle,$
    _extra=extra
endif else begin
  oplot,wns(1:*),spectrum(1:*),$
    line=oplt, _extra=extra
endelse
if keyword_set(fit) then oplot,exp(kfit),exp(yfit)

return
end

;+
 PRO test_power3d,n,iter=iter,index=testindex,_extra=e
;
;  Amplitudes           ~ k^p
;  Square ampl          ~ k^{2p}
;  Power spectrum       ~ k^{2p+2}
;
;  p=-1 should give a flat power spectrum P(k), and p=-2 should give a flat k^2 P(k)
;-
  default,n,16
  default,iter,100
  kk=wavenumbers(n,n,n)
  kk(0,0,0)=1.

  sp1=0.
  sp2=0.
  powerindex=0.
  ft=complexarr(n,n,n)
  for i=1,iter do begin
    ft(0:n/2-1,0:n/2-1,0:n/2-1)=exp(complex(0.,1.)*2.*!pi*randomu(seed,n/2,n/2,n/2))
    ft(0      ,0:n/2-1,0:n/2-1)=sqrt(.5)*ft(0      ,0:n/2-1,0:n/2-1)
    ft(0:n/2-1,0      ,0:n/2-1)=sqrt(.5)*ft(0:n/2-1,0      ,0:n/2-1)
    ft(0:n/2-1,0:n/2-1,0      )=sqrt(.5)*ft(0:n/2-1,0:n/2-1,0      )
    ft=ft/kk^(powerindex/2.+1.)
    ft(0,0,0)=0.
    f=float(fft(/inv,ft))
    power3d,f,sp=p1,wa=k,/nopl,_extra=e
    sp1=sp1+p1
    power3d,f,sp=p2,wa=k,/nopl,/aver,_extra=e
    sp2=sp2+p2
  end
  print,k(0:n/2-1)
  print,sp1(0:n/2-1)/iter
  print,sp2(0:n/2-1)/iter

  sp1=0.
  sp2=0.
  default,testindex,1.75
  for i=1,iter do begin
    ft(0:n/2-1,0:n/2-1,0:n/2-1)=exp(complex(0.,1.)*2.*!pi*randomu(seed,n/2,n/2,n/2))
    ft(0      ,0:n/2-1,0:n/2-1)=sqrt(.5)*ft(0      ,0:n/2-1,0:n/2-1)
    ft(0:n/2-1,0      ,0:n/2-1)=sqrt(.5)*ft(0:n/2-1,0      ,0:n/2-1)
    ft(0:n/2-1,0:n/2-1,0      )=sqrt(.5)*ft(0:n/2-1,0:n/2-1,0      )
    ft=ft/kk^(testindex/2.+1.)
    ft(0,0,0)=0.
    f=float(fft(/inv,ft))
    power3d,f,sp=p1,wa=k,/nopl,_extra=e
    sp1=sp1+p1
    power3d,f,sp=p2,wa=k,/nopl,/aver,_extra=e
    sp2=sp2+p2
  end
  print,testindex
  print,k(0:n/2-1)^testindex*sp1(0:n/2-1)/iter
  print,k(0:n/2-1)^testindex*sp2(0:n/2-1)/iter

END