#This routine creates a flat-field for filter "f"
#and uses it to flatten the object frames for the same filter
#It also saves the  imstats in "FlatStats"
#Added better error checking and logic 2/5/06
#Also gets rid of skys before combining flats...
#and moves the flattened skys at the end
procedure  y4kflat(f)
       string f {prompt='filter'}
begin
       string filt,g
       real x
       int inum,inumo,inumf,i
       filt=f
#Lets deal with any skys before we go looking for flats...
#Identify skys:
delete("skys")
delete("gskys")
delete("bskys")
hselect("Zy*%.fits%%","$I","(filterid='"//filt//"' &&((i_title ?='SKY')||(i_title ?='Sky')||(i_title ?='sky')||(i_title ?='even')||(i_title ?= 'twi')||(i_title ?='Twi')||(i_title ?='morn')))",>"skys")
count("skys") | scan (i)
print("Number of skys"//(i))
if (i>0) {
hedit("@skys","imgtype","OBJECT",up+,ver-,show+)
#Get rid of saturated skies:
list="skys"
while (fscan(list,s1) !=EOF) {
imstat(s1//"[500:1500,500:1500]",fields="midpt",format-,cache+) | scan (x)
if ( (x) < 32000 && (x)>2000)
      files(s1,>>"gskys")
  else
      files(s1,>>"bskys")
}
#Move the bad skys to junk
if (access("bskys")) {
set g="../junk/"
imdelete("g$//@bskys")
imrename("@bskys","g$//@bskys")
}
else
print("No bad skys...")
}
else
print("No skys...")
#OK, now MAKE the flat for filter "f"
#bad flats:
delete("bflats")
#good flats:
delete("gflats")
#all flats
delete("allflats")
hselect("Zy*.fits","$I","imgtype='FLAT'&&filterid='"//filt//"'&&exptime>10",>"allflats")
count("allflats") | scan(inum)
print("inum=",(inum))
if(inum>0) {
#Check to make sure that none of these are saturated...
list="allflats"
while (fscan(list,s1) !=EOF) {
imstat(s1//"[500:1500,500:1500]",fields="midpt",format-,cache+) | scan (x)
if ( x<32000 && x>2000 ) 
    files(s1,>>"gflats")
    else 
    files(s1,>>"bflats")
    }
    }
    else
    print("No flats of type",(filt))
#Get rid of bad flats
if (access("bflats")) {
set g="../junk/"
imdelete("g$//@bflats")
imrename("@bflats","g$//@bflats")
    }
inumf=0
if (access("gflats")) {
count("gflats") | scan (inumf)
    }
#Insist that here be at least 3 good flats:
if(inumf>2) {
imdelete("Flat"//filt//".fits")
#Make the Flat
imcombine("@gflats","Flat"//filt//".fits",combine="average",reject="crreject",scale="median",statsec="[500:1500,500:1500]",rdnoise=12,gain=1.5,>>"../FlatStats")
}
else 
print("There are not enough flats of type ",(filt)," to combine...checking for objs")
#
delete("objs")
#Identify objects
#So, we should be including skys in this...but not the bad skys.
hselect("%Z%%y*.fits","$I","imgtype='OBJECT'&&filterid='"//filt//"'",>"objs")
count("objs") | scan (inumo)
if (inumo>0) {
#Make sure that there is flat-fielding capability
if(inumf<3) {
#There IS no flat field!
#Find out where we are...
!pwd
!pwd >> "../ERRORLOG"
print("Objects but no flat field for filter",(filt))
print("Objects but no flat field for filter",(filt),>>"../ERRORLOG")
}
else {
#Do the flat fielding...
#First, capture the median in the [500:1500,500:1500] section:
imstat("Flat"//filt//".fits[500:1500,500:1500]",fields="midpt",format-,cache+) | scan(x)
delete("Flat"//filt//"n.fits")
imarith("Flat"//filt//".fits","/",(x),"Flat"//filt//"n.fits",verb+)
#Delete the flatttened images, if they exist:
imdelete("F//@objs")
imarith("Z//@objs","/","Flat"//filt//"n.fits","F//@objs",verb+)
#Now move the skys somewhere else...
if (access("gskys")) {
set g="../allskys/"
imdelete("g$F//@gskys")
imrename("F//@gskys","g$F//@gskys")
}
}
}
end