PRO tempernew,date,dtemp,btemp1,btemp2,bt1q,bt2q,bt1a,bt2a,bt1p,bt2p

;Calculation of brightness temperatures
;from a couple of registered frames 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.
;NEW: PLOT Tbr(IR) vs. Tbr(WL), MAP OF TEMPERATURE DIFFERENCES
;	INSTEAD OF TEMPERATURE RATIOS
;INPUT date = 23 or 25 (obligatory)
;OUTPUT: dtemp  - map of temperature differences
;	 btemp1 - map of brightness temperature for WL
;	 btemp2 - map of brightness temperature for IR
;	 bt*    - brightness temperatures of:
;		  1- WL, 2- IR; q- QR, a- AR, p- pore.

on_error,1
print,'OPT.OUTPUT: dt_ima,t8_ima,t16_ima,t8q,t16q,t8a,t16a,t8p,t16p'
if (date ne 23) and (date ne 25) then  $
   message,'Sorry, I do not know how to handle this dataset'

;INPUT:

if date eq 23 then begin
	restore,'besti23a'
	restore,'bestw23a'
	iima=reform(besti23a(*,*,2))
	wima=reform(bestw23a(*,*,2))
endif else begin
	restore,'besti25a'
	restore,'bestw25a'
	iima=reform(besti25a(*,*,2))
	wima=reform(bestw25a(*,*,2))
endelse

;PARAMETERS:

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
ta=(-0.034)

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

;normalization to mean int. of quiet regions

iino=float(iima)
wino=float(wima)
iino=iino/mean(iino(qx1:qx2,qy1:qy2))
wino=wino/mean(wino(qx1:qx2,qy1:qy2))

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

;brightness temperature (for WL)

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

;brightness temperature (for IR)

br=iino*b2		;IR frame 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:  mean   T_bright =',mbt1q,' K'
print,'WL:  median T_bright =',median(bt1q),' K'
print,'IR:  mean   T_bright =',mbt2q,' K'
print,'IR:  median T_bright =',median(bt2q),' K'
refd=mbt2q-mbt1q
print,'mean IR - WL differ. =',refd,' K'
print,' '

;pores

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

;differential images

fact=cwl/cir
dwv=(iino-1.)*fact-(wino-1.)
dwv=smooth(dwv,5)	;smoothing by 5x5 pix boxcar

;active

bt1a=btemp1(where((dwv lt ta) and (wino gt tp)))
bt2a=btemp2(where((dwv lt ta) and (wino gt tp)))

print,'ACTIVE REGION:'
print,'WL:  mean   T_bright =',mean(bt1a),' K'
print,'WL:  median T_bright =',median(bt1a),' K'
print,'IR:  mean   T_bright =',mean(bt2a),' K'
print,'IR:  median T_bright =',median(bt2a),' K'
print,' '

;plotting

loadct,0
window,4,ysiz=700

plot,bt1q,bt2q,psym=3,xra=[5000,6300],yra=[5600,6800],col=0,  $
;plot,btemp1,btemp2,psym=3,xra=[5000,6300],yra=[5600,6800],col=0,  $
	/xst,/yst,back=255,charsiz=1.4,  $
	xtitle='T!Dbr!N (K),  !4k = 0.80 l!3m', $
	ytitle='T!Dbr!N (K),  !4k = 1.55 l!3m'
;oplot,bt1q,bt2q,psym=3,col=30
oplot,bt1a,bt2a,psym=3,col=100	;col=125
oplot,bt1p,bt2p,psym=3,col=0	;col=75
plots,[5600-refd,6800-refd],[5600,6800],col=0

;IR - WL temperature differences map

dtemp=btemp2-btemp1	;difference map
dtemp(0,0)=refd-200	;setting scale
dtemp(1,0)=refd+200

window,2,xsiz=660
OCONT,dtemp,wino,[0,tp]
OCONT,dtemp,wino,[0,tp],/white

window,1,xsiz=388,ysiz=300
tvscl,dtemp

;print,'CTs: for scatter plot 12, for map 3'

END