# Written by Phil Massey, Lowell Obs.  March 21, 2002
# Added support for sround
# Added support for finding the "best" fwhm for fitpsf
# some bug!
#Changed how we did mode April 27, 2005
#Changed how we did min
#Removed gain from list....
        procedure character6 (img,rdnoise,saturate,fname,timestamp)
        string img {prompt="Image"}
        real rdnoise {6., prompt="readnoise"}
        real saturate {30000., prompt="saturation"}
	string fname {"filter", prompt="filter key word"}
	string timestamp {"time-obs", prompt="time key word"}

	string *inlist

        begin

        string s1, s2, s3, s4, s5, s6, s7


real fwhmv,gainv,rdnoisev,saturatev,f1,fd1,fd2,fd3,rms,bot,top,round,best,val,num

# Get parameters
	s1=img
	rdnoisev=rdnoise
	saturatev=saturate
	s6=fname
	s7=timestamp
	imgets(s1,"gain")
	gainv=real(imgets.value)
	hedit(s1,"rdnoise",rdnoisev,add+,verify-,show+)
	hedit(s1,"saturate",saturatev,add+,verify-,show+)
	hedit(s1,"fname",s6,add+,verify-,show+)
	hedit(s1,"timestamp",s7,add+,verify-,show+)
# Now find mode:
	delete("temp")
#REPLACED MODAL DEFINTION
#	imstat(s1,fields="mode",lower=-500,upper=20000,format-,>"temp")
#
	imhist(s1,z1=-500,z2=20000,binwidth=1,nbins=INDEF,listout=yes,>"temp")
	delete("temp2")
	sort("temp",col=2,num+,rev+,>"temp2")
	inlist="temp2"
	i=fscan(inlist,f1)
	hedit(s1,"mode",f1,add+,verify-,show+)
	print "mode=",(f1)
	rms=sqrt(f1*gainv+rdnoisev**2)/gainv
	hedit(s1,"sigma",rms,add+,verify-,show+)
	print "rms=",(rms)
	val=f1-5*rms
        print "Final value for DATAMI=",(val)
	bot=val
	hedit(s1,"DATAMI",bot,add+,veri-,show+)
	top=saturatev
	hedit(s1,"DATAMA",top,add+,veri-,show+)
	delete("tempxy")
	print("running daofind in character")
daofind(s1,output="tempxy",sigma=rms,fwhm=4.0,thresh=30.0,datamin=bot,datamax=top,interact-,verify-,verbose-,filter=s6,obstime=s7)
flpr
	delete("tempround")
	txdump("tempxy","sround","yes",>"tempround")
	delete("temproundans")
	pickround
	inlist="temproundans"
	i=fscan(inlist,round)
	hedit(s1,"sround",round,add+,verify-,show+,update+)
# I'm sorry. This is dumb.  But apparently this is what we have to do in
# order to get a good value.
	print("running fitpsf in character:",(rms),(bot),(top))
# 2.0:
	delete("temprad2.0")
	print("PSF=2.0")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad2.0",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=2.0)
	print("done running fitpsf")
flpr
	delete("rad2.0")
	txdump("temprad2.0",fields="rsigma",expr="rsigma!=INDEF",>"rad2.0")
# 2.5:
	print("PSF=2.5")
	delete("temprad2.5")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad2.5",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=2.5)
flpr
	print("done running fitpsf")
	delete("rad2.5")
	txdump("temprad2.5",fields="rsigma",expr="rsigma!=INDEF",>"rad2.5")

# 3.0:
	print("PSF=3.0")
	delete("temprad3.0")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad3.0",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=3.0)
flpr
	print("done running fitpsf")
	delete("rad3.0")
	txdump("temprad3.0",fields="rsigma",expr="rsigma!=INDEF",>"rad3.0")
# 3.5:
        print("PSF=3.5")
	delete("temprad3.5")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad3.5",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=3.5)
flpr
        print("done running fitpsf")
        delete("rad3.5")
        txdump("temprad3.5",fields="rsigma",expr="rsigma!=INDEF",>"rad3.5")

# 4.0:
	print("PSF=4.0")
	delete("temprad4.0")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad4.0",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=4.0)
flpr
	print("done running fitpsf")
	delete("rad4.0")
	txdump("temprad4.0",fields="rsigma",expr="rsigma!=INDEF",>"rad4.0")
# 4.5:
	print("PSF=4.5")
	delete("temprad4.5")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad4.5",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=4.5)
	print("done running fitpsf")
flpr
	delete("rad4.5")
	txdump("temprad4.5",fields="rsigma",expr="rsigma!=INDEF",>"rad4.5")
# 5.0:
	print("PSF=5.0")
	delete("temprad5.0")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad5.0",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=5.0)
flpr
	print("done running fitpsf")
	delete("rad5.0")
	txdump("temprad5.0",fields="rsigma",expr="rsigma!=INDEF",>"rad5.0")
# 5.5:
	print("PSF=5.5")
	delete("temprad5.5")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad5.5",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=5.5)
flpr
	print("done running fitpsf")
	delete("rad5.5")
	txdump("temprad5.5",fields="rsigma",expr="rsigma!=INDEF",>"rad5.5")
# 6.0:
	print("PSF=6.0")
	delete("temprad6.0")
fitpsf(s1,10,sigma=rms,coords="tempxy",output="temprad6.0",datamin=bot,datamax=top,interact-,verify-,verbose-,cache+,filter=s6,obstime=s7,fwhmpsf=6.0)
flpr
	print("done running fitpsf")
	delete("rad6.0")
	txdump("temprad6.0",fields="rsigma",expr="rsigma!=INDEF",>"rad6.0")

#Now pick one
	delete("tempone")
	pickone
	inlist="tempone"
	i=fscan(inlist,best)
#Now do it one last time for real:
	print("BEST PSF=",(best))
	delete("temprad")
fitpsf(s1,9,sigma=rms,coords="tempxy",output="temprad",datamin=bot,datamax=top,interact-,verify-,verbose-,filter=s6,obstime=s7,fwhmpsf=best)
flpr
	print("done running fitpsf")
	delete("allrad")
	txdump("temprad",fields="rsigma",expr="rsigma!=INDEF",>"allrad")
	delete("tempfwhm")
        radmode
flpr
	inlist="tempfwhm"
	i=fscan(inlist,fwhmv)
	hedit(s1,"fwhm",fwhmv,add+,ver-,show+)
end