; 2007 (c) Shravan Hanasoge
; Email: shravan@solar.stanford.edu
; W. W. Hansen Experimental Physics Laboratory
; Stanford University, Stanford, CA 94305
;
; To compact the data from the code into one giant output
; The final file will be in single precision
; Also to compute the power spectrum
;
; You have to input the horizontal dimensions 
; Leave the 'L' at the end of the number
nx = 256L
ny = nx

; Difference between the indices of subsequent files
steps = 30L
; The start index
start = 0L
; The final index 
finish = 21600L
nt = (finish - start)/steps + 1L

; Cadence in seconds
deltat = 60.0

; Variable name: vx, vy, vz, p, bz (one of these)
varname = 'vz'
; The directory
directory='/tmp13/shravan/cartesian3/'+varname+'_at_'
;
; The final output will be in the same directory
; unless otherwise specified
; 
outputdir = directory

; The variables xlength, ylength etc. are only useful
; if you want to compute the power spectrum
; Dimensions are Mm
xlength = 200.0
ylength = xlength

; Set compute_power_spectrum = 1 if you wish to compute it
; Power spectrum will be placed in outputdir
compute_power_spectrum = 0



filename=STRARR(nt)

J=0
FOR I=start,finish,steps DO BEGIN
IF (I LT 10) THEN BEGIN 
 filename[J]=directory+'0000'+STRTRIM(STRING(i),1)+'_deltat.fits'
ENDIF
IF ((I LT 100) AND (I GE 10)) THEN BEGIN 
 filename[J]=directory+'000'+STRTRIM(STRING(i),1)+'_deltat.fits'
ENDIF
IF ((I LT 1000) AND (I GE 100)) THEN BEGIN 
 filename[J]=directory+'00'+STRTRIM(STRING(i),1)+'_deltat.fits'
ENDIF
IF ((I LT 10000) AND (I GE 1000)) THEN BEGIN 
 filename[J]=directory+'0'+STRTRIM(STRING(i),1)+'_deltat.fits'
ENDIF
IF ((I LT 100000) AND (I GE 10000)) THEN BEGIN 
 filename[J]=directory+STRTRIM(STRING(i),1)+'_deltat.fits'
ENDIF
J = J +1
ENDFOR

a=DBLARR(nx,ny,nt)
FOR I=0,nt-1 DO BEGIN

 temp = readfits(filename[I])
 a[*,*,I] = temp
ENDFOR

writefits,outputdir+varname+'_full.fits',float(a)

if (compute_power_spectrum EQ 1) THEN BEGIN

 a = FFT(a,/DOUBLE)
 dnu = 1.0/(Nt*deltat)
 Rsol=695.9894677 ; in Mm

 dx = xlength/(nx*1.0); in Mm
 dy = ylength/(ny*1.0); in Mm

 dkx = 1.0/(nx*dx)
 dky = 1.0/(ny*dy)

 Nkh = 128L
 khmax = ((nx*dkx*0.5)^2. + (ny*dky*0.5)^2.)^0.5
 khmin = 0.0 
 kh = DINDGEN(Nkh)*(khmax-khmin)/(Nkh*1.0)
 kays = DBLARR(Nx,Ny)
 power = DBLARR(Nkh,Nt/2+1)
 dkh = kh[2]-kh[1] 

 FOR J= 0,Ny/2 DO BEGIN
  FOR I = 0,Nx/2 DO BEGIN
   kays[I,J] = ((dkx*I)^2. + (dky*J)^2.)^0.5
  ENDFOR
  FOR I = Nx/2+1,Nx-1 DO BEGIN
   kays[I,J] = ((dkx*(Nx-I))^2. + (dky*J)^2.)^0.5
  ENDFOR

 ENDFOR

 FOR J= Ny/2+1,Ny-1 DO BEGIN
  FOR I = 0,Nx/2 DO BEGIN
   kays[I,J] = ((dkx*I)^2. + (dky*(Ny-J))^2.)^0.5
  ENDFOR

  FOR I = Nx/2+1,Nx-1 DO BEGIN
   kays[I,J] = ((dkx*(Nx-I))^2. + (dky*(Ny-J))^2.)^0.5
  ENDFOR
 ENDFOR

 power[*,*] = 0.0
 FOR N = 0,Nkh-1 DO BEGIN

  c = nt*(WHERE( ABS(kays-kh(N)) LE dkh))

  FOR I=0,nt/2 DO BEGIN
   POWER[N,I] = MEAN(ABS(a[I + c])^2.)^0.5
  ENDFOR

  FOR I=nt/2+1,nt-1 DO BEGIN
   POWER[N,nt-I] = MEAN(ABS(a[I + c])^2.)^0.5 + power[N,nt-I]
  ENDFOR

 ENDFOR

 nu = DINDGEN(nt/2+1)*dnu*1000.
 loadct,3
 window,retain=2
 contour,(alog(power[*,*])),kh*Rsol,nu,/fill,xtit='!8l!17',ytit='!7m!17',xst=1,yst=1

 set_plot,'ps'
 device,file=outputdir + 'spectrum.eps',/ENCAPSULAT,bits=24,/color
  contour,(alog(power[*,*])),kh*Rsol,nu,/fill,xtit='!8l!17',ytit='!7m!17',xst=1,yst=1
 device,/close
 set_plot,'x'
 read,pause

ENDIF 

END