; program ruball.pro = generic annotated classic version
; init: Sep 18 2010 
; last: Dec 13 2010 
; hist: Sep 2010 SST original from dotdb 2007-04-12 = RR trial
;       Dec 2010 Deil for 2005-07-13-spot = RR EBduds
; note: classic version for DOT data up to spring 2007; Hitachi KPF100 cameras
;       ;;; = Pit Suetterlin (PS) comments in original
;       ; = Rob Rutten (also ;RR) 
;       ; ?? = issue
;       ; @ = todo
;       ; @ adapt? = parameter likely to be adapted for other data sets
;       ; # comment
; does: - rubbersheeting of each cube to a reference made from 3-image
;         sliding time average of the average of gb, bc, rc
;         except ca (standalone)
;       - equalize intensities
;       - enhance power spectrum for all time steps to that of best moment
; todo: @ add ba = barium wavelengths, comparable to ha for KvdL 


; temporary assocs, as in gen-cubes.pro
FUNCTION Readframe, unit, n
COMMON assocsize, sx0, sy0, imsiz
img = intarr(sx0, sy0)
point_lun, unit, 1l*imsiz*n
readu, unit, img
return, img
END

PRO Writeframe, img, unit, n
COMMON assocsize, sx0, sy0, imsiz
tmp = intarr(sx0, sy0)
s = size(img, /dim) < [sx0, sy0]
tmp(0:s(0)-1, 0:s(1)-1) = img(0:s(0)-1, 0:s(1)-1) < 32767  ; cut int overflows
point_lun, unit, 1l*imsiz*n
writeu, unit, tmp
END

; ===================================== start main part

COMMON assocsize, sx0, sy0, imsiz

;RR @ adapt = set choices
plotout=0                   ; @ adapt? get screen output or not?
caflip = 0                  ; @ adapt? have data flipped ca filter or not?
skip_carubb = 0             ; @ adapt? choose whether to rubb ca?
incubedirname = 'Cubesfull'     ; @ adapt?
outcubedirname = 'Cubesfullrubb'   ; @ adapt?

;;;----------  no changes needed past this line? (hopefully)

; set plotvectors value (none or length)
if (plotout eq 0) then plotvectors=0 else plotvectors=3 

; get the cube parameters from the gb cube file name
f = file_search('../'+incubedirname+'/*gb*-cube-*')
f = strsplit(dos_name(f(0)), '-', /extract)  ; dos_name: PS function
r = strsplit(f(5), 'x', /extract) ; 3 dims   ; strsplit: IDL function

base = f(0)    ; example AR0789
date = f(3)    ; example 20050713
dt = '-'+f(6)  ; example -dt30  # the value isn't used in this program

resx = fix(r(0)) ; x size in px
resy = fix(r(1)) ; y size in px

start = 0
num = fix(r(2))  ; nr images
base_new = base

; define file names
names = ['-gb-0000', '-bc-0000', '-rc-0000']
  ; add names and get number of Halpha cubes
tmp = file_search('../'+incubedirname+'/*ha*-cube-*', count=n_ha)
FOR i=0, n_ha-1 DO $
  names = [names, strmid(tmp(i), (strsplit(tmp(i), '-'))(1)-1, 8)]
  ; strmid takes 8 chars after the first hyphen = after object name
  ; so the dir name should NOT contain a hyphen
  ; add names of ca cube(s)
tmp = file_search('../'+incubedirname+'/*ca*-cube-*')
names = [names, strmid(tmp(0), (strsplit(tmp(0), '-'))(1)-1, 8)]
  ; add date and "cube" to names 
names = names+'-'+date+'-cube-'

IF NOT keyword_set(cs) THEN cs = 24 ; adapt? cell size in px both sides
IF NOT keyword_set(sw) THEN sw = 8  ; adapt? step width in px
med = 3
nsum = 3  ; average over 3 temporally?
ns2 = nsum/2  ; 1

nn = n_ha+4 ; number of wavelengths n_ha + bc + gb + rc + ca
sx0 = resx
sy0 = resy
imsiz = 2l*resx*resy
  ;;; input data
u_i = indgen(nn)+1 ; u_i = unit numbers for input starting at 1
FOR i=0, nn-1 DO openr, u_i(i), '../'+incubedirname+'/'+base+names(i)+'*'
  ;;; temporary output data
u_o = indgen(nn)+21  ; u_o = unit nr for output starting at 21
FOR i=0, nn-1 DO openw, u_o(i), 'cube-'+strtrim(i+1, 2)+'.asc'

ref = fltarr(resx, resy)
if (plotout) then show, ref   ;RR show = PS procedure
FOR im=ns2, num-1 DO BEGIN
    ref(*) = 0.
; average over 3 continua, nsum in time, start over half temp interval
; explicit (not calling rubb.pro) because of averaging over 3 cont wavs
    FOR j=im-ns2, ((im+ns2) < (num-1)) DO FOR i=0, 2 DO BEGIN
        tmp = readframe(u_i(i), j+start)
        ref = ref+tmp/avg(tmp)
    ENDFOR
    if (plotout) then begin
      tvscl, ref ; displays the reference image per time step @ out for batch
      xyouts, 20, 20, im, /dev, charsiz=2
      empty
    endif
; remap_image: PS function; it calls PS function shc.pro
    p1 = remap_image(ref, readframe(u_i(0), im+start), cs, sw, /QUINT, $
                     PLOT=plotvectors, /SAME, MEDIAN=med)
    writeframe, p1, u_o(0), im
    p1 = remap_image(ref, readframe(u_i(1), im+start), cs, sw, /QUINT, $
                     PLOT=plotvectors, /SAME, MEDIAN=med)
    writeframe, p1, u_o(1), im
    p1 = remap_image(ref, readframe(u_i(2), im+start), cs, sw, /QUINT, $
                     PLOT=plotvectors, /SAME, MEDIAN=med, MAP=map)
    writeframe, p1, u_o(2), im
    
    FOR i=0, n_ha-1 DO BEGIN
; warp_tri: IDL function
        p1 = warp_tri(map(*, 2), map(*, 3), map(*, 0), map(*, 1), $
                      readframe(u_i(i+3), im+start), /QUINT)
        writeframe, p1, u_o(i+3), im
    ENDFOR
ENDFOR

; now the start
FOR im=ns2-1, 0, -1 DO BEGIN
    ref(*) = 0.
    FOR j=0, ns2 DO FOR i=0, 2 DO BEGIN  ;1 > ns2
        tmp = readframe(u_i(i), j+start)
        ref = ref+tmp/avg(tmp)
    ENDFOR
    if (plotout) then begin ;RR switch
      tvscl, ref
      xyouts, 20, 20, im, /dev, charsiz=2
      empty
    endif
    p1 = remap_image(ref, readframe(u_i(0), im+start), cs, sw, /QUINT, $
                     PLOT=plotvectors, /SAME, MEDIAN=med)
    writeframe, p1, u_o(0), im
    p1 = remap_image(ref, readframe(u_i(1), im+start), cs, sw, /QUINT, $
                     PLOT=plotvectors, /SAME, MEDIAN=med)
    writeframe, p1, u_o(1), im
    p1 = remap_image(ref, readframe(u_i(2), im+start), cs, sw, /QUINT, $
                     PLOT=plotvectors, /SAME, MEDIAN=med, MAP=map)
    writeframe, p1, u_o(2), im
    FOR i=0, n_ha-1 DO BEGIN
        p1 = warp_tri(map(*, 2), map(*, 3), map(*, 0), map(*, 1), $
                      readframe(u_i(i+3), im+start), /QUINT)
        writeframe, p1, u_o(i+3), im
    ENDFOR
ENDFOR

  ;;; now ca - has to be aligned on itself
FOR im=0, num-1 DO writeframe, readframe(u_i(nn-1), im+start), u_o(nn-1), im
ca = assoc(u_o(nn-1), intarr(resx, resy))

IF skip_carubb EQ 1 THEN GOTO, nocarubb

IF caflip EQ 0 THEN BEGIN
     ;;; same wavelength for all images
    rubb_file, ca, num, cell=cs, step=sw, range=4, MEDIAN=med, average=nsum, $
      /quint, plot=4, /same
ENDIF ELSE BEGIN
      ;;; switch wing/core: do separately
    tmp = intarr(resx, resy, num/2)
    FOR i=0, num-1, 2 DO tmp(*, *, i/2) = readframe(u_o(nn-1), i)
    rubb, tmp, cell=cs, step=sw, range=4, MEDIAN=med, average=nsum, $
      plot=4, /same, /quint
    FOR i=0, num-1, 2 DO writeframe, tmp(*, *, i/2), u_o(nn-1), i

    FOR i=1, num-1, 2 DO tmp(*, *, i/2) = readframe(u_o(nn-1), i)
    rubb, tmp, cell=cs, step=sw, range=4, MEDIAN=med, average=nsum, $
      plot=4, /same, /quint
    FOR i=1, num-1, 2 DO writeframe, tmp(*, *, i/2), u_o(nn-1), i
ENDELSE

Nocarubb:

;;;------------
;;; now equalize intensities

av = fltarr(num)
FOR j=0, 2 DO BEGIN
    FOR i=0, num-1 DO av(i) = avg(readframe(u_o(j), i))  ; u_o = output unit
    av = float(avg(av)/av)
    FOR i=0, num-1 DO $
      writeframe, readframe(u_o(j), i)*av(i)+0.5, u_o(j), i
ENDFOR

IF caflip EQ 0 THEN BEGIN
    FOR i=0, num-1 DO av(i) = avg(readframe(u_o(nn-1), i))
    av = float(avg(av)/av)
    FOR i=0, num-1 DO writeframe, $
      readframe(u_o(nn-1), i)*av(i)+0.5, u_o(nn-1), i
ENDIF ELSE BEGIN
    av = fltarr(num/2)
    FOR i=0, num/2-1 DO av(i) = avg(readframe(u_o(nn-1), 2*i))
    av = float(avg(av)/av)
    FOR i=0, num/2-1 DO writeframe, $
      readframe(u_o(nn-1), 2*i)*av(i)+0.5, u_o(nn-1), 2*i

    FOR i=0, num/2-1 DO av(i) = avg(readframe(u_o(nn-1), 2*i+1))
    av = float(avg(av)/av)
    FOR i=0, num/2-1 DO writeframe, $
      readframe(u_o(nn-1), 2*i+1)*av(i)+0.5, u_o(nn-1), 2*i+1
ENDELSE

  ;;; halpha has to be scaled all wavelengths at once

av = fltarr(num)
FOR j=3, nn-2 DO FOR i=0, num-1 DO $
  av(i) = av(i)+avg(readframe(u_o(j), i))
av = float(avg(av)/av)
FOR j=3, nn-2 DO FOR i=0, num-1 DO $
  writeframe, readframe(u_o(j), i)*av(i)+0.5, u_o(j), i

;;; now enhance the cubes

tmp = assoc(u_o(0), intarr(resx, resy)) ; gb
enhance, tmp, ASSOC=num, plotout=plotout     ;RR enhance: PS procedure

tmp = assoc(u_o(1), intarr(resx, resy)) ; bc
enhance, tmp, ASSOC=num, plotout=plotout     ;RR

IF caflip EQ 0 THEN BEGIN
      ;;; Ca:  Only one
    tmp = assoc(u_o(nn-1), intarr(resx, resy))
    enhance, tmp, ASSOC=num, plotout=plotout     ;RR
ENDIF ELSE BEGIN
      ;;; or two positions
    tmp = intarr(resx, resy, num/2)
    FOR i=0, num-1, 2 DO tmp(*, *, i/2) = readframe(u_o(nn-1), i)
    enhance, tmp
    FOR i=0, num-1, 2 DO writeframe, tmp(*, *, i/2), u_o(nn-1), i
    FOR i=1, num-1, 2 DO tmp(*, *, i/2) = readframe(u_o(nn-1), i)
    enhance, tmp
    FOR i=1, num-1, 2 DO writeframe, tmp(*, *, i/2), u_o(nn-1), i
ENDELSE

;;; halpha:  get enhance function from rc, apply to all
; explicit, not calling enhance.pro, because the rc info is needed for
; others

rms = fltarr(num)
FOR i=0, num-1 DO rms(i) = stdev(readframe(u_o(2), i))
  ; stdev obsolete IDL pro but still in SolarSoft/gen (Scott McIntosh)
mx = max(rms, best)
fp = fft(readframe(u_o(2), best), -1)
fpa= abs(fp)  ; @ take out zeroes
fp = fp/fpa
pi = polarint(fpa)   ; polarint: PS function
x = dist(resx, resy)
FOR img=0, num-1 DO BEGIN
    IF img EQ best THEN CONTINUE
    fp1 = fft(readframe(u_o(2), img), -1)
    fpa1 = abs(fp1)
    fp1 = fp1/fpa1
      ;;; careful with zero values!
    IF min(fpa1) EQ 0 THEN fp1(where(fpa1 EQ 0)) = 0
    pi1 = polarint(fpa1)  
    if (plotout) then begin
      wset, 0
      plot_io, pi & oplot, pi1
    endif
    pq = smoother(pi/pi1, 1, 10000)  ; smoother: PS function
    q = interpolate(pq, x, missing=1)
      ;;; average mustn't be changed
    q(0) = 1.
    p1 = float(fft(q*fpa1*fp1, 1))
    writeframe, p1, u_o(2), img
    if (plotout) then begin  ;RR switch
      wset, 1
      tvscl, p1
      xyouts, 10, 10, /dev, string(img+1, num, form="(i3,' of ',i3)")
      empty
    endif
    FOR i=1, n_ha DO BEGIN
        fp1 = fft(readframe(u_o(2+i), img), -1)
        fpa1 = abs(fp1)
        fp1 = fp1/fpa1
        IF min(fpa1) EQ 0 THEN fp1(where(fpa1 EQ 0)) = 0
        p1 = float(fft(q*fpa1*fp1, 1))
        writeframe, p1, u_o(2+i), img
        if (plotout) then begin ;RR switch
          tvscl, p1
          xyouts, 10, 10, /dev, string(img+1, num, form="(i3,' of ',i3)")
          empty
        endif
    ENDFOR
ENDFOR

spawn, 'mkdir -p ../'+outcubedirname

get_lun, unit
FOR im=0, nn-2 DO BEGIN
  ; resstring: PS function
    openw, unit, '../'+outcubedirname+'/'+base_new+names(im)+resstring([resx, resy, num])+dt
    FOR j=0, num-1 DO $
      writeu, unit, readframe(u_o(im), j)
    close, unit
ENDFOR
  ;;; and Ca: two positions
im = nn-1
IF caflip EQ 0 THEN BEGIN
    openw, unit, '../'+outcubedirname+'/'+base_new+names(im)+resstring([resx, resy, num])+dt
    FOR j=0, num-1 DO $
      writeu, unit, readframe(u_o(im), j)
    close, unit
ENDIF ELSE BEGIN
    dt1 = '-dt'+strtrim(2*fix(strmid(dt, 3)), 2)
    openw, unit, '../'+outcubedirname+'/'+base_new+names(im)+resstring([resx, resy, num/2])+dt1+'-even'
    FOR j=0, num-1, 2 DO $
      writeu, unit, readframe(u_o(im), j)
    close, unit
    openw, unit, '../'+outcubedirname+'/'+base_new+names(im)+resstring([resx, resy, num/2])+dt1+'-odd'
    FOR j=1, num-1, 2 DO $
      writeu, unit, readframe(u_o(im), j)
    close, unit
ENDELSE

free_lun, unit
FOR i=0, nn-1 DO close, u_i(i)
FOR i=0, nn-1 DO close, u_o(i)

END