;**************************************************************************
;NAME       : getcalib4 (procedure)
;FUNCTION   : get caliblation data from the file written by mkcaldat.pro
;		    V     <-> B     <--- VBl4.dat
;		    S     <-> Vel   <--- SVel.dat
;MODIFICATION :
;	K.Ichimoto 	'92/10/28
;	K.Ichimoto 	'92/11/21	; extra keyword
;	K.Ichimoto 	'93/02/02	; get from '/data/ftobs/calib/'
;==========================================================================
pro getcalib4,Bl,Vel,dwl=dwl,bmag=bmag,vmax=vmax,smax=smax,extra=extra
; Bl	: longitudinal mag.field in Gauss
; Vel	: Doppler shift in m/s
; dwl   : expected wave length position from line center in mA for V-B relation
; bmag  : expected mag.field strength in G for S-Dop relation
; vmax  : return maximum V without satulation
; smax  : return maximum S without satulation
; extra : extrapolate curve by linear

;dir='/usr/local/lib/idl/lib/flare/'
dir='/data/ftobs/calib/'
file=''
L6337:	
	;--------  VBl4.dat --->   V ==> Bl     
	get_lun,Unit
	openr,Unit,dir+'VBl4.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 ijn 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))	; fpos decrease with i
		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 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

	;--------  SVel.dat --->   S <-> Vel     
	get_lun,Unit
	openr,Unit,dir+'SVel.dat'
	nn=0	&	nb=0
	readf,Unit,format='(I5,I5,F10.3,f10.3,A100)',nn,nb,gam,xai,file
	nbl=fltarr(nb)		; magnetic field strength in Gauss
	smax=fltarr(nb)		; maximum S without satulation
	s=intarr(nn)		; Doppler signel *10000
	readu,Unit,nbl		; magnetic field strength in Gauss, must >0
	readu,Unit,smax
	readu,Unit,s
	if n_params() eq 2 then begin
	if n_elements(bmag) then begin
		bmax=max(abs(nbl))
		if (bmax-bmag)*(-bmax-bmag) gt 0 then begin
			print,'bmag out of range !!'
			return
		endif
		if bmag lt 0 then bmag=-bmag
		i1=max(where(nbl le bmag))	; nbl increase with i
		point_lun,-Unit,pos
		point_lun,Unit,pos+long(nn)*2l*i1
		Vel1=intarr(nn)
		Vel2=intarr(nn)
		readu,Unit,Vel1
		readu,Unit,Vel2
		i2=i1+1
		Vel = fix( ((abs(nbl(i1)-bmag))*Vel2 	$
		   + (abs(nbl(i2)-bmag))*Vel1)/abs(nbl(i1)-nbl(i2)) )
		smax = ((abs(nbl(i1)-bmag))*smax(i2) 	$
		   + (abs(nbl(i2)-bmag))*smax(i1))/abs(nbl(i1)-nbl(i2))
		if keyword_set(extra) then begin
			nn=n_elements(Vel)
			grad=(Vel(nn/2+5)-Vel(nn/2-5))/10.
			nmax=min(where(Vel eq max(Vel)))
			nmax=nmax-100
			Vel(nmax:nn-1)=Vel(nmax)+fix(indgen(nn-nmax)*grad)
			over=where(Vel(nmax:nn-1) lt 0,count)
			if count gt 0 then Vel(nmax+over)=max(Vel)
			nmin=max(where(Vel eq min(Vel)))
			nmin=nmin+100
			Vel(0:nmin)=Vel(nmin)-fix((nmin-indgen(nmin+1))*grad)
			over=where(Vel(0:nmin) gt 0,count)
			if count gt 0 then Vel(over)=min(Vel)
		endif
	endif else begin
		Vel=intarr(nn,nb)	; polarization degree
		readu,Unit,Vel
	endelse
	endif
	close,Unit
	free_lun,unit

	if n_params() eq 0 and n_elements(dwl) eq 0 then dwl=fposl
	if n_params() eq 0 and n_elements(bmag) eq 0 then bmag=nbl

end