PRO temper3,iima,wima,date,btemp1,btemp2

;Calculation of brightness teperatures
;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: btemp1 - maps of brightness temperature for WL
;	 btemp2 - maps of brightness temperature for IR

;PARAMETERS:
on_error,1

l1= 8000.		;wavelength WL (A)
l2=15542.		;wavelength IR (A)
b1=0.196		;abs. intensity at disk center for l1 (WL)
b2=0.043		;abs. intensity at disk center for l2 (IR)

c1=1.1903e20		;Planck function constants
c2=1.4384e8

;masking parameters

tp=0.90 	;WL threshold intensity to include pores
ta=0.028	;threshold for active region in smoothed dif.image

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

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

;brightness temperature (for WL)

br=wino*b1		;WL movie in absolute intensities
btemp1=c2/l1/alog(c1/(l1^5)/br+1)

;brightness temperature (for IR)

br=iino*b2		;IR movie in absolute intensities
btemp2=c2/l2/alog(c1/(l2^5)/br+1)


;MASKING for quiet, active, and pores regions

;quiet

bt1q=btemp1(qx1:qx2,qy1:qy2,*)
bt2q=btemp2(qx1:qx2,qy1:qy2,*)
mbt1q=mean(bt1q)
mbt2q=mean(bt2q)
print,'WL:  T_bright =',mbt1q,' K'
print,'IR:  T_bright =',mbt2q,' K'

;pores

bt1p=btemp1(where((wino lt tp) and (btemp2 lt 6400)))
bt2p=btemp2(where((wino lt tp) and (btemp2 lt 6400)))

;differential images

dwv=wino-2*iino+1.
for i=0,si(3)-1 do dwv(0,0,i)=smooth(dwv(*,*,i),5)
   ;smoothing by 5x5 pix boxcar

;active

bt1a=btemp1(where(dwv gt ta))
bt2a=btemp2(where(dwv gt ta))

;plotting

loadct,0
window,4,ysiz=700

plot,bt2q,bt1q,psym=3,xra=[5600,6800],yra=[5000,6300],col=0,back=255,  $
	/yst,xtitle='IR T_br (K)',ytitle='WL T_br (K)'
oplot,bt2a,bt1a,psym=3,col=100
oplot,bt2p,bt1p,psym=3,col=0
plots,[mbt2q-1000,mbt2q+200],[mbt1q-1000,mbt1q+200],col=100

;histograms
window,2

histy,bt1q,bin=5,min=5600,max=6400
histy,bt1a,bin=5,min=5600,max=6400,/over

print,'WL histogram. Press any key for IR histogram.'
g=get_kbrd(1)

histy,bt2q,bin=5,min=6000,max=6800
histy,bt2a,bin=5,min=6000,max=6800,/over

END