PRO temper,iima,wima,date,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.
;date = 23 or 25
;OUTPUT: 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.

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

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

;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'
print,' '

;pores

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

;differential images
;dwv=wino-2*iino+1.

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,12
window,4,ysiz=700

plot,btemp2,btemp1,psym=3,xra=[5600,6800],yra=[5000,6300],  $
	/yst,col=0,back=255,charsiz=1.4,  $
	xtitle='T!Dbr!N (K),  !4k = 1.6 l!3m', $
	ytitle='T!Dbr!N (K),  !4k = 0.8 l!3m'
;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)',charsiz=2.2
oplot,bt2q,bt1q,psym=3,col=30
oplot,bt2a,bt1a,psym=3,col=125
oplot,bt2p,bt1p,psym=3,col=75
plots,[mbt2q-1000,mbt2q+200],[mbt1q-1000,mbt1q+200],col=0

;histograms
window,2,xsiz=660

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

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

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

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

;WL/IR temperature ratio map
loadct,3

ntemp1=btemp1/mbt1q
ntemp2=btemp2/mbt2q

OCONT,ntemp1/ntemp2,wino,[0,tp]
OCONT,ntemp1/ntemp2,wino,[0,tp]

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

END