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. 

Test Scripts in Python-II (Tues-Fri last week)

Over the next four days I played around with Python's matplotlib's various graphing capabilities.

Tuesday-Wednesday:
Made a file called gauss.py:
-Made a Gaussian function that returns the gaussian at a certain x for a given sigma and mu
-Made a plotting function that plotted the gaussian at a range of x values, using the function described above. This graph definitely had some other nice things thrown in: choice of having a vertical line at mu, x and y labels, title, and resolution specifications.
Graph:

Wednesday-Friday Morning:
Made a file called mc_pi:
-Made a function that estimated the value of pi using the Monte Carlo Method given a number n (it can only go up to 1e8 though, or my computer will complain and crash).
-Created an array of pi values estimated given a certain number n with a specified number of array elements (niter).
-Put these array values into a histogram graph, where the user can specify the number of "bins", n, and niter. Also graphed a gaussian curve over it (sigma=standard deviation, mu=mean). The dotted line is at Pi. 
Graph:

-Finally, I plotted a log graph (x-axis) where pi is estimated with an array of n values, niter times for each value. The error bars are of length sigma=1, and the dots are given by mu.
Graph:

Tim also helped me install matplotlibrc on my computer; however it doesn't seem to work in Python. I still need to try changing all the backends.

Also, on Wednesday I started a new file called utils.py; in here, I'll add random functions that will be useful during this summer/in general so that I won't have to retype them every time I need them. Currently only hosting:
-a distance function that finds distance in 3-space given x,y,z coordinates
-a coordDist function that finds distances in 2-space given two points (tuples)
-a setfig function that has been rather abandoned as Python seems to automatically set the figure number of each graph

Friday
No plotting exercise today! Today I made a file called makecirc.py
-It hosts one function. Given a 2-D data array, it makes an array of the same dimensions and fills them with zeros. Then within a circle with a fixed center and radius specified by the user, it fills the circle with ones. This will definitely be useful when I  have to plot the brightness of stars with respect to a radius (distance from the center of the star).
Here's an example:
Given:


>>>from numpy import*
>>>import makecirc.py
>>>dataArr=zeros((10,10))
>>> dataArr
array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

>>> makecirc.makecirc(dataArr,5,4,3)
[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  0.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  0.  0.  0.]
 [ 0.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  0.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  0.  0.  0.]
 [ 0.  0.  0.  0.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

Tim also helped me install PyFITS on my computer, which is nice as now I open up pretty pictures of stars on my computer (the FITS files). I also learned some command line code as we had so much trouble installing PyFITS on a Windows computer.

Stuff Learned:
-Too much to write on here. Mostly how to use different functions of matplotlib and how making graphics in Python works. I wrote them in my notebook.
Problems: 
-matplotlibrc won't work.
To Do:
-Next week I'll start making a 2-D Gaussian plot.
-Learn more about how to use PyFITS.
-Read assigned research papers.

Test Scripts in Python-I (Monday last week)

Tim was kind enough to provide me with a few Python exercises to familiarize myself with Python's syntax, as I originally programmed in Java.

My first ever Python file contained four simple functions:
-A swap function that swapped two variables.
-A fun 8-ball function that randomly gave you an answer out of an array of six types of answers
-A Fibonacci function that printed Fibonacci numbers up to a maximum specified value
-And a modified Fibonacci numbers that printed a user-selected-character up to a maximum specified time

The rest of the week would be spent on plotting exercises.

Stuff Learned:
-General Python syntax
-How to use tuples
-Python is FUN! (And so much better than Java)

Project Overview

Hi!

Over the next ten weeks (well, nine weeks from now as this post is one week late), I'll be working in Professor John Johnson's Exolab with Tim Morton at Caltech's astrophysics department. The Exolab is interested in the study of extrasolar planets.

PROJECT
Background: Transits and False Positives
An extrasolar planet is commonly detected through studying its transit across a star. Generally, astronomers find this planet by plotting the flux of a star with respect to time and finding the transit signal. If there is a dip in the amount of flux, then something could be blocking the light of the star; hopefully that something is a small planet moving across the star's surface from our point of view. Typically, however, the amount of light blocked is very small. The value is called delta, which varies with (Rplanet/Rstar)^2. For a Jupiter-sized planet, delta ~ 0.01, and for an Earth-sized planet, delta ~10^-4!

However, there are other things that can cause the flux to dip. Binary star systems can also cause the flux to lower, though delta is usually much greater for binary stars in comparison to planets (delta ~ 0.1). Also, for binary star systems, if we plot the flux vs time graph long enough, there should be a second dip as the other companion gets blocked behind the first star. Under certain conditions, however, binary star systems can be misinterpreted as extrasolar planets transiting across stars due to the way their transit signal looks, or the way the light is collected, leading to what we call a false positive. 

Solution (Where I come in)
In order to avoid interpreting false positives as transiting planets, we need to be able to better plot the light from each star as a function of its distance. Thus, for my project, I will be analyzing CCD images taken by the Palomar 60'' telescope fitted with the Robo-AO and plotting the brightness across each CCD image. If there are any anomalies in the light curve, such as a sharp spike in brightness far away from the star, we can conclude that there are other sources contributing to the amount of light in the CCD image. This light can then be subtracted out (along with background noise) to get the light curve of the main star. This technique is called photometry, or the measurement of the star's flux.

In Tim's words, what I'm doing in a nutshell is determining "with a certain degree of confidence" if there are or aren't any bright objects in a patch of sky that can influence the transit signal.

Approach
I will be doing photometry using CCD's/FITS files and writing scripts with Python.

Background: CCDs
A CCD is a device that contains an array of "wells" capable of holding around ~40k photons per well. As a certain number of photons hits the well, an electron will be emitted and registered onto a pixel. These pixels, along with some other information, form an image called a FITS file. The ratio of electrons emitted (in ADU's, or Astronomical Data Units) to the photons absorbed (or counts) is called the Gain of the CCD.

Though a star is considered as a point source of light, due to atmospheric effects, or seeing, the light of the star is smeared out across a circle of pixels determined by the point spread function. The point spread function should theoretically be determined by the resolution (wavelength/diameter of the telescope); however due to seeing it is not. Thus, the FITS files we obtain of the stars will have a bright center surrounded by rings of increasingly dim pixels.

Plotting
We can plot the intensity of the light with respect to radius and obtain a bell curve. Significant deviations from the bell curve can be marked as potential sources of light and subtracted out. However, a trickier problem is finding the brightest an object can be and yet not be detected at five times the noise (which we will call sigma). This has to do with measuring and subtracting the background noise of the sky.

The final plot I will make will be a plot of the magnitude of the star as distance from its center increases. If there are no other significant sources of light, this plot should be smoothly declining.