pro taustar_fit

; v1 - 5jan09 - dhc - makes scatter plot with error bars for 16
;                     lines fit with non-porous wind profile model in zeta Pup
;
; 16 data points and *uneven* error bars are specified at the top of the file

; above comments are from taustar_scatter.pro, on which this program is built
;
; v1 - 7jan09 - dhc - fitting models to the ensemble taustar data
;                     first, fit a constant mode, calculate chisq
;                     next, read in the cmfgen opacity model (take code from
;                     opacity_plot.pro) and use it to generate a wavelength-
;                     dependent model of taustar; treating Mdot as a free parameter,
;                     we can scale this model to fit the ensemble taustar values
;

; PARAMETERS FOR TAUSTAR CALCULATION

Mdot = 8.3e-6
vinf = 2250.
Rstar = 18.6

mdot_msunyr = mdot
Mdot = Mdot*1.99e33/3.15e7
vinf = vinf*1.e5
Rstar = 18.6*6.96e10

print, Mdot, vinf, Rstar

; AXIS RANGES

y_low = 0.0
;y_upp = 4.7
y_upp = 6.

x_low = 4.5 
x_upp = 23.

; READING INPUT FILE WITH TRANSREAD

wave=0.0d
tau=0.0d
tau_err_down=0.0d
tau_err_up=0.0d
Ro=0.0d
Ro_err_down=0.0d
Ro_err_up=0.0d

transread,unit,wave,tau,tau_err_down,tau_err_up,$
              Ro,Ro_err_down,Ro_err_up,count=count,$
              filename='tau_Ro.dat',skiplines=2

;print,unit
;print,count
;print,wave,tau


; FITTING A CONSTANT TAUSTAR MODEL TO THESE 16 DATA POINTS
; USING EITHER THE UPPER OR LOWER ERROR BARS, DEPENDING ON WHETHER
; A POINT LIES ABOVE OR BELOW THE MODEL
;

model = fltarr(16)
chisq_best = 1.e10

for i = 0,200 do begin

chisq = 0.
for j = 0,15 do begin
	model[j] = 0.0 + i/50.
	if (model[j] LT tau[j]) then begin
		chisq = chisq + (model[j]-tau[j])^2/(tau_err_down[j])^2
	endif else begin
		chisq = chisq + (model[j]-tau[j])^2/(tau_err_up[j])^2
	endelse
endfor

;print,model[0],chisq

if (chisq LT chisq_best) then begin
	chisq_best = chisq
	best_tau = model[0]
endif 

endfor

print,' '
print,'chisq_best = ',chisq_best
print,'best_tau = ',best_tau


; READING IN OPACITY MODEL (CMFGEN)

opacity_read,wave_cmfgen,opacity_cmfgen

;print, 'cmfgen first entries =',wave_cmfgen[0],opacity_cmfgen[0]

; SETTING TAUSTAR BASED ON ASSUMED MDOT

taustar_model = opacity_cmfgen*mdot/4./!pi/vinf/Rstar
;print, 'wave_cmfgen,opacity_cmfgen,taustar_model',wave_cmfgen[9100],$
;              opacity_cmfgen[9100],taustar_model[9100]

; note: wavelength grid for cmfgen model is quite non-linear

; TESTING VALUES OF MDOT TO FIND THE BEST-FIT TAUSTAR MODEL

taustar_test = fltarr(10000)
model_mdot = fltarr(51)
chisq_best_mdot = 1.e10

for i = 0,50 do begin

; from 3.5e-6 to 4.5e-6
	mdot_fac = (1./8.3)*(3.5 + i/50.)
;print, 'i, mdot_fac = ',i,mdot_fac,mdot_fac*mdot_msunyr
	model_mdot[i] = mdot_fac*mdot_msunyr
; note: as long as the input mdot at the top is 8.3e-6, this tests mdot_-6 values of 1. to 5.

;	for j = 0,9999 do begin
;		taustar_test[j] = mdot_fac*taustar_model[j]
;	endfor

	taustar_test = mdot_fac*taustar_model

; note: we've got to find the 16 values in model that correspond to the wavelengths where we've
;       got taustar values from the data

; interpolate on wave_cmfgen,taustar_test

	taustar_points = interpol(taustar_test,wave_cmfgen,wave)

;print, tau,taustar_points

; calculate chisq
	chisq = 0.
	for k = 0,15 do begin
		if (taustar_points[k] LT tau[k]) then begin
			chisq = chisq + (taustar_points[k] - tau[k])^2 $
        	                /(tau_err_down[k])^2
		endif else begin
			chisq = chisq + (taustar_points[k] - tau[k])^2 $
        	                /(tau_err_up[k])^2
		endelse
	endfor
	
print,'mdot, chisq = ',model_mdot[i],chisq

if (chisq LT chisq_best_mdot) then begin
	chisq_best_mdot = chisq
	best_mdot = model_mdot[i]
	best_i = i
endif


endfor

print, ' '
print, 'chisq_best_mdot = ', chisq_best_mdot
print, 'best_mdot, ratio to mdot = ',best_mdot,best_mdot/mdot_msunyr
print, 'best_i = ',best_i


xtit = 'Wavelength ('+STRING("305B)+')'
ytit = '!9t!7!d*!n'


; PLOTTING (TO IDL SCREEN WINDOW)

for i=1,2 do begin

if (i EQ 1) then begin
window
!p.font=0
endif else begin
psopen,'taustar_fit'

!p.font=0
device,/times

endelse

psym_circ
plot,wave,tau,ps=8,thick=4,xr=[x_low,x_upp],yr=[y_low,y_upp],xstyle=1,ystyle=1,$
       ytit=ytit,xtit=xtit,chars=1.8,xthick=2,ythick=2,symsize=1.35

errplot,wave,tau-tau_err_down,tau+tau_err_up,thick=2

oplot,[0.,100.],[best_tau,best_tau],linestyle=3,thick=2
oplot,wave_cmfgen,taustar_model,linestyle=2,thick=2
oplot,wave_cmfgen,taustar_model*best_mdot/mdot_msunyr,linestyle=0,thick=2

endfor

if (i EQ 1) then begin
;print, i
endif else begin
psclose
endelse


end