;**************************************************************************
;NAME       : getcalib1 (procedure)
;FUNCTION   : get caliblation data from the file written by mkcaldat.pro
;		    V     <-> B    <--- VBl.dat
;		(Q^2+U^2) <-> Bt   <--- QUBt.dat
;MODIFICATION
;	K.Ichimoto 	'92/10/26
;	K.Ichimoto 	'92/10/28
;	K.Ichimoto 	'92/11/21	; extra keyword
;	K.Ichimoto 	'93/02/02	; get from '/data/ftobs/calib/'
;	K.Ichimoto 	'95/06/01	; caldir keyword
;==========================================================================
pro getcalib1,Bl,Bt,dwl=dwl,vmax=vmax,qmax=qmax,extra=extra,plot=plot, $
	caldir=caldir
; Bl	: longitudinal mag.field in Gauss
; Bt	: tangential mag.field in Gauss
; dwl   : expected wave length position from line center in mA
; vmax  : return maximum V before satulation
; qmax  : return maximum Q before satulation
; extra : extrapolate curve by linear
; plot  : plot the curves

;dir='/usr/local/lib/idl/lib/flare/'
dir='/data/ftobs/calib/'
if keyword_set(caldir) then dir=caldir
file=''
L6303:	
	get_lun,Unit
	openr,Unit,dir+'VBl1.dat'
	readf,Unit,format='(I5,I5,F10.3,f10.3,A100)',nn,npos,gam,xai,file
	fposl=fltarr(npos)	; wave length position from line center in A
	vmax=fltarr(npos)	; maximum V without satulation
	v=intarr(nn)		; magnetic field strength in Gauss
	readu,Unit,fposl	; wavelength position, must <0
	readu,Unit,vmax
	readu,Unit,v
	if n_params() eq 2 then begin
	if n_elements(dwl) then begin
		dwl0=float(dwl)/1000.
		dwlmax=max(abs(fposl))
		if (dwlmax-dwl0)*(-dwlmax-dwl0) gt 0 then begin
			print,'dwl out of range !!'
			return
		endif
		if dwl gt 0 then dwl0=-dwl0
		i1=max(where(fposl ge dwl0))
		point_lun,-Unit,pos
		point_lun,Unit,pos+long(nn)*2l*i1
		Bl1=intarr(nn)
		Bl2=intarr(nn)
		readu,Unit,bl1
		readu,Unit,bl2
		i2=i1+1
		Bl = fix( ((abs(fposl(i1)-dwl0))*Bl2 	$
		   + (abs(fposl(i2)-dwl0))*Bl1)/abs(fposl(i1)-fposl(i2)) )
		vmax = ((abs(fposl(i1)-dwl0))*vmax(i2) 	$
		   + (abs(fposl(i2)-dwl0))*vmax(i1))/abs(fposl(i1)-fposl(i2))
		if dwl gt 0 then bl=-bl
		if keyword_set(extra) then begin
			nn=n_elements(bl)
			grad=(bl(nn/2+5)-bl(nn/2-5))/10.
			nmax=min(where(bl eq max(bl)))
			nmax=nmax-100
			bl(nmax:nn-1)=bl(nmax)+fix(indgen(nn-nmax)*grad)
			over=where(bl(nmax:nn-1) lt 0,count)
			if count gt 0 then bl(nmax+over)=max(bl)
			nmin=max(where(bl eq min(bl)))
			nmin=nmin+100
			bl(0:nmin)=bl(nmin)-fix((nmin-indgen(nmin+1))*grad)
			over=where(bl(0:nmin) gt 0,count)
			if count gt 0 then bl(over)=min(bl)
			print,'calibration curve extrapolated'
		endif
	endif else begin
		Bl=intarr(nn,npos)	; polarization degree
		readu,Unit,Bl
	endelse
	endif
	close,Unit
	free_lun,unit

	get_lun,Unit
	openr,Unit,dir+'QUBt.dat'
	readf,Unit,format='(I5,I5,F10.3,f10.3,A100)',nn,npos,gam,xai,file
	fpost=fltarr(npos)	; wave length position from line center ijn A
	qmax=fltarr(npos)	; maximum Q without satulation
	q=intarr(nn)		; magnetic field strength in Gauss
	Bt=intarr(nn,npos)	; polarization degree
	readu,Unit,fpost	; wavelength position, must <0
	readu,Unit,qmax
	readu,Unit,q
	if n_params() eq 2 then begin
	if n_elements(dwl) then begin
		dwl0=float(dwl)/1000.
		dwlmax=max(abs(fpost))
		if (dwlmax-dwl0)*(-dwlmax-dwl0) gt 0 then begin
			print,'dwl out of range !!'
			return
		endif
		if dwl gt 0 then dwl0=-dwl0
		i1=max(where(fpost ge dwl0))
		point_lun,-Unit,pos
		point_lun,Unit,pos+long(nn)*2l*i1
		Bt1=intarr(nn)
		Bt2=intarr(nn)
		readu,Unit,Bt1
		readu,Unit,Bt2
		i2=i1+1
		Bt = fix( ((abs(fpost(i1)-dwl0))*Bt2 	$
		   + (abs(fpost(i2)-dwl0))*Bt1)/abs(fpost(i1)-fpost(i2)) )
		qmax = ((abs(fposl(i1)-dwl0))*qmax(i2) 	$
		   + (abs(fposl(i2)-dwl0))*qmax(i1))/abs(fposl(i1)-fposl(i2))
		if keyword_set(extra) then begin
			nn=n_elements(bt)
			nmax=min(where(bt eq max(bt)))
			div=( shift(float(bt(0:nmax-10)),-10)	$
				-float(bt(0:nmax-10)) )/10.
			nmax=nmax-100
			grad=min(div(0:nmax))
			bt(nmax:nn-1)=bt(nmax)+fix(indgen(nn-nmax)*grad)
			over=where(bt(nmax:nn-1) lt 0,count)
			if count gt 0 then bt(nmax+over(0):nn-1)=max(bt)
		endif
	endif else begin
		Bt=intarr(nn,npos)	; polarization degree
		readu,Unit,Bt
	endelse
	endif
	close,Unit
	free_lun,unit

	if n_params() eq 0 and n_elements(dwl) eq 0 then dwl=fposl
	if keyword_set(plot) then begin
		title='Calibration Curves for T1 data  : wave length='+	$
			string(dwl,format='(I3)')+' mA'
		plot,(v/100.)>0,(bl/10.)>0,title=title,	$
			xtitle='poralization degree (%)',	$
			ytitle='magnetic field (G)'
		oplot,q/100.,bt/10.
	endif


end