pro norh_polariz,indexa,dataa,indexs,datas,indexpol,polariz,mvalid $
     ,datacri=datacri,tbinput=tbinput

;+
; NAME:
;       norh_polariz
;
; PURPOSE:
;       Derive circular polarization degree
;
; CALLING SEQUENCE:
;       norh_polariz,indexa,dataa,indexs,datas,polariz,mvalid
;
; INPUTS:
;       indexa: NoRH index of 17 GHz R+I images
;       dataa:  2-D or 3-D array of 17 GHz R+I images
;       indexs: NoRH index of 17 GHz R-I images
;       datas:  2-D or 3-D array of 17 GHz R-I images
;
; OUTPUTS:
;      polariz  : 2-D or 3-D array of the circular polarization degree
;                  Tb[R-L]/Tb[R+L]
;      mvalid  : 2-D or 3-D array of valid pixels
;                1 for valid, 0 for non-valid
;
; OPTIONAL INPUTS:
;      datacri   ; criterion to judge the validity. default is 3.e4 K
;                if both dataac datazc > datacri then valid.
;
; OPTIONAL OUTPUTS:
;
; HISTORY:
;       Writtern 2000-7-24 T. Yokoyama
;       Revise 2000-8-24 T. Yokoyama add index
;
;-

if not keyword_set(datacri) then datacri=3.e4
facri=norh_tb2flux(datacri,indexa,/intensity)

if keyword_set(tbinput) then facri=3.e4

sza=size(dataa)
szs=size(datas)

dima=sza(1:2)
dims=szs(1:2)
if ((dims/dima)(0) ne (dims/dima)(1)) then message,'dimension inconsistent'
if (sza(0) eq 3) then if (sza(3) ne szs(3)) then $
   message,'dimension inconsistent'

if (sza(0) eq 3) then nx=sza(3) else nx=1
ix=sza(1)
jx=sza(2)

polariz=fltarr(ix,jx,nx)

;for n=0,nx-1 do begin
;   polariz(*,*,n)=datas(*,*,n)/(dataa(*,*,n)>1)
;endfor

polariz=datas/(dataa>1)

mvalid=bytarr(ix,jx,nx)+1b
whr=where((dataa lt facri),count)
if (count ge 1) then mvalid(whr)=0b

indexpol=indexs
indexpol.norh.polariz='pol'


return
end