Skip to content

Commit

Permalink
fixes
Browse files Browse the repository at this point in the history
fixed starformation output to correctly match columns and columnnames
fixed analysis to corectly work with fitsfiles instead of ascii
fixed plot to work with fits files instead of ascii
added extinction_law files to setup
changed execute to rerun skript
  • Loading branch information
healther committed Apr 8, 2013
1 parent 0c4f11c commit 41bf189
Show file tree
Hide file tree
Showing 5 changed files with 60 additions and 49 deletions.
9 changes: 5 additions & 4 deletions analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from decimal import Decimal
from StringIO import StringIO
import sys
from astropy.io import fits



Expand Down Expand Up @@ -42,17 +43,17 @@ def main(folder, quiet=0):


# Read in
hdulist = fits.open('file.fits')
av = hdulist[1].header['age']
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['corrected_flux %s' % color1]/64130) - np.log10(data[:,c2]/7140))
y = -2.5*(np.log10(data['corrected_flux %s' % color2]/7140))
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))


# efficiency? accuracy?
Expand Down
9 changes: 6 additions & 3 deletions execute.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ def main(quiet = False):
'''
t0 = time() #timing possibility
print(t0)
if quiet:
output_stream = StringIO()
else:
Expand Down Expand Up @@ -49,23 +50,25 @@ def g(x):
sf = [dist.distribution(g, 10000., ages[i]) for i in range(len(ages))]

t1 = time() # finished reading the distributions
print(t1)


# setting up model data
aperas = np.linspace(10000, 50000, 5)
aperas = np.logspace(2, 5, 4)
avs = np.linspace(10.0, 50.0, 5)
l = 0.
parameters = []
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 = .08, \
apera = aperas[j], maxage = ages[k], appendix = "%s_%s_%s_%s" % ('sim',avs[i],aperas[j],ages[k]), quiet=True)
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)
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)
head = ['AV', 'Aperature_size', 'Age']
Expand Down
59 changes: 34 additions & 25 deletions plot.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,11 @@
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits


def main():
aperature = 2
aperature = 1
visual_ex = 1


Expand All @@ -21,52 +22,52 @@ def main():

fig = plt.figure()
ax = fig.gca(projection='3d')
avs = np.linspace(5.0, 50.0, 10)
aperas = np.linspace(15000, 50000, 11)
avs = np.linspace(10.0, 50.0, 5)
aperas = np.linspace(10000, 50000, 5)
ages = np.linspace(500000, 2000000, 11)
numbers = data[:,3].reshape(10, 11, 11)
numbers = data[:,3].reshape(5, 5, 11)
numbers = np.roll(numbers,4,2)


X, Y = np.meshgrid(avs, ages)
Z = 559.*.08/np.transpose(numbers[:,aperature,:])
Z = 559.*.01/np.transpose(numbers[:,aperature,:])
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm)
fig.colorbar(surf)


x = X.reshape(110)
y = Y.reshape(110)
z = Z.reshape(110)
x = X.reshape(55)
y = Y.reshape(55)
z = Z.reshape(55)


ax.scatter(x,y,z)
ax.set_xlabel('av')
ax.set_ylabel('age')
ax.set_zlabel('starformation rate in 0.001 M_sun/year')
ax.set_zlabel('starformation rate in M_sun/year')

plt.savefig('av-age.svg')
plt.savefig('plot/av-age.svg')

fig2 = plt.figure()
ax2 = fig2.gca(projection='3d')


X2, Y2 = np.meshgrid(ages, aperas)
Z2 = numbers[visual_ex,:,:]/559.
Z2 = 559.*.01/numbers[visual_ex,:,:]
surf2 = ax2.plot_surface(X2, Y2, Z2, rstride=1, cstride=1, cmap=cm.coolwarm)
fig2.colorbar(surf2)


x2 = data[:,2].reshape(10, 11, 11)[visual_ex,:,:].reshape(121)
y2 = data[:,1].reshape(10, 11, 11)[visual_ex,:,:].reshape(121)
z2 = 559.*.08/data[:,3].reshape(10, 11, 11)[visual_ex,:,:].reshape(121)
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)
ax2.scatter(x2,y2,z2)
ax2.set_xlabel('age')
ax2.set_ylabel('apera')
ax2.set_zlabel('starformation rate in 0.001 M_sun/year')
ax2.set_zlabel('starformation rate in M_sun/year')



plt.savefig('age-apera.svg')
plt.savefig('plot/age-apera.svg')



Expand All @@ -81,16 +82,24 @@ def cmd(folder, av, apera, age):
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()
#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)))
# av = hdulist[1].header['age']
data = hdulist[1].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))


c1 = headers.index('corrected_flux %s' % color1)
c2 = headers.index('corrected_flux %s' % color2)
#c1 = headers.index('corrected_flux %s' % color1)
#c2 = headers.index('corrected_flux %s' % color2)

x = -2.5*(np.log10(data[:,c1]/64130) - np.log10(data[:,c2]/7140))
y = -2.5*(np.log10(data[:,c2]/7140))
#x = -2.5*(np.log10(data[:,c1]/64130) - np.log10(data[:,c2]/7140))
#y = -2.5*(np.log10(data[:,c2]/7140))
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x,y)
Expand All @@ -105,4 +114,4 @@ def cmd(folder, av, apera, age):
ax.set_ylim(ymin, ymax)
ax.set_title('%s_%s_%s.svg' %(av,apera,age))

plt.savefig('%s_%s_%s.svg' %(av,apera,age))
plt.savefig('plot/%s_%s_%s.svg' %(av,apera,age))
10 changes: 5 additions & 5 deletions setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ def main(quiet=False):
'http://www.mpia-hd.mpg.de/~robitaille/share/andreas/M1.fits',
'http://www.mpia-hd.mpg.de/~robitaille/share/andreas/M2.fits',
'http://www.mpia-hd.mpg.de/~robitaille/share/andreas/M3.fits',
''
'http://caravan.astro.wisc.edu/protostars/files/extinction_law.tar.gz'
]
file_names = [
'models/parameters.fits.gz',
Expand All @@ -58,8 +58,8 @@ def main(quiet=False):
f.close()
print('Downloaded %s from %s' % (file_names[i],urls[i]), file=output_stream)


f = tarfile.open(path, 'r:gz')
try: f.extractall()
finally: f.close()
if not os.path.isfile('modesl/extinction_law.ascii')
f = tarfile.open('models/extinction_law.tar.gz', 'r:gz')
try: f.extractall()
finally: f.close()

22 changes: 10 additions & 12 deletions starformation.py
Original file line number Diff line number Diff line change
Expand Up @@ -157,8 +157,7 @@ def g(x):
# saving data to output: #, log10(age), log10(mass), modelmatch, (flux, flux_error, corrected_flux, corrected_flux_error) for each model
fluxes = np.asarray(fluxes)
newfluxes = np.asarray(newfluxes)
output = np.vstack([np.asarray(output).transpose(), fluxes, newfluxes]).transpose()

output = np.vstack([np.asarray(output).transpose(), matches, fluxes, newfluxes]).transpose()

# creating the output file
#head = ['#', 'age', 'mass', 'model']
Expand All @@ -178,18 +177,17 @@ def g(x):
t.add_column(Column(name='age', data=output[:,1]))
t.add_column(Column(name='mass', data=output[:,2]))
t.add_column(Column(name='model', data=output[:,3]))
#for i in range(len(models)):
#t.add_column(Column(name='flux %s' % models[i], data=output[:,4+i]))
#tmp = 4 + len(models)
for i in range(len(models)):
t.add_column(Column(name='corrected_flux %s' % models[i], data=output[:,4+i]))
t.add_column(Column(name='flux %s' % models[i], data=output[:,4+i]))
for i in range(len(models)):
t.add_column(Column(name='cflux %s' % models[i], data=output[:,4+len(models)+i]))

header = fits.Header()
header['Av'] = A_v
header['star formation rate'] = sfr
header['apperature size'] = apera
header['maximum age of the stars'] = maxage
header['distance of the formation size'] = distance
header['AV'] = A_v
header['SFR'] = sfr
header['APPERA'] = apera
header['MAXAGE'] = maxage
header['DIST'] = distance

fits.writeto('out/%s' % appendix, np.array(t), header, clobber=True)

Expand Down Expand Up @@ -221,4 +219,4 @@ def g(x):
print( "total runtime %f" %(t6-t0), file=output_stream)
print( "finishing script %f" %t6, file=output_stream)

main(sfr = .08)
#main(sfr = .08)

0 comments on commit 41bf189

Please sign in to comment.