#PH4100 Major Poject - Photometry - Star field simulation program (5) #Aaron Andrews #Last edited: 18/10/2016 #('*' indicates newest addition) #This program generates a set number of stars with random positions and max. grey values. #Stars are assumed to have 2d gaussian profiles. *correlation coefficient for x and y is an adjustable variable #Mean value for each pixel (per star) calculated via integration over pixel area #Mean value used to generate random poisson value for each pixel #The output is a 2d array representing an image of a star field - each cell stores a grey value calculated using the cumulative gaussian profiles for all stars generated. #The star field representation generated by this particular program is not accurate. #*Output can be saved as FITS file import math import random as rm import numpy as np import matplotlib import matplotlib.pyplot as plt import matplotlib.image as mpimg import scipy.integrate as si #for integration of gaussian function over pixels from astropy.io import fits from textwrap import wrap #The first two functions are used for calculations of star gaussian profiles def Gaussian(x,y,xStar,yStar,xSig,ySig,corr,peak): #2d gaussian function, returns grey value contribution from a single star for any point in the star field return peak*math.exp(-(1/(2*(1-corr**2)))*( (((x-xStar)**2)/(xSig**2)) + (((y-yStar)**2)/(ySig**2)) - 2*corr*((x-xStar)/xSig)*((y-yStar)/ySig) )) #normalisation: (1/(2*math.pi*xSig/8ySignp.sqrt(1-corr**2))) def IntGaussian(x1,y1,xStar,yStar,xSig,ySig,corr,peak): #returns integrated grey value contribution from a star for any pixel #integrate x and y components of gaussian seperately - define functions for integration first #GaussFuncX=lambda x:math.exp(-(((x-xStar)**2)/(2*xSig**2))) #GaussFuncY=lambda y:math.exp(-(((y-yStar)**2)/(2*ySig**2))) if corr==0: GaussFunc=lambda y,x:math.exp(-(((x-xStar)**2)/(2*xSig**2)))*math.exp(-(((y-yStar)**2)/(2*ySig**2))) #Failsafe for integration - without it, would cause integration to return value of 1 else: GaussFunc=lambda y,x:math.exp(-(1/(2*(1-corr**2)))*( (((x-xStar)**2)/(xSig**2)) + (((y-yStar)**2)/(ySig**2)) - 2*corr*((x-xStar)/xSig)*((y-yStar)/ySig) )) #Take pixel coordinates as centre of pixel - integrate between range of [-0.5,+0.5] pixels #xIntegral,xIntError=si.quad(GaussFuncX,x1-0.5,x1+0.5)#as numerical integration, returns error #yIntegral,yIntError=si.quad(GaussFuncY,y1-0.5,y1+0.5) Integral,IntError=si.dblquad(GaussFunc,x1-0.5,x1+0.5,lambda y:y1-0.5,lambda y:y1+0.5)#2d numerical integration return (peak*Integral) #The following functions can be used to manually create a star field and generate stars with set parameters. #(Optional) Generated image can be converted into array of discrete (poisson) pixel values, and saved as a FITS file. def MakeField(nxPixel,nyPixel): #*Generate unfilled array for star field return np.zeros((nxPixel,nyPixel)) def AddStar(posX,posY,sigX,sigY,corr,maxN,image): #*adds a star to an existing image Image1=image #Loop over all pixels in image pxlX=0 while pxlX