# Difference: StewartBoogertPhotometryStarFieldSimulation (1 vs. 6)

#### Revision 609 Jan 2018 - RebekahChafer1

Line: 1 to 1

 META TOPICPARENT name="StewartBoogertPhotometry2017"
Deleted:
<
<
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```
>
>
import scipy as _sc #from scipy.stats import skewnorm

Changed:
<
<
def StarfieldPositions(CanvasSize, NumStars):
>
>
def StarfieldPositions(CanvasSize, NumStars, Globbed):
"""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. """
Deleted:
<
<
if CanvasSize == "Auto": CanvasSize = (1530,1020)
Max_x = CanvasSize[0] Max_y = CanvasSize[1] A = _np.linspace(1,NumStars, NumStars) X = [] Y = []
>
>
Pos = []

>
>
if Globbed == True: GaussMean = 0.5 GaussSig = 0.2 y_max = 1.6
for i in A:
Changed:
<
<
x = round(_r.random()*Max_x) y = round(_r.random()*Max_y)
>
>
Acceptedx = False while Acceptedx == False: x_Monte = _r.random() y_Monte = _r.random()*y_max if y_Monte<Gaussian(x_Monte, GaussMean, GaussSig): x = round(x_Monte*Max_x)
X.append(x)
>
>
Acceptedx = True else: Acceptedx = False Acceptedy = False while Acceptedy == False: x_Monte = _r.random() y_Monte = _r.random()*y_max if y_Monte<Gaussian(x_Monte, GaussMean, GaussSig): y = round(x_Monte*Max_x)
Y.append(y)
Changed:
<
<
#_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
>
>
Acceptedy = True else: Acceptedy = False else: for i in A: x = round(_r.random()*Max_x) y = round(_r.random()*Max_y) Pos.append((x,y))

return Pos

def MovePixel(Pos, dx, dy): NewPos = (Pos[0] + dx, Pos[1] + dy) return NewPos

def CentretoCorner(Centre, SizeOfBox): return MovePixel(Centre, -SizeOfBox, SizeOfBox)

def MoveRow(Array, CurrentPixel, SizeOfBox, List): Max_x, Max_y = _np.shape(Array) for i in range(0,SizeOfBox*2): CurrentPixel = MovePixel(CurrentPixel,+1, 0) if(CurrentPixel[0]>=0 and CurrentPixel[1]>=0 and CurrentPixel[0]<=Max_x-1 and CurrentPixel[1]<=Max_y-1): List.append(CurrentPixel) return CurrentPixel, List

def ChangeRow(CurrentPixel, SizeOfBox): return MovePixel(CurrentPixel, -2*SizeOfBox , -1)

def GetBoxPixelPos(Array, CentrePos, SizeOfBox): """Returns a list of pixel positions of a box around a centre pixel""" Max_x, Max_y = _np.shape(Array) Box = [] Position = CentretoCorner(CentrePos, SizeOfBox) if(Position[0]>=0 and Position[1]>=0 and Position[0]<=Max_x-1 and Position[1]<=Max_y-1): Box.append(Position) Position, Box = MoveRow(Array,Position, SizeOfBox, Box) for i in range(0,(2*SizeOfBox)-1): Position = ChangeRow(Position, SizeOfBox) if(Position[0]>=0 and Position[1]>=0 and Position[0]<=Max_x-1 and Position[1]<=Max_y-1): Box.append(Position) Position, Box = MoveRow(Array, Position, SizeOfBox, Box) return Box

def phi(x): return 1./(_np.sqrt(2*_np.pi))*_np.exp(-(x)**2/2.)

def PHI(x): return (1./2.)*(1+_sc.special.erf(x/_np.sqrt(2.)))

def SkewNormPdf(x, a, mean, sigma): return (2./sigma)*phi((x-mean)/sigma)*PHI(a*((x-mean)/sigma))

def Gaussian(x, mean, sigma): return 1./(_np.sqrt(2*_m.pi*sigma**2)) * _np.exp(-(x - mean)**2/(2*sigma**2))

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 Gaussian2DWCorr(A, x, y, x_pos, y_pos, sigx, sigy, P, r): """Function that returns the value of a 2 dimensional Gaussian Function""" X = ((x-x_pos)/(2*sigx**2)) Y = ((y-y_pos)/(2*sigy**2)) PCoeff = (1./(2.*(1-P**2))) if r<0.5: return A * _np.exp(-PCoeff*(X**2 + Y**2 + 2*P*X*Y) ) if r>=0.5: return A * _np.exp(-PCoeff*(X**2 + Y**2 - 2*P*X*Y) )

def Poisson(E, n): return ((E**n)/(_m.factorial(n)))*_np.exp(-E)
Changed:
<
<
def montecarloPoisson(E): n = _np.linspace(0,100,101) f = [] for i in n: a = Poisson(E,i) f.append(a)
>
>
def ParetoDist(x): return 1./pow((1+x),2)

def ParetoDistMonteCarlo(): xmax = 4.0 ymax = 1.0

Accepted = False while Accepted == False:
Changed:
<
<
r1 = _r.random() r2 = _r.random() X = round(min(n) + r1 *(max(n)-min(n))) Y = r2*max(f) if Y<Poisson(E,X):
>
>
x = _r.random()*xmax y = _r.random()*ymax if y<ParetoDist(x):
Accepted = True else: Accepted = False
Changed:
<
<
return X
>
>
return x/4.

def bgFunc(A): """A function of the bg, determined by a fit from real image.""" f0,f1,f2 = -4.23943734841e-05, 0.292931970885, 566.077397528 y = f0*pow(A,2) + f1*A + f2 return y

def getBox1(array, x, y ,dim): """x, y = centre of box, dim must be odd."""

if dim%2 == 0: dim-=1

Diag = (dim-1)/2. ULCx = x-Diag ULCy = y-Diag

box = array[ULCx:ULCx+dim,ULCy:ULCy+dim] return box

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

CounterX = -1

>
>
StarPos = StarfieldPositions(CanvasSize, NumStars, globular)
Max_x = CanvasSize[0] Max_y = CanvasSize[1]
Deleted:
<
<
#Img = _np.zeros((Max_x, Max_y))

Changed:
<
<
#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)
>
>
Rows = _np.arange(0,Max_y,1)
ListofColumns = []
Changed:
<
<
for i in Columns:
>
>
for i in Rows:
A = _np.zeros(Max_x) Value = bgFunc(i) A+=Value
Line: 91 to 172
ListofRows = [] for i in NoFlucImg: i = _np.random.poisson(i)
>
>
i = _np.random.normal(i)
ListofRows.append(i) Img = _np.array(ListofRows)

Changed:
<
<
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)
>
>
StarCount = 1 Gausssigma = 0.15 YMAX= 2.7

Gaussmean_P = 0

if shape == True: ShapeRandom = _r.random() P = _r.random()*2 -1 ShapeAccepted = False while ShapeAccepted == False: r_val = _r.random() y_val = r_val*YMAX if y_val<Gaussian(P, Gaussmean_P, Gausssigma): ShapeAccepted = True else: ShapeAccepted = False

StarSigma = 4

for k in StarPos: X = ParetoDistMonteCarlo() Stot = X*BrightStar box = GetBoxPixelPos(Img, k, StarSigma*4) X,Y = _np.shape(Img)

for l in box: i = l[0] j = l[1] if shape == True: p = Gaussian2DWCorr(BrightStar,i,j,k[0], k[1], StarSigma,StarSigma, P, ShapeRandom)*(1./Max_x)*(1./Max_y) else: p = Gaussian2D(BrightStar,i,j,k[0], k[1], StarSigma,StarSigma)*(1./Max_x)*(1./Max_y)

S = p*Stot ExpectedValue = S
Changed:
<
<
#n = _np.random.poisson(ExpectedValue) #n = montecarloPoisson(ExpectedValue) Img[i-1,j-1] += ExpectedValue
>
>
if ExpectedValue<0: ExpectedValue*=-1 n = _np.random.poisson(ExpectedValue) Img[l] += n

print "StarCount:", StarCount StarCount+=1 return Img
Line: 120 to 228
hdulist = fits.HDUList([hdu]) hdulist.writeto(Filename)
Changed:
<
<
def bgFunc(A): f0,f1,f2 = -4.23943734841e-05, 0.292931970885, 566.077397528 y = f0*pow(A,2) + f1*A + f2 return y
>
>
A = MakeAStarField((1250,1500), 100, globular = False, shape = False) MakeFits(A, "3_Stars.fits") _plt.figure() _plt.imshow(A) _plt.colorbar() _plt.show()

#### Revision 510 Nov 2017 - RebekahChafer1

Line: 1 to 1

 META TOPICPARENT name="StewartBoogertPhotometry2017"
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.
Line: 61 to 61
Accepted = False return X
Changed:
<
<
def MakeAStarField(CanvasSize, Stot, bg, NumStars):
>
>
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]
Line: 71 to 72
Max_x = CanvasSize[0] Max_y = CanvasSize[1]
Changed:
<
<
Img = _np.zeros((Max_x, Max_y))
>
>
#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:
Changed:
<
<
CounterX+=1 b = StarPos_y[CounterX]
>
>
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
Changed:
<
<
ExpectedValue = S+bg n = montecarloPoisson(ExpectedValue) Img[i-1,j-1] += n
>
>
ExpectedValue = S #n = _np.random.poisson(ExpectedValue) #n = montecarloPoisson(ExpectedValue) Img[i-1,j-1] += ExpectedValue print "StarCount:", StarCount StarCount+=1
return Img
Line: 94 to 120
hdulist = fits.HDUList([hdu]) hdulist.writeto(Filename)
Changed:
<
<
A = MakeAStarField((500, 500), 150, 25, 5) MakeFits(A,"TrialPoissonStar.fits")
>
>
def bgFunc(A): f0,f1,f2 = -4.23943734841e-05, 0.292931970885, 566.077397528 y = f0*pow(A,2) + f1*A + f2 return y

#### Revision 401 Nov 2017 - RebekahChafer1

Line: 1 to 1

 META TOPICPARENT name="StewartBoogertPhotometry2017"
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.
Line: 6 to 6
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
Line: 39 to 40
"""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)))
Changed:
<
<
def MakeAStar(CanvasSize, Stot, bg, NumStars):
>
>
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, Stot, bg, NumStars):

"""Function that simulates the image of a star.""" StarPos = StarfieldPositions(CanvasSize, NumStars) StarPos_x = StarPos[0]
Line: 50 to 72
Max_x = CanvasSize[0] Max_y = CanvasSize[1] Img = _np.zeros((Max_x, Max_y))
Deleted:
<
<
Img+=bg
for a in StarPos_x: CounterX+=1
Line: 61 to 82
for j in Y: p = Gaussian2D(1,i,j,a, b, 5,5) S = p*Stot
Changed:
<
<
Img[i-1,j-1] += S
>
>
ExpectedValue = S+bg n = montecarloPoisson(ExpectedValue) Img[i-1,j-1] += n
return Img
>
>
def MakeFits(Array, Filename):
>
>
"""Converts array into a fits image"""
hdu = fits.PrimaryHDU(Array) hdulist = fits.HDUList([hdu]) hdulist.writeto(Filename)
Changed:
<
<
A = MakeAStar((1530,1020), 500, 50, 20) MakeFits(A, "50Stars.fits") _plt.imshow(A) _plt.show()
>
>
A = MakeAStarField((500, 500), 150, 25, 5) MakeFits(A,"TrialPoissonStar.fits")

#### Revision 325 Oct 2017 - RebekahChafer1

Line: 1 to 1

 META TOPICPARENT name="StewartBoogertPhotometry2017"
Changed:
<
<
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 functon implemented for one star
>
>
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```
Line: 22 to 22
Y = []

for i in A:

Changed:
<
<
x = _r.random()*Max_x y = _r.random()*Max_y
>
>
x = round(_r.random()*Max_x) y = round(_r.random()*Max_y)
X.append(x) Y.append(y)
Changed:
<
<
_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])
>
>
#_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)))

Changed:
<
<
def MakeAStar(CanvasSize, StarPos, Stot, bg):
>
>
def MakeAStar(CanvasSize, Stot, bg, NumStars):
"""Function that simulates the image of a star."""
>
>
StarPos = StarfieldPositions(CanvasSize, NumStars) StarPos_x = StarPos[0] StarPos_y = StarPos[1]

CounterX = -1

Max_x = CanvasSize[0] Max_y = CanvasSize[1]
Deleted:
<
<
StarPos_x = StarPos[0] StarPos_y = StarPos[0]
Img = _np.zeros((Max_x, Max_y))
Changed:
<
<
Stot = 500 bg = 50
>
>
Img+=bg

>
>
for a in StarPos_x: CounterX+=1 b = StarPos_y[CounterX]
X = _np.linspace(0,Max_x, Max_x+1) Y = _np.linspace(0,Max_y, Max_y+1) for i in X: for j in Y:
Changed:
<
<
p = Gaussian2D(1,i,j,StarPos_x, StarPos_y, 5,5)
>
>
p = Gaussian2D(1,i,j,a, b, 5,5)
S = p*Stot
Changed:
<
<
v = S+bg Img[i-1,j-1] = v
>
>
Img[i-1,j-1] += S
return Img

def MakeFits(Array, Filename):

Line: 65 to 69
hdulist = fits.HDUList([hdu]) hdulist.writeto(Filename)
>
>
A = MakeAStar((1530,1020), 500, 50, 20) MakeFits(A, "50Stars.fits") _plt.imshow(A) _plt.show()

- - BekiChafer - 22 Oct 2017</verbatim>

#### Revision 225 Oct 2017 - RebekahChafer1

Line: 1 to 1

 META TOPICPARENT name="StewartBoogertPhotometry2017"
Changed:
<
<
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.
>
>
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 functon implemented for one star

```import numpy as _np
import random as _r
import matplotlib.pyplot as _plt```
Changed:
<
<
#Variable Canvas size and number of stars
>
>
from astropy.io import fits
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)
Line: 20 to 26
y = _r.random()*Max_y X.append(x) Y.append(y)
>
>
_plt.figure(facecolor = "w")
_plt.scatter(X,Y, marker = '*')
Changed:
<
<
_plt.title("Position of stars in simulated image")
>
>
_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
Changed:
<
<
Canvas = (100,100)
>
>
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 MakeAStar(CanvasSize, StarPos, Stot, bg): """Function that simulates the image of a star."""

Max_x = CanvasSize[0] Max_y = CanvasSize[1] StarPos_x = StarPos[0] StarPos_y = StarPos[0] Img = _np.zeros((Max_x, Max_y)) Stot = 500 bg = 50

X = _np.linspace(0,Max_x, Max_x+1) Y = _np.linspace(0,Max_y, Max_y+1) for i in X: for j in Y: p = Gaussian2D(1,i,j,StarPos_x, StarPos_y, 5,5) S = p*Stot v = S+bg Img[i-1,j-1] = v return Img

def MakeFits(Array, Filename): hdu = fits.PrimaryHDU(Array) hdulist = fits.HDUList([hdu]) hdulist.writeto(Filename)

Deleted:
<
<
StarfieldPositions(Canvas,100) _plt.show()

- - BekiChafer - 22 Oct 2017</verbatim>

#### Revision 122 Oct 2017 - RebekahChafer1

Line: 1 to 1
>
>
 META TOPICPARENT name="StewartBoogertPhotometry2017"
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.
```import numpy as _np
import random as _r
import matplotlib.pyplot as _plt

#Variable Canvas size and number of stars

def StarfieldPositions(CanvasSize, NumStars):
Max_x = CanvasSize[0]
Max_y = CanvasSize[1]
A = _np.linspace(1,NumStars, NumStars)
X = []
Y = []

for i in A:
x = _r.random()*Max_x
y = _r.random()*Max_y
X.append(x)
Y.append(y)
_plt.scatter(X,Y, marker = '*')
_plt.title("Position of stars in simulated image")
_plt.xlabel("x position")
_plt.ylabel("y position")
_plt.xlim(0,CanvasSize[0])
_plt.ylim(0,CanvasSize[1])
return X, Y

Canvas = (100,100)

StarfieldPositions(Canvas,100)
_plt.show()
```

- - BekiChafer - 22 Oct 2017</verbatim>

Copyright © 2008-2022 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