from astwro.sampledata import fits_image
frame = fits_image()

"""
Explanation: Deriving a Point-Spread Function in a Crowded Field
following Appendix III of Peter Stetson's User's Manual for DAOPHOT II
Using pydaophot form astwro python package
All italic text here have been taken from Stetson's manual.
The only input file for this procedure is a FITS file containing reference frame image. Here we use sample FITS form astwro package (NGC6871 I filter 20s frame). Below we get filepath for this image, as well as create instances of  Daophot and Allstar classes - wrappers around daophot and allstar respectively.
One should also provide daophot.opt, photo.opt and allstar.opt in apropiriete constructors. Here default, build in, sample, opt files are used.
End of explanation
"""


from astwro.pydaophot import Daophot, Allstar
dp = Daophot(image=frame)
al = Allstar(dir=dp.dir)

"""
Explanation: Daophot object creates temporary working directory (runner directory), which is passed to Allstar constructor to share.
End of explanation
"""


res = dp.FInd(frames_av=1, frames_sum=1)

"""
Explanation: Daophot got FITS file in construction, which will be automatically ATTACHed. 
(1) Run FIND on your frame
Daophot FIND parameters Number of frames averaged, summed are defaulted to 1,1, below are provided for clarity.
End of explanation
"""


print ("{} pixels analysed, sky estimate {}, {} stars found.".format(res.pixels, res.sky, res.stars))

"""
Explanation: Check some results returned by FIND, every method for daophot command returns results object.
End of explanation
"""


!ls -lt $dp.dir

"""
Explanation: Also, take a look into runner directory
End of explanation
"""


res = dp.PHotometry(apertures=[8], IS=35, OS=50)

"""
Explanation: We see symlinks to input image and opt files, and i.coo - result of FIND
(2) Run PHOTOMETRY on your frame
Below we run photometry, providing explicitly radius of aperture A1 and IS, OS sky radiuses.
End of explanation
"""


stars = res.photometry_starlist

"""
Explanation: List of stars generated by daophot commands, can be easily get as astwro.starlist.Starlist being essentially pandas.DataFrame:
End of explanation
"""


stars.sort_values('mag_err').iloc[:10]

"""
Explanation: Let's check 10 stars with least A1 error (mag_err column). (pandas style)
End of explanation
"""


sorted_stars = stars.sort_values('mag')
sorted_stars.renumber()

"""
Explanation: (3) SORT the output from PHOTOMETRY
in order of increasing apparent magnitude  decreasing
stellar brightness   with the renumbering feature.  This step is optional  but it can be more convenient than not.
SORT command of daophor is not implemented (yet) in pydaohot. But we do sorting by ourself.
End of explanation
"""


dp.write_starlist(sorted_stars, 'i.ap')

!head -n20 $dp.PHotometry_result.photometry_file

dp.PHotometry_result.photometry_file

"""
Explanation: Here we write sorted list back info photometry file at default name (overwriting existing one), because it's convenient to use default files in next commands.
End of explanation
"""


pick_res = dp.PIck(faintest_mag=20, number_of_stars_to_pick=40)

"""
Explanation: (4) PICK to generate a set of likely PSF stars
How many stars you want to use is a function of the degree of variation you expect  and the frequency with which stars are contaminated by cosmic rays or neighbor stars. [...]
End of explanation
"""


ls $dp.dir

"""
Explanation: If no error reported, symlink to image file (renamed to i.fits), and all daophot output files (i.*) are  in the working directory of runner:
End of explanation
"""


dp.set_options('VARIABLE PSF', 2)
psf_res = dp.PSf()

"""
Explanation: One may examine and improve i.lst list of PSF stars. Or use astwro.tools.gapick.py to obtain list of PSF stars optimised by genetic algorithm.
*(5) Run PSF *
tell it the name of your complete  (sorted  renumbered)  aperture photometry  file,  the name of the  file with the list of PSF stars,  and the name of the disk  file you want the point spread function stored in (the default should be  fine) [...]
If the frame is crowded  it is probably worth your while to generate the  first PSF with the  "VARIABLE PSF"  option set to -1 ---  pure analytic PSF.  That way,  the companions will not generate ghosts in the model PSF that will come back to haunt you later.  You should also have specified a reasonably generous  fitting radius ---  these stars have been preselected to be as isolated as possible  and you want the best  fits you can get.  But remember to avoid letting neighbor stars intrude within one  fitting radius of the center of any PSF star.
For illustration we will set VARIABLE PSF option, before PSf()
End of explanation
"""


alls_res = al.ALlstar(image_file=frame, stars=psf_res.nei_file, subtracted_image_file='is.fits')

"""
Explanation: (6) Run GROUP and NSTAR or ALLSTAR on your NEI  file
If your PSF stars have many neighbors this may take some minutes of real time.  Please be patient  or submit it as a batch job and perform steps       on your next frame while you wait.
We use allstar. (GROUP and NSTAR command are not implemented in current version of pydaophot). We use prepared above Allstar object: al operating on the same runner dir that dp.
As parameter we set input image (we haven't do that on constructor), and nei file produced by PSf(). We do not remember name i.psf so use psf_res.nei_file property. 
Finally we order allstar to produce subtracted FITS .
End of explanation
"""


print (alls_res.get_buffer())

"""
Explanation: All result objects, has get_buffer() method, useful to lookup unparsed daophot or allstar output:
End of explanation
"""


sub_img = alls_res.subtracted_image_file

"""
Explanation: *(8) EXIT from DAOPHOT and send this new picture to the image display *
Examine each of the PSF stars and its environs.  Have all of the PSF stars subtracted out more or less cleanly,  or should some of them be rejected from further use as PSF stars?  (If so  use a text editor to delete these stars from the  LST  file.)   Have the neighbors mostly disappeared,  or have they left behind big zits?  Have you uncovered any faint companions that FIND missed?[...] 
The absolute path to subtracted file (like for most output files) is available as result's property:
End of explanation
"""


from astwro.starlist.ds9 import write_ds9_regions
reg_file_path = dp.file_from_runner_dir('lst.reg')
write_ds9_regions(pick_res.picked_starlist, reg_file_path)

# One can run ds9 directly from notebook:
!ds9 $sub_img -regions $reg_file_path 

"""
Explanation: We can also generate region file for psf stars:
End of explanation
"""


ls $al.dir

sub_res = dp.SUbstar(subtract=alls_res.profile_photometry_file, leave_in=pick_res.picked_stars_file)

"""
Explanation: (9) Back in DAOPHOT II  ATTACH the original picture and run SUBSTAR
specifying the  file created in step (6) or in step (8f) as the stars to subtract,  and the stars in the  LST  file as the stars to keep.
Lookup into runner dir:
End of explanation
"""


for i in range(3):
    print ("Iteration {}: Allstar chi: {}".format(i, alls_res.als_stars.chi.mean()))
    dp.image = 'is.fits'
    respsf = dp.PSf()
    print ("Iteration {}: PSF chi: {}".format(i, respsf.chi))
    alls_res = al.ALlstar(image_file=frame, stars='i.nei')
    dp.image = frame
    dp.SUbstar(subtract='i.als', leave_in='i.lst')
print ("Final:       Allstar chi: {}".format(alls_res.als_stars.chi.mean()))

alls_res.als_stars

"""
Explanation: You have now created a new picture which has the PSF stars still in it  but from which the known neighbors of these PSF stars have been mostly removed
(10) ATTACH the new star subtracted frame and repeat step (5) to derive a new point spread function
(11+...) Run GROUP   NSTAR or ALLSTAR
End of explanation
"""


!ds9 $dp.SUbstar_result.subtracted_image_file -regions $reg_file_path 

"""
Explanation: Check last image with subtracted PSF stars neighbours.
End of explanation
"""


dp.image = 'is.fits'
psf_res = dp.PSf()
print ("PSF file: {}".format(psf_res.psf_file))

"""
Explanation: Once you have produced a frame in which the PSF stars and their neighbors all subtract out cleanly, one more time through PSF should produce a point-spread function you can be proud of.
End of explanation
"""


