diff --git a/scopesim_templates/__init__.py b/scopesim_templates/__init__.py index 4aad28a..76a97c4 100644 --- a/scopesim_templates/__init__.py +++ b/scopesim_templates/__init__.py @@ -17,6 +17,7 @@ from .misc.misc import source_from_image, source_from_file, source_from_cube from .stellar.stellar import star, stars, star_grid, star_field from .stellar.clusters import cluster +from .stellar.galactic_centre import galactic_centre from .extragalactic.galaxies import galaxy, galaxy3d, elliptical from .calibration.calibration import empty_sky from .micado import darkness, flatlamp diff --git a/scopesim_templates/stellar/J_ApJ_837_30_table3.dat.fits b/scopesim_templates/stellar/J_ApJ_837_30_table3.dat.fits new file mode 100644 index 0000000..06f9f04 --- /dev/null +++ b/scopesim_templates/stellar/J_ApJ_837_30_table3.dat.fits @@ -0,0 +1 @@ +SIMPLE = T / Standard FITS Format BITPIX = 8 / Character data NAXIS = 0 / No Image --- just extension(s) EXTEND = T / There are standard extensions ORIGIN = 'CDS ' / File generated at CDS, Strasbourg, France (tofits, Version 3.4) e-mail: question@simbad.u-strasbg.fr COMMENT ARG='-m -1 /ftp/cats/J/ApJ/837/30/./table3.dat' LONGSTRN= 'OGIP 1.0' / Long string convention (&/CONTINUE) may be usedDATE = '2026-05-14' / Written on 2026-05-14:10:34:39 (GMT) by: www-data@cdsarc.cds.unistra.fr CDS-CAT = 'J/ApJ/837/30' / Catalogue designation in CDS nomenclature COMMENT 25yrs monitoring of stellar orbits in the GC (Gillessen+, 2017) COMMENT Title: An update on monitoring stellar orbits in the Galactic Center. AUTHOR = 'Gillessen S., Plewa P.M., Eisenhauer F., Sari R., Waisberg I., Habi&'CONTINUE 'bi M.,&' CONTINUE 'Pfuhl O., George E., Dexter J., von Fellenberg S., Ott T., Genzel R&'CONTINUE '. ' REFERENC= 'Astrophys. J., 837, 30-30 (2017)' BIBCODE = '2017ApJ...837...30G' / 19-digit SIMBAD/NED/ADS BibCode COMMENT ADC_Keywords: Milky Way ; Radial velocities ; Photometry, infrared COMMENT Keywords: astrometry; black hole physics; Galaxy: center; Galaxy: fundamental parameters; techniques: high angular resolution COMMENT Abstract: Using 25 years of data from uninterrupted monitoring of stellar orbits in the Galactic Center (GC), we present an update of the main results from this unique data set: a measurement of mass and distance to SgrA*. Our progress is not only due to the eight-year increase in time base, but also to the improved definition of the coordinate system. The star S2 continues to yield the best constraints on the mass of and distance to Sgr A*; the statistical errors of 0.13x10^6^M_{sun}_ and 0.12kpc have halved compared to the previous study. The S2 orbit fit is robust and does not need any prior information. Using coordinate system priors, the star S1 also yields tight constraints on mass and distance. For a combined orbit fit, we use 17 stars, which yields our current best estimates for mass and distance: M=4.28+/-0.10|_stat._+/-0.21|_sys_x10^6^M_{sun}_ and R_0_=8.32+/-0.07|_stat._+/-0.14|_sys_kpc. These numbers are in agreement with the recent determination of R_0_ from the statistical cluster parallax. The positions of the mass, of the near-infrared flares from Sgr A*, and of the radio source Sgr A* agree to within 1mas. In total, we have determined orbits for 40 stars so far, a sample which consists of 32 stars with randomly oriented orbits and a thermal eccentricity distribution, plus eight stars that we can explicitly show are members of the clockwise disk of young stars, and which have lower-eccentricity orbits. COMMENT Description: This work is an update and improvement over our previous orbital study (Gillessen+ 2009ApJ...692.1075G). Since the previous work we have added eight more years of data, using the adaptive optics (AO) imager NACO on the VLT and the AO-assisted integral field spectrograph SINFONI. This extends our time base from 17 to 25 years for the imaging and from five to 13 years for the spectroscopy. We add (in the best cases) 78 epochs of imaging and 19 epochs of spectroscopy. Our imaging data set contains an interruption. During 2014 and in spring 2015, NACO was not available at the VLT, resulting in significant gaps in our time series. The first high-resolution imaging data of the GC region were obtained in 1992 with the SHARP camera at ESO's 3.5m NTT in Chile. The first AO imaging data available to us of the GC region were obtained in 2002 with the Naos-Conica (NACO) system mounted at the fourth unit telescope Yepun of the VLT (see Gillessen+ 2009ApJ...692.1075G) for more details. COMMENT File Summary: ----------------------------------------------------------------------- FileName Lrecl Records Explanations ----------------------------------------------------------------------- ReadMe 80 . This file table3.dat 118 40 Orbital parameters of the 40 stars for which we were able to determine orbits table5.dat 557 216 *Orbital data used in this work -----------------------------------------------------------------------COMMENT Note on table5.dat: The astrometric data presented here will change in future publications on monitoring stellar orbits, since they refer to the current definition of the reference system, which will improve in the future. Using these data requires some care concerning the coordinate system. The position and the velocity of the origin of the coordinate system in which these data are given do not correspond directly to the central point mass, but only to our best estimate where the radio source Sgr A* is located in the coordinate system (Plewa+ 2015MNRAS.453.3234P ; ({alpha} ,{delta})=(0,0)+/-(0.2,0.2)mas at our reference epoch 2009.0). Even Keplerian orbits therefore are expected to yield non-closing orbit figures, see figure 2 and table 1. Non-closing orbits do thus not mean automatically that post-Newtonian effects have been detected, or that the mass distribution is extended. Due to the uncertainty in the reference frame, it is also non-trivial to add to these data other astrometric measurements, see Gillessen+ (2009ApJ...692.1075G). Other authors wishing to use the data presented here are invited to contact the authors of this work at COMMENT See also: J/ApJ/697/1741 : Warped disks of YSOs in Galactic center (Bartko+, 2009) J/ApJ/707/L114 : The orbit of S2 star around Sgr A* (Gillessen+, 2009) J/ApJ/737/73 : Infrared extinction toward the Galactic Centre (Fritz+, 2011) J/ApJ/783/131 : Kinematic of stars in Galactic center (Yelda+, 2014) J/A+A/570/A2 : Spectra & RVs of nuclear stars (Feldmeier+, 2014) J/ApJ/821/44 : Star motions in the nuclear cluster of the MW (Fritz+, 2016) J/ApJ/830/17 : GC secondary IR astrometric standards (Boehle+, 2016) COMMENT Unused Described file: table3.dat COMMENT Unused Described file: table5.dat COMMENT Global notes: COMMENT History: From electronic version of the journal =======================================================================HISTORY (End) Prepared by [AAS], Emmanuelle Perret [CDS] 11-Oct-2017 END XTENSION= 'TABLE ' / Ascii Table Extension BITPIX = 8 / Character data NAXIS = 2 / Simple 2-D matrix NAXIS1 = 118 / Number of bytes per record NAXIS2 = 40 / Number of records PCOUNT = 0 / Get rid of random parameters GCOUNT = 1 / Only one group (isn't it obvious?) TFIELDS = 19 / Number of data fields (columns) EXTNAME = 'table3.dat' / Orbital parameters of the 40 stars for which we were able to determine orbits =======================================================================TBCOL1 = 1 / ============== Start column +0 TFORM1 = 'A4 ' / Fortran Format TTYPE1 = 'Star ' / Star Name TBCOL2 = 6 / ============== Start column +5 TFORM2 = 'A1 ' / Fortran Format TTYPE2 = 'f_Star ' / [a] Flag on Object S111 (1) TASET2 = 'a ' / Allowed character set TBCOL3 = 8 / ============== Start column +7 TUNIT3 = 'arcsec ' / Unit: second of arc TFORM3 = 'F8.4 ' / Fortran Format TDISP3 = 'F8.4 ' / Display Format for Binary Tables TTYPE3 = 'a ' / [-12.3/4.6] Semi-major axis of the orbit, from Sag A* TAMIN3 = -12.3000 / Allowed minimal value TAMAX3 = 4.6000 / Allowed maximal value TBCOL4 = 17 / ============== Start column +16 TUNIT4 = 'arcsec ' / Unit: second of arc TFORM4 = 'F6.4 ' / Fortran Format TDISP4 = 'F6.4 ' / Display Format for Binary Tables TTYPE4 = 'e_a ' / [0.0002/8.4] Uncertainty in a TAMIN4 = 0.0002 / Allowed minimal value TAMAX4 = 8.4000 / Allowed maximal value TBCOL5 = 24 / ============== Start column +23 TFORM5 = 'F6.4 ' / Fortran Format TDISP5 = 'F6.4 ' / Display Format for Binary Tables TTYPE5 = 'e ' / [0.1/1.1] Eccentricity TAMIN5 = 0.1000 / Allowed minimal value TAMAX5 = 1.1000 / Allowed maximal value TBCOL6 = 31 / ============== Start column +30 TFORM6 = 'F6.4 ' / Fortran Format TDISP6 = 'F6.4 ' / Display Format for Binary Tables TTYPE6 = 'e_e ' / [0.0003/0.3] Uncertainty in e TAMIN6 = 0.0003 / Allowed minimal value TAMAX6 = 0.3000 / Allowed maximal value TBCOL7 = 38 / ============== Start column +37 TUNIT7 = 'deg ' / Unit: degree TFORM7 = 'F6.2 ' / Fortran Format TDISP7 = 'F6.2 ' / Display Format for Binary Tables TTYPE7 = 'i ' / [24.7/171.1] Inclination TAMIN7 = 24.70 / Allowed minimal value TAMAX7 = 171.10 / Allowed maximal value TBCOL8 = 45 / ============== Start column +44 TUNIT8 = 'deg ' / Unit: degree TFORM8 = 'F4.2 ' / Fortran Format TDISP8 = 'F4.2 ' / Display Format for Binary Tables TTYPE8 = 'e_i ' / [0.06/8.3] Uncertainty in i TAMIN8 = 0.06 / Allowed minimal value TAMAX8 = 8.30 / Allowed maximal value TBCOL9 = 50 / ============== Start column +49 TUNIT9 = 'deg ' / Unit: degree TFORM9 = 'F6.2 ' / Fortran Format TDISP9 = 'F6.2 ' / Display Format for Binary Tables TTYPE9 = 'Omega ' / [7.9/345] Position angle of the ascending node TAMIN9 = 7.90 / Allowed minimal value TAMAX9 = 345.00 / Allowed maximal value TBCOL10 = 57 / ============== Start column +56 TUNIT10 = 'deg ' / Unit: degree TFORM10 = 'F5.2 ' / Fortran Format TDISP10 = 'F5.2 ' / Display Format for Binary Tables TTYPE10 = 'e_Omega ' / [0.07/19] Uncertainty in {Omega} TAMIN10 = 0.07 / Allowed minimal value TAMAX10 = 19.00 / Allowed maximal value TBCOL11 = 63 / ============== Start column +62 TUNIT11 = 'deg ' / Unit: degree TFORM11 = 'F6.2 ' / Fortran Format TDISP11 = 'F6.2 ' / Display Format for Binary Tables TTYPE11 = 'w ' / [17.9/357] Longitude of the pericenter TAMIN11 = 17.90 / Allowed minimal value TAMAX11 = 357.00 / Allowed maximal value TBCOL12 = 70 / ============== Start column +69 TUNIT12 = 'deg ' / Unit: degree TFORM12 = 'F5.2 ' / Fortran Format TDISP12 = 'F5.2 ' / Display Format for Binary Tables TTYPE12 = 'e_w ' / [0.07/24] Uncertainty in w TAMIN12 = 0.07 / Allowed minimal value TAMAX12 = 24.00 / Allowed maximal value TBCOL13 = 76 / ============== Start column +75 TUNIT13 = 'yr ' / Unit: year TFORM13 = 'F7.2 ' / Fortran Format TDISP13 = 'F7.2 ' / Display Format for Binary Tables TTYPE13 = 'Tp ' / [611/2132] Epoch of pericenter passage, fractional year (G1) TAMIN13 = 611.00 / Allowed minimal value TAMAX13 = 2132.00 / Allowed maximal value TBCOL14 = 84 / ============== Start column +83 TUNIT14 = 'yr ' / Unit: year TFORM14 = 'F6.2 ' / Fortran Format TDISP14 = 'F6.2 ' / Display Format for Binary Tables TTYPE14 = 'e_Tp ' / [0.01/154] Uncertainty in Tp TAMIN14 = 0.01 / Allowed minimal value TAMAX14 = 154.00 / Allowed maximal value TBCOL15 = 91 / ============== Start column +90 TUNIT15 = 'yr ' / Unit: year TFORM15 = 'F7.2 ' / Fortran Format TDISP15 = 'F7.2 ' / Display Format for Binary Tables TTYPE15 = 'Per ' / [12.8/3580]? Orbital period (G1) TAMIN15 = 12.80 / Allowed minimal value TAMAX15 = 3580.00 / Allowed maximal value TBNUL15 = ' ' / NULL (undefined) value TBCOL16 = 99 / ============== Start column +98 TUNIT16 = 'yr ' / Unit: year TFORM16 = 'F7.2 ' / Fortran Format TDISP16 = 'F7.2 ' / Display Format for Binary Tables TTYPE16 = 'e_Per ' / [0.02/2550]? Uncertainty in Per TAMIN16 = 0.02 / Allowed minimal value TAMAX16 = 2550.00 / Allowed maximal value TBNUL16 = ' ' / NULL (undefined) value TBCOL17 = 107 / ============== Start column +106 TFORM17 = 'A1 ' / Fortran Format TTYPE17 = 'SpT ' / [el] Spectral type: e=early-type stars, l=late-type stars TASET17 = 'el ' / Allowed character set TBCOL18 = 109 / ============== Start column +108 TUNIT18 = 'mag ' / Unit: magnitude TFORM18 = 'F5.2 ' / Fortran Format TDISP18 = 'F5.2 ' / Display Format for Binary Tables TTYPE18 = 'Kmag ' / [10/18] K-band magnitude TAMIN18 = 10.00 / Allowed minimal value TAMAX18 = 18.00 / Allowed maximal value TBCOL19 = 115 / ============== Start column +114 TFORM19 = 'F4.2 ' / Fortran Format TDISP19 = 'F4.2 ' / Display Format for Binary Tables TTYPE19 = 'r ' / [0.9/3.4] Global rescaling factor TAMIN19 = 0.90 / Allowed minimal value TAMAX19 = 3.40 / Allowed maximal value =======================================================================COMMENT Note (1): a = Object S111 formally has a negative semi major axis, indicative for a hyperbolic orbit with e>1 and no solution for an orbital period. =======================================================================END S1 0.595 0.024 0.556 0.018 119.14 0.21 342.04 0.32 122.3 1.4 2001.80 0.15 166.0 5.8 e 14.7 1.75S2 0.1255 0.0009 0.8839 0.0019 134.18 0.40 226.94 0.60 65.51 0.57 2002.33 0.01 16.00 0.02 e 13.95 1.13S4 0.3570 0.0037 0.3905 0.0059 80.33 0.08 258.84 0.07 290.8 1.5 1957.4 1.2 77.0 1.0 e 14.4 1.25S6 0.6574 0.0006 0.8400 0.0003 87.24 0.06 85.07 0.12 116.23 0.07 2108.61 0.03 192.0 0.17 e 15.4 1.58S8 0.4047 0.0014 0.8031 0.0075 74.37 0.30 315.43 0.19 346.70 0.41 1983.64 0.24 92.9 0.41 e 14.5 1.18S9 0.2724 0.0041 0.644 0.020 82.41 0.24 156.60 0.10 150.6 1.0 1976.71 0.92 51.3 0.70 e 15.1 1.65S12 0.2987 0.0018 0.8883 0.0017 33.56 0.49 230.1 1.8 317.9 1.5 1995.59 0.04 58.9 0.22 e 15.5 2.37S13 0.2641 0.0016 0.4250 0.0023 24.70 0.48 74.5 1.7 245.2 2.4 2004.86 0.04 49.00 0.14 e 15.8 3.25S14 0.2863 0.0036 0.9761 0.0037 100.59 0.87 226.38 0.64 334.59 0.87 2000.12 0.06 55.3 0.48 e 15.7 2.16S17 0.3559 0.0096 0.397 0.011 96.83 0.11 191.62 0.21 326.0 1.9 1991.19 0.41 76.6 1.0 l 15.3 3.00S18 0.2379 0.0015 0.471 0.012 110.67 0.18 49.11 0.18 349.46 0.66 1993.86 0.16 41.9 0.18 e 16.7 2.28S19 0.520 0.094 0.750 0.043 71.96 0.35 344.60 0.62 155.2 2.3 2005.39 0.16 135. 14. e 16. 2.57S21 0.2190 0.0017 0.764 0.014 58.8 1.0 259.64 0.62 166.4 1.1 2027.40 0.17 37.00 0.28 l 16.9 1.60S22 1.31 0.28 0.449 0.088 105.76 0.95 291.7 1.4 95. 20. 1996.9 10.2 540. 63. e 16.6 2.78S23 0.253 0.012 0.56 0.14 48.0 7.1 249. 13. 39.0 6.7 2024.7 3.7 45.8 1.6 e 17.8 2.08S24 0.944 0.048 0.8970 0.0049 103.67 0.42 7.93 0.37 290. 15. 2024.50 0.03 331. 16. l 15.6 1.54S29 0.428 0.019 0.728 0.052 105.8 1.7 161.96 0.80 346.5 5.9 2025.96 0.94 101.0 2.0 e 16.7 3.32S31 0.449 0.010 0.5497 0.0025 109.03 0.27 137.16 0.30 308.0 3.0 2018.07 0.14 108. 1.2 e 15.7 3.16S33 0.657 0.026 0.608 0.064 60.5 2.5 100.1 5.5 303.7 1.6 1928. 12. 192.0 5.2 e 16. 2.21S38 0.1416 0.0002 0.8201 0.0007 171.1 2.1 101.06 0.24 17.99 0.25 2003.19 0.01 19.2 0.02 l 17. 2.48S39 0.370 0.015 0.9236 0.0021 89.36 0.73 159.03 0.10 23.3 3.8 2000.06 0.06 81.1 1.5 16.8 3.27S42 0.95 0.18 0.567 0.083 67.16 0.66 196.14 0.75 35.8 3.2 2008.24 0.75 335. 58. e 17.5 1.65S54 1.20 0.87 0.893 0.078 62.2 1.4 288.35 0.70 140.8 2.3 2004.46 0.07 477. 199. e 17.5 2.60S55 0.1078 0.0010 0.7209 0.0077 150.1 2.2 325.5 4.0 331.5 3.9 2009.34 0.04 12.80 0.11 17.5 1.61S60 0.3877 0.0070 0.7179 0.0051 126.87 0.30 170.54 0.85 29.37 0.29 2023.89 0.09 87.1 1.4 e 16.3 1.65S66 1.502 0.095 0.128 0.043 128.5 1.6 92.3 3.2 134. 17. 1771. 38. 664. 37. e 14.8 1.70S67 1.126 0.026 0.293 0.057 136.0 1.1 96.5 6.4 213.5 1.6 1705. 22. 431. 10. e 12.1 1.43S71 0.973 0.040 0.899 0.013 74.0 1.3 35.16 0.86 337.8 4.9 1695. 21. 346. 11. e 16.1 1.87S83 1.49 0.19 0.365 0.075 127.2 1.4 87.7 1.2 203.6 6.0 2046.8 6.3 656. 69. e 13.6 1.82S85 4.6 3.30 0.78 0.15 84.78 0.29 107.36 0.43 156.3 6.8 1930.2 9.8 3580. 2550. l 15.6 1.50S87 2.74 0.16 0.224 0.027 119.54 0.87 106.32 0.99 336.1 7.7 611. 154. 1640. 105. e 13.6 1.38S89 1.081 0.055 0.639 0.038 87.61 0.16 238.99 0.18 126.4 4.0 1783. 26. 406. 27. l 15.3 1.16S91 1.917 0.089 0.303 0.034 114.49 0.32 105.35 0.74 356.4 1.6 1108. 69. 958. 50. e 12.2 1.33S96 1.499 0.057 0.174 0.022 126.36 0.96 115.66 0.59 233.6 2.4 1646. 16. 662. 29. e 10. 1.31S97 2.32 0.46 0.35 0.11 113.0 1.3 113.2 1.4 28. 14. 2132. 29. 1270. 309. e 10.3 1.22S111 a -12.3 8.4 1.092 0.064 102.68 0.40 52.34 0.75 132.4 3.3 1947.7 4.5 l 13.8 0.97S145 1.12 0.18 0.50 0.25 83.7 1.6 263.92 0.94 185. 16. 1808. 58. 426. 71. l 17.5 1.46S175 0.414 0.039 0.9867 0.0018 88.53 0.60 326.83 0.78 68.52 0.40 2009.51 0.01 96.2 5.0 e 17.5 2.72R34 1.81 0.15 0.641 0.098 136.0 8.3 330. 19. 57.0 8.0 1522. 52. 877. 83. e 14. 1.35R44 3.9 1.4 0.27 0.27 131.0 5.2 80.5 7.1 217. 24. 1963. 85. 2730. 1350. e 14. 1.11 \ No newline at end of file diff --git a/scopesim_templates/stellar/__init__.py b/scopesim_templates/stellar/__init__.py index 8437542..46ba9f9 100644 --- a/scopesim_templates/stellar/__init__.py +++ b/scopesim_templates/stellar/__init__.py @@ -3,3 +3,4 @@ from .clusters import * from .stellar import * +from .galactic_centre import galactic_centre diff --git a/scopesim_templates/stellar/_gillessen.py b/scopesim_templates/stellar/_gillessen.py new file mode 100644 index 0000000..dbf4efa --- /dev/null +++ b/scopesim_templates/stellar/_gillessen.py @@ -0,0 +1,186 @@ +"""Evaluate Gillessen+2017 Galactic Centre stellar orbits at arbitrary epochs. + +Loads the catalogue (CDS J/ApJ/837/30 table 3, 40 stars orbiting Sgr A*) and +returns, for any time, each star's projected sky position (d_ra, d_dec in +arcsec, east-positive / north-positive) and line-of-sight radial velocity +(km/s, positive = receding), alongside catalogue Kmag and spectral type. + +Public entry point: ``stars_at_time(time)``. + +Vendored from the standalone ``gillessen`` script collection. Catalogue +reference: Gillessen et al. 2017, ApJ 837, 30 — CDS table J/ApJ/837/30. +The companion FITS catalogue ``J_ApJ_837_30_table3.dat.fits`` is shipped +alongside this module. +""" + +from __future__ import annotations + +import warnings +from datetime import datetime +from numbers import Real +from pathlib import Path + +import numpy as np +from astropy import units as u +from astropy.table import QTable, Table +from astropy.time import Time + +DEFAULT_CATALOGUE = Path(__file__).resolve().parent / "J_ApJ_837_30_table3.dat.fits" +DEFAULT_R0 = 8.32 * u.kpc # Gillessen+2017 + + +def _load_catalogue(path: Path | str = DEFAULT_CATALOGUE) -> Table: + return Table.read(path) + + +def _to_jyear(t) -> float: + """Convert input time to Julian year (J2000-based, matches catalogue Tp). + + Accepts ``astropy.time.Time``, ``datetime``, ISO string, or numeric MJD. + """ + if isinstance(t, Time): + return float(t.jyear) + if isinstance(t, datetime): + return float(Time(t).jyear) + if isinstance(t, Real) and not isinstance(t, bool): + return float(Time(float(t), format="mjd").jyear) + return float(Time(t).jyear) + + +def _solve_kepler(M, e, tol: float = 1e-12, max_iter: int = 50) -> np.ndarray: + """Solve Kepler's equation ``M = E - e sin E`` for ``E``. + + Vectorised Newton iteration; ``M`` is wrapped to ``[-pi, pi)``. + """ + M = np.asarray(M, dtype=float) + e = np.asarray(e, dtype=float) + M = np.mod(M + np.pi, 2.0 * np.pi) - np.pi + E = M + e * np.sin(M) + for _ in range(max_iter): + f = E - e * np.sin(E) - M + fp = 1.0 - e * np.cos(E) + dE = f / fp + E = E - dE + if np.all(np.abs(dE) < tol): + break + return E + + +def _orbital_state(a, e, n, E): + """Position and velocity in the orbital plane. + + Pericenter on +x. Returns ``(x_orb, y_orb, vx_orb, vy_orb)`` in + ``(a-units, a-units / time-units-of-n)``. + """ + cosE = np.cos(E) + sinE = np.sin(E) + sqrt1me2 = np.sqrt(1.0 - e * e) + x_orb = a * (cosE - e) + y_orb = a * sqrt1me2 * sinE + Edot = n / (1.0 - e * cosE) + vx_orb = -a * sinE * Edot + vy_orb = a * sqrt1me2 * cosE * Edot + return x_orb, y_orb, vx_orb, vy_orb + + +def _thiele_innes(Omega, inc, omega): + """Thiele-Innes constants (A, B, F, G, C, H), all angles in radians. + + Mapping (Murray & Dermott 1999): + d_dec (north) = A*x + F*y + d_ra (east) = B*x + G*y + z (away from observer) = C*x + H*y + """ + cO, sO = np.cos(Omega), np.sin(Omega) + ci, si = np.cos(inc), np.sin(inc) + cw, sw = np.cos(omega), np.sin(omega) + A = cO * cw - sO * sw * ci + B = sO * cw + cO * sw * ci + F = -cO * sw - sO * cw * ci + G = -sO * sw + cO * cw * ci + C = sw * si + H = cw * si + return A, B, F, G, C, H + + +def _arcsec_per_yr_to_km_per_s(v_arcsec_per_yr: np.ndarray, R0: u.Quantity) -> u.Quantity: + """Convert angular speed at distance ``R0`` to linear km/s.""" + return (v_arcsec_per_yr * (u.arcsec / u.yr) * R0).to( + u.km / u.s, equivalencies=u.dimensionless_angles() + ) + + +def stars_at_time(time, catalogue: Path | str | None = None, R0: u.Quantity = DEFAULT_R0) -> QTable: + """Predict positions and radial velocities of all catalogue stars at ``time``. + + Parameters + ---------- + time : astropy.time.Time, datetime, ISO string, or float MJD + Epoch at which to evaluate the orbits. + catalogue : path-like, optional + Path to ``J_ApJ_837_30_table3.dat.fits``. Defaults to the file + shipped next to this module. + R0 : astropy.units.Quantity, optional + Distance to Sgr A*. Default 8.32 kpc (Gillessen+2017). Used only to + convert orbital angular velocity to line-of-sight km/s. + + Returns + ------- + astropy.table.QTable + Columns: ``Star`` (str), ``d_ra`` [arcsec, east-positive], + ``d_dec`` [arcsec, north-positive], ``RV`` [km/s, positive=receding], + ``Kmag`` [mag], ``SpT`` (str: ``e`` early / ``l`` late). + Hyperbolic / open-orbit rows (e.g. S111) are skipped with a warning. + """ + cat = _load_catalogue(catalogue if catalogue is not None else DEFAULT_CATALOGUE) + t_jyear = _to_jyear(time) + + a = np.asarray(cat["a"], dtype=float) + e = np.asarray(cat["e"], dtype=float) + per = np.asarray(cat["Per"], dtype=float) + valid = (a > 0) & (e >= 0) & (e < 1) & (per > 0) + if not valid.all(): + skipped = [ + (s.decode().strip() if isinstance(s, bytes) else str(s).strip()) + for s in np.asarray(cat["Star"])[~valid] + ] + warnings.warn( + f"Skipping hyperbolic / open-orbit rows: {skipped}", + UserWarning, + stacklevel=2, + ) + + cat = cat[valid] + a = a[valid] + e = e[valid] + per = per[valid] + inc = np.deg2rad(np.asarray(cat["i"], dtype=float)) + Omega = np.deg2rad(np.asarray(cat["Omega"], dtype=float)) + omega = np.deg2rad(np.asarray(cat["w"], dtype=float)) + Tp = np.asarray(cat["Tp"], dtype=float) + + n = 2.0 * np.pi / per # rad / Julian yr + M = n * (t_jyear - Tp) + E = _solve_kepler(M, e) + x_orb, y_orb, vx_orb, vy_orb = _orbital_state(a, e, n, E) + A, B, F, G, C, H = _thiele_innes(Omega, inc, omega) + + d_dec = A * x_orb + F * y_orb + d_ra = B * x_orb + G * y_orb + vz_arcsec_per_yr = C * vx_orb + H * vy_orb + rv = _arcsec_per_yr_to_km_per_s(vz_arcsec_per_yr, R0) + + def _decode(x): + return x.decode().strip() if isinstance(x, bytes) else str(x).strip() + + star_names = [_decode(s) for s in np.asarray(cat["Star"])] + spt = [_decode(s) for s in np.asarray(cat["SpT"])] + + out = QTable() + out["Star"] = star_names + out["d_ra"] = d_ra * u.arcsec + out["d_dec"] = d_dec * u.arcsec + out["RV"] = rv + out["Kmag"] = np.asarray(cat["Kmag"], dtype=float) * u.mag + out["SpT"] = spt + return out diff --git a/scopesim_templates/stellar/galactic_centre.py b/scopesim_templates/stellar/galactic_centre.py new file mode 100644 index 0000000..9507bb9 --- /dev/null +++ b/scopesim_templates/stellar/galactic_centre.py @@ -0,0 +1,120 @@ +# -*- coding: utf-8 -*- +"""Galactic Centre star-field template (Gillessen+2017 orbits).""" + +import warnings + +import numpy as np +from astropy import units as u +from astropy.table import Table + +import pyckles +from spextra import Spextrum + +from ..rc import Source +from ..utils.general_utils import add_function_call_str +from . import _gillessen +from . import stars_utils as su + + +__all__ = ["galactic_centre"] + + +# Sgr A* (J2000) +RA_SGR_A_STAR = 266.41684 # deg +DEC_SGR_A_STAR = -29.00781 # deg + + +@add_function_call_str +def galactic_centre(time, + filter_name="Paranal/NACO.Ks", + early_type_spec="B0V", + late_type_spec="M8III"): + """ + Source of the Sgr A* star field at a user-specified epoch. + + Stars and orbital elements from Gillessen et al. 2017 + (CDS J/ApJ/837/30, vendored). Per-star spectra are RV-shifted to the + line-of-sight velocity at ``time`` and scaled to each star's catalogue + K-band magnitude. + + Parameters + ---------- + time : astropy.time.Time, datetime, ISO string, or float MJD + Epoch at which to evaluate the orbits. + filter_name : str + SVO filter ID used to scale the spectra to ``Kmag``. Defaults to + ``"Paranal/NACO.Ks"`` to match the photometric system of the + Gillessen catalogue. + early_type_spec, late_type_spec : str + Spectral types assigned to early-type (``SpT == "e"``) and + late-type (``SpT == "l"``) stars respectively. + + Returns + ------- + scopesim.Source + + Examples + -------- + >>> from astropy.time import Time + >>> from scopesim_templates.stellar import galactic_centre + >>> src = galactic_centre(Time("2024-01-01")) + """ + tbl = _gillessen.stars_at_time(time) + n_stars = len(tbl) + + # Map SpT codes to user-supplied spectral types. + spec_types = [] + for code in tbl["SpT"]: + code_norm = str(code).strip().lower() + if code_norm == "e": + spec_types.append(early_type_spec) + elif code_norm == "l": + spec_types.append(late_type_spec) + else: + warnings.warn( + f"galactic_centre: unknown SpT code {code!r}; " + f"defaulting to early-type {early_type_spec!r}", + UserWarning, + stacklevel=2, + ) + spec_types.append(early_type_spec) + + # Base spectrum per unique user-supplied spectral type, scaled to 0 mag. + unique_user_types = [early_type_spec, late_type_spec] + pickles_lib = pyckles.SpectralLibrary("pickles", return_style="synphot") + cat_types = su.nearest_spec_type(unique_user_types, pickles_lib.table) + base = { + orig: Spextrum(modelclass=pickles_lib[cat]).scale_to_magnitude( + filter_curve=filter_name, amplitude=0 * u.mag) + for orig, cat in zip(unique_user_types, cat_types) + } + + # Per-star spectra: each gets its own RV-shifted copy of the base. + rv_kms = tbl["RV"].to(u.km / u.s).value + spectra = [base[spec_types[i]].redshift(vel=float(rv_kms[i])) + for i in range(n_stars)] + + # Weights bake in the K-band magnitude (base is at 0 mag). + weight = 10 ** (-0.4 * tbl["Kmag"].to(u.mag).value) + + # Identity ref: each star points at its own (unique) spectrum. + ref = np.arange(n_stars, dtype=int) + + field = Table( + names=["x", "y", "ref", "weight", "spec_types"], + data=[tbl["d_ra"].to(u.arcsec), + tbl["d_dec"].to(u.arcsec), + ref, + weight, + spec_types], + ) + + src = Source(spectra=spectra, table=field) + src.meta.update({ + "object": "Galactic Centre (Gillessen+2017)", + "epoch": str(time), + "ra": RA_SGR_A_STAR, + "dec": DEC_SGR_A_STAR, + "filter_name": filter_name, + }) + return src diff --git a/scopesim_templates/tests/test_stellar/test_galactic_centre.py b/scopesim_templates/tests/test_stellar/test_galactic_centre.py new file mode 100644 index 0000000..036bdb1 --- /dev/null +++ b/scopesim_templates/tests/test_stellar/test_galactic_centre.py @@ -0,0 +1,110 @@ +# -*- coding: utf-8 -*- + +import warnings + +import pytest +import numpy as np +from numpy import testing as npt +from astropy import units as u +from astropy.table import QTable +from astropy.time import Time + +from scopesim_templates.rc import Source +from scopesim_templates.stellar import galactic_centre +from scopesim_templates.stellar import _gillessen + + +EPOCH = Time("2024-01-01") + + +@pytest.mark.webtest +class TestGalacticCentre: + + def test_returns_source(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) # hyperbolic-orbit skip + src = galactic_centre(EPOCH) + assert isinstance(src, Source) + field = src.fields[0].field + assert len(field) >= 30 # 40 catalogue rows minus a few hyperbolic + + def test_spectra_one_per_star(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + src = galactic_centre(EPOCH) + n = len(src.fields[0].field) + assert len(src.spectra) == n + # ref column is the identity mapping into spectra + npt.assert_array_equal( + src.fields[0].field["ref"], np.arange(n, dtype=int) + ) + + def test_two_distinct_spectral_types_used(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + src = galactic_centre(EPOCH) + assert set(src.fields[0].field["spec_types"]) == {"B0V", "M8III"} + + def test_rv_applied_per_star(self): + """Two stars of the same class but different RV must have shifted spectra.""" + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + src = galactic_centre(EPOCH) + + field = src.fields[0].field + # Need two early-type stars with appreciably different RVs. + # We can reconstruct RVs from the gillessen output for the same epoch. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + ref_tbl = _gillessen.stars_at_time(EPOCH) + early_mask = np.array([s.strip().lower() == "e" for s in ref_tbl["SpT"]]) + early_idx = np.where(early_mask)[0] + rv = ref_tbl["RV"].to(u.km / u.s).value[early_idx] + # pick the pair with the largest RV spread + i, j = early_idx[np.argmax(rv)], early_idx[np.argmin(rv)] + sp_i, sp_j = src.spectra[int(i)], src.spectra[int(j)] + # Sample both spectra on a shared wavelength grid; nonzero difference + # proves the RV shift took effect. + wave = np.linspace(2.0e4, 2.4e4, 200) # K-band, Angstroms + flux_i = sp_i(wave) + flux_j = sp_j(wave) + assert not np.allclose(flux_i, flux_j, rtol=1e-6, atol=0.0) + + def test_weights_match_kmag(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + src = galactic_centre(EPOCH) + ref_tbl = _gillessen.stars_at_time(EPOCH) + expected = 10 ** (-0.4 * ref_tbl["Kmag"].to(u.mag).value) + npt.assert_allclose( + src.fields[0].field["weight"], expected, rtol=1e-6 + ) + + def test_positions_are_arcsec(self): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + src = galactic_centre(EPOCH) + field = src.fields[0].field + # Positions are in arcsec offsets from Sgr A* — expect |x|,|y| < a few " + x = np.asarray(field["x"]) + y = np.asarray(field["y"]) + assert np.max(np.abs(x)) < 30.0 + assert np.max(np.abs(y)) < 30.0 + + def test_unknown_spt_warns_and_defaults_to_early(self, monkeypatch): + """Inject a star with an unrecognised SpT and confirm the warning + fallback.""" + # Real run for one row, then mutate. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", UserWarning) + real_tbl = _gillessen.stars_at_time(EPOCH) + fake = QTable(real_tbl[:1]) + fake["SpT"] = ["?"] + + def fake_stars_at_time(*args, **kwargs): + return fake + + monkeypatch.setattr(_gillessen, "stars_at_time", fake_stars_at_time) + + with pytest.warns(UserWarning, match="unknown SpT"): + src = galactic_centre(EPOCH) + assert src.fields[0].field["spec_types"][0] == "B0V"