PRO figscaq,imov,wmov,date,num

;Plots a scatter plot of WL versus IR intensities of quiet,
;active, and pore regions in IR and WL frame no. 'num'.
;imov - 3D array with IR movie
;wmov - 3d array with WL movie
;date - 23 or 25

on_error,1
if (date ne 23) and (date ne 25) then  $
   message,'Sorry, I do not know how to handle this dataset'
dat=strtrim(date,2)
nus=strtrim(num,2)

si=size(imov)		;orig. size
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

;te=0.95	;IR threshold intensity to exclude pores
tp=0.88 	;WL threshold intensity to include pores
ta=0.028	;threshold for active region in smoothed dif.image
imini=0.8	;limits for plotting
wmini=0.6
imaxi=1.1
wmaxi=1.2

;extract images, remove apodized edges and convert to float

imai=float(imov(10:si(1)-11,10:si(2)-11,num))
imaw=float(wmov(10:si(1)-11,10:si(2)-11,num))

;normalization to mean int. of quiet regions

imai=imai/mean(imai(qx1:qx2,qy1:qy2))
imaw=imaw/mean(imaw(qx1:qx2,qy1:qy2))

;quiet regions

mqi=imai(qx1:qx2,qy1:qy2)
mqw=imaw(qx1:qx2,qy1:qy2)

;mask for pores

bmpor=imaw lt tp

;differential images

dwv=imaw-2*imai+1.
dwv=smooth(dwv,5)	;smoothing by 5x5 pix boxcar

;excluding mask for pores
;bmp=imai gt te
;mask for active regions, excluding pores
;bma=(dwv gt ta)*bmp

;mask for active regions, including some pores pixels

bma=dwv gt ta

;active regions

mai=imai*bma
maw=imaw*bma

;pore regions

mpi=imai*bmpor
mpw=imaw*bmpor

;plotting

loadct,3
window,4,ysiz=760

plot,mqi,mqw,psym=3,xra=[imini,imaxi],yra=[wmini,wmaxi],col=0,back=255,  $
	xtitle='IR Intensity',ytitle='WL Intensity'
oplot,mai,maw,psym=3,col=100
;oplot,mai,maw,psym=1,symsiz=0.5,col=100
oplot,mpi,mpw,psym=1,symsiz=0.5,col=0

;PS plotting

;loadct,0
;set_plot,'ps'
;device,file='scaq'+dat+nus+'.ps',yoffset=5,ysize=20
;
;plot,xia,(yia1+yia2)/2,psym=10,linestyle=1,thick=3,/xst,  $
;	title='INFRARED  1.5542 um',xtitle='Intensity',  $
;	ytitle='Normalized frequency',charsiz=1
;oplot,xiq,(yiq1+yiq2)/2,psym=10,linestyle=0,thick=1
;
;plot,xwa,(ywa1+ywa2)/2,psym=10,linestyle=1,thick=3,/xst,  $
;	title='WHITE LIGHT  0.8000 um',xtitle='Intensity',  $
;	ytitle='Normalized frequency',charsiz=1
;oplot,xwq,(ywq1+ywq2)/2,psym=10,linestyle=0,thick=1
;
;device,/close
;set_plot,'x'

END