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

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.
    """
    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:
            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)
                    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)

def ParetoDist(x):
    return 1./pow((1+x),2)
    
def ParetoDistMonteCarlo():
    xmax = 4.0
    ymax = 1.0
    Accepted = False
    while Accepted==False:
        x = _r.random()*xmax
        y = _r.random()*ymax   
        if y<ParetoDist(x):
            Accepted = True
        else:
            Accepted = False
    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

def MakeAStarField(CanvasSize, NumStars, globular, shape):
    """Function that simulates the image of a star."""
    BrightStar = 65000
    StarPos = StarfieldPositions(CanvasSize, NumStars, globular) 

    Max_x = CanvasSize[0]
    Max_y = CanvasSize[1]
    
    Rows = _np.arange(0,Max_y,1)
    ListofColumns = []
    for i in Rows:
        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)
        i = _np.random.normal(i)
        ListofRows.append(i)
    Img = _np.array(ListofRows)
     
    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
            if ExpectedValue<0:
                ExpectedValue*=-1
            n = _np.random.poisson(ExpectedValue)
            Img[l] += n
            
        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)     

A = MakeAStarField((1250,1500), 100, globular = False, shape = False)
MakeFits(A, "3_Stars.fits")
_plt.figure()
_plt.imshow(A)
_plt.colorbar()
_plt.show()

- - BekiChafer - 22 Oct 2017</verbatim>