improved statistics and plot

changed execute to a constant number of stars instead of a constant sfr
added better plotting ability each 2d-plane of the data grid is
accessible
added labels and discriptions to the plots

execute.py

@@ -66,7 +66,11 @@ def g(x):
         return sfr
     g = np.vectorize(g)
     ages = np.logspace(5,7,7)
-    sf = [dist.distribution(g, 1000., ages[i]) for i in range(len(ages))]
+    #sf = [dist.distribution(g, 1000., ages[i]) for i in range(len(ages))]
+    sf = []
+    for i in range(len(ages)):
+        sfr = 150000*mf.mean()/(ages[i]-1000.)           # using sfr so that we expect ~150000 stars
+        sf.append( dist.distribution(g, 1000., ages[i]))
 
     t1 = time()                 # finished reading the distributions
     print(t1,file=output_stream)

plot.py

@@ -9,7 +9,7 @@
 from astropy.io import fits
 
 
-def main(aperature = 1, visual_ex = 0):
+def main(aperature = 1, visual_ex = 0, age=0):
     '''main(aperature = 1, visual_ex = 0)
 
 Parameters
@@ -59,7 +59,7 @@ def main(aperature = 1, visual_ex = 0):
     ax.set_xlabel('av')
     ax.set_ylabel('log(age)')
     ax.set_zlabel('starformation rate in M_sun/year')
-    ax.set_title('Starformation as a function of\n visual extinction Av and age at apperaturesize=1kpc')
+    ax.set_title('Starformation as a function of\n visual extinction Av and age at apperaturesize=%s' % aperas[aperature])
 
     plt.savefig('plot/3dav-age.svg')
 
@@ -80,24 +80,46 @@ def main(aperature = 1, visual_ex = 0):
     ax2.set_xlabel('log(age)')
     ax2.set_ylabel('log(apera)')
     ax2.set_zlabel('starformation rate in M_sun/year')
-    ax2.set_title('Starformation as a function of\napperaturesize and age at Av=10')
+    ax2.set_title('Starformation as a function of\napperaturesize and age at Av=%s' % avs[visual_ex])
 
     plt.savefig('plot/3dage-apera.png')
 
 
+    fig3 = plt.figure()
+    ax3 = fig3.gca(projection='3d')
+    ax3.view_init(45,-135)
+
+    X3, Y3 = np.meshgrid(np.log10(aperas), avs)
+    Z3 = 559.*sfr/numbers[:,:,age]
+    surf3 = ax3.plot_surface(X3, Y3, Z3, rstride=1, cstride=1, cmap=cm.coolwarm)
+    cbar3 = fig3.colorbar(surf3)
+    cbar3.set_label('starformation rate in M_sun/year')
+
+    x3 = X3.reshape(numbers.shape[0]*numbers.shape[1])
+    y3 = Y3.reshape(numbers.shape[0]*numbers.shape[1])
+    z3 = Z3.reshape(numbers.shape[0]*numbers.shape[1])   
+    ax3.scatter(x3,y3,z3)
+    ax3.set_xlabel('log(apera)')
+    ax3.set_ylabel('AV')
+    ax3.set_zlabel('starformation rate in M_sun/year')
+    ax3.set_title('Starformation as a function of\napperaturesize and AV at age=%s years' % ages[age])
+
+    plt.savefig('plot/3dav-apera.png')
+
+
 
 
 def cmd(folder, av, apera, age, color1 = "I4", color2 = "M1", corrected=True, old=False, color='mass'):
     '''cmd(folder, av, apera, age, color1 = "I4", color2 = "M1") - creates a CMD
 
 Parameters
 ----------
-folder   String
-av       float
-apera    float
-age      float
-color1   String
-color2   String
+folder   String  folder in which the datafile is to be found
+av       float   value for the visual extinction as in the filename
+apera    float   value of the aperture size as in the filename 
+age      float   value of the age as in the filename
+color1   String  band filter to be used for the first color 
+color2   String  band filter to be used for the second color
 
 Returns
 ----------
@@ -146,7 +168,7 @@ def cmd(folder, av, apera, age, color1 = "I4", color2 = "M1", corrected=True, ol
     plt.savefig('plot/%s_%s_%s.png' %(av,apera,age))
 
 
-def plot_2d(aperature = 1, visual_ex = 0):
+def plot_2d(aperature = 1, visual_ex = 0, age = 0):
     '''plot_2d(aperature = 1, visual_ex = 0) - creates two 2d contour plots
 
 Parameters
@@ -160,22 +182,6 @@ def plot_2d(aperature = 1, visual_ex = 0):
 of the data in 'out/__expected_number'
 '''
 
-    #f = open('out/__head', 'r')
-    #headers = f.readline().strip().split(',')[1:]
-    #data = np.loadtxt(f)
-    #f.close()
-    #avs = np.unique(data[:,headers.index('AV')])
-    #aperas = np.unique(data[:,headers.index('Aperature_size')])
-    #ages = np.unique(data[:,headers.index('Age')])
-
-    #hdulist = fits.open('%s/%s' %('out','sim_%03d_%06d_%09d'%(avs[0],aperas[0],ages[0])))
-    #sfr = hdulist[2].data['SFR'][0]
-
-    #f = open('out/__expected_number', 'r')
-    #headers = f.readline().strip().split(',')[1:]
-    #data = np.loadtxt(f)
-    #f.close()
-
     avs, aperas, ages, sfr, headers, data = init()
 
     numbers = np.clip(data[:,3].reshape(len(avs), len(aperas), len(ages)),1.,1000000.) #prevent inf if no stars are detected
@@ -192,24 +198,26 @@ def plot_2d(aperature = 1, visual_ex = 0):
 
     cs = ax.contourf(X,Y,Z)
     #cs1 = ax.contour(X,Y,Z,[0.01,.011,.012,.013,.014,.015],colors=['red','red','red','red','red','red'])
-    cs1 = ax.contour(X,Y,Z,levels=[0.01,.011,.012,.013,.014,.015],cmap=plt.cm.jet)
+    cs1 = ax.contour(X,Y,Z,levels=[0.1,.125,.15,.175,.2], colors=['red','red','red','red','red'])
+    ax.clabel(cs1, fontsize=10, inline=1)
     cbar = plt.colorbar(cs)
 
     ax.scatter(x,y)
     ax.set_xlabel('av')
     ax.set_ylabel('log(age)')
     cbar.set_label('starformation rate in M_sun/year')
-    ax.set_title('Starformation as a function of visual extinction Av\nand age at apperaturesize=1kpc')
+    ax.set_title('Starformation as a function of visual extinction Av\nand age at apperaturesize=%s' % aperas[aperature])
 
-    plt.savefig('plot/2dav-age.png')
+    plt.savefig('plot/2dav-age-%s.png' % aperas[aperature])
 
     fig2 = plt.figure()
     ax2 = fig2.gca()
 
     X2, Y2 = np.meshgrid(np.log10(ages), np.log10(aperas))
     Z2 = 559.*sfr/numbers[visual_ex,:,:]
     cs2 = ax2.contourf(X2, Y2, Z2)
-    cs21 = ax2.contour(X2,Y2,Z2,[0.01,.011,.012,.013,.014,.015])
+    cs21 = ax2.contour(X2,Y2,Z2,levels=[0.1,.125,.15,.175,.2], colors=['red','red','red','red','red'])
+    ax2.clabel(cs21, fontsize=10, inline=1)
     cbar2 = plt.colorbar(cs2)
 
     x2 = np.log10(data[:,2].reshape(len(avs), len(aperas), len(ages))[visual_ex,:,:].reshape(numbers.shape[1]*numbers.shape[2]))
@@ -218,10 +226,32 @@ def plot_2d(aperature = 1, visual_ex = 0):
     ax2.set_xlabel('log(age)')
     ax2.set_ylabel('log(apera)')
     cbar2.set_label('starformation rate in M_sun/year')
-    ax2.set_title('Starformation as a function of\n apperaturesize and age at Av=10')
+    ax2.set_title('Starformation as a function of\n apperaturesize and age at Av=%s' % avs[visual_ex])
+
+    plt.savefig('plot/2dage-apera-%s.png' % avs[visual_ex])
+
 
-    plt.savefig('plot/2dage-apera.png')
+    fig3 = plt.figure()
+    ax3 = fig3.gca()
 
+    X3, Y3 = np.meshgrid(np.log10(aperas), avs)
+    Z3 = 559.*sfr/numbers[:,:,age]
+    cs3 = ax3.contourf(X3, Y3, Z3)
+    cs31 = ax3.contour(X3, Y3, Z3,levels=[0.1,.125,.15,.175,.2], colors=['red','red','red','red','red'])
+    ax3.clabel(cs31, fontsize=10, inline=1)
+
+    cbar3 = fig3.colorbar(cs3)
+    cbar3.set_label('starformation rate in M_sun/year')
+
+    x3 = X3.reshape(numbers.shape[0]*numbers.shape[1])
+    y3 = Y3.reshape(numbers.shape[0]*numbers.shape[1])
+    ax3.scatter(x3,y3)
+    ax3.set_xlabel('log(apera)')
+    ax3.set_ylabel('AV')
+    ax3.set_title('Starformation as a function of\napperaturesize and AV at age=%s years' % ages[age])
+    cbar3.set_label('starformation rate in M_sun/year')
+
+    plt.savefig('plot/2dav-apera-%s.png' % ages[age])
 
 def init():
     '''init() - aquires the base data
@@ -246,23 +276,39 @@ def init():
     return avs, aperas, ages, sfr, headers, data
 
 
-def histogram(folder, av, apera, age, color1 = "I4", color2 = "M1"):
+def histogram(folder, av, apera, age):
+    '''histogram(folder, av, apera, age, color1 = "I4", color2 = "M1")
+
+Parameters
+----------
+folder   String  folder in which the datafile is to be found
+av       float   value for the visual extinction as in the filename
+apera    float   value of the aperture size as in the filename 
+age      float   value of the age as in the filename
+
+Returns
+----------
+A histogram of the massdistribution in logspace with 20 bins from -.5 to 1.5 M_sun
+of all the sampled stars and restricted to the selected stars 
+'''
 
     hdulist = fits.open('%s/%s' %(folder,'sim_%03d_%06d_%09d'%(av,apera,age)))
     data = hdulist[1].data
 
-
-
     fig = plt.figure()
     ax = fig.add_subplot(111)
 
 
-    pdf, bins, patches = ax.hist(data['mass'], range=(-.5,1.5))
+    pdf, bins, patches = ax.hist(data['mass'], range=(-.5,1.5), label='all', bins=20)
     dat = data['mass'][data['sel'] == 1]
-    pdf, bins, patches = ax.hist(dat, range=(-.5,1.5))
+    pdf, bins, patches = ax.hist(dat, range=(-.5,1.5), label='selected', bins=20)
+    ax.legend()
     #ax.set_yscale('log')
+    ax.set_ylabel('number of stars')
+    ax.set_xlabel('log(mass/M_sun)')
+    ax.set_title('mass distribution simulated\nwith AV=%s, age=%s years, apertursize=%s AU' % (av,age,apera))
 
-    plt.savefig('plot/hist.png')
+    plt.savefig('plot/hist-%s-%s-%s.png' % (av,apera,age))
 
 
 

starformation.py

@@ -85,7 +85,7 @@ def f(x):               # Kouper IMF
         def g(x):
             return sfr
         g = np.vectorize(g)
-        starformationhistory = dist.distribution(g, 10000., maxage)
+        starformationhistory = dist.distribution(g, 1000., maxage)