1 """
2 Functions for Red Giant Seismology.
3
4 Contents:
5 =========
6
7 This library contains several functions which allow to identify modes
8 in the power-spectra of red giants. For the description of the pressure modes,
9 the formulism of Mosser et al. (2011A&A...525L...9M) for the 'Universal
10 Oscillation Pattern' is available.
11
12 For the description of the frequencies and rotational splitting of the mixed
13 dipole modes, the formalism of Mosser et al. (2012A&A...540A.143M) and (2012arXiv1209.3336M)
14 has been scripted. To apply this method, one needs to know the value of the true
15 period spacing. The modulation of the mixing between pressure- and gravity-dipole
16 modes is described through a Lorentzian profile. This approach assumes a constant
17 Period Spacing for all radial orders. Comparing the modelled frequencies with the
18 real frequencies, found in an spectrum, will work well for the several radial
19 orders but deviations can be found in other regions of the spectrum.
20 A similar Lorentzian description is valid for rotational splitting.
21
22
23 Contents:
24 ---------
25 * purePressureMode: calculates the position of the pure pressure pressure
26 mode of a spherical degree \ell and radial order.
27
28 * universalPattern: creates an array with the frequencies of all pure pressure
29 modes with \ell=0,1,2 and 3, whereby the index i selectes all m
30 odes of the same spherical degree:
31 e.g.: asymptoticRelation[0] = all \ell=0, asymptoticRelation[2] = all \ell=2,
32 asymptoticRelation[-1] gives an array with the radial orders (as float).
33
34 * asymptoticRelation: returns an array with the frequencies of the pure
35 pressure modes of the pure pressure modes
36
37 * rotational splitting: calculates the expected rotational splitting
38 of a dipole mode based on the degree of mixing
39
40
41 Example: KIC 4448777
42 --------------------
43 This star is a H-shell burning red giant with rotational splitting.
44 The powerspectrum can be found here: ...
45
46 >>> import ivs.asteroseismology.redgiantfreqs as rg
47 >>> centralRadial= 226.868 # [muHz], Frequency of the central radial mode
48 >>> largeSep = 17.000 # [muHz], large Frequency separation in powerspectrum
49 >>> smallSep = 2.201 # [muHz], small Frequency separation in powerspectrum
50 >>>
51 >>> DeltaPg = 89.9 # [seconds], true period spacing of dipole modes
52 >>> dfMax = 0.45 # [muHz], normalized rotational splitting (i.e. |freq_{m=0} - freq_{m=+/-1}|
53 >>>
54 >>> qCoupling= 0.15; lambdaParam = 0.5; beParam = 0.08 #Default values
55
56 First model the frequencies of the pure pressure modes.
57
58 >>> frequenciesUniversalPattern = rg.universalPattern(largeSep, smallSep, centralRadial)
59
60 Next, we calculate the frequencies of the mixed dipole modes:
61
62 >>> mixedModesPattern = rg.asymptoticRelation(centralRadial,largeSep, DeltaPg, qCoupling)
63
64 >>> rotationalSplittingModulation = rg.symmetricRotationalModulation(dfMax, largeSep, frequenciesUniversalPattern, mixedModesPattern, lambdaParam, beParam)
65
66
67
68
69 """
70
71 import sys, math, os
72 import scipy as scipy
73 import pylab as pl
74 import numpy as np
75 import numpy.random as npr
76 import random as rd
77 from time import gmtime, strftime
78 import glob as glob
79 from scipy import stats
80 import logging
81
82
83
84 format = ''
85 datefmt = ''
86
87 logging.basicConfig(level=logging.INFO,format=format,datefmt=datefmt)
88 logger = logging.getLogger('IVS.RG')
89
90
91
93 """
94 Calculates and returns the frequency of a single pure pressure mode for a spherical (\ell=0,1,2, and 3)[Nota Bene: 3 to be implemented].
95
96 smallSep can be set to an abitrary value if not calculating \ell=2 modes.
97
98 @param nnn: radial order of the mode with respect to the radial order of the central radial mode.
99 @type nnn: float or integer
100 @param largeSep: large Frequency separation in powerspectrum, in muHz
101 @type largeSep: float
102 @param smallSep: small Frequency separation in powerspectrum, in muHz
103 @type smallSep: float
104 @param centralRadial: Frequency of the central radial mode, in muHz
105 @type centralRadial: float
106 @param epsilon: echelle phase shift of the central Radial mode: epsilon = centralRadial / largeSep - centralRadial // largeSep
107 @type epsilon: float
108 @param ell: spherical degree of a mode
109 @type ell: float or integer
110 @return frequency in muHz
111 @rtype float array
112 """
113
114 ell = np.float(ell)
115 nnn = np.float(nnn)
116 smallSep /= largeSep
117 dominantRadialOrderPressure = np.floor(centralRadial / largeSep)
118 alphaCorr = 0.015 * largeSep**(-0.32)
119
120 if ell == 0.: constantDEll = 0.0
121 if ell == 1.: constantDEll = -(ell/2.+0.025)
122 if ell == 2.: constantDEll = smallSep
123 if ell == 3.: constantDEll = 0.0
124
125 purePressureModeEll = largeSep * ( (dominantRadialOrderPressure + nnn) + epsilon - constantDEll + alphaCorr/2.*(nnn)**2.)
126
127 return purePressureModeEll
128
129
130
131
132
133
134 -def universalPattern(largeSep, smallSep, centralRadial, numberOfRadialOrders = 3):
135
136 """
137 Calculates and returns an array containing frequencies the frequencies of the pure pressure modes
138 of the degrees spherical degrees \ell= 0,1,2 and 3 for a defined range of radial orders.
139
140 When slicing, the index number is equal to the spherical degree \ell, universalPattern[0] = all \ell=0, universalPattern[2] = all \ell=2,
141 The number of the radial order is given as output[-1]. If desired, the array can also be printed to the screen.
142
143 @param largeSep: large Frequency separation in powerspectrum, in muHz
144 @type largeSep: float
145 @param smallSep: small Frequency separation in powerspectrum, in muHz
146 @type smallSep: float
147 @param centralRadial: Frequency of the central radial mode, in muHz
148 @type centralRadial: float
149 @param numberOfRadialOrders: for how many radial orders above and below the Central Radial Mode should the frequencies be calculated.
150 @type numberOfRadialOrders: float or integer
151
152 @return array with frequencies in muHz, [-1] indicates the radial order.
153 @rtype float array
154 """
155
156
157 epsilon = centralRadial / largeSep - centralRadial // largeSep
158 radialOrderRange = np.arange(-numberOfRadialOrders,numberOfRadialOrders+1,dtype='float')
159 dominantRadialOrderPressure = np.floor(centralRadial/largeSep)
160 logger.info('dominantRadialOrderPressure {0}'.format(dominantRadialOrderPressure))
161 univPatternPerOrder = []
162
163 for nnn in radialOrderRange:
164 radial = purePressureMode(nnn,centralRadial, largeSep, smallSep, epsilon, 0.)
165 dipole = purePressureMode(nnn,centralRadial, largeSep, smallSep, epsilon, 1.)
166 quadrupole = purePressureMode(nnn,centralRadial, largeSep, smallSep, epsilon, 2.)
167 septupole = purePressureMode(nnn,centralRadial, largeSep, smallSep, epsilon, 3.)
168
169 univPatternPerOrder.append([radial,dipole,quadrupole,septupole,nnn+dominantRadialOrderPressure])
170
171 univPatternPerOrder = np.array(univPatternPerOrder, dtype='float')
172
173
174 logger.info('radial\t\tdipole\t\tquadrupole\tseptupole\tRadial Order (Pressure)')
175 logger.info('--------------------------------------------------------------------------------------')
176 for modeValue in univPatternPerOrder: logger.info('{0:.3f}\t\t{1:.3f}\t\t{2:.3f}\t\t{3:.3f}\t\t{4}'.format(modeValue[0],modeValue[1],modeValue[2],modeValue[3],modeValue[4]))
177 logger.info('--------------------------------------------------------------------------------------')
178
179 return univPatternPerOrder.T
180
181
182
183
184
185 -def asymptoticRelation(centralRadial, largeSep, DeltaPg, qCoupling= 0.15,
186 approximationTreshold=0.001, numberOfRadialOrders = 4):
187 """ Calculates and returns the frequencies mixed dipole modes in the powerspectrum of a red-giant star.
188
189 This function the modulation of the mixing between pressure- and gravity-dipole and derives the frequencies of the mixed modes from it.
190
191 @param centralRadial: Frequency of the central radial mode, in muHz
192 @type centralRadial: float
193 @param largeSep: large Frequency separation in powerspectrum, in muHz
194 @type largeSep: float
195 @param smallSep: small Frequency separation in powerspectrum, in muHz
196 @type smallSep: float
197 @param DeltaPg: value of the true period spacing
198 @type DeltaPg: float
199 @param qCoupling: coupling factor
200 @type qCoupling: float
201
202 @param approximationTreshold: prefiltering of the solution. if modes seem to be missing, increase this parameter
203 @type approximationTreshold:float
204 @param numberOfRadialOrders: for how many radial orders above and below the Central Radial Mode should the frequencies be calculated.
205 @type numberOfRadialOrders: float or integer
206 @return frequencies of the m=0 component of the mixed dipole modes
207 @rtype float array
208 """
209
210 DeltaPg /= 10.**6.
211 epsilon = centralRadial / largeSep - centralRadial // largeSep
212 radialOrderRange = np.arange(-numberOfRadialOrders,numberOfRadialOrders+1.,dtype='float')
213 dominantRadialOrderPressure = np.floor(centralRadial/largeSep)
214 univPatternPerOrder = []; positionMixedModes = []
215 rotationalSplitting = []
216 logger.info('calculating dipole mixed modes for {0} radial ordres ...'.format(len(radialOrderRange)))
217
218
219 for nnn in radialOrderRange:
220 valueDipole = purePressureMode(nnn,centralRadial, largeSep, 2., epsilon, 1.)
221 freqRange = np.linspace((centralRadial+largeSep*(nnn-0.25)), (centralRadial+largeSep*(nnn+1.25)), 100000)
222 equalityLine = freqRange/largeSep
223 arcTerm = largeSep / scipy.pi * np.arctan( qCoupling * np.tan( scipy.pi / (DeltaPg*freqRange) ) )
224
225
226
227 residualsToEquealityLine = np.array((freqRange/largeSep,(valueDipole+arcTerm)/largeSep-equalityLine))
228
229
230 filteredResidualsToEquealityLine = residualsToEquealityLine[:,(residualsToEquealityLine[1] < approximationTreshold/1.)]
231
232 """ # plotting the method to solve the implicit formula
233 pl.figure(num=2)
234 pl.plot(freqRange/largeSep,equalityLine,'r', alpha=0.5)
235 pl.plot(valueDipole/largeSep,valueDipole/largeSep,'mo')
236 pl.plot(freqRange/largeSep,(valueDipole+arcTerm)/largeSep,'b', alpha=0.5)
237 pl.plot(residualsToEquealityLine[0],residualsToEquealityLine[1],'k', alpha=0.5)
238 pl.plot(filteredResidualsToEquealityLine[0],filteredResidualsToEquealityLine[1])
239 #"""
240
241
242 for jjj in np.arange(0,len(filteredResidualsToEquealityLine[0])-1):
243 if filteredResidualsToEquealityLine[1,jjj] >= 0. and filteredResidualsToEquealityLine[1,jjj+1] <= 0.:
244 modePosition = scipy.average(filteredResidualsToEquealityLine[0,jjj:jjj+2])*largeSep
245 rotationalSplitting = 0.
246
247 positionMixedModes.append(modePosition)
248
249 dipoleMixedModes = np.array(positionMixedModes, dtype='float').T
250
251 return dipoleMixedModes
252
253
254
255
256
258 """ Calculates and returns rotational splitting
259
260 This function the modulation of the mixing between pressure- and gravity-dipole and derives the frequencies of the mixed modes from it.
261
262 @param dfMax: the largest rotational splitting measured for a given star
263 @type dfMax: float
264 @param largeSep: large Frequency separation in powerspectrum, in muHz
265 @type largeSep: float
266
267 @param frequenciesUniversalPattern: output from rg.universalPattern
268 @type frequenciesUniversalPattern: array, float
269
270 @param mixedModesPattern: output from rg.asymptoticRelation
271 @type mixedModesPattern: array, float
272
273 @param lambdaParam: parameter for fitting, Default: lambdaParam = 0.5
274 @type lambdaParam: float
275
276 @param beParam: parameter for fitting, Default: beParam = 0.08
277 @type beParam: float
278
279 @return array with the rotational splitting
280 @rtype array, float
281 """
282
283
284
285 radialModes = frequenciesUniversalPattern[0]
286 pureDipolePressure = frequenciesUniversalPattern[1]
287 mixedDipoleMode = mixedModesPattern
288
289 rotationalCentreMode = []; rotationalSplitting = []
290 for nnn in np.arange(0,len(radialModes.T)-1):
291 pureDipole = pureDipolePressure[ (pureDipolePressure > radialModes[nnn]) & (pureDipolePressure < radialModes[nnn+1])]
292 mixedDipoles= mixedDipoleMode[ (mixedDipoleMode > radialModes[nnn]) & (mixedDipoleMode < radialModes[nnn+1])]
293
294 for mode in mixedDipoles:
295 modulationProfile = 1.- lambdaParam / (1. + ( (mode - pureDipole[0]) / (beParam*largeSep) )**2. )
296 rotationalCentreMode.append(mode)
297 rotationalSplitting.append(modulationProfile)
298
299 rotationalCentreMode = np.array(rotationalCentreMode)
300 rotationalSplitting = np.array(rotationalSplitting).T*dfMax
301
302 rotationalSplitting = np.vstack((rotationalCentreMode,rotationalSplitting))
303 rotationalSplitting = np.array(rotationalSplitting, dtype='float')
304
305 return rotationalSplitting
306