;**************************************************************************
; NAME       : quv2b (procedure)
; PURPOSE    : convert pol.deg. to mag.strength (Gauss) for T1
; CALLING SEQUENCE :
;	quv2b, q, u, v, bx, by, bl, Bl_cal=Bl_cal, Bt_cal=Bt_cal, dwl=dwl,
;		scatter=scatter,kdpfact=kdpfact
; INPUTS :
;     q,u,v     -- pol. degree *10000 (q,u,v already devided by img)
; KEYWORD PARAMETERS :
;     Bl_cal    -- calibration function of Bl,    p(*10000) => B(G*10)
;     Bt_cal    -- calibration function of Bt,    p(*10000) => B(G*10)
;     dwl       -- wave length position of Lyot filter from line center (mA)
;			if dwl is not specified, -80 mA is assumed
; OUTPUTS :
;     bx,by,bl  -- magnetic field (transversal & long.) in Gauss*10
; 
; MODIFICATION HISTORY :
;	K.Ichimoto, 1992/10
;	K.Ichimoto, 1992/10/25	rename to f1calib
;	K.Ichimoto, 1992/10/28	change order of Bl and Bx,By
;	K.Ichimoto, 1992/11/11	if v is scalar then return, rename of f1calib
;	K.I. '92/11/19	change v sign
;	K.I. '92/11/21	getcalib1,/extra
;	K.I. '92/05/18	q,u,v devided by polfact=2
;==========================================================================
;
pro quv2b,q,u,v,bx,by,bl,Bl_cal=Bl_cal,Bt_cal=Bt_cal,dwl=dwl

METHOD=1
blmodif=0	; if 1, make modification of Bl magnitude to barance with Bt
polfact=1	; pol.degree = q,u,v /10000./polfact

case METHOD of

    0:	begin	; ===== simple linear relation
	lfact=1000.*10./(0.1*10000.)	;   Bl=1000 G    <-->  V=0.1
	tfact2=4e5*100./(0.025*10000.)	; Bq^2=4e5 G^2  <-->  Q=0.025
	bl=fix(v*lfact)
	pl=sqrt(float(q)^2+float(u)^2)
	bt=fix(sqrt(pl*tfact2))
	ii=where( q eq 0 and u eq 0 , count )
	if count ne 0 then q(ii)=1
	;th=fix(atan(u,q)/2.*1000.)
	th=atan(u,q)/2.
	bx=fix(bt*cos(th))
	by=-fix(bt*sin(th))
	end

    1:	begin	; ===== use theoretical relation
	if not keyword_set(dwl) then dwl=-80
	print,'quv2b: '
	if not keyword_set(Bl_cal) or not keyword_set(Bl_cal) then begin
		getcalib1,Bl_cal,Bt_cal,dwl=dwl,/extra
		print,'calibration table loaded,  dwl=',dwl,' mA'
	endif
	nn2=n_elements(Bl_cal)/2
	print,'getting Bx, By & Bl ...'
	bl=Bl_cal( (-v/polfact+nn2) >0<(n_elements(Bl_cal)-1) )
	pl=float(q/polfact)^2+float(u/polfact)^2
	pl=fix(sqrt(pl>0))
	bt=Bt_cal( pl >0<(n_elements(Bt_cal)-1) )
	ii=where( q eq 0 and u eq 0 , count )
	if count ne 0 then q(ii)=1
	th=atan(u,q)/2.
	bx=fix(bt*cos(th))
	by=-fix(bt*sin(th))
	; ----- get potential field -----
	s=size(bl)
	nx=s(1)	&	ny=s(2)
	if nx eq 1 or ny eq 1 then return
	print,'getting potential field...'
	pfact=nx/64
	cff1n,rebin(bl,nx/pfact,ny/pfact),bxp,byp
	bxp=rebin(bxp,nx,ny)
	byp=rebin(byp,nx,ny)
	; ----- correct the Bl magnitude ------
	if blmodif eq 1 then begin
		btp=fix(sqrt(float(bxp)^2+float(byp)^2))
		costh=(float(bx)*bxp+float(by)*byp)/float(bt>1)/float(btp>1)
		btmax=max(bt)
		ii=where( (abs(costh) gt 0.8) and (bt gt btmax*0.7) )
		blfact=total(float(bt(ii))/btp(ii))
		blfact=blfact/n_elements(ii)
		print,blfact,'   is multiplied on Bl'
		bl=fix(bl*blfact)
	endif
	; ----- vectorize -----
	print,'setting Bt direction...'
	flip = where( bx*bxp + by*byp lt 0.0, count )
	if count gt 0 then begin
		bx( flip ) = -bx( flip )
		by( flip ) = -by( flip )
	end
	end
endcase

end