PRO STRIPCOR,sph,corsph,skipstrip=SKIPSTRIP
;
; Subroutine for DIG_MAIN2,3
; Correction for horizontal structures and asymmetric intensity distribution
; in rectified (by DISKCOR) spectroheliograms.
;
; INPUT:  SPH - rectified spectroheliogram
; OUTPUT: CORSPH - corrected spectroheliogram
; KEYWORD: SKIPSTRIP - only the asymmetric intensity distribution is corrected
;
; Copyright M. Sobotka and M. Klvana, ASU AVCR Ondrejov, Czech Rep., 2009,2014

sph=sph>1

; sph centre
si=size(sph)
cenx=si(1)/2+1
ceny=si(2)/2+1

; cartesian coordinates
x1=findgen(si(1))
xx=fltarr(si(1),si(2))
for j=0,si(2)-1 do xx(*,j)=x1
y1=findgen(si(2))
yy=fltarr(si(1),si(2))
for i=0,si(1)-1 do yy(i,*)=y1
xx=xx-cenx
yy=yy-ceny

; distance array
rr=sqrt(xx^2+yy^2)

; intensities averaged in bins with increasing distance from the centre
;----------------------------------------------------------------------
bin=1   ;binsize parameter

h=histogram(rr,binsize=bin,location=r1,reverse_indices=ris)
nh=n_elements(h)

cli1=fltarr(nh)     ;array to store avg. intensities depending on rr
cld1=fltarr(nh)     ;array to store stdev intensities depending on rr
cli2=fltarr(nh)     ;array to store avg. int in +/- stdev range

for i=20,nh-21 do begin     ;averaging intensities in each rr bin

    chlindex=(ris(ris(i):ris(i+1)-1))   ;list of indices for the bin
    mom=MOMENT(sph(chlindex))
    cli1(i)=mom(0)
    cld1(i)=sqrt(mom(1))
        ;take only intensities within +/- cld1 (= stdev)
    w=where(sph(chlindex) gt cli1(i)-cld1(i) and sph(chlindex) lt cli1(i)+cld1(i),c)
    if c gt 0 then newindex=chlindex(w) else MESSAGE,'Not enough points'
    cli2(i)=total(sph(newindex))/c

endfor

    ; refill the beginning with a constant
cli2(0:19)=total(cli2(20:49)/30.)
    ; abscissa values
raxis=float(r1)+bin/2.
    ; symmetrization of a range 0:400 pix
raxis2=fltarr(801)
cli22=fltarr(801)
raxis2(0:399)=reverse(raxis(1:400))*(-1.)
raxis2(400:*)=raxis(0:400)
cli22(0:399)=reverse(cli2(1:400))
cli22(400:*)=cli2(0:400)

    ; fit a polynomial
coef=POLY_FIT(raxis2,cli22,6)   ;,yfit=cli23)

; intensity reference surface IRF
;--------------------------------
irf=float(sph)       ;most of solar disc will be replaced by IRF
w400=where(rr lt 400.)  ;pixels with distances rr < 400 from centre
r=rr(w400)              ;array of radial distances
; calculation of reference intensity in the pixels with rr < 400
irf(w400)=coef(0)+coef(1)*r+coef(2)*r^2+coef(3)*r^3+coef(4)*r^4+coef(5)*r^5+coef(6)*r^6

; CLV and tilt removal
;------------------------------
spir=sph/irf    ; CLV removal
xq1=cenx-281    ; square subfield (282=399*sqrt(2))
xq2=cenx+281
yq1=ceny-281
yq2=ceny+281
plan=sfit(spir(xq1:xq2,yq1:yq2),1,kx=kxy) ;tilt calculation
xx=xx+cenx
yy=yy+ceny
plan=kxy(0,0)+kxy(1,0)*yy+kxy(0,1)*xx+kxy(1,1)*xx*yy  ;tilt reference plane
sprt=spir/plan  	; flat spectroheliogram with removed CLV and tilt

corsph=sph/plan		;corsph is the centered spectroheliogram

if keyword_set(skipstrip) then goto,fin ;*********************************************

;adaptation of IRF intensity to spectroheliogram (new, 7.10.2014)
mecs=MOMENT(corsph(w400))
mirf=MOMENT(irf(w400))
ftor=mecs(0)/mirf(0)	; ratio of averages
irf=irf*ftor
;print,'IRF intensity factor',ftor

;active regions localization and removal
;----------------------------------------
mom=MOMENT(sprt(w400))
mea=mom(0)
std=sqrt(mom(1))
std=std>0.03       ;tunable - minimum stdev

wout=where(rr gt 400.,c)
if c gt 0 then sprt(wout)=mea   ; set outer region to mean
wk1=where(sprt lt (mea-2.5*std),ck1)    ;tunable - stdev coefficient 2.5
wk2=where(sprt gt (mea+1.0*std),ck2)    ;tunable - stdev coefficient 2.5

if ck1 gt 0 then corsph(wk1)=irf(wk1)   ;active regions removed from corsph
if ck2 gt 0 then corsph(wk2)=irf(wk2)

;completing outer part of the IRF
;--------------------------------
;symmetrization of the outer ring of irf by averaging 8 images rotated by 45 deg
if si(1) gt si(2) then irf1=irf((si(1)-si(2))/2:(si(1)+si(2))/2-1,*) $
  else irf1=irf(*,(si(2)-si(1))/2:(si(2)+si(1))/2-1)
irf2=irf1/8.
for i=1,7 do irf2=irf2+rot(irf1,i*45)/8.
if si(1) gt si(2) then irf((si(1)-si(2))/2:(si(1)+si(2))/2-1,*)=irf2  $
  else irf(*,(si(2)-si(1))/2:(si(2)+si(1))/2-1)=irf2

;removing intensity jump
w401=where(rr ge 400.)
wi1=where(rr gt 398. and rr lt 400.,ci1)
wi2=where(rr ge 400. and rr lt 402.,ci2)
mi1=total(irf(wi1))/ci1
mi2=total(irf(wi2))/ci2
irf(w401)=irf(w401)*mi1/mi2

; correction for horizontal strips
;---------------------------------

;1D average profiles and 1D strip profile
scuk=total(corsph(cenx-100:cenx+100,*),1)/201.  ;avg. column of corsph with strips
rcuk=total(irf(cenx-100:cenx+100,*),1)/201.     ;avg. column of ref.surface (no strips)
strip=(scuk/rcuk)*(total(rcuk)/total(scuk))     ;normalized strip profile
strip=strip>0.9
strip=strip<1.1

; scaling strips profile for disk edges
x=x1-cenx
scf=coef(0)+coef(1)*x+coef(2)*x^2+coef(3)*x^3+coef(4)*x^4+coef(5)*x^5+coef(6)*x^6
scf=scf/max(scf)

;cleaning horizontal strips to obtain the final corrected spectroheliogram
corsph1=sph/plan          ;centering intensity in original spectroheliogram
for i=0,si(1)-1 do corsph(i,*)=corsph1(i,*)/(strip*scf(i)+1-scf(i))

fin:
END