#PH4100 Major Poject - Photometry - Star field simulation program (2)
#Aaron Andrews
#Last edited: 05/10/2016

#This program generates a set number of stars with random positions and max. grey values.
#Stars are assumed to have 2d gaussian profiles.
#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

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 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

    #Loop over all stars
    n=0
    while n<nStars:
        #Random value generation (eventually will be replaced with retrieval from database of stars)
        PosX[n]=rm.random()*nxPixel #generates random value between 0 and 1, multiplies by relevent range
        PosY[n]=rm.random()*nyPixel
        Peak[n]=rm.random()*MaxPixelValue
        
        #Image1[int(PosX[n]),int(PosY[n])]=1
        #Loop over all pixels in image
        pxlX=0
        while pxlX<nxPixel:
            pxlY=0
            while pxlY<nyPixel:
                #For every pixel, a grey value as determined by the gaussian profile of each star is calculated, then added to the existing grey value of the pixel.
                Image1[pxlX,pxlY]=Image1[pxlX,pxlY]+int(Gaussian(pxlX,pxlY,PosX[n],PosY[n],3,3,Peak[n])) #Note: arbitrary values (3 pixels) have been input for sigma.
                pxlY=pxlY+1
            pxlX=pxlX+1
        
        n=n+1
    return Image1
    
    

#Defining example variables
ImageX=640 #pixels, x axis
ImageY=480 #pixels, y axis
MaxGrey=65535

Nstars=50 #number of stars to generate

SimStarField(ImageX,ImageY,MaxGrey,Nstars)