Notebook

TRIPPy examples¶

The first thing to do is import all the necessary packages. Note that this notebook assumes you have the optional packages installed, as well as SExtractor available on your command line.

In [34]:

%matplotlib inline
import numpy as num, astropy.io.fits as pyf,pylab as pyl
from trippy import psf, pill, psfStarChooser
from trippy import scamp,MCMCfit
import scipy as sci
from os import path

NOTE: proper use of psfStarChooser requires plot interaction. So for this tutorial you'd best comment out the first line, %matplotlib inline. But for my web presentation, I leave inline.

The function trim catalog is a convenience function to simply return only those sources that are well enough isolated for PSF generation. We may fold this into psfStarChooser in the future.

In [35]:

def trimCatalog(cat):
    good=[]
    for i in range(len(cat['XWIN_IMAGE'])):
        try:
            a=int(cat['XWIN_IMAGE'][i])
            b=int(cat['YWIN_IMAGE'][i])
            m=num.max(data[b-4:b+5,a-4:a+5])
        except: pass
        dist=num.sort(((cat['XWIN_IMAGE']-cat['XWIN_IMAGE'][i])**2+(cat['YWIN_IMAGE']-cat['YWIN_IMAGE'][i])**2)**0.5)
        d=dist[1]
        if abs(cat['XWIN_IMAGE'][i]-812)<5500 and  abs(cat['YWIN_IMAGE'][i]-4004)<5000 and cat['FLAGS'][i]==0 and d>30 and m<70000:
            good.append(i)
    good=num.array(good)
    outcat={}
    for i in cat:
        outcat[i]=cat[i][good]
    return outcat

Get the image this tutorial assumes you have. If wget fails then you are likely on a mac, and should just download it manually

In [36]:

inputFile='Polonskaya.fits'
if not path.isfile(inputFile):
    os.system('wget http://www.canfar.phys.uvic.ca/vospace/nodes/fraserw/Polonskaya.fits?view=data')
else:
    print "We already have the file."

We already have the file.

First load the fits image and get out the header, data, and exposure time.

In [37]:

with pyf.open(inputFile) as han:
    data=han[0].data
    header=han[0].header
    EXPTIME=header['EXPTIME']
    

Next run sextractor on the images, and use trimCatalog to create a trimmed down list of isolated sources. Finally, find the source closest to 811, 4005 which is the bright asteroid, 2006 Polonskaya. Also, set the rate and angle of motion. These were found from JPL horizons. The 1 degree increase is to account for the slight rotation of the image.

In [38]:

scamp.makeParFiles.writeSex('example.sex',
                    minArea=3.,
                    threshold=5.,
                    zpt=27.8,
                    aperture=20.,
                    min_radius=2.0,
                          catalogType='FITS_LDAC',
                          saturate=55000)
scamp.makeParFiles.writeConv()
scamp.makeParFiles.writeParam(numAps=1) #numAps is thenumber of apertures that you want to use. Here we use 1

scamp.runSex('example.sex', inputFile ,options={'CATALOG_NAME':'example.cat'})
catalog=trimCatalog(scamp.getCatalog('example.cat',paramFile='def.param'))
dist=((catalog['XWIN_IMAGE']-811)**2+(catalog['YWIN_IMAGE']-4005)**2)**0.5
args=num.argsort(dist)
xt=catalog['XWIN_IMAGE'][args][0]
yt=catalog['YWIN_IMAGE'][args][0]

rate=18.4588 # "/hr
angle=31.11+1.1 # degrees counter clockwise from horizontal, right

sex Polonskaya.fits -c example.sex -CATALOG_NAME example.cat
> 
----- SExtractor 2.19.5 started on 2016-01-06 at 15:48:57 with 1 thread

> Setting catalog parameters
> Reading detection filter
> Initializing catalog
> Looking for Polonskaya.fits
----- Measuring from: Polonskaya.fits
      "CDE1303G0" / no ext. header / 2112x4644 / 32 bits (floats)
Detection+Measurement image: > Setting up background maps
> Setting up background map at line:   64
> Setting up background map at line:  128
> Setting up background map at line:  192
> Setting up background map at line:  256
> Setting up background map at line:  320
> Setting up background map at line:  384
> Setting up background map at line:  448
> Setting up background map at line:  512
> Setting up background map at line:  576
> Setting up background map at line:  640
> Setting up background map at line:  704
> Setting up background map at line:  768
> Setting up background map at line:  832
> Setting up background map at line:  896
> Setting up background map at line:  960
> Setting up background map at line: 1024
> Setting up background map at line: 1088
> Setting up background map at line: 1152
> Setting up background map at line: 1216
> Setting up background map at line: 1280
> Setting up background map at line: 1344
> Setting up background map at line: 1408
> Setting up background map at line: 1472
> Setting up background map at line: 1536
> Setting up background map at line: 1600
> Setting up background map at line: 1664
> Setting up background map at line: 1728
> Setting up background map at line: 1792
> Setting up background map at line: 1856
> Setting up background map at line: 1920
> Setting up background map at line: 1984
> Setting up background map at line: 2048
> Setting up background map at line: 2112
> Setting up background map at line: 2176
> Setting up background map at line: 2240
> Setting up background map at line: 2304
> Setting up background map at line: 2368
> Setting up background map at line: 2432
> Setting up background map at line: 2496
> Setting up background map at line: 2560
> Setting up background map at line: 2624
> Setting up background map at line: 2688
> Setting up background map at line: 2752
> Setting up background map at line: 2816
> Setting up background map at line: 2880
> Setting up background map at line: 2944
> Setting up background map at line: 3008
> Setting up background map at line: 3072
> Setting up background map at line: 3136
> Setting up background map at line: 3200
> Setting up background map at line: 3264
> Setting up background map at line: 3328
> Setting up background map at line: 3392
> Setting up background map at line: 3456
> Setting up background map at line: 3520
> Setting up background map at line: 3584
> Setting up background map at line: 3648
> Setting up background map at line: 3712
> Setting up background map at line: 3776
> Setting up background map at line: 3840
> Setting up background map at line: 3904
> Setting up background map at line: 3968
> Setting up background map at line: 4032
> Setting up background map at line: 4096
> Setting up background map at line: 4160
> Setting up background map at line: 4224
> Setting up background map at line: 4288
> Setting up background map at line: 4352
> Setting up background map at line: 4416
> Setting up background map at line: 4480
> Setting up background map at line: 4544
> Setting up background map at line: 4608
> Filtering background map(s)
> Computing background d-map
> Computing background-noise d-map
(M+D) Background: 1223.9     RMS: 28.4164    / Threshold: 142.082    
> Scanning image
> Line:   25  Objects:       11 detected /        0 sextracted
> Line:   50  Objects:       26 detected /        0 sextracted
> Line:   75  Objects:       35 detected /        0 sextracted
> Line:  100  Objects:       46 detected /        0 sextracted
> Line:  125  Objects:       54 detected /        0 sextracted
> Line:  150  Objects:       64 detected /        0 sextracted
> Line:  175  Objects:       71 detected /        0 sextracted
> Line:  200  Objects:       85 detected /        0 sextracted
> Line:  225  Objects:       93 detected /        0 sextracted
> Line:  250  Objects:      100 detected /        0 sextracted
> Line:  275  Objects:      114 detected /        0 sextracted
> Line:  300  Objects:      120 detected /        0 sextracted
> Line:  325  Objects:      128 detected /        0 sextracted
> Line:  350  Objects:      133 detected /        0 sextracted
> Line:  375  Objects:      138 detected /        0 sextracted
> Line:  400  Objects:      140 detected /        0 sextracted
> Line:  425  Objects:      151 detected /        0 sextracted
> Line:  450  Objects:      158 detected /        0 sextracted
> Line:  475  Objects:      173 detected /        0 sextracted
> Line:  500  Objects:      181 detected /        0 sextracted
> Line:  525  Objects:      195 detected /        0 sextracted
> Line:  550  Objects:      203 detected /        0 sextracted
> Line:  575  Objects:      209 detected /        0 sextracted
> Line:  600  Objects:      215 detected /        0 sextracted
> Line:  625  Objects:      226 detected /        0 sextracted
> Line:  650  Objects:      230 detected /        0 sextracted
> Line:  675  Objects:      238 detected /        0 sextracted
> Line:  700  Objects:      245 detected /        0 sextracted
> Line:  725  Objects:      253 detected /        0 sextracted
> Line:  750  Objects:      259 detected /        0 sextracted
> Line:  775  Objects:      269 detected /        0 sextracted
> Line:  800  Objects:      275 detected /        0 sextracted
> Line:  825  Objects:      280 detected /        0 sextracted
> Line:  850  Objects:      286 detected /        0 sextracted
> Line:  875  Objects:      300 detected /        0 sextracted
> Line:  900  Objects:      306 detected /        0 sextracted
> Line:  925  Objects:      309 detected /        0 sextracted
> Line:  950  Objects:      314 detected /        0 sextracted
> Line:  975  Objects:      324 detected /        0 sextracted
> Line: 1000  Objects:      329 detected /        0 sextracted
> Line: 1022  Objects:      339 detected /        0 sextracted
> Line: 1025  Objects:      340 detected /        8 sextracted
> Line: 1050  Objects:      350 detected /       23 sextracted
> Line: 1075  Objects:      356 detected /       30 sextracted
> Line: 1100  Objects:      362 detected /       46 sextracted
> Line: 1125  Objects:      374 detected /       57 sextracted
> Line: 1150  Objects:      386 detected /       66 sextracted
> Line: 1175  Objects:      397 detected /       74 sextracted
> Line: 1200  Objects:      404 detected /       81 sextracted
> Line: 1225  Objects:      411 detected /       90 sextracted
> Line: 1250  Objects:      420 detected /       99 sextracted

> WARNING: Pixel stack overflow at position 2032,1256

> Line: 1275  Objects:      426 detected /      111 sextracted
> Line: 1300  Objects:      432 detected /      118 sextracted
> Line: 1325  Objects:      439 detected /      124 sextracted
> Line: 1350  Objects:      443 detected /      130 sextracted
> Line: 1375  Objects:      449 detected /      135 sextracted
> Line: 1400  Objects:      458 detected /      138 sextracted
> Line: 1425  Objects:      469 detected /      147 sextracted
> Line: 1450  Objects:      475 detected /      155 sextracted
> Line: 1475  Objects:      484 detected /      173 sextracted
> Line: 1500  Objects:      492 detected /      181 sextracted
> Line: 1525  Objects:      495 detected /      192 sextracted
> Line: 1550  Objects:      503 detected /      202 sextracted
> Line: 1575  Objects:      506 detected /      209 sextracted
> Line: 1600  Objects:      512 detected /      214 sextracted
> Line: 1625  Objects:      524 detected /      225 sextracted
> Line: 1650  Objects:      528 detected /      229 sextracted
> Line: 1675  Objects:      533 detected /      234 sextracted

> WARNING: Pixel stack overflow at position 18,1686

> Line: 1700  Objects:      541 detected /      244 sextracted
> Line: 1725  Objects:      543 detected /      251 sextracted
> Line: 1750  Objects:      551 detected /      256 sextracted
> Line: 1775  Objects:      557 detected /      260 sextracted
> Line: 1800  Objects:      568 detected /      260 sextracted
> Line: 1825  Objects:      572 detected /      269 sextracted
> Line: 1850  Objects:      577 detected /      272 sextracted
> Line: 1875  Objects:      584 detected /      285 sextracted
> Line: 1900  Objects:      592 detected /      294 sextracted
> Line: 1925  Objects:      596 detected /      300 sextracted
> Line: 1950  Objects:      601 detected /      302 sextracted
> Line: 1975  Objects:      606 detected /      310 sextracted
> Line: 2000  Objects:      608 detected /      317 sextracted
> Line: 2025  Objects:      619 detected /      328 sextracted
> Line: 2050  Objects:      627 detected /      335 sextracted
> Line: 2075  Objects:      633 detected /      340 sextracted
> Line: 2100  Objects:      641 detected /      347 sextracted
> Line: 2125  Objects:      647 detected /      357 sextracted
> Line: 2150  Objects:      652 detected /      368 sextracted
> Line: 2175  Objects:      658 detected /      379 sextracted
> Line: 2200  Objects:      666 detected /      386 sextracted
> Line: 2225  Objects:      679 detected /      395 sextracted
> Line: 2247  Objects:      686 detected /      400 sextracted

> WARNING: Pixel stack overflow at position 65,2250

> Line: 2250  Objects:      687 detected /      401 sextracted
> Line: 2275  Objects:      694 detected /      409 sextracted
> Line: 2300  Objects:      698 detected /      415 sextracted
> Line: 2325  Objects:      708 detected /      421 sextracted
> Line: 2350  Objects:      715 detected /      426 sextracted
> Line: 2375  Objects:      727 detected /      432 sextracted
> Line: 2400  Objects:      739 detected /      443 sextracted
> Line: 2425  Objects:      750 detected /      452 sextracted
> Line: 2450  Objects:      758 detected /      460 sextracted
> Line: 2475  Objects:      764 detected /      467 sextracted
> Line: 2500  Objects:      772 detected /      475 sextracted
> Line: 2525  Objects:      781 detected /      480 sextracted
> Line: 2550  Objects:      788 detected /      486 sextracted
> Line: 2575  Objects:      793 detected /      488 sextracted
> Line: 2600  Objects:      797 detected /      494 sextracted

> WARNING: Pixel stack overflow at position 74,2619

> Line: 2625  Objects:      807 detected /      503 sextracted
> Line: 2650  Objects:      815 detected /      509 sextracted
> Line: 2675  Objects:      827 detected /      514 sextracted
> Line: 2700  Objects:      833 detected /      520 sextracted
> Line: 2725  Objects:      838 detected /      522 sextracted
> Line: 2750  Objects:      845 detected /      533 sextracted
> Line: 2775  Objects:      862 detected /      535 sextracted
> Line: 2800  Objects:      872 detected /      546 sextracted
> Line: 2825  Objects:      881 detected /      552 sextracted
> Line: 2850  Objects:      889 detected /      558 sextracted
> Line: 2875  Objects:      906 detected /      562 sextracted
> Line: 2900  Objects:      911 detected /      573 sextracted
> Line: 2925  Objects:      917 detected /      576 sextracted
> Line: 2950  Objects:      921 detected /      579 sextracted
> Line: 2975  Objects:      926 detected /      583 sextracted
> Line: 3000  Objects:      934 detected /      587 sextracted
> Line: 3025  Objects:      939 detected /      599 sextracted
> Line: 3050  Objects:      951 detected /      603 sextracted
> Line: 3075  Objects:      958 detected /      611 sextracted

> WARNING: Pixel stack overflow at position 13,3079

> Line: 3100  Objects:      966 detected /      619 sextracted
> Line: 3125  Objects:      972 detected /      624 sextracted
> Line: 3150  Objects:      979 detected /      629 sextracted
> Line: 3175  Objects:      984 detected /      634 sextracted
> Line: 3200  Objects:      992 detected /      646 sextracted
> Line: 3225  Objects:     1002 detected /      656 sextracted
> Line: 3250  Objects:     1007 detected /      665 sextracted
> Line: 3275  Objects:     1014 detected /      673 sextracted
> Line: 3300  Objects:     1024 detected /      677 sextracted
> Line: 3325  Objects:     1029 detected /      684 sextracted
> Line: 3350  Objects:     1040 detected /      693 sextracted
> Line: 3375  Objects:     1046 detected /      702 sextracted
> Line: 3400  Objects:     1051 detected /      714 sextracted
> Line: 3425  Objects:     1059 detected /      728 sextracted
> Line: 3450  Objects:     1069 detected /      734 sextracted
> Line: 3475  Objects:     1075 detected /      741 sextracted
> Line: 3500  Objects:     1083 detected /      751 sextracted
> Line: 3525  Objects:     1087 detected /      761 sextracted
> Line: 3550  Objects:     1091 detected /      767 sextracted

> WARNING: Pixel stack overflow at position 65,3574

> Line: 3575  Objects:     1101 detected /      772 sextracted
> Line: 3600  Objects:     1108 detected /      775 sextracted
> Line: 3625  Objects:     1115 detected /      785 sextracted
> Line: 3650  Objects:     1124 detected /      792 sextracted
> Line: 3670  Objects:     1131 detected /      800 sextracted
> Line: 3675  Objects:     1132 detected /      804 sextracted
> Line: 3700  Objects:     1142 detected /      811 sextracted
> Line: 3725  Objects:     1153 detected /      817 sextracted
> Line: 3750  Objects:     1164 detected /      825 sextracted
> Line: 3775  Objects:     1171 detected /      838 sextracted
> Line: 3800  Objects:     1181 detected /      847 sextracted
> Line: 3825  Objects:     1185 detected /      856 sextracted
> Line: 3850  Objects:     1194 detected /      862 sextracted
> Line: 3875  Objects:     1204 detected /      875 sextracted
> Line: 3900  Objects:     1209 detected /      882 sextracted
> Line: 3925  Objects:     1221 detected /      890 sextracted
> Line: 3950  Objects:     1230 detected /      894 sextracted
> Line: 3975  Objects:     1242 detected /      899 sextracted
> Line: 4000  Objects:     1252 detected /      906 sextracted
> Line: 4025  Objects:     1269 detected /      909 sextracted
> Line: 4050  Objects:     1277 detected /      919 sextracted
> Line: 4075  Objects:     1287 detected /      928 sextracted
> Line: 4100  Objects:     1300 detected /      938 sextracted
> Line: 4125  Objects:     1306 detected /      944 sextracted
> Line: 4150  Objects:     1309 detected /      950 sextracted
> Line: 4175  Objects:     1314 detected /      956 sextracted
> Line: 4200  Objects:     1329 detected /      965 sextracted

> WARNING: Pixel stack overflow at position 47,4215

> Line: 4225  Objects:     1338 detected /      972 sextracted
> Line: 4250  Objects:     1349 detected /      980 sextracted
> Line: 4275  Objects:     1355 detected /      987 sextracted
> Line: 4300  Objects:     1363 detected /      994 sextracted
> Line: 4325  Objects:     1373 detected /     1000 sextracted
> Line: 4350  Objects:     1379 detected /     1008 sextracted
> Line: 4375  Objects:     1383 detected /     1017 sextracted
> Line: 4400  Objects:     1393 detected /     1022 sextracted
> Line: 4425  Objects:     1396 detected /     1028 sextracted
> Line: 4450  Objects:     1407 detected /     1037 sextracted
> Line: 4475  Objects:     1414 detected /     1045 sextracted
> Line: 4500  Objects:     1426 detected /     1053 sextracted
> Line: 4525  Objects:     1435 detected /     1057 sextracted
> Line: 4550  Objects:     1440 detected /     1062 sextracted
> Line: 4575  Objects:     1445 detected /     1072 sextracted
> Line: 4600  Objects:     1450 detected /     1079 sextracted
> Line: 4625  Objects:     1453 detected /     1094 sextracted
> Line: 4644  Objects:     1453 detected /     1200 sextracted
      Objects: detected 1453     / sextracted 1381            

> Closing files
> 
> All done (in 2.3 s: 2038.7 lines/s , 606.3 detections/s)

Now use psfStarChooser to select the PSF stars. The first and second parameters to starChooser are the fitting box width in pixels, and the SNR minimum required for a star to be considered as a potential PSF star. This will pop-up a multipanel window. Top left: histogram of fit chi values. Top right: chi vs. FWHM for each fitted source. Middle right: histogram of FWHM. Bottom right: image display of the currently selected source. Bottom left: Radial profiles of all sources displayed in the top right scatter plot.

The point of this window is to select only good stars for PSF generation, done by zooming to the good sources, and rejecting those that are bad.

Use the zoom tool to select the region containing the stars. In this image, that's a cluser at FWHM~3.5 pixels.

Left and right clicks will select a source, now surrounded by a diamond, displaying the radial profile bottom left, and the actual image bottom right.

Right click will oscillate between accepted source and rejected source (blue and red respectively).

When the window is closed, only those sources shown as blue points, and within the zoom of the top right plot will be used to generate the PSF.

In [39]:

starChooser=psfStarChooser.starChooser(data,catalog['XWIN_IMAGE'],catalog['YWIN_IMAGE'],catalog['FLUX_AUTO'],catalog['FLUXERR_AUTO'])
(goodFits,goodMeds,goodSTDs)=starChooser(30,200)
print goodFits
print goodMeds

Fitting stars with moffat profiles...
 1657.51   157.63  15.16  2.45  17.64
 1009.42   363.71  2.84  2.56  3.18
  251.90   684.70  2.86  2.59  3.44
 1211.61   936.20  3.16  1.55  4.89
 1587.22   945.88  2.67  2.35  3.18
 1081.55   914.77  2.78  2.49  3.42
 1315.03  1023.28  2.84  2.57  3.19
 1652.57  1014.63  2.77  2.53  3.13
  383.82  1238.34  2.89  2.62  3.18
 1241.75  1286.83  3.56  2.26  4.30
  510.65  1902.63  20.19  3.00  20.81
 1106.80  2196.58  2.90  2.62  3.43
  855.31  2362.30  2.99  2.69  3.50
 1606.23  2700.56  2.76  2.51  3.16
  433.79  2761.05  2.93  2.61  3.45
 1458.17  3023.19  7.59  2.55  8.71
 1071.70  3291.51  2.90  2.62  3.42
 1633.35  3827.17  2.96  2.68  3.49
  357.50  3746.18  2.90  2.62  3.21
 1269.99  4429.10  2.86  2.55  3.47
 1483.94  4477.91  2.94  2.66  3.45
  820.38  4243.44  2.84  2.55  3.22
 1919.63  4017.65  2.93  2.67  3.42
  294.44  4080.45  2.65  2.32  3.22
  812.20  4004.29  5.67  2.55  6.40
 1635.10  3908.77  3.27  1.37  5.47

[[  1.76365860e+01   8.47076869e+01   1.51633490e+01   2.44672509e+00
    1.65751340e+03   1.57628646e+02   1.25583203e+03]
 [  3.18047133e+00   4.16146488e+01   2.84091207e+00   2.55914980e+00
    1.00941845e+03   3.63707502e+02   1.21809766e+03]
 [  3.44006198e+00   3.55634993e+01   2.85614090e+00   2.59008509e+00
    2.51904915e+02   6.84704324e+02   1.20371142e+03]
 [  4.89032838e+00   9.40025485e+01   3.15819079e+00   1.54545918e+00
    1.21161443e+03   9.36202565e+02   1.23887108e+03]
 [  3.17931758e+00   2.33440738e+02   2.66682869e+00   2.34832967e+00
    1.58722404e+03   9.45882915e+02   1.23580535e+03]
 [  3.42376880e+00   3.96245660e+01   2.77502737e+00   2.48925830e+00
    1.08155135e+03   9.14770754e+02   1.22461080e+03]
 [  3.19035275e+00   1.55863559e+02   2.83935119e+00   2.57341680e+00
    1.31503491e+03   1.02328459e+03   1.22805218e+03]
 [  3.13232261e+00   6.28573733e+01   2.77334395e+00   2.52547706e+00
    1.65257487e+03   1.01462975e+03   1.22662066e+03]
 [  3.18206829e+00   4.81066235e+01   2.89152259e+00   2.62417586e+00
    3.83816005e+02   1.23834251e+03   1.21280239e+03]
 [  4.29756344e+00   3.90522734e+01   3.56057212e+00   2.26028407e+00
    1.24174842e+03   1.28682566e+03   1.22804821e+03]
 [  2.08119995e+01   6.09206642e+01   2.01894709e+01   2.99718748e+00
    5.10649808e+02   1.90262656e+03   1.24476172e+03]
 [  3.42732966e+00   6.70109832e+01   2.90382391e+00   2.62133781e+00
    1.10680319e+03   2.19658165e+03   1.23492188e+03]
 [  3.49637218e+00   2.07430028e+02   2.99001717e+00   2.69432889e+00
    8.55306404e+02   2.36229528e+03   1.22864870e+03]
 [  3.16308084e+00   3.69086644e+01   2.76211131e+00   2.51194441e+00
    1.60622802e+03   2.70056419e+03   1.23157137e+03]
 [  3.45338931e+00   9.95061041e+01   2.92882439e+00   2.61482931e+00
    4.33790608e+02   2.76105321e+03   1.22493726e+03]
 [  8.70547099e+00   6.66627137e+01   7.59389955e+00   2.54821532e+00
    1.45816523e+03   3.02318655e+03   1.23625909e+03]
 [  3.41717646e+00   4.06792693e+01   2.89628631e+00   2.62390841e+00
    1.07169935e+03   3.29150670e+03   1.22649637e+03]
 [  3.49169241e+00   1.96305421e+02   2.95667581e+00   2.67809536e+00
    1.63334565e+03   3.82717334e+03   1.23647403e+03]
 [  3.20871784e+00   1.67596546e+02   2.89667888e+00   2.62019203e+00
    3.57495177e+02   3.74617700e+03   1.22044441e+03]
 [  3.46742090e+00   3.82549479e+01   2.86143047e+00   2.55208340e+00
    1.26998994e+03   4.42910156e+03   1.23267272e+03]
 [  3.45375334e+00   2.05480576e+02   2.94295360e+00   2.66163930e+00
    1.48393614e+03   4.47791465e+03   1.23492188e+03]
 [  3.21660368e+00   6.98053060e+01   2.84064504e+00   2.55278356e+00
    8.20382356e+02   4.24344218e+03   1.22624835e+03]
 [  3.42134715e+00   4.09623031e+01   2.92679915e+00   2.67208634e+00
    1.91962918e+03   4.01765066e+03   1.23524358e+03]
 [  3.22003210e+00   3.62346976e+01   2.65367009e+00   2.31600015e+00
    2.94444052e+02   4.08044641e+03   1.21702298e+03]
 [  6.39871703e+00   6.97399542e+02   5.66740946e+00   2.55032674e+00
    8.12198315e+02   4.00429471e+03   1.23492188e+03]
 [  5.47002098e+00   4.95513808e+01   3.27088219e+00   1.36811572e+00
    1.63509662e+03   3.90877086e+03   1.23492188e+03]]
[    4.16519518    63.16116785     3.00716584     2.50293744  1110.51644068
   524.205913    1228.48436776]

Generate the PSF. We want a 61 pixel wide PSF, adopt a repFactor of 10, and use the mean star fits chosen above.

Repfactors of 5 and 10 have been tested thoroughly. Larger is pointless, smaller is inaccurate. 5 is faster than 10, 10 is more accurate than 5.

The PSF has to be wide/tall enough to handle the trailing length and the seeing disk. For Polonskaya, the larger is trailing, at ~19"/hr*480s/3600/0.185"/pix = 14 pixels. Choose something a few times larger. Also, stick with odd width PSFs, as the even ones have some funny centroid stuff that I haven't fully sorted out.

The full PSF is created with instantiation, and running both genLookupTable and genPSF.

In [40]:

goodPSF=psf.modelPSF(num.arange(61), num.arange(61), alpha=goodMeds[2],beta=goodMeds[3],repFact=10)
fwhm=goodPSF.FWHM() ###this is the pure moffat FWHM
goodPSF.genLookupTable(data,goodFits[:,4]-0.5,goodFits[:,5]-0.5,verbose=False)
goodPSF.genPSF() ###
fwhm=goodPSF.FWHM() ###this is the FWHM with lookuptable included

Now generate the TSF, which we call the line/long PSF interchangeably through the code...

Rate is in units of length/time and pixScale is in units of length/pixel, time and length are in units of your choice. Sanity suggests arcseconds and hours. Then rate in "/hr and pixScale in "/pix. Angle is in degrees counter clockwise from horizontal between +-90 degrees.

This can be rerun to create a TSF with different rate/angle of motion, though keep in mind that the psf class only contains one longPSF (one rate/angle) at any given time.

In [41]:

goodPSF.line(rate,angle,EXPTIME/3600.,pixScale=0.185,useLookupTable=True)

Using the lookup table when generating the long PSF.

Now calculate aperture corrections for the PSF and TSF. Store for values of r=1.4*FWHM.

Note that the precision of the aperture correction depends lightly on the sampling from the compute functions. 10 is generally enough to preserve 1% precision in the .roundAperCorr() and lineAperCorr() functions which use linear interpolation to get the value one actually desires.

NOTE: Set useLookupTable=False if one wants to calculate from the moffat profile alone. Generally, not accuarate for small apertures however.

In [42]:

goodPSF.computeRoundAperCorrFromMoffat(psf.extent(0.8*fwhm,4*fwhm,10),display=False,
                                                          displayAperture=False,
                                                          useLookupTable=True)
roundAperCorr=goodPSF.roundAperCorr(1.4*fwhm)

goodPSF.computeLineAperCorrFromMoffat(psf.extent(0.1*fwhm,4*fwhm,10),
                                                         l=(EXPTIME/3600.)*rate/0.185,a=angle,display=False,displayAperture=False)
lineAperCorr=goodPSF.lineAperCorr(1.4*fwhm)
print lineAperCorr,roundAperCorr

    0.345475009184 7720.71598514 -9.71914394189
    0.520503078829 10692.3814537 -10.0726861098
    0.784205652705 15717.1770591 -10.4909363642
    1.18150791177 22837.1850609 -10.896606428
    1.78009548994 31174.6839711 -11.2345051489
    2.68194560675 40388.3547581 -11.5156404058
    4.04070021985 49060.5354354 -11.7268307092
    6.08784094116 57572.9756252 -11.9005466907
    9.17212495567 66541.6090376 -12.0577332458
    13.8190003674 75552.2796207 -12.1956189318
0.4012456164 0.16589487044

Store the PSF.

In [43]:

goodPSF.psfStore('psf.fits')

We could save most of the above by restoring a previously constructed PSF by the following commented out code.

In [44]:

#goodPSF=psf.modelPSF(restore='psf.fits')
#goodPSF.fitted=False
#goodPSF.line(rate,angle,EXPTIME/3600.,pixScale=0.185,useLookupTable=True)

Now let's do some pill aperture photometry. Instantiate the class, then call the object you created to get photometry of Polonskaya. Again assume repFact=10.

pillPhot assumes iraf coordinates, unlike psf which assumes numpy coordinates. Hence the +1+1 in the below.

enableBGselection=True will cause a popup display of the source, in which one can zoom to a section with no background source.

mode="smart" is the default background selection mode. See the bgFinder documentation for other acceptable modes. Really though, you'd be silly to select anything else.

display=True to see the image subsection

r is the radius of the pill, l is the length, a is the angle. Sky radius is the radius of a larger pill aperture. The pixels in this larger aperture, but outside the smaller aperture are ignored. Anything outside the larger pill, but inside +-width is used for background estimation.

Trimbghighpix is mostly made not important if mode=smart. But if you want to use a mean or median for some reason, then this value is used to reject pixels with values trimBGhighPix standard deviations above the mean of the cutout.

In [45]:

phot=pill.pillPhot(data,repFact=10)
#get photometry, assume ZPT=26.0
#enableBGselection=True allows you to zoom in on a good background region in the aperture display window
#trimBGhighPix is a sigma cut to get rid of the cosmic rays. They get marked as blue in the display window
#background is selected inside the box and outside the skyRadius value
#mode is th background mode selection. Options are median, mean, histMode (JJ's jjkmode technique), fraserMode (ask me about it), gaussFit, and "smart". Smart does a gaussian fit first, and if the gaussian fit value is discrepant compared to the expectation from the background std, it resorts to the fraserMode. "smart" seems quite robust to nearby bright sources

####make sure to use IRAF coordinates not numpy/sextractor coordinates
phot(xt+1,yt+1,radius=fwhm*1.4,l=(EXPTIME/3600.)*rate/0.185,a=angle,
             skyRadius=4*fwhm,width=6*fwhm,
              zpt=26.0,exptime=EXPTIME,enableBGSelection=True,display=True,
              mode="smart",trimBGHighPix=3.)

Current background value: 1235.816

the SNR function calculates the SNR of the aperture. verbose=True puts some nice terminal output in your face. These values can be accessed with their internal names.

In [24]:

phot.SNR(verbose=True)

#get those values
print phot.magnitude
print phot.dmagnitude
print phot.sourceFlux
print phot.snr
print phot.bg

   SNR: 985.066057126
   Flux: 769368.115411
   Background: 1235.78197211
   Background STD: 29.6530471375
   Num Pixels : 185.63

17.9879508358
0.00110219634196
769368.115411
985.066057126
1235.78197211

Finally, let's do some PSF source subtraction. This is only possible with emcee and sextractor installed.

First get the cutout. This makes everything faster later. Also, remove the background, just because.

In [25]:

Data=data[int(yt)-200:int(yt)+200,int(xt)-200:int(xt)+200]-phot.bg

Now instantiate the MCMCfitter class, and then perform the fit. Verbose=False will not put anything to terminal. Setting to true will dump the result of each step. Only good idea if you insist on seeing what's happening. Do you trust black boxes?

Set useLinePSF to True if you are fitting a trailed source, False if a point source.

Set useErrorMap to True if you care to use an estimate of the poisson noise in each pixel during your fit. This produces honest confidence ranges.

I personally like nWalkers=nBurn=nStep=40. To get a reasonable fit however, that's overkill. But to get the best... your mileage will vary.

In [31]:

fitter=MCMCfit.MCMCfitter(goodPSF,Data)
fitter.fitWithModelPSF(200+xt-int(xt)-1,200+yt-int(yt)-1, m_in=1000.,fitWidth=10, nWalkers=40, nBurn=40, nStep=40, bg=phot.bg, useLinePSF=True, verbose=False,useErrorMap=False)

Now get the fits results, including best fit and confidence region using the input value. 0.67 for 1-sigma is shown

In [32]:

(fitPars,fitRange)=fitter.fitResults(0.67)
print fitPars
print fitRange

Best point: [  199.16410475   199.31572032  1260.68351167 -1726.16344823]
[  199.16410475   199.31572032  1260.68351167 -1726.16344823]
[[199.18231786306143, 199.21709676807941], [199.30792212714854, 199.33408252794339], [1258.5966662233327, 1260.5158660021466]]

Finally, lets produce the model best fit image, and perform a subtraction. Plant will plant a fake source with the given input x,y,amplitude into the input data. If returnModel=True, then no source is planted, but the model image that would have been planted is returned.

remove will do the opposite of plant given input data (it actually just calls plant).

In [46]:

modelImage=goodPSF.plant(fitPars[0],fitPars[1],fitPars[2],Data,addNoise=False,useLinePSF=True,returnModel=True)
pyl.imshow(modelImage)
pyl.show()

removed=goodPSF.remove(fitPars[0],fitPars[1],fitPars[2],Data,useLinePSF=True)
pyl.imshow(removed)
pyl.show()

In [ ]: