Starfield Simulation program will be used to test star finding algorithm.
22/10/17 - Random position function of stars with variable canvas size and number of stars.
25/10/17 - Gaussian point spread function implemented for multiple stars, however, algorithm takes too long, need to make more efficient.
import numpy as _np
import random as _r
import matplotlib.pyplot as _plt
from astropy.io import fits
import math as _m

def StarfieldPositions(CanvasSize, NumStars):
    """A function, given a canvas size and number of stars, returns random 
        positions of stars in that canvas. 
        CanvasSize must be tuple of ints or "Auto", which makes it (1530, 1020), the size of the typical FITs image.
        NumStars must be an int.
    """
    if CanvasSize == "Auto":
        CanvasSize = (1530,1020)
    Max_x = CanvasSize[0]
    Max_y = CanvasSize[1]
    A = _np.linspace(1,NumStars, NumStars)
    X = []
    Y = []
    
    for i in A:
        x = round(_r.random()*Max_x)
        y = round(_r.random()*Max_y)
        X.append(x)
        Y.append(y)
    #_plt.figure(facecolor = "w")
    #_plt.scatter(X,Y, marker = '*')
    #_plt.title("Positions "+str(Max_x)+ " by "+str(Max_y)+" pixels with "+str(NumStars)+" simulated stars.")
    #_plt.xlabel("x position")
    #_plt.ylabel("y position")
    #_plt.xlim(0,CanvasSize[0])
    #_plt.ylim(0,CanvasSize[1])
    return X, Y

def Gaussian2D(A, x, y, x_pos, y_pos, sigx, sigy):
    """Function that returns the value of a 2 dimensional Gaussian Function"""
    return A*_np.exp(-((x-x_pos)**2/(2*sigx**2)+(y-y_pos)**2/(2*sigy**2)))

def Poisson(E, n):
    return ((E**n)/(_m.factorial(n)))*_np.exp(-E)

def montecarloPoisson(E):
    n = _np.linspace(0,100,101)
    f = []
    for i in n:
        a = Poisson(E,i)
        f.append(a)
    Accepted = False
    while Accepted == False:
        r1 = _r.random()
        r2 = _r.random()
        X = round(min(n) + r1 *(max(n)-min(n)))
        Y = r2*max(f)
        if Y<Poisson(E,X):
            Accepted = True
        else:
            Accepted = False 
    return X

def MakeAStarField(CanvasSize, NumStars):
    """Function that simulates the image of a star."""
    BrightStar = 65000
    StarPos = StarfieldPositions(CanvasSize, NumStars) 
    StarPos_x = StarPos[0]
    StarPos_y = StarPos[1]
    
    CounterX = -1

    Max_x = CanvasSize[0]
    Max_y = CanvasSize[1]
    #Img = _np.zeros((Max_x, Max_y))
    
    
    #background fluctuation dependent on position on image (bias to left of image)
    #CanvasSize = (3352,2352)
    #Max_x = CanvasSize[0]
    #Max_y = CanvasSize[1]
    Columns = _np.arange(0,Max_x,1)
    ListofColumns = []
    for i in Columns:
        A = _np.zeros(Max_x)
        Value = bgFunc(i)
        A+=Value
        ListofColumns.append(A)

    NoFlucImg = _np.rot90(_np.rot90(_np.rot90(_np.array(ListofColumns))))
    ListofRows = []
    for i in NoFlucImg:
        i = _np.random.poisson(i)
        ListofRows.append(i)
    Img = _np.array(ListofRows)
        
    StarCount = 0    
    for a in StarPos_x:
        Stot = _r.random()*BrightStar
        X = _np.linspace(0,Max_x, Max_x+1)
        Y = _np.linspace(0,Max_y, Max_y+1)
        CounterX+=1
        b = StarPos_y[CounterX]
        for i in X:
            for j in Y:
                p = Gaussian2D(1,i,j,a, b, 5,5)
                S = p*Stot
                ExpectedValue = S
                #n = _np.random.poisson(ExpectedValue)
                #n = montecarloPoisson(ExpectedValue)
                Img[i-1,j-1] += ExpectedValue
        print "StarCount:", StarCount
        StarCount+=1
    return Img


def MakeFits(Array, Filename):
    """Converts array into a fits image"""
    hdu = fits.PrimaryHDU(Array)   
    hdulist = fits.HDUList([hdu])
    hdulist.writeto(Filename)     

def bgFunc(A):
    f0,f1,f2 = -4.23943734841e-05, 0.292931970885, 566.077397528
    y = f0*pow(A,2) + f1*A + f2
    return y

- - BekiChafer - 22 Oct 2017</verbatim>

Physics WebpagesRHUL WebpagesCampus Connect • Royal Holloway, University of London, Egham, Surrey TW20 0EX; Tel/Fax +44 (0)1784 434455/437520

Topic revision: r5 - 10 Nov 2017 - BekiChafer

 
This site is powered by the TWiki collaboration platformCopyright © 2008-2017 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding RHUL Physics Department TWiki? Send feedback