PRO darkfacula,iima,wima,date,dfac

;Visualization of dark faculae
;from couples of registered frames
;stored in movies iima (IR) and wima (WL).
;FOV must be close to the disk center.
;Frames must have the same size (not checked).
;Apodized edges in these frames have to be already removed.
;date = 23 or 25
;OUTPUT: dfac - movie of dark faculae

;PARAMETERS:
on_error,1

tp=0.90 	;WL threshold intensity to include pores

if (date ne 23) and (date ne 25) then  $
   message,'Sorry, I do not know how to handle this dataset'

if date eq 23 then begin
	qx1=246		; subfield of quiet region
	qx2=379
	qy1=182
	qy2=291
endif else begin
	qx1=304		; subfield of quiet region
	qx2=363
	qy1=135
	qy2=274
endelse

si=size(iima)

;normalization to mean int. of quiet regions

iino=float(iima)
wino=float(wima)
for i=0,si(3)-1 do begin
  iino(0,0,i)=iino(*,*,i)/mean(iino(qx1:qx2,qy1:qy2,i))
  wino(0,0,i)=wino(*,*,i)/mean(wino(qx1:qx2,qy1:qy2,i))
endfor

cwl=stdev(wino(qx1:qx2,qy1:qy2,*))
cir=stdev(iino(qx1:qx2,qy1:qy2,*))
print,'QUIET REGION:'
print,'WL:  contrast =',cwl
print,'IR:  contrast =',cir
fact=cwl/cir
print,'WL/IR: factor =',fact

dfac=(iino-1.)*fact-(wino-1.)

displ:
print,'Which frame do you want to display?  0 - ',strtrim(si(3)-1,2)
read,i
if i gt si(3)-1 then goto,fin
window,0,xsiz=660
OCONT,dfac(*,*,i),wino(*,*,i),[0,tp],/white
plots,[qx1,qx2,qx2,qx1,qx1],[qy1,qy1,qy2,qy2,qy1],color=255
goto,displ

fin:
END