docstrings

added docstrings to analysis, distribution, execute, setup and
starformation
fixed a bug in analysis (wrong detection area due to magnetudes)
fixed a bug in plot due to inf-values if no stars are detected

analysis.py

@@ -11,9 +11,19 @@
 
 
 def main(folder, quiet=0):
-    '''This script gives the number of "observed" stars from the sampled datafiles in "folder"/ 
+    '''analysis.main(folder, quiet=0)
+
+This script gives the number of "observed" stars from the sampled datafiles in "folder"/ 
 according to the selection criteria from Yusef-Zadeh et al
 
+input:
+folder  String   Specifiecs the folder where the files are
+quiet   boolean  =1 surpresses all standard output
+
+output:
+returns a file named __expected_number in which contains a numpy array of the simulation parameters
+number, A_v, Aperature_size, Age and the expected 'detected' number
+The first line of the file is an ','-seperated head of the contained informations
 '''
 
     if quiet:
@@ -36,21 +46,12 @@ def main(folder, quiet=0):
 #only using files created by the automated simulation
         if fil.startswith('sim_') and not 'settings' in fil.encode("ascii"):
             print ("%s/%s" % (folder,fil.encode("ascii")), file=output_stream)
-            #f = open("%s/%s" % (folder,fil.encode("ascii")), 'r')
-            #headers = f.readline().strip().split(',')
-            #data = np.loadtxt(f)
-            #f.close()
 
 
       # Read in
             hdulist = fits.open('%s/%s' %(folder,fil))
-#            av = hdulist[1].header['age']
             data = hdulist[1].data
 
-      ##selecting the relevant columns
-            #c1 = headers.index('corrected_flux %s' % color1)
-            #c2 = headers.index('corrected_flux %s' % color2)
-
       #calculating magnitudes from fluxes and converting to CMD-data
             x = -2.5*(np.log10(data['cflux %s' % color1]/64130) - np.log10(data['cflux %s' % color2]/7140))
             y = -2.5*(np.log10(data['cflux %s' % color2]/7140))
@@ -59,9 +60,10 @@ def main(folder, quiet=0):
           # efficiency? accuracy?
             n=0
       #selecting "observed" stars
-            for i in range(len(x)):
-                if (y[i] < -10./3. * (x[i]-1.) + 10.) and (max_mag < y[i] < min_mag):
-                    n = n+1
+            #for i in range(len(x)):
+                #if (y[i] > -10./3. * (x[i]-1.) + 10.) and (max_mag < y[i] < min_mag):
+                    #n = n+1
+            n = sum(np.logical_and( (y > -10./3. * (x-1.) + 10.), np.logical_and(max_mag < y, y < min_mag)))
             tmp, av, apera, age = fil.split('_')
             out.append([Decimal(av), Decimal(apera), Decimal(age), n])
 

distribution.py

@@ -9,8 +9,27 @@
 
 
 class distribution( object ):
-    '''to be added
+    '''distribution(func, a, b)
 
+A distribution object represents a probability distribution, the provided function is taken as the
+probability density function and an cumulative density function and its inverse a calculated. a and
+b are taken as the lower and the upper bound respectively.
+func need not to be normalized
+
+input:
+func    function  specifies the probability density function of the distribution
+a       float     specifies the lower bound of the distribution
+b       float     specifies the upper bound of the distribution
+
+members:
+pdf()          returns the probability density function, takes a float as an argument
+cdf()          returns the cumulative probability density function, takes a float as an argument
+ppf()          returns the inverse cumulative probability funciton, takes a float as an argument
+norm()         returns the normalisation factor (integral over pdf)
+sample(int n)  returns an array of n sampled values
+mean()         returns the mean value of the distribution
+_lowerbound    lowerbound of the distribution
+_upperbound    upperbound of the distribution
 '''
     def __init__(self, func = lambda x: 1, a = 0., b = 1000.):
         self._pdf = func

execute.py

@@ -11,7 +11,9 @@
 
 
 def main(quiet = False):
-    ''' This script produces a grid of expected numbers of stars according to the selection 
+    '''execute.main(quiet = False)
+
+This script produces a grid of expected numbers of stars according to the selection 
     criteria of Yusef-Zedah et al. 2009, 702,178-225 The Astrophysical Journal.
     The grid is in av for visual extinction, apera for aperature size and age for the maxage 
     of the starformation size
@@ -25,7 +27,7 @@ def main(quiet = False):
     else:
         output_stream = sys.stdout
 
-    sfr = .08
+    sfr = .01
     # star mass function
     def f(x):               # Kouper IMF
     #http://adsabs.harvard.edu/abs/2001MNRAS.322..231K
@@ -61,16 +63,17 @@ def g(x):
     for i in range(len(avs)):
         for j in range(len(aperas)):
             for k in range(len(ages)):
-                starformation.main(massfunction = mf, starformationhistory = sf[k], A_v = avs[i], sfr = .01, \
-                    apera = aperas[j], maxage = ages[k], appendix = "%s_%03d_%06d_%07d" % ('sim',avs[i],aperas[j],ages[k]), quiet=True)
+                starformation.main(massfunction = mf, starformationhistory = sf[k], \
+                    A_v = avs[i], sfr = sfr, apera = aperas[j], maxage = ages[k], \
+                    appendix = "%s_%03d_%06d_%07d" % ('sim',avs[i],aperas[j],ages[k]), quiet=True)
                 print(avs[i],aperas[j],ages[k], l/len(avs)/len(aperas)/len(ages), file=output_stream)
                 l = l+1
                 parameters.append([avs[i],aperas[j],ages[k]])
 
     t2 = time()                 # end of simulation
     print(t2, t1, t2-t1)
 
-    print ('number of simulations run: %s' %k , file=output_stream)  
+    print ('number of simulations run: %s' %l , file=output_stream)  
     head = ['AV', 'Aperature_size', 'Age']
     f = open('out/__head', 'w')
     f.write( ','.join(head)+'\n' )
@@ -79,7 +82,7 @@ def g(x):
 
     t3 = time()                 # end of saving data
 
-    analysis.main('out', quiet=True)
+    analysis.main('out')
     print ('analysis complete' , file=output_stream)  
 
     t4 = time()                 # end of analysing data

plot.py

@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-#test commentary
+from __future__ import print_function, division
 from mpl_toolkits.mplot3d import Axes3D
 from matplotlib import cm
 from matplotlib.mlab import griddata
@@ -11,7 +11,7 @@
 
 def main():
     aperature = 1
-    visual_ex = 1
+    visual_ex = 0
 
 
     f = open('out/__expected_number', 'r')
@@ -23,14 +23,13 @@ def main():
     fig = plt.figure()
     ax = fig.gca(projection='3d')
     avs = np.linspace(10.0, 50.0, 5)
-    aperas = np.linspace(10000, 50000, 5)
+    aperas = np.logspace(2, 5, 4)
     ages = np.linspace(500000, 2000000, 11)
-    numbers = data[:,3].reshape(5, 5, 11)
-    numbers = np.roll(numbers,4,2)
+    numbers = np.clip(data[:,3].reshape(5, 4, 11),5.,100000.) #prevent inf if no stars are detected
 
 
     X, Y = np.meshgrid(avs, ages)
-    Z = 559.*.01/np.transpose(numbers[:,aperature,:])
+    Z = 559.*.08/np.transpose(numbers[:,aperature,:])
     surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm)
     fig.colorbar(surf)
 
@@ -44,30 +43,34 @@ def main():
     ax.set_xlabel('av')
     ax.set_ylabel('age')
     ax.set_zlabel('starformation rate in M_sun/year')
+    #ax.set_zlim(0.,1.)
+    ax.set_title('Starformation as a function of\n visual extinction Av and age at apperaturesize=1kpc')
 
-    plt.savefig('plot/av-age.svg')
+    #plt.savefig('plot/av-age.svg')
 
     fig2 = plt.figure()
     ax2 = fig2.gca(projection='3d')
 
 
-    X2, Y2 = np.meshgrid(ages, aperas)
-    Z2 = 559.*.01/numbers[visual_ex,:,:]
+    X2, Y2 = np.meshgrid(ages, np.log10(aperas))
+    Z2 = 559.*.08/numbers[visual_ex,:,:]
     surf2 = ax2.plot_surface(X2, Y2, Z2, rstride=1, cstride=1, cmap=cm.coolwarm)
     fig2.colorbar(surf2)
 
 
-    x2 = data[:,2].reshape(5, 5, 11)[visual_ex,:,:].reshape(55)
-    y2 = data[:,1].reshape(5, 5, 11)[visual_ex,:,:].reshape(55)
-    z2 = 559.*.01/data[:,3].reshape(5, 5, 11)[visual_ex,:,:].reshape(55)
+    x2 = data[:,2].reshape(5, 4, 11)[visual_ex,:,:].reshape(44)
+    y2 = np.log10(data[:,1].reshape(5, 4, 11)[visual_ex,:,:].reshape(44))
+    z2 = 559.*.08/data[:,3].reshape(5, 4, 11)[visual_ex,:,:].reshape(44)
     ax2.scatter(x2,y2,z2)
     ax2.set_xlabel('age')
-    ax2.set_ylabel('apera')
+    ax2.set_ylabel('log(apera)')
     ax2.set_zlabel('starformation rate in M_sun/year')
+   # ax2.set_zlim(0.,1.)
+    ax2.set_title('Starformation as a function of\n visual extinction apperaturesize and age at Av=10')
 
 
 
-    plt.savefig('plot/age-apera.svg')
+   # plt.savefig('plot/age-apera.png')
 
 
 
@@ -78,16 +81,16 @@ def cmd(folder, av, apera, age):
     xmin = -1.
     xmax = 8.
     ymin = 15.
-    ymax = 0.
+    ymax = -1.
     selectmax = 0.
     selectmin = 8.
 
     #f = open("%s/sim_%s_%s_%s" % (folder,av,apera,age), 'r')
     #headers = f.readline().strip().split(',')
     #data = np.loadtxt(f)
     #f.close()
-    
-    hdulist = fits.open('%s/%s' %(folder,'sim_%s_%s_%s'%(av,apera,age)))
+
+    hdulist = fits.open('%s/%s' %(folder,'sim_%03d_%06d_%07d'%(av,apera,age)))
 #            av = hdulist[1].header['age']
     data = hdulist[1].data
 

setup.py

@@ -9,8 +9,17 @@
 
 
 def main(quiet=False):
-    '''Pulling all necessary files for using the starformation script
+    '''setup.main(quiet=False)
 
+Pulling all necessary files for using the starformation script
+
+This script pulls the fits-files for the radiation models from the MPIA-folder of Thomas
+Robitaille, if the files have been moved feel free to contact [email protected]
+
+It also pulls the extinction law from 
+http://caravan.astro.wisc.edu/protostars/info.php?topic=sedfitter_results 
+if you have problems with them please contact
+http://caravan.astro.wisc.edu/protostars/contact.php
 '''
     if quiet:
         output_stream = StringIO()

starformation.py

@@ -11,23 +11,40 @@
 from astropy.io import fits
 
 
-def main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, apera = 24000, maxage = 2000000., distance = 8.0, appendix='default', quiet=0, precise=0):
-    '''Creates a sample of stars
+def main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, apera = 24000,\
+ maxage = 2000000., distance = 8.0, appendix='default', quiet=0, precise=0):
+    '''main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, apera = 24000,\
+ maxage = 2000000., distance = 8.0, appendix='default', quiet=0, precise=0)
+
+Creates a sample of stars
 
 input:
-A_v     float   value for the visual extinction 
-sfr     float   Star formation rate in M_sun/year, is assumed to be constant over maxage
-apera   float   used aperature size for selecting the fluxes of the protostars
-maxage  float   age of the star formation site, sfr is assumed to be constant
-appendix String sets the outputfilename, default is the starting time (via time.time())
-quiet   boolean if true (=1) surpresses all output
-precise boolean if true (=1) sample single star till expected mass reached, else calculate
-                expected number (viaand sample as an array
+massfunction            distribution
+starformation history   distribution
+A_v       float     value for the visual extinction 
+sfr       float     Star formation rate in M_sun/year, is assumed to be constant over maxage
+apera     float     used aperature size for selecting the fluxes of the protostars
+maxage    float     age of the star formation site, sfr is assumed to be constant
+distance  float     distance to the simulated starformation site
+appendix  String    sets the outputfilename, default is the starting time (via time.time())
+quiet     boolean   if true (=1) surpresses all standard output
+precise   boolean   if true (=1) sample single star till expected mass reached, else calculate
+                    expected number and sample as an array
+
+The distributions must provide an object which has the following members:
+    float    cdf(float x)   returns the integrated distribution up to x, is used to calculate
+                            the expected mass
+    float    _upperbound    returns the upper limit of the distribution, is used to calculate
+                            the expected mass
+    float[]  sample(int n)  returns an array of n floats, sampled from the distribution
+    float    mean()         returns the mean value of the distribution
+
 
 output:
-returns two files in the folder 'out/' the _settings file contains the used values of 
-          A_v, sfr, apera, maxage, number of sampled stars, their cumulated mass and
-          the expected mass
+returns a fits file in the out-folder, either using the appendix as filename or the time of the
+      starting of the script in order to prevent overwriting existing files
+      In the header of the fits-file are the values: A_v, sfr, apera, maxage and distance recorded
+      In the data part are the age, mass, modelnumber and the uncorrected and corrected fluxes
 '''
 
     if quiet: