;=======================================================================
; NAME :	calparm1 (procedure)
; PURPOSE :
;	Determine the calibration parameters, wl and scatter, for T1 data.
; CATEGORY :
;	idl/lib/flare
; CALLING SEQUENCE :
;	calparm1,img,q,u,v,scatter,wl,kdpfact=kdpfact
; INPUTS :
;	img     : image (*,*) or i,q,u,v as (*,*,4)
;       q,u,v   : polarization degree*10000
; OPTIONAL INPUT PARAMETERS : 
;	none
; KEYWORD PARAMETERS :
; OUTPUTS :
;	none
; COMMON BLOCKS : 	none
; SIDE EFFECTS : 	none
; RESTRICTIONS :	none
; PROCEDURE :
;	* If imin/imax < 0.17, it assumes that the data is limb image 
;	  including sky and scatter=imin, Bl(max) or Bt(max) is saturated.
;	* If imin/imax > 0.17, it assumes that 
;		scatter=(1-0.7)*imin*imax/(imax-0.7*imin), 
;	  Bl(max) or Bt(max) is saturated.
; MODIFICATION HISTORY :
;	K.I. '92/10/28	change order of Bl & bx,by
;	K.I. '92/11/20	change algorithm
;----------------------------------------------------------------------
pro calparm1,img,q,u,v,scatter,wlv,wlp,kdpfact=kdpfact,vmaxp=vmaxp,pmaxp=pmaxp

wlrange=[-30.,-130.]	; wavelength range allowed in mA
			; satulation level at -120 mA is
			;	Bl = 2620 G
			;	Bt = 2560 G

if not keyword_set(kdpfact) then kdpfact=1.

getcalib1,vmax=vmax,qmax=pmax,dwl=dwl
rr=vmax/pmax
nn=n_elements(vmax)
n1=min(where(dwl*1000. le wlrange(0)))
n2=max(where(dwl*1000. ge wlrange(1)))

;plot,dwl,vmax,title='saturated level of V & Q',xtitle='dwl (mA)'
;oplot,dwl,pmax
;plot,dwl,rr,title='Vmax/Qmax',xtitle='dwl (mA)'

s=size(v)	&	nx=s(1)	&	ny=s(2)
imax=max(img)
imin=min(img>0,pos)
r=float(imin)/float(imax)
p=fix(sqrt(float(q)^2+float(u)^2))

if r lt 0.17 then begin	; <===  limb image including sky
	print,'This is limb image....'
	scatter=imin 
	ii=where(img gt imax*0.17)
	imgd=img(ii)
	vmax_o=max(abs(v(ii)),pos)
	iv=imgd(pos)
	y=ii(pos)/nx	&	x=ii(pos)-y*nx	&	vmaxp=[x,y]
	pmax_o=max(p(ii),pos)
	ip=imgd(pos)
	y=ii(pos)/nx	&	x=ii(pos)-y*nx	&	pmaxp=[x,y]
endif else begin	; <===  disk image without sky
	print,'This is disk image....'
	sctfact=0.7
	scatter=fix((1.-sctfact)*imin*imax/(imax-sctfact*imin))
	vmax_o=max(abs(v),pos)
	iv=img(pos)
	y=pos/nx	&	x=pos-y*nx	&	vmaxp=[x,y]
	pmax_o=max(p,pos)
	ip=img(pos)
	y=pos/nx	&	x=pos-y*nx	&	pmaxp=[x,y]
endelse

vmax_o=vmax_o*float(iv)/(iv-scatter)
pmax_o=pmax_o*float(ip)/(ip-scatter)
zv=zero(vmax_o-vmax(n1:n2)*10000.)
zp=zero(pmax_o-pmax(n1:n2)*10000.)
case zv of
    -1: wlv=dwl(n1)*1000.
    -2: wlv=dwl(n2)*1000.
    else: wlv=(dwl(n1+fix(zv))*(1.-zv+fix(zv))	$
		+dwl(n1+fix(zv)+1)*(zv-fix(zv)))*1000.
endcase
case zp of
    -1: wlp=dwl(n1)*1000.
    -2: wlp=dwl(n2)*1000.
    else: wlp=(dwl(n1+fix(zp))*(1.-zp+fix(zp))	$
		+dwl(n1+fix(zp)+1)*(zp-fix(zp)))*1000.
endcase
	

print,format='("vmax=",i5,"  at (",i3,",",i3,")  I=",f5.3,"   dwl=",f6.1)', $
	vmax_o,vmaxp(0),vmaxp(1),float(iv)/imax,wlv
print,format='("pmax=",i5,"  at (",i3,",",i3,")  I=",f5.3,"   dwl=",f6.1)', $
	pmax_o,pmaxp(0),pmaxp(1),float(ip)/imax,wlp

;wl=min([wlv,wlp])
wl=wlv
rsct=float(scatter)/imax*100.
print,scatter,rsct,wl,	$
	format='("scatter=",I5," (",I3,"% of imax )    dwl=",f6.1," mA")'

end