; THIS script MUST be run from same dir as galaxy being debanded
;
; To run just type:
;
; IDL> .RUN fixband.pro     ;if prompted, to proceed type:
; IDL> .continue 
function fixband2,m,ybad,width,x1l,x2l,x1r,x2r,xmin,xmax,thresh,diag=diag
ml=fixband(m,ybad,width,x1l,x2l,thresh,band_l,diag=diag)
mr=fixband(m,ybad,width,x1r,x2r,thresh,band_r,diag=diag)

sz=size(m)
mf=m
for i=0,sz(2)-1 do begin &$
    mf(xmin:x2l,i)=m(xmin:x2l,i)-band_l(i) &$
    mf(x1r:xmax,i)=m(x1r:xmax,i)-band_r(i) &$
    mf(x2l:x1r,i)=m(x2l:x1r,i)-interpol([band_l(i),band_r(i)],[x2l,x1r],x2l+indgen(x1r-x2l+1)) &$
endfor
return,mf
end

function fixband,m,ybad,width,x1,x2,thresh,band,diag=diag
sz=size(m)
prof=fltarr(sz(2))
for i=0,sz(2)-1 do prof(i)=median(m(x1:x2,i))

x=indgen(sz(2))
lbase=where(abs(x-ybad) gt 2*width and abs(x-ybad) lt 7*width)
cc=poly_fit(x(lbase),prof(lbase),1)
base=poly(x,cc)
band=prof-base
lok=where(abs(x-ybad) gt 2.*width)
band(lok)=0.

if keyword_set(diag) then begin
    window,2
    !p.multi=[0,1,2]
    plot,x,prof,xr=ybad+width*8*[-1,1],/ys,psym=10
    oplot,x,base
    plot,x,band,xr=ybad+width*8*[-1,1],psym=10
    !p.multi=0
endif

mnew=m
for i=0,sz(2)-1 do mnew(x1:x2,i)=m(x1:x2,i)-band(i)
return,mnew

end

; #1
print,'reading the original IRAC image...'
m=readfits('mosaic.fits',h)

; #2
print,'displaying original image...'
atv,m,min=0,max=4


sz=size(m)

; #3
; left top galaxy fix region
print,'left top galaxy...'
x1=127 & x2=379 & ybad=316 & width=3
x1l=x1 & x2l=x2
mf_l=fixband(m,ybad,width,x1,x2,0.05,band_l)

; #4
; right top galaxy fix
print,'right top galaxy...'
x1=380 & x2=621 & ybad=316 & width=3
x1r=x1 & x2r=x2
mf_r=fixband(m,ybad,width,x1,x2,0.05,band_r)


; #5
; define top fix/replace region
xmin=127 & xmax=621
; combine top left and right fixes
x1l=127 & x2l=379 & ybad=316 & width=3 & x1r=380 & x2r=621 
mnew=fixband2(m,ybad,width,x1l,x2l,x1r,x2r,xmin,xmax,0.05)


; #8
;; left bottom galaxy fix
;print,'left bottom galaxy...'
;x1=456 & x2=640 & ybad=332 & width=12
;x1l=x1 & x2l=x2
;mf_lb=fixband(m,ybad,width,x1,x2,0.05,band_l)

;; right bottom galaxy fix
;print,'left top galaxy...'
;x1=641 & x2=815 & ybad=332 & width=12
;x1r=x1 & x2r=x2
;mf_rb=fixband(m,ybad,width,x1,x2,0.05,band_r)

;; define new fix/replace region
;xmin=456 & xmax=815
;x1l=456 & x2l=640 & ybad=332 & width=12 & x1r=641 & x2r=815 
;mnew2=fixband2(mnew,ybad,width,x1l,x2l,x1r,x2r,xmin,xmax,0.05)


; #6
writefits,'arp242ch4_db.fits',mnew,h

stop

; #7
atv,mnew,min=0,max=4

end