; program ruball.pro = annotated version for 2007-04-12/QuietSun
; init: Sep 18 2010 SST 
; last: Sep 22 2010 SST
; note: ;;; = Pit Suetterlin (PS) comments in original version 2007-04-12
;       ;   = Rob Rutten (RR) comments annotated version Sep 2010 
;       ; @ = todo
;       ; @ adapt? = parameter likely to be adapted for other data sets

; rubball does:
; - rubbersheeting of each cube to a reference made from 3-image
;   sliding time average of the average of gb, bc, rc
; 

; 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

caflip = 1
skip_carubb = 0

f = file_search('../Cubes/*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)
date = f(3)
dt = '-'+f(6)

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

tmp = file_search('../Cubes/*ha*-cube-*', count=n_ha)

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

names = ['-gb-0000', '-bc-0000', '-rc-0000']
FOR i=0, n_ha-1 DO $
  names = [names, strmid(tmp(i), (strsplit(tmp(i), '-'))(1)-1, 8)]
tmp = file_search('../Cubes/*ca*-cube-*')
names = [names, strmid(tmp(0), (strsplit(tmp(0), '-'))(1)-1, 8)]
names = names+'-'+date+'-cube-'

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

;;;----------  no neccessary changes past this line (hopefully)
; @ add ba = barium wavelengths, comparable to ha for KvdL 

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), '../Cubes/'+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)
show, ref   ; 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
    tvscl, ref ; displays the reference image per time step @ out for batch
    xyouts, 20, 20, im, /dev, charsiz=2
    empty
; remap_image: PS function; it calls PS function shc.pro
    p1 = remap_image(ref, readframe(u_i(0), im+start), cs, sw, /QUINT, $
                     PLOT=3, /SAME, MEDIAN=med)
    writeframe, p1, u_o(0), im
    p1 = remap_image(ref, readframe(u_i(1), im+start), cs, sw, /QUINT, $
                     PLOT=3, /SAME, MEDIAN=med)
    writeframe, p1, u_o(1), im
    p1 = remap_image(ref, readframe(u_i(2), im+start), cs, sw, /QUINT, $
                     PLOT=3, /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
    tvscl, ref
    xyouts, 20, 20, im, /dev, charsiz=2
    empty
    p1 = remap_image(ref, readframe(u_i(0), im+start), cs, sw, /QUINT, $
                     PLOT=3, /SAME, MEDIAN=med)
    writeframe, p1, u_o(0), im
    p1 = remap_image(ref, readframe(u_i(1), im+start), cs, sw, /QUINT, $
                     PLOT=3, /SAME, MEDIAN=med)
    writeframe, p1, u_o(1), im
    p1 = remap_image(ref, readframe(u_i(2), im+start), cs, sw, /QUINT, $
                     PLOT=3, /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     ; enhance: PS procedure

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

IF caflip EQ 0 THEN BEGIN
      ;;; Ca:  Only one
    tmp = assoc(u_o(nn-1), intarr(resx, resy))
    enhance, tmp, ASSOC=num
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)  
    wset, 0
    plot_io, pi & oplot, pi1
    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
    wset, 1
    tvscl, p1
    xyouts, 10, 10, /dev, string(img+1, num, form="(i3,' of ',i3)")
    empty
    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
        tvscl, p1
        xyouts, 10, 10, /dev, string(img+1, num, form="(i3,' of ',i3)")
        empty
    ENDFOR
ENDFOR

spawn, '[ -d ../Cubes1 ] || mkdir ../Cubes1'

get_lun, unit
FOR im=0, nn-2 DO BEGIN
  ; resstring: PS function
    openw, unit, '../Cubes1/'+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, '../Cubes1/'+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, '../Cubes1/'+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, '../Cubes1/'+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