diff --git a/analysis.py b/analysis.py index 89c07cd..71f2d5c 100755 --- a/analysis.py +++ b/analysis.py @@ -6,6 +6,7 @@ from decimal import Decimal from StringIO import StringIO import sys +from astropy.io import fits @@ -42,8 +43,8 @@ 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 @@ -51,8 +52,8 @@ def main(folder, quiet=0): #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? diff --git a/execute.py b/execute.py index 4d99b91..e1320f0 100644 --- a/execute.py +++ b/execute.py @@ -19,6 +19,7 @@ def main(quiet = False): ''' t0 = time() #timing possibility + print(t0) if quiet: output_stream = StringIO() else: @@ -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'] diff --git a/plot.py b/plot.py index 6f930f1..5a10ab8 100644 --- a/plot.py +++ b/plot.py @@ -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 @@ -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') @@ -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) @@ -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)) diff --git a/setup.py b/setup.py index 8e6f7f9..253c865 100644 --- a/setup.py +++ b/setup.py @@ -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', @@ -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() diff --git a/starformation.py b/starformation.py index 3f30d3c..81de9c7 100755 --- a/starformation.py +++ b/starformation.py @@ -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'] @@ -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) @@ -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) \ No newline at end of file +#main(sfr = .08) \ No newline at end of file