real commit

analysis.py

@@ -7,6 +7,7 @@
 from StringIO import StringIO
 import sys
 from astropy.io import fits
+from astropy.table import Table, Column
 
 
 
@@ -19,7 +20,7 @@ def main(folder, quiet=0):
 Parameters
 ----------
 folder  String   Specifiecs the folder where the files are
-quiet   boolean  =1 surpresses all standard output
+quiet   boolean  =1 suppresses all standard output
 
 Returns
 -------
@@ -58,8 +59,19 @@ def main(folder, quiet=0):
             x = -2.5*(np.log10(data['c%s' % color1]/64130) - np.log10(data['c%s' % color2]/7140))
             y = -2.5*(np.log10(data['c%s' % color2]/7140))
 
-            n = sum(np.logical_and( (y > -10./3. * (x-1.) + 10.), np.logical_and(max_mag < y, y < min_mag)))
+
+            sel = np.logical_and( (y > -10./3. * (x-1.) + 10.), np.logical_and(max_mag < y, y < min_mag))
+            n = sum(sel)
+            t = Table(hdulist[1].data)
+            if 'sel' in t.columns:
+                t.remove_column('sel')
+            t.add_column(Column(name='sel', data=sel.astype('int')))
+
+            hdulist[1].data = np.array(t)
             tmp, av, apera, age = fil.split('_')
+            #hdulist = fits.open('%s/%s' %(folder,fil), 'write')
+            #hdulist.writeto('%s/%s' %(folder,fil), clobber=True)
+            fits.update('%s/%s' %(folder,fil), np.array(t), ext = 1, clobber=True)
             out.append([Decimal(av), Decimal(apera), Decimal(age), n])
 
     #writing obtained data to "folder/__expected_number"

execute.py

@@ -20,7 +20,7 @@ def main(quiet = False):
 
 Parameters
 ----------
-quiet   boolean  if true surpresses all standard output
+quiet   boolean  if true suppresses all standard output
 
 
 Returns

plot.py

@@ -232,3 +232,31 @@ def init():
 
     return avs, aperas, ages, sfr, headers, data
 
+
+def histogram(folder, av, apera, age, color1 = "I4", color2 = "M1"):
+
+    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))
+    dat = data['mass'][data['sel'] == 1]
+    pdf, bins, patches = ax.hist(dat, range=(-.5,1.5))
+    #ax.set_yscale('log')
+
+
+
+
+
+
+
+
+
+
+
+

setup.py

@@ -23,7 +23,7 @@ def main(quiet=False):
 
 Parameter
 ---------
-quiet  boolean  if true surpresses all standard output
+quiet  boolean  if true suppresses all standard output
 '''
     if quiet:
         output_stream = StringIO()

starformation.py

@@ -28,7 +28,7 @@ def main(massfunction = 0, starformationhistory = 0, A_v = 10.0, sfr = .01, aper
 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
+quiet     boolean   if true (=1) suppresses all standard output
 precise   boolean   if true (=1) sample single star till expected mass reached, else calculate
                     expected number and sample as an array
 
@@ -178,18 +178,6 @@ def g(x):
     newfluxes = np.asarray(newfluxes)
     output = np.vstack([np.asarray(output).transpose(), matches, fluxes, newfluxes]).transpose()
 
-    # creating the output file
-    #head = ['#', 'age', 'mass', 'model']
-    #for mod in models:
-        #head.append('flux %s' % mod)
-    #for mod in models:
-        #head.append('corrected_flux %s' % mod)
-    #f = open('out/%s' % appendix, 'w')
-    #f.write( ','.join(head)+'\n' )
-    ##np.savetxt(f, output)
-    #np.savetxt(f, output)
-    #f.close()
-
     # create table
     # data table
     t = Table()
@@ -200,6 +188,7 @@ def g(x):
         t.add_column(Column(name='%s' % models[i], data=output[:,4+i]))
     for i in range(len(models)):
         t.add_column(Column(name='c%s' % models[i], data=output[:,4+len(models)+i]))
+   # t.add_column(Column(name='sel', data=zeros(len(models[0]))))
     # head table
     header = Table()
     header.add_column(Column(name='AV', data = [A_v]))