1 """
2 Compute Least-Square-Deconvolutions of spectra.
3
4 Section 1. Single stars
5 =======================
6
7 Section 1.1 Without Tikhonov regularization
8 -------------------------------------------
9
10 Generate some spectra of a single star:
11
12 >>> velos,V,S,masks = __generate_test_spectra(1,binary=False,noise=0.01)
13
14 Compute the LSD profiles in certain radial velocity range:
15
16 >>> rvs = np.linspace(-10,10,100)
17 >>> Z,cc = lsd(velos,V,S,rvs,masks)
18
19 Plot both the input spectrum and the LSD and CCF profiles:
20
21 >>> p = pl.figure()
22 >>> p = pl.subplot(121)
23 >>> p = pl.plot(velos,V[:,0],'k-',label='Observation')
24 >>> p = pl.vlines(masks[0][0],1,1-np.array(masks[0][1]),color='r',label='Mask')
25 >>> p = pl.legend(loc='best')
26 >>> p = pl.subplot(122)
27 >>> p = pl.plot(rvs,Z[0][0],'k-',label='LSD')
28 >>> p = pl.plot(rvs,cc[0][0],'r-',label='CCF')
29 >>> p = pl.legend(loc='best')
30
31 ]]include figure]]ivs_spectra_lsd01.png]
32
33 Section 1.2 With Tikhonov regularization
34 ----------------------------------------
35
36 Do the same as above, but for more noise and different Tikhonov regularizations.
37
38 >>> velos,V,S,masks = __generate_test_spectra(1,binary=False,noise=0.1)
39
40 >>> rvs = np.linspace(-10,10,100)
41 >>> lambdas = [0.,0.1,0.5,1.0]
42 >>> output = [lsd(velos,V,S,rvs,masks,Lambda=lam) for lam in lambdas]
43
44 We plot both the input spectrum and the LSD and CCF profiles:
45
46 >>> p = pl.figure()
47 >>> p = pl.subplot(121)
48 >>> p = pl.plot(velos,V[:,0],'k-',label='Observation')
49 >>> p = pl.vlines(masks[0][0],1,1-np.array(masks[0][1]),color='r',label='Mask')
50 >>> p = pl.legend(loc='best')
51 >>> p = pl.subplot(122)
52 >>> for lam,(Z,cc) in zip(lambdas,output):
53 ... p = pl.plot(rvs,Z[0][0],'-',label='$\Lambda$=%.1f'%(lam),lw=2)
54 >>> p = pl.legend(loc='best')
55
56 ]]include figure]]ivs_spectra_lsd02.png]
57
58 Section 2. Binary stars
59 =======================
60
61 Generate some spectra of a binary star:
62
63 >>> velos,V,S,masks = __generate_test_spectra(1,binary=True,noise=0.01)
64
65 Compute the LSD profiles in certain radial velocity range, first using only
66 one line mask, then both:
67
68 >>> rvs = np.linspace(-10,10,100)
69 >>> Z1,cc1 = lsd(velos,V,S,rvs,masks[:1])
70 >>> Z2,cc2 = lsd(velos,V,S,rvs,masks)
71
72 Plot both the spectrum and the LSD and CCF profiles. Note that the CCF in the
73 binary case is exactly the same as in the single star case (see implementation).
74 First, we plot the LSD profile under the assumption of a single star. Second,
75 we plot the LSD profile when taking binarity into account.
76
77 >>> p = pl.figure()
78 >>> p = pl.subplot(121)
79 >>> p = pl.plot(velos,V[:,0],'k-',label='Observation')
80 >>> p = pl.vlines(masks[0][0],1,1-np.array(masks[0][1]),color='r',label='Mask 1',lw=2)
81 >>> p = pl.legend(loc='lower right')
82 >>> p = pl.subplot(122)
83 >>> p = pl.plot(rvs,Z1[0][0],'k-',label='LSD 1')
84 >>> p = pl.plot(rvs,cc1[0][0],'r-',label='CCF 1')
85 >>> p = pl.legend(loc='best')
86
87 ]]include figure]]ivs_spectra_lsd03.png]
88
89 >>> p = pl.figure()
90 >>> p = pl.subplot(121)
91 >>> p = pl.plot(velos,V[:,0],'k-',label='Observation')
92 >>> p = pl.vlines(masks[0][0],1,1-np.array(masks[0][1]),color='r',label='Mask 1',lw=2)
93 >>> p = pl.vlines(masks[1][0],1,1-np.array(masks[1][1]),color='b',label='Mask 2',lw=2)
94 >>> p = pl.legend(loc='lower right')
95 >>> p = pl.subplot(122)
96 >>> p = pl.plot(rvs,Z2[0][0],'r-',label='LSD 1',lw=2)
97 >>> p = pl.plot(rvs,Z2[0][1],'b-',label='LSD 2',lw=2)
98 >>> p = pl.legend(loc='best')
99
100 ]]include figure]]ivs_spectra_lsd04.png]
101
102 """
103 import pylab as pl
104 import numpy as np
105 import numpy.linalg as la
106 from ivs.sigproc import evaluate
107 import itertools
108
109 -def lsd(velos,V,S,rvs,masks,Lambda=0.):
110 """
111 Compute LSD profiles and cross correlation functions.
112
113 Possibility to include Tikhonov regularization to clean up the profiles,
114 when setting C{Lambda>0}.
115
116 Possibility to include multiprofile LSD. Parameter C{masks} should be a list
117 of C{(centers,weights)} (if you give only one mask, give
118 C{masks=[(centers,weights)]}.
119
120 See Donati, 1997 for the original paper and Kochukhov, 2010 for extensions.
121
122 @parameter velos: velocity vector of observations
123 @type velos: array of length N_spec
124 @parameter V: observation array
125 @type V: N_obs x N_spec array
126 @parameter S: weights of individual pixels
127 @type S: array of length N_spec
128 @parameter rvs: radial velocity vector to compute the profile on
129 @type rvs: array of length N_rv
130 @parameter Lambda: Tikhonov regularization parameter
131 @parameter masks: list of tuples (center velocities, weights)
132 @type masks: list (length N_mask) of tuples of 1D arrays
133 @type Lambda: float
134 @return: LSD profile, CCF of shape (N_obs x (N_rv.N_mask))
135 @rtype: 2D array, 2D array
136 """
137
138 m,n = len(rvs),len(velos)
139 Nspec = V.shape[1]
140 Nmask = len(masks)
141 V = np.matrix(V)-1
142
143
144 S = np.matrix(np.diag(S))
145
146 M = np.matrix(np.zeros((n,m*len(masks))))
147 for N,(line_centers,weights) in enumerate(masks):
148 for l,lc in enumerate(line_centers):
149 for i in range(n):
150 for j in range(m-1):
151 vi = velos[i]-lc
152 if not (rvs[j]<vi<rvs[j+1]): continue
153 M[i,j+ N*m] = weights[l]*(rvs[j+1]-vi)/(rvs[j+1]-rvs[j])
154 M[i,j+1+N*m] = weights[l]*(vi -rvs[j])/(rvs[j+1]-rvs[j])
155 M = np.matrix(M)
156
157 if Lambda:
158 R = np.matrix(np.zeros((m*Nmask,m*Nmask)))
159 for i in range(1,m-1):
160 R[i,i] = 2
161 R[i-1,i] = -1
162 R[i+1,i] = -1
163 R[0,0] = 1
164 R[1,0] = -1
165 R[-1,-1] = 1
166 R[-2,-1] = -1
167
168 X = M.T*(S**2)
169 XM = X*M
170 if Lambda:
171 XM = XM+Lambda*R
172 cc = X*V
173
174
175
176
177 Z,res,rank,s = la.lstsq(XM,cc)
178
179 Z = np.array(Z.T)
180 cc = np.array(cc.T)
181
182 Z_ = []
183 C_ = []
184 for i in range(len(Z)):
185 Z_.append([])
186 C_.append([])
187 for N in range(Nmask):
188 Z_[-1].append(Z[i][N*m:(N+1)*m])
189 C_[-1].append(cc[i][N*m:(N+1)*m])
190
191 return Z_,C_
192
194 spec_length = 1000
195 velo_length = 100
196 obs_velo1 = [-3.14,-3.5,-4,3.25][:Nspec]
197 obs_velo2 = [+3.14,+3.5,+4,-3.25][:Nspec]
198 m,n = velo_length,spec_length
199
200 velos = np.linspace(-30,30,spec_length)
201 line_centers1 = [-11,3,9]
202 weights1 = [0.5,0.1,0.3]
203 line_centers2 = [-5,2,10]
204 weights2 = [0.3,0.4,0.1]
205
206
207 S = np.ones(spec_length)
208
209 V = [np.zeros(spec_length) for i in range(Nspec)]
210 masks = [(line_centers1,weights1)]
211 for i,prof_velo in enumerate(obs_velo1):
212 for line_center,weight in zip(line_centers1,weights1):
213 V[i] += evaluate.gauss(velos,[weight,line_center-prof_velo,1.])
214 V[i] += np.random.normal(size=n,scale=noise)
215 if binary:
216 masks.append((line_centers2,weights2))
217 for i,prof_velo in enumerate(obs_velo2):
218 for line_center,weight in zip(line_centers2,weights2):
219 V[i] += evaluate.gauss(velos,[weight,line_center-prof_velo,1.])
220 V = 1-np.array(V)
221 V = np.matrix(V).T
222
223 return velos,V,S,masks
224
225 if __name__=="__main__":
226 import doctest
227 doctest.testmod()
228 pl.show()
229