#! /usr/bin/env python # chip's program to do photometry on a series and compute polarization # vers 2017 Oct 13 # Finds two brightest stars in each image (the O and E pair), does photometry on them # and uses the counts to compute Q,U and errors using the Tinbergen method. # Input as command line argument like A*.fits and it will use the four HWP angles # (00,22,45,67) to read in each and perform computations. # invoke as python polphot.py files*.fit # or use a single file for a list # python polphot.py files.list import argparse import os.path from os import makedirs from astropy.io import fits import numpy as np # need these for photometry from photutils import CircularAperture from photutils import aperture_photometry from photutils import CircularAnnulus from photutils.utils import calc_total_error from photutils import DAOStarFinder from astropy.stats import sigma_clipped_stats FWHM=35. # FWHM of star PSF in pixels ; set this correctly!! myphotradius=0.8*FWHM # photom radius in pixels myannulusradius=myphotradius+0.1*FWHM # inner sky annulus for photom gain=2.5 # unverified rdn=5 # uncertain?? not used anyway PAoffset=-81.5 # converts instrumental PA to real PAs for WIRO def dophot(filename): hdu = fits.open(filename) data=hdu[0].data mean, median, std = sigma_clipped_stats(data[300:400,300:400], sigma=3.0, iters=5) daofind = DAOStarFinder(fwhm=FWHM, threshold=20.*std) sources = daofind(data - median) sindices=(sources.argsort(['peak'])) # sort and pick out two brightest sources rs=sindices[::-1] # print(rs) # print(sources) xcen=sources['xcentroid'] ycen=sources['ycentroid'] if ycen[rs[0]] < ycen[rs[1]] : ycen1=ycen[rs[0]] # top half of array xcen1=xcen[rs[0]] # ycen2=ycen[rs[1]] # bottom half of array xcen2=xcen[rs[1]] # if ycen[rs[0]] > ycen[rs[1]] : ycen1=ycen[rs[1]] # top half of array xcen1=xcen[rs[1]] # ycen2=ycen[rs[0]] # top half of array xcen2=xcen[rs[0]] # # print(xcen1,ycen1,xcen2,ycen2) # now do photom on two brightest positions = [(xcen1, ycen1), (xcen2, ycen2)] apertures = CircularAperture(positions, r=myphotradius) annulus_apertures = CircularAnnulus(positions, r_in=myannulusradius, r_out=myannulusradius+8) apers = [apertures, annulus_apertures] datae=data/gain # convert data to electrons data_error = np.sqrt(np.absolute(datae)) # estimate data error based on photons only ## take abs value!!! phot_table = aperture_photometry(datae, apers,error=data_error) bkg_mean = phot_table['aperture_sum_1'] / annulus_apertures.area() bkg_sum = bkg_mean * apertures.area() final_sum = phot_table['aperture_sum_0'] - bkg_sum phot_table['residual_aperture_sum'] = final_sum # print(phot_table['residual_aperture_sum']) # print(phot_table) return phot_table['xcenter'].data,phot_table['ycenter'].data,phot_table['residual_aperture_sum'].data, phot_table['aperture_sum_err_0'].data def main(): parser = argparse.ArgumentParser(description='Perform photometry of images') parser.add_argument('filenames',metavar='filenames',nargs='+',help='List of files to photometer') args = parser.parse_args() if len(args.filenames) == 1: # one input filename implies that this file is a list (since why would you try to stack one file?) with open(args.filenames[0],'r') as f: filenames = [x.strip() for x in f.readlines()] args.filenames = filenames myfilenames=args.filenames # print(myfilenames) # loop over filenames nfiles=np.size(myfilenames) # print(nfiles) xc2=np.zeros(4) # create centroid arrays to hold each of 4 positions yc2=np.zeros(4) # xc1=np.zeros(4) # create centroid arrays to hold each of 4 positions yc1=np.zeros(4) # farray1=np.zeros(4) # create arrays to hold each of 4 pol measurements farray2=np.zeros(4) earray1=np.zeros(4) # create arrays to hold each of 4 pol error measurements earray2=np.zeros(4) print('make sure input file list contains multiples of four, e.g., _00/22/45/67') print(' Q U P% Perr PA PAerr') for i in range(0,int(np.size(myfilenames)),4) : # iterate over all images for j in range(0,4) : # loop over sach seq of four images at 0,22.5,45,67.5 # print(myfilenames[i+j]) xcentroid,ycentroid,flux,error=dophot(myfilenames[i+j]) xc1[j]=xcentroid[0] yc1[j]=ycentroid[0] xc2[j]=xcentroid[1] yc2[j]=ycentroid[1] # print('xcen1 ycen1 xcen2 ycen2') # print(xc1[j], yc1[j],xc2[j], yc2[j] ) farray1[j]=flux[0] farray2[j]=flux[1] earray1[j]=error[0] earray2[j]=error[1] i=i+4 q=100*((farray1[0]+farray2[2])-(farray1[2]+farray2[0]))/(farray1[0]+farray1[2]+farray2[0]+farray2[2]) # gives us Q u=100*((farray1[1]+farray2[3])-(farray1[3]+farray2[1]))/(farray1[1]+farray1[3]+farray2[1]+farray2[3]) # gives us Q eq=100*np.sqrt(earray1[0]**2+earray1[2]**2+earray2[0]**2+earray2[2]**2 )/ (farray1[0]+farray1[2]+farray2[0]+farray2[2]) eu=100*np.sqrt(earray1[1]**2+earray1[3]**2+earray2[1]**2+earray2[3]**2 )/ (farray1[1]+farray1[3]+farray2[1]+farray2[3]) # totals Traaditional method! # P=np.sqrt(q**2+u**2) # Perr=(eq+eu)/np.sqrt(2) # PA=0.5 * 57.3*np.arctan(u/q) # if u < 0 and q < 0 : #3rd quadrant # PA=PA+90 # if u < 0 and q > 0 :#4th quadrant # PA=PA+180 # if u > 0 and q < 0 : # 2nd quadrant # PA=PA+90 # PA=PA+PAoffset # if PA < 0 : # PA=PA+180. # PAerr=57.3*np.sqrt((u**2*eq**2 + q**2*eu**2))/(u**2+q**2) # print('Q U P% Perr PA PAerr') # print("%6.3f"% q ,"%6.3f"% u,"%6.3f"% P,"%6.3f"% Perr,"%6.1f"% PA,"%6.1f"% PAerr) # Tinbergen method...better! Rqsq=(farray1[0]/farray2[0])/(farray1[2]/farray2[2]) Rusq=(farray1[1]/farray2[1])/(farray1[3]/farray2[3]) Tq=100*(np.sqrt(Rqsq)-1.) / (np.sqrt(Rqsq)+1.) Tu=100*(np.sqrt(Rusq)-1.) / (np.sqrt(Rusq)+1.) TP=np.sqrt(Tq**2+Tu**2) TPA=0.5 * 57.3*np.arctan(u/q) if Tu < 0 and Tq < 0 : #3rd quadrant TPA=TPA+90 if Tu < 0 and Tq > 0 :#4th quadrant TPA=TPA+180 if Tu > 0 and Tq < 0 : # 2nd quadrant TPA=TPA+90 TPA=TPA+PAoffset Teq=100*earray1[0]/np.sqrt(Rqsq)/farray1[0] /2 # Teu=100*earray1[1]/np.sqrt(Rusq)/farray1[1] /2 # TPerr=(Teq+Teu)/np.sqrt(2) TPAerr=57.3*np.sqrt((Tu**2*Teq**2 + Tq**2*Teu**2))/(Tu**2+Tq**2) if TPA < 0 : TPA=TPA+180. print("%s"% myfilenames[i-1] , "%6.3f"% Tq ,"%6.3f"% Tu,"%6.3f"% TP,"%6.3f"% TPerr,"%6.1f"% TPA,"%6.1f"% TPAerr ) if __name__ == '__main__': main()