Hey Look I found Point Sources! ...Oh wait that's a little bit too many.

The first prototype of findptsources(...) is now done. How does it work?
Great! In fact, too great. Though findptsources(...) lists a lot of tuples as coordinates, there seems to be way too many points. I suppose at some point I didn't filter out the cosmic rays well enough and received some false results.

Anyways, here's some methods I finished today:

def potptsources(self,dataarr,darkpix=0.95): Finds potential point sources and returns a list of tuples of "coordinates"

def findpeak(self, datalist): Returns a list of points brighter than their surrounding points. I actually rewrote this function and completed it without using recursion. This function still needs to be modified, as I haven't really told it what to do if it gets an index out of bound error. 


def findptsources(self, dataarr, darkpix=0.95): Finds the point sources, taking into account background radiation and cosmic rays that might skew results. 


This function is really simple; it uses the two functions above to first find a list of potential sources, and then find the brightest pixels in that list. 


For example:

>>> from startup import*  #imports numpy, phot, etc...

>>> im=phot.Image(f)

>>> im.findptsources(im.data)

[(21, 5), (22, 114), (22, 115), (22, 116), (25, 132), (25, 133), (40, 82), (40, 83), (53, 24), (53, 25), (72, 88), (72, 89), (72, 90), (72, 91), (72, 92), (72, 93), (72, 94), (76, 90), (76, 91), (76, 92), (76, 93), (76, 94), (80, 69), (100, 17), (100, 18), (100, 19), (100, 20), (100, 21), (100, 22), (128, 50), (128, 51), (128, 52)]

Unfortunately, this seems like WAY too many coordinates for potential point sources. Clearly my code is not picky enough. One good thing to note is that the brightest pixel in the whole picture, at position (72,94), is picked up. However, it seems to also have picked up some false data. One possibility of why this happened is that I only checked the pixels adjacent to a "point source" to check if they were "bright" enough. However, a cosmic ray could have "spilled" electrons into neighboring wells, as Professor Johnson kindly pointed out. I will probably have to check pixels farther away or check for a Point Spread Function (PSF) characteristic of stars. I shall continue to work on this next week. 

-----
Learned:
-Little bit more about lists
To-Do:
-Fix up my code so it doesn't give me an overwhelming number of point sources.

Finding Point Sources

Thursdays are the best. We get free donuts/bagels. :D

I continued to modify my phot.py file while munching on a delicious bagel.
Today my main focus was on creating a function that would plot potential point sources. However, in order to do that, I had to create a lot of little sub-methods.

First though, I made two more plotting methods:

def fluxcol(self,col,cx,title="Flux v. Radial Distance from Object Center",fig=1): Plots flux v. radial distance from object along a specified column

Assuming we have an im Image object: (im.xc, im.yc are the 'coordinates' of the brightest pixel in the im.data array)
>>> fluxcol(im.yc,im.xc)

def fluxrow(self,row,cy,title="Flux v. Radial Distance from Object Center",fig=1): Plots flux v. radial distance from object along a specified row

Assuming we have an im Image object:
>>> fluxcol(im.xc,im.yc)

Here are the sub-methods I worked on today:

def bg(self,darkpix=0.95): Finds all the pixels that are in the background, where the number of background pixels is specified by the percentage of darkpix

def get_bgsigma(self,darkpix=0.95): Gets background sigma (excludes bright objects)


def get_bgmean(self,darkpix=0.95): Gets background mean value (excludes bright objects)


def potptsources(self,dataarr,darkpix=0.95): Finds a list of potential point sources

The goal is to eventually use all of these methods in a method called findptsources(). My rough outline of what to do goes somewhat along these lines:
-Use bg() and get_bgsigma() to find the average background variation across the image.
-Then we can use potptsources(im.data) to find the potential point sources (the method will return an array of potential point sources). The method will take all pixels that have a brightness value greater than 5*sigma of the average background brightness and return their coordinates.
-However, we cannot just assume that all bright points are point sources - occasionally cosmic rays might strike one of the CCD wells and produce a shower of electrons, resulting in a very high measurement for that certain well. However, this will be a well with a single incredibly bright pixel surrounded by pixels that have the same brightness as the background. Thus we have to create another function that checks the adjacent pixels for values > 5*sigma+im.get_bgmean(darkpix). The function will keep only those pixels that satisfy these requirements.
-Finally, we will write a recursive function that starts at a pixel, and if it finds a brighter value in one of the adjacent pixels, moves to that pixels and again looks for brighter adjacent pixels. Eventually we will get a few pixels that are brighter than all of their respectively surrounding pixels. These will be the point sources we are looking for.

-----

Learned:
-the histogram plotting function actually can take a range of values for bin as well. For example, if we set a = arange(0,10,0.01), the function will create bins spanning 0-10 that differ by 0.01. Basically, when the function sees an array for bins, it will only make bars at those values in the array. Very neat!
-sorting algorithms (a.sort(), sort(a), etc)
Problems:
-Figuring out the recursive methods/how to code the sub-methods
To-Do:
-Finish the findptsources() code.

Classy!

Created my first class Image in phot.py:
It's far from done, but I'll post up some stuff I've written so far.

Here's the "constructor". It takes a fits file and converts it into an Image instance.
class Image:
   """Class to be used with a FITS file image"""
   def __init__(self,filename):
       self.filename=filename
       self.data,self.hdr=f.getdata(self.filename,header=True)
       self.object=self.hdr['TARGNAME']


Methods:
def header(self): prints the header of the FITS file


def phot(self,xc,yc,rad=5,ann_in=8,ann_out=10): basically a remake of the aper_phot function I made yesterday. The data argument is simply replaced with self.data (data of the FITS file)

def show(self): displays the fits file image on a graph


def histcol(self,col,bins=50,tite="",fig=1): takes one column (all rows in that column), and creates a histogram of all the values found in the elements in that column


def histrow(self,row,bins=50,tite="",fig=1): takes one row (all columns in that row), and creates a histogram of all the values found in the elements in that row

------
Some examples: (assuming I created an instance called im)
The original FITS file's data looks like this:
>>>im.show()

>>> a=im.phot(xc=94,yc=71)
>>> im.show

>>> im.histcol(71)

>>> p.show()

-----

I also created a nice file called startup.py so that I wouldn't have to import all the files by hand every time I opened up the command prompt. 

Learned:

-Stuff about displaying an array as an image

-A little bit more about how classes work (bound, unbound methods, etc)

-a pretty cool numpy.argmax(array) function I haven't quite used yet

Problems:

-Wow setfig(fig,title) is still causing problems =[

-How to manipulate/control what colors arrays display

To-Do:

-make a findptsource(self) function that finds possible point sources depending on the array values in self.data

Aperture Photometry

Hello!

More coding today:

aper_phot.py:
-It calculates the flux of a star (or other object) centered at a pixel, taking into an estimated background flux. It makes use of the makecirc.py file I created earlier.

aper_phot(data,xc,yc,rad=5,ann_in=8,ann_out=10)
data: data array that you want to calculate the flux from
xc,yc: center x and y coordinates of pixel (essentially an array element)
rad: aperture radius
ann_in, ann_out: inner and outer annulus radius, respectively


Annulus refers to a ring around the object that ideally contains only background light from the sky. The aperture is the "size" of the object.

For example:
>>> from numpy import*
>>> import aper_phot
>>> b=arange(49).reshape(7,7)
>>> aper_phot.aper_phot(b,3,3,2,4,5)
array([[  0.,   0.,   0.,   0.,   0.,   0.,   0.],
       [  0.,   0.,   0., -14.,   0.,   0.,   0.],
       [  0.,   0.,  -8.,  -7.,  -6.,   0.,   0.],
       [  0.,  -2.,  -1.,   0.,   1.,   2.,   0.],
       [  0.,   0.,   6.,   7.,   8.,   0.,   0.],
       [  0.,   0.,   0.,  14.,   0.,   0.,   0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.]])


I'm also starting to use emacs as my main Python code editor.


Learned:
-Roughly how aperture photometry works
-Basic emac commands.
Problems:
-setfig still doesn't work correctly. I'm leaving it as it is for now.
To-Do:
-My next big ... subproject is to create a class file that handles images (fits files). I'm saving this class in aper_phot.py, and generally whatever function I need that has to do with aperture photometry will go into this file. 

Contour Plotting

I'm doing a lot of graphing this summer, so I'm still doing many kinds of Python exercises. Today I made a contour map plotter function of a 2-D Gaussian Distribution. I'm pretty sure this code will be very useful in the future, when I start plotting the brightness of each pixel in CCD images.

The graph plots a Gaussian function defined by these equations:

G(x,y)=A0+A1 exp(-U/2)

Where:

U=(x'/a)² + (y'/b)²

and

x'=(x-h)cos T - (y-k)sin T

y'=(x-h)sin T + (y-k)cos T


A0 = constant offset
A1 = amplitude
a = s width of Gaussian in the X direction
b = s width of Gaussian in the Y direction
h = X centroid
k = Y centroid
T = Theta, the rotation of the ellipse from the X axis in radians.

You can pick between two different types of contour maps: Solid or not solid. 

Solid looks like this:

While not solid looks like this:

As you can see, you can also set a bunch of different arguments, including the x and y range, resolution, the angle of rotation, x-axis sigma, y-axis sigma, amplitude, offset, etc.... There is also a handy little column on the right indicating the corresponding values of the colors. 

One very important thing to note is that the solid contour map takes ages to plot. You would be better off not using it. By default, the argument solid is set to false, so that a line contour map is plotted. 

------

Aside from that, today I created my blog and caught up on some posting I should have done last week. I also managed to finish reading one of the research papers I was supposed to. 

Learned:

Contour mapping on Python. Also that it's best not to use fancy solid rainbow colors and stick to lines. 

*I should probably mention that the exercises I do come from this page: 

http://astro.berkeley.edu/~johnjohn/idl_course/

I just use Python instead of IDL to code.