FAOM API documentation
foam.gmode_rotation_scaling
Calculate g-mode period spacing patterns in the asymptotic regime using the TAR.
1"""Calculate g-mode period spacing patterns in the asymptotic regime using the TAR.""" 2 3# Class from T. Van Reeth 4import sys 5 6import astropy.units as u 7import numpy as np 8 9 10################################################################################ 11class Asymptotic(object): 12 """ 13 A python class to calculate g-mode period spacing patterns in the asymptotic regime using the TAR. 14 """ 15 16 def __init__(self, gyre_dir, kval=0, mval=1, nmin=1, nmax=150): 17 """ 18 Setting up the environment/asymptotic object to calculate g-mode period spacing patterns. 19 20 Parameters 21 ---------- 22 gyre_dir: string 23 The GYRE installation directory. 24 kval: int 25 Part of the mode identification of the pulsation pattern. For gravito-inertial 26 modes with an equivalent in a non-rotating star, k = l - |m|. Default value = 0. 27 mval: int 28 Part of the mode identification of the pulsation pattern: azimuthal order. 29 mval > 0 for prograde modes. Default value = 1. 30 nmin: int 31 The minimum radial order that we will (attempt to) calculate. Default value = 1. 32 nmax: int 33 The maximum radial order that we will (attempt to) calculate. Default value = 150. 34 """ 35 36 self.kval = int(kval) 37 self.mval = int(mval) 38 self.nvals = np.arange(nmin, nmax + 0.1, 1.0) 39 40 self.lam_fun = self._retrieve_laplacegrid(gyre_dir) 41 self.spin, self.lam, self.spinsqlam = self._sample_laplacegrid() 42 43 ################################################################################ 44 def _retrieve_laplacegrid(self, gyre_dir): 45 """ 46 Retrieving the function lambda(nu) given in GYRE. 47 48 Parameters 49 ---------- 50 self: asymptotic object 51 gyre_dir: string 52 The GYRE installation directory. 53 54 Returns 55 ---------- 56 lam_fun: function 57 A function to calculate lambda, given spin parameter values as input. 58 """ 59 60 if self.kval >= 0: 61 kstr = f"+{self.kval}" 62 else: 63 kstr = f"{self.kval}" 64 if self.mval >= 0: 65 mstr = f"+{self.mval}" 66 else: 67 mstr = f"{self.mval}" 68 69 infile = f"{gyre_dir}/data/tar/tar_fit.m{mstr}.k{kstr}.h5" 70 71 sys.path.append(gyre_dir + "/src/tar/") 72 import gyre_cheb_fit 73 import gyre_tar_fit 74 75 tf = gyre_tar_fit.TarFit.load(infile) 76 lam_fun = np.vectorize(tf.lam) 77 78 return lam_fun 79 80 ################################################################################ 81 def _sample_laplacegrid(self, spinmax=1000.0, spindensity=1.0): 82 """ 83 Sampling the function lambda(nu) that was set up in _retrieve_laplacegrid(). 84 This subroutine includes a har-coded custom sampling density function for the spin parameter. 85 86 Parameters 87 ---------- 88 self: asymptotic object 89 spinmax: float 90 The maximum spin parameter value for which lambda eigenvalues will be retrieved 91 (following the Laplace tidal equation). Default value = 1000. 92 spindensity: float 93 A scaling factor for the sampling density function. The default value (=1) results 94 in 20000 data points for the spin parameter range [0, 100]. 95 96 Returns 97 ---------- 98 spin: numpy array, dtype=float 99 The spin parameter values for which lambda eigenvalues are returned 100 lam: numpy array, dtype=float 101 The lambda eigenvalues corresponding to the spin parameter values in the array 'spin' 102 spinsqlam: numpy array, dtype=float 103 = spin * sqrt(lam) 104 """ 105 106 if (self.kval >= 0) & (self.mval != 0): 107 spinmin = -0.1 108 # Relatively "ad hoc" optimal sampling (based on "experience") 109 nspinmin = round( 110 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmin)) / np.log10(101.0)) ** 0.5 111 ) 112 nspinmax = round( 113 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmax)) / np.log10(101.0)) ** 0.5 114 ) 115 if (spinmin < 0.0) & (spinmax <= 0.0): 116 spin = ( 117 -( 118 10.0 119 ** ( 120 np.linspace( 121 np.log10(1.0 + abs(spinmin)) ** 0.5, 122 np.log10(1.0 + abs(spinmax)) ** 0.5, 123 int(nspinmin - nspinmax), 124 ) 125 ** 2.0 126 ) 127 ) 128 + 1.0 129 ) 130 elif (spinmin >= 0.0) & (spinmax > 0.0): 131 spin = ( 132 10.0 133 ** ( 134 np.linspace( 135 np.log10(1.0 + spinmin) ** 0.5, 136 np.log10(1.0 + spinmax) ** 0.5, 137 int(nspinmax - nspinmin), 138 ) 139 ** 2.0 140 ) 141 - 1.0 142 ) 143 else: 144 spinneg = ( 145 -( 146 10.0 147 ** ( 148 np.linspace(np.log10(1.0 + abs(spinmin)) ** 0.5, 0.0, int(nspinmin)) 149 ** 2.0 150 ) 151 ) 152 + 1.0 153 ) 154 spinpos = ( 155 10.0 ** (np.linspace(0.0, np.log10(1.0 + spinmax) ** 0.5, int(nspinmax)) ** 2.0) 156 - 1.0 157 ) 158 spin = np.unique(np.hstack((spinneg, spinpos))) 159 160 elif (self.kval >= 0) & (self.mval == 0): 161 spinmin = -spinmax 162 # Relatively "ad hoc" optimal sampling (based on "experience") 163 nspinmin = round( 164 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmin)) / np.log10(101.0)) ** 0.5 165 ) 166 nspinmax = round( 167 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmax)) / np.log10(101.0)) ** 0.5 168 ) 169 if (spinmin < 0.0) & (spinmax <= 0.0): 170 spin = ( 171 -( 172 10.0 173 ** ( 174 np.linspace( 175 np.log10(1.0 + abs(spinmin)) ** 0.5, 176 np.log10(1.0 + abs(spinmax)) ** 0.5, 177 int(nspinmin - nspinmax), 178 ) 179 ** 2.0 180 ) 181 ) 182 + 1.0 183 ) 184 elif (spinmin >= 0.0) & (spinmax > 0.0): 185 spin = ( 186 10.0 187 ** ( 188 np.linspace( 189 np.log10(1.0 + spinmin) ** 0.5, 190 np.log10(1.0 + spinmax) ** 0.5, 191 int(nspinmax - nspinmin), 192 ) 193 ** 2.0 194 ) 195 - 1.0 196 ) 197 else: 198 spinneg = ( 199 -( 200 10.0 201 ** ( 202 np.linspace(np.log10(1.0 + abs(spinmin)) ** 0.5, 0.0, int(nspinmin)) 203 ** 2.0 204 ) 205 ) 206 + 1.0 207 ) 208 spinpos = ( 209 10.0 ** (np.linspace(0.0, np.log10(1.0 + spinmax) ** 0.5, int(nspinmax)) ** 2.0) 210 - 1.0 211 ) 212 spin = np.unique(np.hstack((spinneg, spinpos))) 213 else: 214 spinmin = float( 215 (abs(self.mval) + abs(self.kval)) * (abs(self.mval) + abs(self.kval) - 1) 216 ) / abs(self.mval) 217 # Relatively "ad hoc" optimal sampling (based on "experience") 218 nspinmax = round( 219 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmax)) / np.log10(101.0)) ** 0.5 220 ) 221 222 spinpos = ( 223 10.0 224 ** ( 225 np.linspace( 226 np.log10(1.0) ** 0.5, np.log10(1.0 + spinmax) ** 0.5, int(nspinmax) 227 )[1:] 228 ** 2.0 229 ) 230 - 1.0 231 + spinmin 232 ) 233 spinneg = np.linspace(0.0, spinmin, 100) 234 spin = np.unique(np.hstack((spinneg, spinpos))) 235 236 # Obtaining the eigenvalues lambda 237 lam = self.lam_fun(spin) 238 239 # Limiting ourselves to the part of the arrays that exists... 240 lam_exists = np.isfinite(lam) 241 spin = spin[lam_exists] 242 lam = lam[lam_exists] 243 244 # Calculating the derivatives 245 dlamdnu = np.gradient(lam) / np.gradient(spin) 246 d2lamdnu2 = np.gradient(dlamdnu) / np.gradient(spin) 247 248 # A required array to determine the near-core rotation rate 249 spinsqlam = spin * np.sqrt(lam) 250 251 return spin, lam, spinsqlam 252 253 ################################################################################ 254 def update_laplacegrid(self, spinmax=10000.0, spindensity=1.0): 255 """ 256 Recalculating the spin parameter range and the corresponding lambda eigenvalues included within the asymptotic object. 257 258 Parameters 259 ---------- 260 self: asymptotic object 261 spinmax: float 262 The maximum spin parameter value for which lambda eigenvalues will be retrieved 263 (following the Laplace tidal equation). Default value = 1000. 264 spindensity: float 265 A scaling factor for the sampling density function. The default value (=1) results 266 in 20000 data points for the spin parameter range [0, 100]. 267 """ 268 269 self.spin, self.lam, self.spinsqlam = self._sample_laplacegrid( 270 spinmax=spinmax, spindensity=spindensity 271 ) 272 273 ################################################################################ 274 def uniform_pattern(self, frot, Pi0, alpha_g=0.5, unit="days"): 275 """ 276 Calculates the asymptotic period spacing pattern for a star with uniform rotation, following Eq.4 in Van Reeth et al. 2016, A&A, 593, A120. 277 278 Parameters 279 ---------- 280 self: asymptotic object 281 frot: astropy quantity (type frequency) 282 the rotation frequency of the star 283 Pi0: astropy quantity (type time) 284 the buoyancy radius of the star 285 alpha_g: float 286 phase term, dependent on the boundary conditions of the g-mode cavity. Typically, alpha_g = 1/2. 287 unit: string 288 The preferred (astropy) unit of the calculated g-mode pattern. The options are 'days', 'seconds', 289 'cycle_per_day', 'muHz', 'nHz', 'Hz'. Default = 'days'. 290 291 Returns 292 ---------- 293 pattern: array, dtype = astropy quantity matching the preferred 'unit'. 294 """ 295 296 ### Verifying the input 297 298 # Make sure that the rotation frequency and the buoyancy radius have appropriate units 299 assert ( 300 frot.unit.physical_type == "frequency" 301 ), "The provided rotation frequency does not have a frequency unit (as provided in astropy.units)." 302 assert ( 303 Pi0.unit.physical_type == "time" 304 ), "The provided buoyancy radius does not have a time unit (as provided in astropy.units)." 305 306 # Can we handle the preferred output unit? 307 allowed_output_units = { 308 "days": u.day, 309 "seconds": u.s, 310 "cycle_per_day": u.day**-1.0, 311 "muHz": u.uHz, 312 "nHz": u.nHz, 313 "Hz": u.Hz, 314 } 315 assert ( 316 unit in allowed_output_units.keys() 317 ), f"Please ensure that the requested output unit is one of: {', '.join(str(x_unit) for x_unit in allowed_output_units.keys())}." 318 astrounit = allowed_output_units[unit] 319 320 ### Alright, start the calculation! 321 322 # Safety check: the following routine only works if frot != 0 323 if frot != 0.0 * frot.unit: 324 # Calculating the pattern -- in spin parameter 325 basefun_lhs = 2.0 * float(frot * Pi0) * (self.nvals + alpha_g) 326 basefun_rhs = self.spinsqlam 327 selected_spin = np.interp(basefun_lhs, basefun_rhs, self.spin) 328 329 # Setting the values which lie outside of the computable range to "nan" 330 out_of_range = ~np.r_[ 331 (basefun_lhs >= np.amin(basefun_rhs)) & (basefun_lhs <= np.amax(basefun_rhs)) 332 ] 333 selected_spin[out_of_range] = np.nan 334 335 # Converting to frequencies -- in the unit of frot 336 puls_freq = (2.0 * frot / selected_spin) + (self.mval * frot) 337 338 # if frot == 0, the Taylor expansion cannot be used to model differential rotation: 339 # everything reverts to the case of non-rotation. 340 else: 341 if np.amin(self.spin) < 0.0: 342 # Pulsation frequencies -- in the unit of 1/Pi0 343 puls_freq = np.sqrt(np.interp(0.0, self.spin, self.lam)) / ( 344 Pi0 * (self.nvals + alpha_g) 345 ) 346 else: 347 # These are likely r-modes -- not present in a non-rotating star! 348 puls_freq = self.nvals * np.nan 349 350 # Converting to the preferred output unit 351 if astrounit.physical_type == "frequency": 352 pattern = puls_freq.to(astrounit) 353 else: 354 pattern = (1.0 / puls_freq).to(astrounit) 355 356 return pattern 357 358 ################################################################################ 359 def scale_pattern(self, pattern_in, frot_in, frot_out): 360 """ 361 Scale an input g-mode pulsation pattern according to a different rotation rate 362 363 Parameters 364 ---------- 365 self: asymptotic object 366 pattern_in: astropy quantity array 367 the input g-mode pulsation pattern (associated with the same mode identification as the asymptotic object) 368 frot_in: astropy quantity (frequency) 369 the cyclic rotation frequency associated with the input g-mode pulsation pattern 370 frot_out: astropy quantity (frequency) 371 the cyclic rotation frequency to which we want to scale the input g-mode pulsation pattern 372 373 Returns 374 ---------- 375 pattern_out: astropy quantity array 376 the output (scaled) g-mode pulsation pattern, with the same unit as pattern_in 377 """ 378 379 # checking if the input pattern has been given as frequencies or periods 380 if pattern_in.unit.physical_type == "frequency": 381 in_freq_co = pattern_in.copy() - self.mval * frot_in 382 else: 383 in_freq_co = (1.0 / pattern_in) - self.mval * frot_in 384 385 spin_in = 2.0 * frot_in.value / in_freq_co.to(frot_in.unit).value 386 spin_sqlam = np.interp(spin_in, self.spin, self.spinsqlam) 387 spin_out = np.interp( 388 spin_sqlam * frot_out.value / frot_in.to(frot_out.unit).value, self.spinsqlam, self.spin 389 ) 390 out_freq_inert = (2.0 * frot_out / spin_out) + self.mval * frot_out 391 392 if pattern_in.unit.physical_type == "frequency": 393 pattern_out = out_freq_inert.to(pattern_in.unit) 394 else: 395 out_per_inert = 1.0 / out_freq_inert 396 pattern_out = out_per_inert.to(pattern_in.unit) 397 398 return pattern_out
class
Asymptotic:
12class Asymptotic(object): 13 """ 14 A python class to calculate g-mode period spacing patterns in the asymptotic regime using the TAR. 15 """ 16 17 def __init__(self, gyre_dir, kval=0, mval=1, nmin=1, nmax=150): 18 """ 19 Setting up the environment/asymptotic object to calculate g-mode period spacing patterns. 20 21 Parameters 22 ---------- 23 gyre_dir: string 24 The GYRE installation directory. 25 kval: int 26 Part of the mode identification of the pulsation pattern. For gravito-inertial 27 modes with an equivalent in a non-rotating star, k = l - |m|. Default value = 0. 28 mval: int 29 Part of the mode identification of the pulsation pattern: azimuthal order. 30 mval > 0 for prograde modes. Default value = 1. 31 nmin: int 32 The minimum radial order that we will (attempt to) calculate. Default value = 1. 33 nmax: int 34 The maximum radial order that we will (attempt to) calculate. Default value = 150. 35 """ 36 37 self.kval = int(kval) 38 self.mval = int(mval) 39 self.nvals = np.arange(nmin, nmax + 0.1, 1.0) 40 41 self.lam_fun = self._retrieve_laplacegrid(gyre_dir) 42 self.spin, self.lam, self.spinsqlam = self._sample_laplacegrid() 43 44 ################################################################################ 45 def _retrieve_laplacegrid(self, gyre_dir): 46 """ 47 Retrieving the function lambda(nu) given in GYRE. 48 49 Parameters 50 ---------- 51 self: asymptotic object 52 gyre_dir: string 53 The GYRE installation directory. 54 55 Returns 56 ---------- 57 lam_fun: function 58 A function to calculate lambda, given spin parameter values as input. 59 """ 60 61 if self.kval >= 0: 62 kstr = f"+{self.kval}" 63 else: 64 kstr = f"{self.kval}" 65 if self.mval >= 0: 66 mstr = f"+{self.mval}" 67 else: 68 mstr = f"{self.mval}" 69 70 infile = f"{gyre_dir}/data/tar/tar_fit.m{mstr}.k{kstr}.h5" 71 72 sys.path.append(gyre_dir + "/src/tar/") 73 import gyre_cheb_fit 74 import gyre_tar_fit 75 76 tf = gyre_tar_fit.TarFit.load(infile) 77 lam_fun = np.vectorize(tf.lam) 78 79 return lam_fun 80 81 ################################################################################ 82 def _sample_laplacegrid(self, spinmax=1000.0, spindensity=1.0): 83 """ 84 Sampling the function lambda(nu) that was set up in _retrieve_laplacegrid(). 85 This subroutine includes a har-coded custom sampling density function for the spin parameter. 86 87 Parameters 88 ---------- 89 self: asymptotic object 90 spinmax: float 91 The maximum spin parameter value for which lambda eigenvalues will be retrieved 92 (following the Laplace tidal equation). Default value = 1000. 93 spindensity: float 94 A scaling factor for the sampling density function. The default value (=1) results 95 in 20000 data points for the spin parameter range [0, 100]. 96 97 Returns 98 ---------- 99 spin: numpy array, dtype=float 100 The spin parameter values for which lambda eigenvalues are returned 101 lam: numpy array, dtype=float 102 The lambda eigenvalues corresponding to the spin parameter values in the array 'spin' 103 spinsqlam: numpy array, dtype=float 104 = spin * sqrt(lam) 105 """ 106 107 if (self.kval >= 0) & (self.mval != 0): 108 spinmin = -0.1 109 # Relatively "ad hoc" optimal sampling (based on "experience") 110 nspinmin = round( 111 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmin)) / np.log10(101.0)) ** 0.5 112 ) 113 nspinmax = round( 114 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmax)) / np.log10(101.0)) ** 0.5 115 ) 116 if (spinmin < 0.0) & (spinmax <= 0.0): 117 spin = ( 118 -( 119 10.0 120 ** ( 121 np.linspace( 122 np.log10(1.0 + abs(spinmin)) ** 0.5, 123 np.log10(1.0 + abs(spinmax)) ** 0.5, 124 int(nspinmin - nspinmax), 125 ) 126 ** 2.0 127 ) 128 ) 129 + 1.0 130 ) 131 elif (spinmin >= 0.0) & (spinmax > 0.0): 132 spin = ( 133 10.0 134 ** ( 135 np.linspace( 136 np.log10(1.0 + spinmin) ** 0.5, 137 np.log10(1.0 + spinmax) ** 0.5, 138 int(nspinmax - nspinmin), 139 ) 140 ** 2.0 141 ) 142 - 1.0 143 ) 144 else: 145 spinneg = ( 146 -( 147 10.0 148 ** ( 149 np.linspace(np.log10(1.0 + abs(spinmin)) ** 0.5, 0.0, int(nspinmin)) 150 ** 2.0 151 ) 152 ) 153 + 1.0 154 ) 155 spinpos = ( 156 10.0 ** (np.linspace(0.0, np.log10(1.0 + spinmax) ** 0.5, int(nspinmax)) ** 2.0) 157 - 1.0 158 ) 159 spin = np.unique(np.hstack((spinneg, spinpos))) 160 161 elif (self.kval >= 0) & (self.mval == 0): 162 spinmin = -spinmax 163 # Relatively "ad hoc" optimal sampling (based on "experience") 164 nspinmin = round( 165 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmin)) / np.log10(101.0)) ** 0.5 166 ) 167 nspinmax = round( 168 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmax)) / np.log10(101.0)) ** 0.5 169 ) 170 if (spinmin < 0.0) & (spinmax <= 0.0): 171 spin = ( 172 -( 173 10.0 174 ** ( 175 np.linspace( 176 np.log10(1.0 + abs(spinmin)) ** 0.5, 177 np.log10(1.0 + abs(spinmax)) ** 0.5, 178 int(nspinmin - nspinmax), 179 ) 180 ** 2.0 181 ) 182 ) 183 + 1.0 184 ) 185 elif (spinmin >= 0.0) & (spinmax > 0.0): 186 spin = ( 187 10.0 188 ** ( 189 np.linspace( 190 np.log10(1.0 + spinmin) ** 0.5, 191 np.log10(1.0 + spinmax) ** 0.5, 192 int(nspinmax - nspinmin), 193 ) 194 ** 2.0 195 ) 196 - 1.0 197 ) 198 else: 199 spinneg = ( 200 -( 201 10.0 202 ** ( 203 np.linspace(np.log10(1.0 + abs(spinmin)) ** 0.5, 0.0, int(nspinmin)) 204 ** 2.0 205 ) 206 ) 207 + 1.0 208 ) 209 spinpos = ( 210 10.0 ** (np.linspace(0.0, np.log10(1.0 + spinmax) ** 0.5, int(nspinmax)) ** 2.0) 211 - 1.0 212 ) 213 spin = np.unique(np.hstack((spinneg, spinpos))) 214 else: 215 spinmin = float( 216 (abs(self.mval) + abs(self.kval)) * (abs(self.mval) + abs(self.kval) - 1) 217 ) / abs(self.mval) 218 # Relatively "ad hoc" optimal sampling (based on "experience") 219 nspinmax = round( 220 spindensity * 20000.0 * (np.log10(1.0 + abs(spinmax)) / np.log10(101.0)) ** 0.5 221 ) 222 223 spinpos = ( 224 10.0 225 ** ( 226 np.linspace( 227 np.log10(1.0) ** 0.5, np.log10(1.0 + spinmax) ** 0.5, int(nspinmax) 228 )[1:] 229 ** 2.0 230 ) 231 - 1.0 232 + spinmin 233 ) 234 spinneg = np.linspace(0.0, spinmin, 100) 235 spin = np.unique(np.hstack((spinneg, spinpos))) 236 237 # Obtaining the eigenvalues lambda 238 lam = self.lam_fun(spin) 239 240 # Limiting ourselves to the part of the arrays that exists... 241 lam_exists = np.isfinite(lam) 242 spin = spin[lam_exists] 243 lam = lam[lam_exists] 244 245 # Calculating the derivatives 246 dlamdnu = np.gradient(lam) / np.gradient(spin) 247 d2lamdnu2 = np.gradient(dlamdnu) / np.gradient(spin) 248 249 # A required array to determine the near-core rotation rate 250 spinsqlam = spin * np.sqrt(lam) 251 252 return spin, lam, spinsqlam 253 254 ################################################################################ 255 def update_laplacegrid(self, spinmax=10000.0, spindensity=1.0): 256 """ 257 Recalculating the spin parameter range and the corresponding lambda eigenvalues included within the asymptotic object. 258 259 Parameters 260 ---------- 261 self: asymptotic object 262 spinmax: float 263 The maximum spin parameter value for which lambda eigenvalues will be retrieved 264 (following the Laplace tidal equation). Default value = 1000. 265 spindensity: float 266 A scaling factor for the sampling density function. The default value (=1) results 267 in 20000 data points for the spin parameter range [0, 100]. 268 """ 269 270 self.spin, self.lam, self.spinsqlam = self._sample_laplacegrid( 271 spinmax=spinmax, spindensity=spindensity 272 ) 273 274 ################################################################################ 275 def uniform_pattern(self, frot, Pi0, alpha_g=0.5, unit="days"): 276 """ 277 Calculates the asymptotic period spacing pattern for a star with uniform rotation, following Eq.4 in Van Reeth et al. 2016, A&A, 593, A120. 278 279 Parameters 280 ---------- 281 self: asymptotic object 282 frot: astropy quantity (type frequency) 283 the rotation frequency of the star 284 Pi0: astropy quantity (type time) 285 the buoyancy radius of the star 286 alpha_g: float 287 phase term, dependent on the boundary conditions of the g-mode cavity. Typically, alpha_g = 1/2. 288 unit: string 289 The preferred (astropy) unit of the calculated g-mode pattern. The options are 'days', 'seconds', 290 'cycle_per_day', 'muHz', 'nHz', 'Hz'. Default = 'days'. 291 292 Returns 293 ---------- 294 pattern: array, dtype = astropy quantity matching the preferred 'unit'. 295 """ 296 297 ### Verifying the input 298 299 # Make sure that the rotation frequency and the buoyancy radius have appropriate units 300 assert ( 301 frot.unit.physical_type == "frequency" 302 ), "The provided rotation frequency does not have a frequency unit (as provided in astropy.units)." 303 assert ( 304 Pi0.unit.physical_type == "time" 305 ), "The provided buoyancy radius does not have a time unit (as provided in astropy.units)." 306 307 # Can we handle the preferred output unit? 308 allowed_output_units = { 309 "days": u.day, 310 "seconds": u.s, 311 "cycle_per_day": u.day**-1.0, 312 "muHz": u.uHz, 313 "nHz": u.nHz, 314 "Hz": u.Hz, 315 } 316 assert ( 317 unit in allowed_output_units.keys() 318 ), f"Please ensure that the requested output unit is one of: {', '.join(str(x_unit) for x_unit in allowed_output_units.keys())}." 319 astrounit = allowed_output_units[unit] 320 321 ### Alright, start the calculation! 322 323 # Safety check: the following routine only works if frot != 0 324 if frot != 0.0 * frot.unit: 325 # Calculating the pattern -- in spin parameter 326 basefun_lhs = 2.0 * float(frot * Pi0) * (self.nvals + alpha_g) 327 basefun_rhs = self.spinsqlam 328 selected_spin = np.interp(basefun_lhs, basefun_rhs, self.spin) 329 330 # Setting the values which lie outside of the computable range to "nan" 331 out_of_range = ~np.r_[ 332 (basefun_lhs >= np.amin(basefun_rhs)) & (basefun_lhs <= np.amax(basefun_rhs)) 333 ] 334 selected_spin[out_of_range] = np.nan 335 336 # Converting to frequencies -- in the unit of frot 337 puls_freq = (2.0 * frot / selected_spin) + (self.mval * frot) 338 339 # if frot == 0, the Taylor expansion cannot be used to model differential rotation: 340 # everything reverts to the case of non-rotation. 341 else: 342 if np.amin(self.spin) < 0.0: 343 # Pulsation frequencies -- in the unit of 1/Pi0 344 puls_freq = np.sqrt(np.interp(0.0, self.spin, self.lam)) / ( 345 Pi0 * (self.nvals + alpha_g) 346 ) 347 else: 348 # These are likely r-modes -- not present in a non-rotating star! 349 puls_freq = self.nvals * np.nan 350 351 # Converting to the preferred output unit 352 if astrounit.physical_type == "frequency": 353 pattern = puls_freq.to(astrounit) 354 else: 355 pattern = (1.0 / puls_freq).to(astrounit) 356 357 return pattern 358 359 ################################################################################ 360 def scale_pattern(self, pattern_in, frot_in, frot_out): 361 """ 362 Scale an input g-mode pulsation pattern according to a different rotation rate 363 364 Parameters 365 ---------- 366 self: asymptotic object 367 pattern_in: astropy quantity array 368 the input g-mode pulsation pattern (associated with the same mode identification as the asymptotic object) 369 frot_in: astropy quantity (frequency) 370 the cyclic rotation frequency associated with the input g-mode pulsation pattern 371 frot_out: astropy quantity (frequency) 372 the cyclic rotation frequency to which we want to scale the input g-mode pulsation pattern 373 374 Returns 375 ---------- 376 pattern_out: astropy quantity array 377 the output (scaled) g-mode pulsation pattern, with the same unit as pattern_in 378 """ 379 380 # checking if the input pattern has been given as frequencies or periods 381 if pattern_in.unit.physical_type == "frequency": 382 in_freq_co = pattern_in.copy() - self.mval * frot_in 383 else: 384 in_freq_co = (1.0 / pattern_in) - self.mval * frot_in 385 386 spin_in = 2.0 * frot_in.value / in_freq_co.to(frot_in.unit).value 387 spin_sqlam = np.interp(spin_in, self.spin, self.spinsqlam) 388 spin_out = np.interp( 389 spin_sqlam * frot_out.value / frot_in.to(frot_out.unit).value, self.spinsqlam, self.spin 390 ) 391 out_freq_inert = (2.0 * frot_out / spin_out) + self.mval * frot_out 392 393 if pattern_in.unit.physical_type == "frequency": 394 pattern_out = out_freq_inert.to(pattern_in.unit) 395 else: 396 out_per_inert = 1.0 / out_freq_inert 397 pattern_out = out_per_inert.to(pattern_in.unit) 398 399 return pattern_out
A python class to calculate g-mode period spacing patterns in the asymptotic regime using the TAR.
Asymptotic(gyre_dir, kval=0, mval=1, nmin=1, nmax=150)
17 def __init__(self, gyre_dir, kval=0, mval=1, nmin=1, nmax=150): 18 """ 19 Setting up the environment/asymptotic object to calculate g-mode period spacing patterns. 20 21 Parameters 22 ---------- 23 gyre_dir: string 24 The GYRE installation directory. 25 kval: int 26 Part of the mode identification of the pulsation pattern. For gravito-inertial 27 modes with an equivalent in a non-rotating star, k = l - |m|. Default value = 0. 28 mval: int 29 Part of the mode identification of the pulsation pattern: azimuthal order. 30 mval > 0 for prograde modes. Default value = 1. 31 nmin: int 32 The minimum radial order that we will (attempt to) calculate. Default value = 1. 33 nmax: int 34 The maximum radial order that we will (attempt to) calculate. Default value = 150. 35 """ 36 37 self.kval = int(kval) 38 self.mval = int(mval) 39 self.nvals = np.arange(nmin, nmax + 0.1, 1.0) 40 41 self.lam_fun = self._retrieve_laplacegrid(gyre_dir) 42 self.spin, self.lam, self.spinsqlam = self._sample_laplacegrid()
Setting up the environment/asymptotic object to calculate g-mode period spacing patterns.
Parameters
- gyre_dir (string): The GYRE installation directory.
- kval (int): Part of the mode identification of the pulsation pattern. For gravito-inertial modes with an equivalent in a non-rotating star, k = l - |m|. Default value = 0.
- mval (int): Part of the mode identification of the pulsation pattern: azimuthal order. mval > 0 for prograde modes. Default value = 1.
- nmin (int): The minimum radial order that we will (attempt to) calculate. Default value = 1.
- nmax (int): The maximum radial order that we will (attempt to) calculate. Default value = 150.
def
update_laplacegrid(self, spinmax=10000.0, spindensity=1.0):
255 def update_laplacegrid(self, spinmax=10000.0, spindensity=1.0): 256 """ 257 Recalculating the spin parameter range and the corresponding lambda eigenvalues included within the asymptotic object. 258 259 Parameters 260 ---------- 261 self: asymptotic object 262 spinmax: float 263 The maximum spin parameter value for which lambda eigenvalues will be retrieved 264 (following the Laplace tidal equation). Default value = 1000. 265 spindensity: float 266 A scaling factor for the sampling density function. The default value (=1) results 267 in 20000 data points for the spin parameter range [0, 100]. 268 """ 269 270 self.spin, self.lam, self.spinsqlam = self._sample_laplacegrid( 271 spinmax=spinmax, spindensity=spindensity 272 )
Recalculating the spin parameter range and the corresponding lambda eigenvalues included within the asymptotic object.
Parameters
self (asymptotic object):
spinmax (float): The maximum spin parameter value for which lambda eigenvalues will be retrieved (following the Laplace tidal equation). Default value = 1000.
- spindensity (float): A scaling factor for the sampling density function. The default value (=1) results in 20000 data points for the spin parameter range [0, 100].
def
uniform_pattern(self, frot, Pi0, alpha_g=0.5, unit='days'):
275 def uniform_pattern(self, frot, Pi0, alpha_g=0.5, unit="days"): 276 """ 277 Calculates the asymptotic period spacing pattern for a star with uniform rotation, following Eq.4 in Van Reeth et al. 2016, A&A, 593, A120. 278 279 Parameters 280 ---------- 281 self: asymptotic object 282 frot: astropy quantity (type frequency) 283 the rotation frequency of the star 284 Pi0: astropy quantity (type time) 285 the buoyancy radius of the star 286 alpha_g: float 287 phase term, dependent on the boundary conditions of the g-mode cavity. Typically, alpha_g = 1/2. 288 unit: string 289 The preferred (astropy) unit of the calculated g-mode pattern. The options are 'days', 'seconds', 290 'cycle_per_day', 'muHz', 'nHz', 'Hz'. Default = 'days'. 291 292 Returns 293 ---------- 294 pattern: array, dtype = astropy quantity matching the preferred 'unit'. 295 """ 296 297 ### Verifying the input 298 299 # Make sure that the rotation frequency and the buoyancy radius have appropriate units 300 assert ( 301 frot.unit.physical_type == "frequency" 302 ), "The provided rotation frequency does not have a frequency unit (as provided in astropy.units)." 303 assert ( 304 Pi0.unit.physical_type == "time" 305 ), "The provided buoyancy radius does not have a time unit (as provided in astropy.units)." 306 307 # Can we handle the preferred output unit? 308 allowed_output_units = { 309 "days": u.day, 310 "seconds": u.s, 311 "cycle_per_day": u.day**-1.0, 312 "muHz": u.uHz, 313 "nHz": u.nHz, 314 "Hz": u.Hz, 315 } 316 assert ( 317 unit in allowed_output_units.keys() 318 ), f"Please ensure that the requested output unit is one of: {', '.join(str(x_unit) for x_unit in allowed_output_units.keys())}." 319 astrounit = allowed_output_units[unit] 320 321 ### Alright, start the calculation! 322 323 # Safety check: the following routine only works if frot != 0 324 if frot != 0.0 * frot.unit: 325 # Calculating the pattern -- in spin parameter 326 basefun_lhs = 2.0 * float(frot * Pi0) * (self.nvals + alpha_g) 327 basefun_rhs = self.spinsqlam 328 selected_spin = np.interp(basefun_lhs, basefun_rhs, self.spin) 329 330 # Setting the values which lie outside of the computable range to "nan" 331 out_of_range = ~np.r_[ 332 (basefun_lhs >= np.amin(basefun_rhs)) & (basefun_lhs <= np.amax(basefun_rhs)) 333 ] 334 selected_spin[out_of_range] = np.nan 335 336 # Converting to frequencies -- in the unit of frot 337 puls_freq = (2.0 * frot / selected_spin) + (self.mval * frot) 338 339 # if frot == 0, the Taylor expansion cannot be used to model differential rotation: 340 # everything reverts to the case of non-rotation. 341 else: 342 if np.amin(self.spin) < 0.0: 343 # Pulsation frequencies -- in the unit of 1/Pi0 344 puls_freq = np.sqrt(np.interp(0.0, self.spin, self.lam)) / ( 345 Pi0 * (self.nvals + alpha_g) 346 ) 347 else: 348 # These are likely r-modes -- not present in a non-rotating star! 349 puls_freq = self.nvals * np.nan 350 351 # Converting to the preferred output unit 352 if astrounit.physical_type == "frequency": 353 pattern = puls_freq.to(astrounit) 354 else: 355 pattern = (1.0 / puls_freq).to(astrounit) 356 357 return pattern
Calculates the asymptotic period spacing pattern for a star with uniform rotation, following Eq.4 in Van Reeth et al. 2016, A&A, 593, A120.
Parameters
self (asymptotic object):
frot (astropy quantity (type frequency)): the rotation frequency of the star
- Pi0 (astropy quantity (type time)): the buoyancy radius of the star
- alpha_g (float): phase term, dependent on the boundary conditions of the g-mode cavity. Typically, alpha_g = 1/2.
- unit (string): The preferred (astropy) unit of the calculated g-mode pattern. The options are 'days', 'seconds', 'cycle_per_day', 'muHz', 'nHz', 'Hz'. Default = 'days'.
Returns
- pattern (array, dtype = astropy quantity matching the preferred 'unit'.):
def
scale_pattern(self, pattern_in, frot_in, frot_out):
360 def scale_pattern(self, pattern_in, frot_in, frot_out): 361 """ 362 Scale an input g-mode pulsation pattern according to a different rotation rate 363 364 Parameters 365 ---------- 366 self: asymptotic object 367 pattern_in: astropy quantity array 368 the input g-mode pulsation pattern (associated with the same mode identification as the asymptotic object) 369 frot_in: astropy quantity (frequency) 370 the cyclic rotation frequency associated with the input g-mode pulsation pattern 371 frot_out: astropy quantity (frequency) 372 the cyclic rotation frequency to which we want to scale the input g-mode pulsation pattern 373 374 Returns 375 ---------- 376 pattern_out: astropy quantity array 377 the output (scaled) g-mode pulsation pattern, with the same unit as pattern_in 378 """ 379 380 # checking if the input pattern has been given as frequencies or periods 381 if pattern_in.unit.physical_type == "frequency": 382 in_freq_co = pattern_in.copy() - self.mval * frot_in 383 else: 384 in_freq_co = (1.0 / pattern_in) - self.mval * frot_in 385 386 spin_in = 2.0 * frot_in.value / in_freq_co.to(frot_in.unit).value 387 spin_sqlam = np.interp(spin_in, self.spin, self.spinsqlam) 388 spin_out = np.interp( 389 spin_sqlam * frot_out.value / frot_in.to(frot_out.unit).value, self.spinsqlam, self.spin 390 ) 391 out_freq_inert = (2.0 * frot_out / spin_out) + self.mval * frot_out 392 393 if pattern_in.unit.physical_type == "frequency": 394 pattern_out = out_freq_inert.to(pattern_in.unit) 395 else: 396 out_per_inert = 1.0 / out_freq_inert 397 pattern_out = out_per_inert.to(pattern_in.unit) 398 399 return pattern_out
Scale an input g-mode pulsation pattern according to a different rotation rate
Parameters
self (asymptotic object):
pattern_in (astropy quantity array): the input g-mode pulsation pattern (associated with the same mode identification as the asymptotic object)
- frot_in (astropy quantity (frequency)): the cyclic rotation frequency associated with the input g-mode pulsation pattern
- frot_out (astropy quantity (frequency)): the cyclic rotation frequency to which we want to scale the input g-mode pulsation pattern
Returns
- pattern_out (astropy quantity array): the output (scaled) g-mode pulsation pattern, with the same unit as pattern_in