#PH4100 Major Poject - Photometry - Star field simulation program (3) #Aaron Andrews #Last edited: 011/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. #*Mean value for each pixel (per star) calculated via integration over pixel area #* #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. 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 textwrap import wrap def Gaussian(x,y,xStar,yStar,xSig,ySig,peak): #2d gaussian function, returns grey value contribution from a single star for any point in the star field return peak*math.exp(-((((x-xStar)**2)/(2*xSig**2))+(((y-yStar)**2)/(2*ySig**2)))) def IntGaussian(x,y,xStar,yStar,xSig,ySig,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))) #Take pixel coordinates as centre of pixel - integrate between range of [-0.5,+0.5] pixels xIntegral,xIntError=si.quad(GaussFuncX,x-0.5,x+0.5)#as numerical integration, returns error yIntegral,yIntError=si.quad(GaussFuncY,y-0.5,y+0.5) return (peak*xIntegral*yIntegral) def SimStarField(nxPixel,nyPixel,MaxPixelValue,nStars): #Create empty arrays Image1=np.zeros((nxPixel,nyPixel)) #array to be output once filled - each cell contains grey values (no. of photons detected) PosX=np.zeros(nStars) #x coordinates PosY=np.zeros(nStars) #y coordinates Peak=np.zeros(nStars) #peak grey value for each star sigX=3 #standard deviations, set to arbitrary values (for now) sigY=3 #Loop over all stars n=0 while n