pro set_vth_bpow, vth_bpow_par, xs, ys

;+
; Name: pro set_vth_bpow
;
; Usage: set_vth_bpow, vth_bpow_par [, xs, ys]
;
; from 5 points selected graphically, 2 for the thermal,
; and 3 for the broken power law,
; below the break above the break, and below the break,
; return the standard thermal + broken power-law parameterization
; used in f_vth_bpow:
;	vth_bpow_par = [ a0, a1, a2, a3, a4, a5 ] 
;	a0 - em_49, emission measure units of 10^49
;	a1 - KT, plasma temperature in keV
;	a2 - normalization at Epivot, normally 50 keV
;	a3 - negative power-law index, a1>0, at energies less than the break
;	a4 - break energy
;	a5 - negative power-law index for energies greater than the break
;	NB   Only change thermal parameters if the first two points are within the plot box!
; routines called
;	point
;
; Inputs - XS, YS each arrays of three points chosen graphically
;
;	If XS or YS is undefined, then point is called to select the
;	three points.
; Outputs - Function returns the standard parameters for a thermal + 
;	broken power-law     - See F_VTH_BPOW
;-

@function_com

if n_elements(xs) ne 5 or n_elements(ys) ne 5 then begin
	print,'Select the starting values for a thermal/broken power-law'
	print,'	graphically'
        print,'Pick a point in the Low-energy, thermal range.
	PRINT,'SELECT THE FIRST TWO POINTS TO THE LEFT OF THE PLOT BOX TO 
	PRINT,'LEAVE THE INITIAL THERMAL PARAMETERS,',VTH_BPOW_PAR(0:1),' UNCHANGED.'
       	point,npoints=1, x0,y0
	print, x0,y0
	print,'Pick the transition from the thermal to the power law spectrum'
       	point,npoints=1, x1,y1
	print, x1,y1
	print,'Pick a point at an energy less than the break energy.'
	point,npoints=1, x2,y2
	print, x2,y2
	print,'Pick the break energy.'
	point,npoints=1, x3,y3
	print, x3,y3
	print,'Pick a point at an energy greater than the break energy.'
	point,npoints=1, x4,y4
	print, x4,y4
	xs = [x0,x1,x2,x3,x4]
	ys = [y0,y1,y2,y3,y4]
endif

alpha = 0.4        ; alpha is from Gaunt factor. This value is an approx
xmin = min(crange('x') )
wmore = where( xs(0:1) gt xmin, nmore)

;Only change thermal parameters if the first two points are within the plot box!
if nmore eq 2 then begin
	a1 = (xs(1)-xs(0)) / ( alog(ys(0)/ys(1)) + (1+ alpha) * alog(xs(0)/xs(1)) )
	temp0 = ys(0) * xs(0)^(1+alpha) * exp((xs(0)-Epivot)/a1)
	a0 = temp0 * a1^0.1 * exp(Epivot/a1)/1.3e3 * 1e-4  ;1e-4 is to get units right
endif else begin
	a1 = vth_bpow_par(1)
	a0 = vth_bpow_par(0)
endelse

a4 = xs(3)
a3 = alog( ys(2)/ys(3) ) / alog( a4/xs(2) )
a2 = ys(3) / (Epivot/xs(3))^a3
a5 = alog( ys(4)/ys(3) ) / alog( a4/xs(4) )

vth_bpow_par = [a0,a1,a2,a3,a4,a5]
print, 'Vth_bpow parameters are: '
print, vth_bpow_par

end