PRO scatter,aur,lam
;
; Determination of scattered light in sunspots/pores
; in the central zone of the solar disk by interactive graphical
; fitting the computed aureole profile to the measured one.
; Measured aureole profile must extend from 25" inside
; the limb to 30" outside the limb.
;
; Method:  Martinez Pillet 1992, Solar Phys. 140, 207
;
; INPUTS: aur - measured aureola profile
;	  lam - working wavelength
;	  pre-calculated aureole contributions in save files
; OUTPUT: aureole spread parameters B(Lorentz), b1 and b2 (Gauss)
;	  and corresponding weights M, m1, m2; stray light parameter c;
;	  stray-light contamination 1-c (pores) and (1-c)/2 (spots)
;
; Subroutines called: clv.pro, deriv (IDL)
;
; Used for correction of SST 2004 observations.
; Michal, 13.7.2004
;
on_error,1
; -----------------------------

; MEASURED AUREOLE

scale=0.041		;arcsec per pixel
;----------------------------------------

; Determine limb position - maximum gradient

gra=DERIV(aur)
ma=max(abs(gra),pli)	;pli=position of the limb (pixels)

; Find I10 - intensity 25" inside the limb, for intensity scaling

d10=round(25/scale)		;25" in pixels

if gra(pli) lt 0 then begin	;disk is on the left
  p10=pli-d10
  if p10 lt 0 then message,'ERROR - disk < 25" '
endif else begin		;disk is on the right
  aur=reverse(aur)		; - reverse aur to have it on the left
  gra=DERIV(aur)
  ma=max(abs(gra),pli)
  p10=pli-d10
  if p10 lt 0 then message,'ERROR - disk < 25" '
endelse			;p10 - position 25" inside the limb -
			; it corresponds to 935" from disk center
i10=aur(p10)

au=aur(p10:*)	; remove pixels before p10 (p10 set to zero)
n=n_elements(au)
xau=findgen(n)*scale-25.	;abscissa for plotting, zero at limb

; Intensity normalization to disk_center = 1

au=au*CLV(25,lam)/i10
; -------------------------------------

; COMPUTED AUREOLE

xac=findgen(61)-30.	;plotting scale, step 1" in computed aureole
window,0,ys=640
loadct,3
plot_io,xau,au,/xst

; Restoring pre-calculated aureole contributions

restore,'phi_val_'+strtrim(fix(lam),2)	;l10 and g* arrays

; Asking for parameters

new1:
read,'  Gaussian b1 (3,5,10,15,20): ',b1
case b1 of
	 3: pg1=g3
	 5: pg1=g5
	10: pg1=g10
	15: pg1=g15
	20: pg1=g20
	else: goto, new1
endcase

new2:
read,'  Gaussian b2 (40,50,60,70,80,90): ',b2
case b2 of
	40: pg2=g40
	50: pg2=g50
	60: pg2=g60
	70: pg2=g70
	80: pg2=g80
	90: pg2=g90
	else: goto,new2
endcase

new3:
read,'  Gaussian weights m1 and m2: ',m1,m2
mlo=1.-m1-m2
if mlo lt 0. then begin
	print,'*** WARNING: sum of weights > 1, repeat!
	goto,new3
endif

new4:
read,'  Scaling factor c0: ',c0

phi=c0*(mlo*l10(1:*)+m1*pg1(1:*)+m2*pg2(1:*))	; computed aureole
						;from 930" to 990"

oplot,xac(31:*),phi(31:*),col=210	; aureole from 1" out of the limb 

print,'Change scaling c0 (y/n) ?'
if get_kbrd(1) eq 'y' then begin
	oplot,xac(31:*),phi(31:*),col=150
	goto,new4
endif

print,'Change all parameters (y/n) ?'
if get_kbrd(1) eq 'y' then begin
	oplot,xac(31:*),phi(31:*),col=150
	goto,new1
endif

c=1.-c0*(mlo*l10(0)+m1*pg1(0)+m2*pg2(0))
print,'Stray light parameter c =',c
print,'Contamination in pore:',1.-c
print,'Contamination in spot:',(1.-c)/2

for i=0,30 do phi(i)=phi(i)+c*CLV(30-i,lam)   ;adding disk intensity
oplot,xac,phi,col=210	;final plot from -25" to 30"

window,2,ys=640		;linear-scale plot
plot,xau,au,/xst
oplot,xac,phi,col=210

END