/
myutil.py~
122 lines (109 loc) · 3.68 KB
/
myutil.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
#!/usr/bin/python
from numpy import *
from pyraf import iraf
from iraf import stsdas,hst_calib,ttools,synphot
import my_pysynphot as S
import arrayspec as A
import os
pc_cm = 3.0857e18
def closeint(x):
if fmod(x,1)>=0.5:
return int(ceil(x))
else:
return int(floor(x))
def closepar(a,b): ### function to search the closest element to b in list a
a = sort(a)
i = searchsorted(a,b)
if i == len(a): ### b is larger than any element in a
b = a[-1]
elif i == 0: ### b is smaller than any element in a
b = a[0]
else:
if fabs(b-a[i-1])<=fabs(b-a[i]):
b = a[i-1]
else: b = a[i]
return b
def searchpar(Z,tau,rcy,age,ext): ### function that returns the closest parameter value
taumin = agemin = 70
taumax = agemax = 98
taustep = 2
rcymax = 35
rcystep = 5
Zlist = arange(2,8)
taulist = arange(taumin,taumax,taustep)
rcylist = arange(0,rcymax,rcystep)
extlist = arange(0,60,10)
Z = closepar(Zlist,Z)
rcy = closepar(rcylist,rcy)
tau = closepar(taulist,tau)
a = arange(agemin,tau-1,1)
b = arange(tau,agemax,2)
agelist = concatenate((a,b)) ### agelist be constructed after tau is selected!!
age = closepar(agelist,age)
ext = closepar(extlist,ext) ### discrete extinction
#if ext > 0.5: ext = 0.5 ### continuous extinction
#elif ext < 0.0: ext = 0.0
return (Z,tau,rcy,age,ext)
def convertfit(outfile,cdfile,infile):
iraf.tcreate(outfile,cdfile,infile)
def renorm_mag(spec,band,normval,magform):
# input spec as pysynphot spectrum object
# band as pysynphot bandpass object
# normval is the value of abmag that one wants to normalize to
# will return a renormalized spectrum object AND normalization factor
if magform == 'abmag':
M = A.ABmag
elif magform == 'vegamag':
M = A.Vegamag
else:
raise('unsupported magnitude system')
m2 = M(spec,band)
factor = 10**(-0.4*(normval - m2))
sp = spec*factor
return sp, factor
def renorm_absmag(spec,band,normval,magform,parsec=pc_cm):
# normalize the spectrum to a certain absolute magnitude
# placed at 10 pc
if magform == 'abmag':
M = A.ABmag
elif magform == 'vegamag':
M = A.Vegamag
else:
raise('unsupported magnitude system')
m2 = M(spec,band)
absnormmag = 2.5*log10(4*pi*(10*parsec)**2)
mag = m2 + absnormmag # mag of spec at 10pc
normmag = normval - mag
normfactor = 10**(-0.4*normmag)
sp = spec*normfactor
return normmag, normfactor
def trapezoidIntegration(x,y):
""" Input x and y are numpy arrays"""
npoints = x.size
indices = arange(npoints)[:-1]
deltas = x[indices+1] - x[indices]
integrand = 0.5*(y[indices+1] + y[indices])*deltas
return integrand.sum()
def chisquare(model_array, data_array, error_array):
'''Calculates chi-square given arrays of model output, data points
and uncertainties of data points. If multiple sources of error exist,
must add them (in quadrature if neccessary) beforehand'''
if len(data_array) != len(model_array):
raise IndexError, 'data and error length do not match'
temp = (model_array - data_array)/error_array
temp = temp**2
chisq = sum(temp)
return chisq
def ffindall(strarray, dir=os.getcwd()):
'''Finds and returns a list of file names in directory dir that
contains ALL the strings in str'''
flist = []
for f in os.listdir(dir):
n = 1
for str in strarray:
if str not in f:
n = 0
break
if n == 1:
flist.append(f)
return flist