log xspec_temp.log
set xs_return_result 1
# Fit X-ray spectral annuli in X-spec
# additional functionality to tie parameters together such as metallicity in annuli pairs to give the same fit


#Filename inputs

#####Check for input file
##put "Enter filename with comma separated data obsid, etc."
##gets stdin filename


# Inputs:
# List of Spectra - currently constructs based on prefix and suffix 
# Redshift
# nH in 10^22 cm^2
# Tx guess
# Fe guess
# * combine bin and unbinned and set stat and binning parameters


#set check 0  # if there was a typo doesn't crash
#while {$check<1} {
#puts "Enter Obs ID: "
#gets stdin obsid
#puts "Enter number of spectra to fit: "
#gets stdin numspec
#puts "Enter file prefix (e.g. ann): "
#gets stdin prefix
#puts "Enter file suffix (e.g. .ai.fits): "
#gets stdin suffix
##############Enter parameters
#puts "Enter redshift: "
#gets stdin red 
#puts "Enter nH column in units of 10^22 cm^2: "
#gets stdin nH
#puts "Enter Tx guess: "
#gets stdin Tx
#puts "Enter Abundance guess: "
#gets stdin Fe
#puts "The obsid is $obsid"
#puts "The number of spectra to be fit is $numspec"
#puts "The first file to be fit is $prefix\_$obsid\_0$suffix"
#puts "If this first file looks correct enter 1, otherwise enter 0"
#gets stdin check
#} 

#################

#Fit Xspec data
statistic cstat
weight standard
method simplex
setplot energy
setplot xlog on
setplot ylog on
# Keep going until fit converges.
query yes
# Open the file to put the results in.
set fileid [open fit_result.dat w]
cosmo 70 0 0.7
#Format is nH, kT, nH(fixed at 1.0), Abundance, redshift, switch, and normalize
abund grsa

set prefix ann
set obsid 11724
set suffix .pi
set nH 0.010
set Tx 7.1
set Fe 0.3
set red 0.213
set Tsoft 0.4
set Fesoft 0.3
set numspec 9

set groups "1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40"
#load all the annuli first so regions can be tied across
set y 1
setplot device $prefix\_$obsid\_bin-all.ps/ps
for {set x 0} {$x < $numspec} { incr x } {
	data $y:$y $prefix\_$obsid\_$x$suffix	
	set curr_group [string range $groups 0 [expr {2 * $x}]]
	puts "$curr_group"
	setplot group $curr_group
	plot data
	incr y
}

# focus on 0.7-7.0 keV range only need to constrain spectra, not groups	
ignore bad	
ignore 1-$numspec:**-0.7 7.0-** 

# input initial values for all the spectra phabs(mekal_Tx+mekal_softXraybg)
# this setup automatically ties all parameters across spectra
# total of 13 parameters
model phabs(mekal+mekal) & $nH & $Tx & 1.0 & $Fe,0 & $red & 0 & 0.01 & $Tsoft & 1.0 & $Fesoft,0 & $red & 0 & 0.01
# freeze: nH phabs, Tsoft. Rest are already frozen
freeze 1 
freeze 8
freeze 10
#untie: all Tx, all norm, central Fe, pair Fe after that
set paramlen 13
set Fe1 17
set Fe2 30
set norm1 20
set norm2 26
set Tx1 15
for {set x 0} {$x < [expr {$numspec - 1}]} { incr x } {
	untie $Tx1 $norm1 $norm2
	set norm1 [expr {$norm1 + $paramlen}]
	set norm2 [expr {$norm2 + $paramlen}]
	set Tx1 [expr {$Tx1 + $paramlen}]
	
	if {[expr {$x % 2}] == 0} {
		newpar $Fe2=$Fe1
		set Fe1 [expr {$Fe1 + $paramlen}]
		set Fe2 [expr {$Fe2 + $paramlen}]
	}
	
}
puts "Waiting here before fit. Press return to continue."
gets stdin empty

# fit the spectra
fit 10000 0.001
fit
# compute errors on all free parameters	
# output data to files

for {set x 0} {$x < $numspec} { incr x } {
	error [expr {1 + $x * $paramlen}]-[expr {13 + $x * $paramlen}]
	puts "#### Hey Spectra # $x"
	tclout param [expr {2 + $x * $paramlen}]
	puts "#### Hey tx $xspec_tclout"
	tclout param [expr {4 + $x * $paramlen}]
	puts "#### Hey fe $xspec_tclout"
	tclout param [expr {5 + $x * $paramlen}]
	puts "#### Hey z $xspec_tclout"
	tclout stat
	puts "#### Hey stat $xspec_tclout"
	tclout dof
	puts "Hey dof $xspec_tclout"
	setplot energy
	setplot xlog on
	setplot ylog on
	setplot 
	setplot device $prefix\_$obsid\_$x\_bin.ps/ps
	plot data resid
}
close $fileid
log none
file copy -force xspec_temp.log xspec_tx_fit_$obsid\_bin.log
file delete xspec_temp.log
exit