-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathstarformation.py
executable file
·196 lines (158 loc) · 7.79 KB
/
starformation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
# -*- coding: utf-8 -*-
from __future__ import print_function
from time import time
import distribution as dist
import numpy as np
import matplotlib.pyplot as plt
import scipy.spatial
from astropy.io import fits
from StringIO import StringIO
import sys
def main(A_v = 10.0, sfr = .001, apera = 24000, maxage = 2000000., 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
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
'''
if quiet:
output_stream = StringIO()
else:
output_stream = sys.stdout
t0 = time()
if appendix=='default': # making sure not to overwrite former output
appendix=t0 # by using the starting time as an unique id
#parameter settings
k_v = 211.4 # opacity in v_band in cm^2/g
# wavelength of the corresponding filterband in microns
wavelength = [1.235, 1.662, 2.159, 3.550, 4.493, 5.731, 7.872, 23.68, 71.42, 155.9]
models = ['2H', '2J', '2K', 'I1', 'I2', 'I3', 'I4', 'M1', 'M2', 'M3']
# star mass function
def f(x): # Kouper IMF
#http://adsabs.harvard.edu/abs/2001MNRAS.322..231K
if x<.01:
return 0
elif x < .08:
return 62.46192*x**-.3
elif x < .5:
return 1.413352*x**-1.8
elif x < 1.:
return x**-2.3 # also value 2.7 given eq.(6) or eq (2)
else:
return x**-2.3
f = np.vectorize(f)
mf = dist.distribution(f, .1, 50.)
# star formation history
def g(x):
return sfr
g = np.vectorize(g)
sf = dist.distribution(g, 10000., maxage)
# todo - change sample an expected number of stars (expmass/averagemass)
cumass = 0. #sampled mass
exmass = sf.cdf()(sf._upperbound)-sf.cdf()(sf._lowerbound) #expected mass formed
stars = [] #storing for the sample
n = 0
t1 = time() # startup completed
if precise == True:
while cumass < exmass:
mass, age = mf.sample(1)[0], sf.sample(1)[0]
cumass = cumass + mass
stars.append([n, age, mass])
if n % 10000 == 0:
print (n, cumass, file=output_stream) #reporting progress
n = n+1
else:
n = int(exmass/.49)
mass, age = mf.sample(n), sf.sample(n)
cumass = np.sum(mass)
stars = [[i, age[i], mass[i]] for i in range(n)]
print ('number of sampled stars: %s' %n , file=output_stream)
print ('mass of sampled stars: %s' % cumass , file=output_stream)
print ('mean mass: %s' % (cumass/n), file=output_stream)
print ('expected mass of stars: %s' % exmass , file=output_stream)
t2 = time() # sampleing completed
# python code for model contact
#initial parameters
model = [ fits.open('models/%s.fits' % mod) for mod in models ] # fits-data for the model
param = fits.open('models/parameters.fits.gz') # modelparameter
app_num = [ np.interp(apera, model[i][2].data.field(0), range(model[i][2].data.field(0).size)) for i in range(len(models)) ]
# to do:
# check for interpolation of aperature size
# sampling viewing angle
angle = [np.random.random_integers(9) for i in range(len(stars))]
#reading model grid
grid = [param[1].data[10*i][1:3] for i in range(param[1].data.size /10) ] #model data 0: name, 1: age, 2: mass ; only one per model
#converting to logspace
stars = np.asarray(stars)
grid = np.log10(grid)
stars[:,1:] = np.log10(stars[:,1:])
output = stars.tolist() #creating output
#normalizing for nearest neighbor search
grid[0,:] = grid[0,:]/(grid[0,:].max() - grid[0,:].min())
grid[1,:] = grid[1,:]/(grid[1,:].max() - grid[1,:].min())
stars[1,:] = stars[1,:]/(grid[0,:].max() - grid[0,:].min())
stars[2,:] = stars[2,:]/(grid[1,:].max() - grid[1,:].min())
t3 = time() #model data load complete
tree = scipy.spatial.cKDTree(grid,leafsize=10) #search tree
matches = [tree.query(star[1:] , k=1)[1] for star in stars] #saves matches with (dist, index)
t4 = time() #matching sample to data complete
# extracting fluxes
fluxes = [0 for j in range(len(models)) ]
for j in range(len(models)):
# fluxes[j] = [ [model[j][1].data[10*matches[i] + angle[i]][1][app_num[j]], model[j][1].data[10*matches[i] + angle[i]][2][app_num[j]] ] for i in range(len(matches)) ]
fluxes[j] = [ model[j][1].data[10*matches[i] + angle[i]][1][app_num[j]] for i in range(len(matches)) ]
# applying extinction
extinction = np.loadtxt('models/extinction_law.ascii')
k_lambda = np.interp(wavelength, extinction[:,0], extinction[:,1])
correctionfactor = 10.**(-.4 * A_v * k_lambda / k_v)
newfluxes = [0 for j in range(len(models)) ]
for j in range(len(models)):
newfluxes[j] = np.asarray(fluxes[j]) * correctionfactor[j]
t5 = time() #extracting fluxes complete
# saving data to output: #, log10(age), log10(mass), modelmatch, (flux, flux_error, corrected_flux, corrected_flux_error) for each model
for i in range(len(output)):
output[i].append(matches[i])
for j in range(len(models)):
output[i].append(fluxes[j][i])
output[i].append(newfluxes[j][i])
# output[i].append(fluxes[j][i][0]) possible to use if fluxerrors
# output[i].append(fluxes[j][i][1]) are necessary
# output[i].append(newfluxes[j][i][0])
# output[i].append(newfluxes[j][i][1])
# creating the output file
head = ['#', 'age', 'mass', 'model']
for mod in models:
head.append('flux %s' % mod)
head.append('corrected_flux %s' % mod)
f = open('out/%s' % appendix, 'w')
f.write( ','.join(head)+'\n' )
np.savetxt(f, output)
f.close()
# creating the settings file
f = open('out/%s_settings' % appendix, 'w')
settings = '%s #visual extinction A_v \n' % A_v
settings = settings + '%s #star formation rate sfr \n' % sfr
settings = settings + '%s #star formation time time \n' % maxage
settings = settings + '%s #aperature size apera\n' % apera
settings = settings + '%s #number of sampled stars \n' % len(stars)
settings = settings + '%s #cumulated sampled mass \n' % cumass
settings = settings + '%s #expected mass \n' % exmass
f.write(settings)
f.close()
t6 = time() #saving complete
# timing possibility for optimization efforts
print( 'starting script at %f' %(t0), file=output_stream)
print( 'initializing %f' %(t1-t0), file=output_stream)
print( "sampleing %f" %(t2-t1), file=output_stream)
print( "model data load %f" %(t3-t2), file=output_stream)
print( "matching model %f" %(t4-t3), file=output_stream)
print( "extracting fluxes %f" %(t5-t4), file=output_stream)
print( "saving %f" %(t6-t5), file=output_stream)
print( "________________________", file=output_stream)
print( "total runtime %f" %(t6-t0), file=output_stream)
print( "finishing script %f" %t6, file=output_stream)