-
Notifications
You must be signed in to change notification settings - Fork 0
/
runaway_functionsv2.py
2133 lines (1812 loc) · 111 KB
/
runaway_functionsv2.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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import os
import astropy
from astropy.table import Table, Column, QTable, join,vstack
import yaml
# import pandas # Renamed for clarity
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe
from IPython.display import IFrame
from matplotlib.widgets import CheckButtons
from astroquery.vizier import Vizier
from scipy.stats import norm
from astropy.visualization import astropy_mpl_style, quantity_support
from astropy.coordinates import AltAz, EarthLocation, SkyCoord
from astropy.time import Time
import astropy.units as u
from astropy.coordinates import get_sun, get_body
from IPython.display import Math
from astropy.wcs import WCS
from astroquery.skyview import SkyView
import os
# from runaway_functions import *
# import pandas as pd
import shutil
from regions import CircleSkyRegion
from selenium import webdriver
from selenium.webdriver.chrome.service import Service
from selenium.webdriver.common.by import By
from selenium.webdriver.support.ui import Select
from selenium.webdriver.chrome.options import Options
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from astropy.io import fits
import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
import matplotlib.patheffects as path_effects
import matplotlib.patches as patches
from astropy.visualization.wcsaxes import add_scalebar
from astropy.coordinates import Angle
from psrqpy import QueryATNF
from astroquery.simbad import Simbad
import astropy.coordinates as coord
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy import units as u
from IPython.display import display, Math
import yaml
from astropy.stats import sigma_clip
matplotlib.rcParams.update({'font.size': 12})
def simbad(keyword:str):
return IFrame(f"https://simbad.cds.unistra.fr/simbad/sim-basic?Ident={keyword}+&submit=SIMBAD+search", 1400,500)
cluster_list = Table.read('Clusters_from_dias_and_a99', format='ascii.ecsv')
cluster_list_mod = Table.read('Clusters_from_dias_and_a99_mod', format='ascii.ecsv')
def read_yaml_file(file_path):
'''
Read the configuration file for the rest of the code.
This contains the various parameters for the code to run.
'''
with open(file_path, 'r') as yaml_file:
config = yaml.safe_load(yaml_file)
return config
config = read_yaml_file('config.yaml')
class Cluster:
from typing import List
def __init__(self, cluster_name,version='dr3', add_members: List[int]=[], Pmemb = config['memb_prob'], PlxQuality = config['parallax_quality_threshold']):
'''
add_members: Gaia sourceIDs of new members to be added to the cluster, only works with the default `version='dr3'`
'''
# Assuming 'cluster_list' is a predefined list containing your Astropy table
cluster_list = Table.read('Clusters_from_dias_and_a99', format='ascii.ecsv')
cluster_list_mod = Table.read('Clusters_from_dias_and_a99_mod', format='ascii.ecsv')
if cluster_list[cluster_list['Cluster'] == cluster_name][0] == cluster_list_mod[cluster_list_mod['Cluster'] == cluster_name][0]:
cluster_list = cluster_list
print('Original Data')
else:
cluster_list = cluster_list_mod
print('Modified Data')
self.cluster_table = cluster_list[cluster_list['Cluster'] == cluster_name]
self.all = self.cluster_table[0]
self.name = self.cluster_table['Cluster'][0]
ti = (Time('J2000')+1*u.Myr)
self.coordinates = SkyCoord(ra=self.cluster_table['RA_ICRS'],dec=self.cluster_table['DE_ICRS'],pm_ra_cosdec=self.cluster_table['pmRA'],pm_dec=self.cluster_table['pmDE'],obstime=ti)[0]
self.diameter = self.cluster_table['Diameter'][0]*u.arcmin
self.diameter_dias = (self.cluster_table['r50'][0]*u.deg*2).to(u.arcmin)
self.dias_distance = self.cluster_table['Dist'][0]*u.pc
self.Pmemb = Pmemb
self.PlxQuality = PlxQuality
if len(self.cluster_table) == 0:
raise ValueError(f"No data found for cluster '{cluster_name}' in the cluster list.")
if version == 'dr3':
d = self.dias_members(memb_prob=Pmemb)
dm = d[d['Gmag']<config['gmag_lim']]['Source','RAdeg','DEdeg','Gmag','e_Gmag','Plx','e_Plx','Pmemb']
sources = np.array(dm['Source'])
sources = np.insert(sources,0,np.array(add_members))
sources_int_list = sources.astype(np.int64)
sr=self.stars_in_region(no_rmRArmDE=True)
mruwe = (sr['RUWE']<config['ruwe_lim'])
madd = np.isin(sr['Source'],np.array(add_members))
sr = sr[mruwe | madd]
mask = np.isin(sr['Source'], sources_int_list)
# print(sources_int_list)
dr3dias = sr[mask]
print(f"{len(dr3dias)}/{len(d)} members matched in dr3 search region,{len(d)-len(dm)} have > {config['gmag_lim']} Gmag")
# dr3dias.add_row()
dr3dias.sort('Gmag')
mask_plx = [value >= PlxQuality for value in (dr3dias['Plx']/dr3dias['e_Plx'])]
dr3dias = dr3dias[mask_plx]
print(f"{len(dr3dias)}/{len(dr3dias)+(len(mask_plx)-np.count_nonzero(mask_plx))} members have good parallax")
# print(dr3dias['rgeo'].mean(),dr3dias['rgeo'].std())
self.all['Dist'] = dr3dias['rgeo'].mean()
self.all['e_Dist'] = dr3dias['rgeo'].std()
self.all['pmRA'] = dr3dias['pmRA'].mean()
self.all['e_pmRA'] = dr3dias['pmRA'].std()
self.all['pmDE'] = dr3dias['pmDE'].mean()
self.all['e_pmDE'] = dr3dias['pmDE'].std()
self.all['Dist'] = dr3dias['rgeo'].mean()
self.all['Plx'] = dr3dias['Plx'].mean()
self.all['e_Plx'] = dr3dias['Plx'].std()
self.all['N'] = len(dr3dias)
elif version == 'dr2':
dr3dias = self.dias_members(memb_prob=Pmemb)
mask_plx = [value >= PlxQuality for value in (dr3dias['Plx']/dr3dias['e_Plx'])]
dr3dias = dr3dias[mask_plx]
self.N = self.cluster_table['N'][0] #no. of members
elif version == 'dr2raw':
dr3dias = self.dias_members(memb_prob=0)
self.N = self.cluster_table['N'][0] #no. of members
# Extracting values from the table and setting attributes
if not os.path.exists(f'./Clusters/{self.name}'):
os.mkdir(f'./Clusters/{self.name}')
self.distance = self.all['Dist']*u.pc
self.RV = self.cluster_table['RV'][0]*u.km/u.s
self.Av = self.cluster_table['Av'][0]
self.logage = self.cluster_table['logage'][0]
self.age = 10**(self.cluster_table['logage'][0])
self.FeH = self.cluster_table['__Fe_H_'][0]
self.members = dr3dias #gets dias_members with the default settings for memb_prob and parallax_threshold_quality
#to change the memb_prob and parallax_threshold_quality filters for member selection, use something like dias_members(0.6,11)
self.N = len(dr3dias)
self.runaways = self.read_table('runaways')
self.runaways_all = self.read_table('runaways_all')
def changeParam(self,change: tuple):
param,new_value = change
_old_value = cluster_list[cluster_list['Cluster'] == self.name][0][param]
cluster_list_mod[param][cluster_list_mod['Cluster'] == self.name] = new_value
cluster_list_mod.write('Clusters_from_dias_and_a99_mod', format='ascii.ecsv', overwrite=True)
print(f'Changed {param} from {_old_value} --> {new_value}')
def restoreParam(self, param: str):
# Read the original and modified tables
cluster_list_original = Table.read('Clusters_from_dias_and_a99', format='ascii.ecsv')
cluster_list_modified = Table.read('Clusters_from_dias_and_a99_mod', format='ascii.ecsv')
# Get the name of the cluster
cluster_name = self.name
if param == 'all':
# Iterate over all columns in the original table
for column in cluster_list_original.colnames:
# Get the original value for the current column
original_value = cluster_list_original[cluster_list_original['Cluster'] == cluster_name][0][column]
# Get the current value for the current column in the modified table (for logging)
current_value = cluster_list_modified[cluster_list_modified['Cluster'] == cluster_name][0][column]
# Restore the original value in the modified table
cluster_list_modified[column][cluster_list_modified['Cluster'] == cluster_name] = original_value
if current_value != original_value:
print(f'Restored {column}: current {current_value} --> {original_value}')
else:
# Handle the case for a single parameter as before
original_value = cluster_list_original[cluster_list_original['Cluster'] == cluster_name][0][param]
current_value = cluster_list_modified[cluster_list_modified['Cluster'] == cluster_name][0][param]
cluster_list_modified[param][cluster_list_modified['Cluster'] == cluster_name] = original_value
print(f'Restored {param}: current {current_value} --> {original_value}')
# Save the modified table
cluster_list_modified.write('Clusters_from_dias_and_a99_mod', format='ascii.ecsv', overwrite=True)
def dias_members(self,memb_prob=0):
dias_members = Table.read(f'./Clusters_Dias/{self.name}.dat', format='ascii.tab')
# dias_members = Table.read(f'Clusters/{self.name}.dat', format='ascii.tab')
dias_members.remove_rows([0,1])
# Get the column names in the table
column_names = dias_members.colnames
column_names.remove('Source')
# Convert all columns from string to float, replacing non-float values with None
for column_name in column_names:
column_data = dias_members[column_name]
converted_data = []
for entry in column_data:
try:
converted_data.append(float(entry))
except ValueError:
converted_data.append(None)
dias_members[column_name] = converted_data
dias_members.sort('Source')
mask1 = [value >= memb_prob for value in dias_members['Pmemb']]
# mask2 = [value >= parallax_quality_threshold for value in (dias_members['Plx']/dias_members['e_Plx'])]
# final_mask = (np.array(mask1) & np.array(mask2))
final_mask = (np.array(mask1))
dias_members = dias_members[final_mask]
#create a filter to get rid of the stars which don't have a Bp-Rp value in the .dat table
mask_bprp = [value is not None for value in dias_members["BP-RP"] ]
plottable_members = dias_members[mask_bprp]
print(len(dias_members)-len(plottable_members),"memebrs do not have BP-RP from Dias.")
return plottable_members
def clean(self,what='except_downloads'):
print(f"Deleting: {what}")
folder_path = f'./Clusters/{self.name}'
# List all files in the folder
files = os.listdir(folder_path)
# Iterate over each file and delete it
for file_name in files:
file_path = os.path.join(folder_path, file_name)
if what=='everything' and os.path.isfile(file_path):
os.remove(file_path)
print(f'{file_path} removed')
if what=='runaways' and os.path.isfile(file_path) and 'runaway' in file_path:
print(f'{file_path} removed')
os.remove(file_path)
if what=='except_downloads' and os.path.isfile(file_path) and not (('compare_data' in file_path) or ('stars_from' in file_path) or ('extra' in file_path)):
print(f'{file_path} removed')
os.remove(file_path)
def calculate_search_arcmin(self, extra=config['Cluster']['search_extent'], output=False):
"""
Calculate the search arcminute for the cluster.
Parameters:
extra (float, optional): Extent of search around cluster in pc. Default is taken from config file `search_extent`.
output (bool, optional): Whether to display the output of some useful parameters. Default is False.
Returns:
Quantity: The search arcminute value.
"""
theta = self.diameter / 2 # radius of the cluster in arcmin
D = self.dias_distance # Assuming you have a 'distance' attribute in your Cluster class
r = np.tan(theta) * D
# display(r)
search_arcmin = np.arctan((r + extra * u.pc) / D)
search_arcmin = search_arcmin.to(u.arcminute)
r, search_arcmin = r.round(3), search_arcmin.round(3)
if output:
display(Math(r'\mathrm{Dist.}' + f'= {D}' + r'\mathrm{\ Ang. \ Radius}' + f'= {theta}' +
r'\mathrm{\ Phy. \ Radius}' + f'= {r}'))
display(Math(r'\mathrm{Extra \ } '+f'= {extra} pc'))
display(Math(r'\mathrm{\ search \ arcmin}' + f'= {search_arcmin}'))
return search_arcmin
def plot_search_region(self, extra=config['Cluster']['search_extent'],display=True,**kwargs):
"""
Plots and saves the fits file for the region around the given cluster.
Parameters:
- extra (float): Additional extent for the search region (default is from config['Cluster']['search_extent']).
Returns:
None
"""
search_arcmin = self.calculate_search_arcmin(extra=extra)
# Define the file path
fits_file_path = f'./Clusters/{self.name}/{self.name}_extra{extra}pc.fits'
# Check if the file already exists
if os.path.exists(fits_file_path):
print(f'fits image exists in {self.name} folder')
# File exists, no need to download, use the existing file
images = [fits.open(fits_file_path)]
# Extract the WCS information
wcs = WCS(images[0][0].header)
else:
# File doesn't exist, get the image data from SkyView
images = SkyView.get_images(position=self.coordinates,
survey=config['Plots']['skyview_survey'],
radius=2*search_arcmin,
**kwargs)
# Extract the WCS information
wcs = WCS(images[0][0].header)
hdu = fits.PrimaryHDU(data=images[0][0].data, header=images[0][0].header)
hdulist = fits.HDUList([hdu])
# Save the fits file
hdulist.writeto(fits_file_path, overwrite=True)
if not display:
return None
# Plot the image
print(f'Plotting a {extra} pc region around cluster.')
# print(wcs.wcs.ctype) #projection is TAN by default for skyview images
fig, ax = plt.subplots(figsize=(6,6), subplot_kw={'projection': wcs})
ax.imshow(images[0][0].data, cmap='gray')
data = images[0][0].data
# Add labels and title
ax.set_xlabel('Right Ascension (degrees)')
ax.set_ylabel('Declination (degrees)')
ax.set_title(f"{self.name} with {extra}pc search region")
ax.grid(color='lightgrey',ls='dotted')
# Add a circle representing the cluster radius
circle_cluster = plt.Circle((data.shape[1] / 2, data.shape[0] / 2), radius=(data.shape[0] / 2)*self.diameter/(2*search_arcmin),
edgecolor='red', facecolor='none',ls='dashed',label=f'Cluster Diameter = {self.diameter}')
circle_search_region = plt.Circle((data.shape[1] / 2, data.shape[0] / 2), radius=(data.shape[0] / 2)*search_arcmin/search_arcmin,
edgecolor='green', facecolor='none',ls='dashed',label=f'Search Region Diameter = {2*search_arcmin}')
ax.add_artist(circle_cluster)
ax.add_artist(circle_search_region)
_ = search_arcmin.round(-1)/2 #round to the nearest 5 (by round first to nearest 10 and then divide by 2
scalebar_length = ((_.to(u.rad))*(self.distance.to(u.m)).to(u.pc)).round(2)
__ = self.distance.value/1000
add_scalebar(ax, _,color="yellow",label=f"or {scalebar_length.value}pc (at dist {__:.2f}kpc)",size_vertical=0.5)
x_min, x_max = ax.get_xlim()
y_min, y_max = ax.get_ylim()
ax.annotate(f"{_:.2f}", xy=(0.83*x_max,0.08*y_max), color='yellow', ha='center',fontweight='bold')
ax.legend()
images[0].close() # Close the opened FITS file
# Show the plot
plt.show()
def stars_in_region(self,allTables=False,no_rmRArmDE=False):
"""
Finds the stars around the cluster within a conical frustum.
Saves the three tables in the cluster folder.
Qurries VizieR only if tables are not already present in the cluster folder.
Parameters:
allTables (bool, optional): Default is false, and returns 1 argument - stars_in_region
which is the stars_fromDR3 and stars_fromDR3_dis joined.
Returns:
stars_in_region (astropy table)
or (stars_in_region, stars_fromDR3, stars_fromDR3_dis) if allTables is True
"""
if not os.path.exists(f'./Clusters/{self.name}'):
os.mkdir(f'./Clusters/{self.name}')
stars_fromDR3_path = f'./Clusters/{self.name}/{self.name}_stars_fromDR3.tsv'
stars_fromDR3_dis_path = f'./Clusters/{self.name}/{self.name}_stars_fromDR3_dis.tsv'
stars_in_region_path = f'./Clusters/{self.name}/{self.name}_stars_in_region.tsv'
fs_path = f'./Clusters/{self.name}/{self.name}_fs.tsv'
nsfs_path = f'./Clusters/{self.name}/{self.name}_nsfs.tsv'
if allTables and (os.path.exists(stars_in_region_path)) and os.path.exists(stars_fromDR3_path) and os.path.exists(stars_fromDR3_dis_path) and os.path.exists(fs_path) and os.path.exists(nsfs_path):
#checks for generate_tables() if all the tables are already present.
stars_in_region = Table.read(stars_in_region_path,format='ascii.ecsv')
stars_fromDR3 = Table.read(stars_fromDR3_path,format='ascii.ecsv')
stars_fromDR3_dis = Table.read(stars_fromDR3_dis_path,format='ascii.ecsv')
fs = Table.read(fs_path,format='ascii.ecsv')
nsfs = Table.read(nsfs_path,format='ascii.ecsv')
print("All tables present (stars_in_region,stars_fromDR3,stars_fromDR3_dis,fs,nsfs)")
return stars_in_region,stars_fromDR3,stars_fromDR3_dis,fs,nsfs
if os.path.exists(stars_in_region_path) and (not allTables):
stars_in_region = Table.read(stars_in_region_path,format='ascii.ecsv')
# print(f'{stars_in_region_path} exists in {self.name} with {len(stars_in_region)} stars')
return stars_in_region
if os.path.exists(stars_fromDR3_path):
stars_fromDR3 = Table.read(stars_fromDR3_path,format='ascii.ecsv')
print(f'Table stars_fromDR3 exists in {self.name} with {len(stars_fromDR3)} stars, skipping VizieR query')
else:
print(f'Querying VizieR {config["cat_stars"]} \n at: {self.coordinates.ra}, {self.coordinates.dec} \n within: {self.calculate_search_arcmin()}')
print(f'Filters: \n {config["cat_stars_filters"]}')
stars_fromDR3 = Vizier(columns=["*","+_r"],row_limit = -1).query_region(self.coordinates,
radius=self.calculate_search_arcmin(),
catalog=config['cat_stars'],
column_filters=config['cat_stars_filters'])[0]
print(f'{len(stars_fromDR3)} sources found')
if os.path.exists(stars_fromDR3_dis_path):
stars_fromDR3_dis = Table.read(stars_fromDR3_dis_path,format='ascii.ecsv')
print(f'Table stars_fromDR3_dis exists in {self.name} with {len(stars_fromDR3_dis)} stars, skipping VizieR query')
else:
print(f'Querying VizieR {config["cat_star_distances"]} \n at: {self.coordinates.ra}, {self.coordinates.dec} \n within: {self.calculate_search_arcmin()}')
_l = (1-config['distance_tolerance'])*self.dias_distance.value
_u = (1+config['distance_tolerance'])*self.dias_distance.value
_filter = {'rgeo':f'>{_l} && <{_u}'}
print(f'Filters: \n {_filter}')
stars_fromDR3_dis = Vizier(columns=["*","+_r"],row_limit = -1).query_region(self.coordinates,
radius=self.calculate_search_arcmin(),
catalog=config['cat_star_distances'],
column_filters=_filter)[0]
print(f'{len(stars_fromDR3_dis)} sources found')
stars_in_region = join(stars_fromDR3,stars_fromDR3_dis,keys='Source',join_type='inner')
# relative to the cluster stuff
if not no_rmRArmDE:
dias_members= self.members
rmRA = stars_in_region['pmRA']-dias_members['pmRA'].mean()
rmDE = stars_in_region['pmDE']-dias_members['pmDE'].mean()
e_rmRA = stars_in_region['e_pmRA']+np.sqrt((dias_members['e_pmRA']**2).mean())
e_rmDE = stars_in_region['e_pmDE']+np.sqrt((dias_members['e_pmDE']**2).mean())
e_rmRA = stars_in_region['e_pmRA']+self.all['e_pmRA']
e_rmDE = stars_in_region['e_pmDE']+self.all['e_pmDE']
µ_pec = np.sqrt(rmRA**2+rmDE**2)
D = stars_in_region['rgeo']/1000
stars_in_region.add_column(µ_pec, name='µ_pec', index=0)
stars_in_region.add_column(µ_pec*D*4.74,name='v_pec',index=0)
stars_in_region.add_column(rmRA,name='rmRA',index=5)
stars_in_region.add_column(e_rmRA,name='e_rmRA',index=6)
stars_in_region.add_column(rmDE,name='rmDE',index=7)
stars_in_region.add_column(e_rmDE,name='e_rmDE',index=8)
stars_in_region['v_pec'].unit = u.km / u.s
stars_in_region.sort('Gmag')
stars_in_region = stars_in_region['RA_ICRS_1','DE_ICRS_1','e_RA_ICRS','e_DE_ICRS','_r_1',
'HIP','TYC2','Source','rgeo','Plx','e_Plx',
'v_pec','µ_pec','rmRA','e_rmRA','rmDE','e_rmDE','pmRA','pmDE','e_pmRA','e_pmDE',
'RUWE','Teff','logg','Gmag','BP-RP','BPmag','RPmag','RV','e_RV',
'b_rgeo','B_rgeo','FG','e_FG','FBP','e_FBP','FRP','e_FRP','RAVE5','RAVE6']
# Adding row for uncertainties in Gmag and BPmag and RPmag
# values for Gaia G, G_BP, G_RP zero point uncertainties
sigmaG_0 = 0.0027553202
sigmaGBP_0 = 0.0027901700
sigmaGRP_0 = 0.0037793818
stars_in_region['e_Gmag'] = np.sqrt((-2.5/np.log(10)*stars_in_region['e_FG']/stars_in_region['FG'])**2 + sigmaG_0**2)
stars_in_region['e_BPmag'] = np.sqrt((-2.5/np.log(10)*stars_in_region['e_FBP']/stars_in_region['FBP'])**2 + sigmaGBP_0**2)
stars_in_region['e_RPmag'] = np.sqrt((-2.5/np.log(10)*stars_in_region['e_FRP']/stars_in_region['FRP'])**2 + sigmaGRP_0**2)
stars_in_region['e_BP-RP'] = stars_in_region['e_BPmag']+stars_in_region['e_RPmag']
fs = stars_in_region[stars_in_region['v_pec'] >= config['v_runaway']]
nsfs = stars_in_region[np.array(stars_in_region['v_pec'] >= config['v_walkaway']) & np.array(stars_in_region['v_pec'] < config['v_runaway'])]
#save the tables
fs.write(f'./Clusters/{self.name}/{self.name}_fs.tsv',format='ascii.ecsv',overwrite=True)
nsfs.write(f'./Clusters/{self.name}/{self.name}_nsfs.tsv',format='ascii.ecsv',overwrite=True)
stars_in_region.write(f'./Clusters/{self.name}/{self.name}_stars_in_region.tsv',format='ascii.ecsv',overwrite=True)
stars_fromDR3.write(f'./Clusters/{self.name}/{self.name}_stars_fromDR3.tsv',format='ascii.ecsv',overwrite=True)
stars_fromDR3_dis.write(f'./Clusters/{self.name}/{self.name}_stars_fromDR3_dis.tsv',format='ascii.ecsv',overwrite=True)
if allTables and not no_rmRArmDE:
print(f'{len(stars_in_region)} stars in the region')
return stars_in_region,stars_fromDR3,stars_fromDR3_dis,fs,nsfs
else:
print(f'{len(stars_in_region)} stars in the region')
return stars_in_region
# def refine_stars_in_region(stars_in_region,v_dispersion=config['v_dispersion'],sigma_clip=config['sigma_clip']):
# mask_vdisp = [value < v_dispersion for value in stars_in_region['v_pec']]
def generate_tables(self):
self.stars_in_region(allTables=True)
def read_table(self, *tables):
"""
Read one or more tables associated with the cluster.
See the particular cluster folder for the tables present.
Run genetate_tables() to generate the tables.
Parameters:
- tables (str): Variable number of table names to be read. The name of the table is the text after the last "_".
For eg., "NGC_4103_psrs.tsv" is called "psrs"
Returns:
- table (astropy.table.table.Table or list): If a single table is provided, returns the corresponding table.
If multiple tables are provided, returns a list of tables.
If any of the specified tables do not exist, prints a message and skips them.
Example:
cluster = Cluster('YourClusterName')
member_table = cluster.read_table('members')
stars_in_region,stars_fromDR3,stars_fromDR3_dis,fs,nsfs = cluster.read_table('stars_in_region','stars_fromDR3','stars_fromDR3_dis','fs','nsfs')
"""
tl = []
for table in tables:
file_path = f'./Clusters/{self.name}/{self.name}_{table}.tsv'
# Check if the file exists before trying to read it
if os.path.exists(file_path):
tl.append(Table.read(file_path, format='ascii.ecsv'))
else:
print(f"File {file_path} not present.")
return tl[0] if len(tl) == 1 else tl
def plot_theoretical_isochrone(self):
fig = plt.figure(figsize=(8,8))
BP_RP_theo,Gmag_theo = theoretical_isochrone(self)
ax = fig.add_subplot()
ax.plot(BP_RP_theo,Gmag_theo,label='Theoretical Isochrone')
ax.set_xlabel(r"$G_{BP}-G_{RP}$ (mag)")
ax.set_ylabel(r"$G$ (mag)")
ax.set_title(f"CMD for {self.name}")
ax.invert_yaxis()
ax.legend()
def latex_text(self):
ra, dec = self.coordinates.ra, self.coordinates.dec
ra_str, dec_str = ra.to_string(format='latex')[1:-1], dec.to_string(format='latex')[1:-1]
dist_str = str(self.all['Dist'])+r"\pm"+str(self.all['e_Dist'])
print(rf"{self.name.replace('_', ' ')} is located at $\alpha = {ra_str}, \delta = {dec_str}$ at a distance of ${dist_str}$ pc.")
def latex_table_kinematics(self):
t_dict = {
# 'Name':[self.name.replace('_',' ')],
r'$\alpha (^\circ)$': [f"{(self.coordinates.ra):.2f}".replace(' deg','')],
r'$\delta (^\circ)$': [f"{(self.coordinates.dec):.2f}".replace(' deg','')],
"Diam. (')":[self.diameter.value],
"N":[len(self.members)],
r"$\mu_{\alpha}^*$ (mas/yr)":[f"{self.all['pmRA']:.2f}"+"$\pm$"+f"{self.all['e_pmRA']:.2f}"],
r"$\mu_{\delta}$ (mas/yr)":[f"{self.all['pmDE']:.2f}"+"$\pm$"+f"{self.all['e_pmDE']:.2f}"],
r"$v_R$ (km/s)":[f"{self.all['RV']:.2f}"+"$\pm$"+f"{self.all['e_RV']:.2f}"] if not isinstance(self.all['RV'],np.ma.core.MaskedConstant) else ['N/A'],
r"$d$ (pc)":[f"{self.all['Dist']:.0f}"+"$\pm$"+f"{self.all['e_Dist']:.0f}"],
}
latexdict={'tabletype':'table*','tablealign':'h',
'header_start':r'\hline','data_start':r'\hline','data_end': r'\hline',
'caption':f'Kinematic parameters of {self.name.replace("_"," ")}',
'preamble':'\label{tab:'+f'{self.name}-kinematics'+'}'}
latexdict['tablehead'] = r'head'
# Convert the dictionary to an Astropy table
astropy_table = Table(t_dict)
astropy.io.ascii.write(astropy_table, format='latex',output='text.txt',overwrite=True,
latexdict=latexdict)
with open('text.txt', 'r+') as f:
lines = f.readlines()
lines[1], lines[2] = lines[2], lines[1]
f.seek(0)
f.writelines(lines)
print(''.join(lines))
os.remove('text.txt')
def latex_table_members(self,n=None):
"""
Prints a LaTeX table for the cluster members of the cluster object. The members are the same as what you
get from `cluster.members`. Table is sorted according to the Gmag of the stars.
Parameters:
- n (int): number of cluster members to be returned
Example usage:
cluster = Cluster('IC_2395')
latex_table_members(cluster,n=10)
"""
table = self.members
table.sort('Gmag')
table = table[:n]
# Create an empty list to store the formatted distances
formatted_distances = []
formatted_parallaxes = []
formatted_gmags = []
formatted_bprps = []
# Iterate over each row in the table
for row in table:
distance = row['rgeo']
plx = row['Plx']
e_plx = row['e_Plx']
upper_error = row['B_rgeo']-row['rgeo']
lower_error = -row['b_rgeo']+row['rgeo']
gmag = row['Gmag']
e_gmag = row['e_Gmag']
bprp = row['BP-RP']
e_bprp = row['e_BP-RP']
# Format strings
formatted_distance = f"${distance:.0f}^{{+{upper_error:.0f}}}_{{-{lower_error:.0f}}}$"
formatted_parallax = f"${plx:.4f}\pm{e_plx:.4f}$"
formatted_gmag = f"${gmag:.3f}\pm{e_gmag:.3f}$"
formatted_bprp = f"${bprp:.3f}\pm{e_bprp:.3f}$"
# Append the formatted distance to the list
formatted_distances.append(formatted_distance)
formatted_parallaxes.append(formatted_parallax)
formatted_gmags.append(formatted_gmag)
formatted_bprps.append(formatted_bprp)
# Add a new column to the table with the formatted distances
table['formatted_distance'] = formatted_distances
table['formatted_parallax'] = formatted_parallaxes
table['formatted_gmag'] = formatted_gmags
table['formatted_bprp'] = formatted_bprps
t_dict = {
'Gaia DR3 Source':table['Source'],
r"$r_{\text{geo}}$ (pc)":table['formatted_distance'],
r"$\pi$ (mas)":table['formatted_parallax'],
r"$G$ (mag)":table['formatted_gmag'],
r"$G_{\text{BP}}-G_{\text{RP}}$ (mag)":table['formatted_bprp']
# r"$v_R$ (km/s)":[f"{self.all['RV']:.2f}"+"$\pm$"+f"{self.all['e_RV']:.2f}"] if not isinstance(self.all['RV'],np.ma.core.MaskedConstant) else ['N/A'],
# r"$d$ (pc)":[f"{self.all['Dist']:.0f}"+"$\pm$"+f"{self.all['e_Dist']:.0f}"],
}
latexdict={'tabletype':'table*','tablealign':'h',
'header_start':r'\hline','data_start':r'\hline','data_end': r'\hline',
'caption':f'Selected members of {self.name.replace("_"," ")}',
'preamble':'\label{tab:'+f'{self.name}-members'+'}'}
latexdict['tablehead'] = r'head'
# Convert the dictionary to an Astropy table
astropy_table = Table(t_dict)
astropy.io.ascii.write(astropy_table, format='latex',output='text.txt',overwrite=True,
latexdict=latexdict)
with open('text.txt', 'r+') as f:
lines = f.readlines()
lines[1], lines[2] = lines[2], lines[1]
f.seek(0)
f.writelines(lines)
print(''.join(lines))
os.remove('text.txt')
def find_cluster(stars_in_region,refCluster,sigma = config['sigma_clip'],plx_quality = config['plx_quality']):
"""
Takes in the stars in a region and sigma_clip their pmRA and pmDE separately to keep values within
a specified sigma. Then further filters stars which have good parallax quality (default>10) and are
within the cluster radius (from dias).
Example Usage:
cluster = Cluster('Berkeley_97')
my_cluster = find_cluster(cluster.stars_in_region(),refCluster=cluster.name)
"""
within_r = Cluster(refCluster).diameter.value/2
ms = stars_in_region[np.array(~sigma_clip(stars_in_region['pmRA'],sigma=sigma).mask)&np.array(~sigma_clip(stars_in_region['pmDE'],sigma=sigma).mask)]
my_stars = ms[(np.array(ms['_r_1']<within_r) & np.array(ms['Plx']/ms['e_Plx']>plx_quality))]
print('mean(v_pec): ',my_stars['v_pec'].mean())
print('std(v_pec): ',my_stars['v_pec'].std())
print('v_pec_max: ',my_stars['v_pec'].max())
print(f"{len(my_stars)} kinematic members")
return my_stars
def cluster_v_dispersion(cluster_name,memb_prob=config['memb_prob'],parallax_quality_threshold=config['parallax_quality_threshold'],plot=False):
"""
Calculate the peculiar velocities of cluster members and analyze velocity dispersion.
Parameters:
- cluster_name (str): Name of the cluster.
- memb_prob (float): Membership probability threshold for selecting cluster members.
- parallax_quality_threshold (float): Parallax quality threshold for selecting cluster members.
- plot (bool): If True, plot a histogram with a bell curve representing the velocity dispersion.
Returns:
- v_pec (numpy.ndarray): Array of peculiar velocities for cluster members.
- mean_velocity (float): Mean peculiar velocity of the cluster members.
- std_deviation (float): Standard deviation of peculiar velocities.
Example:
v_pec, mean_velocity, std_deviation = cluster_v_dispersion('Berkeley_97', memb_prob=0.8, plot=True)
"""
cluster = Cluster(cluster_name)
cluster_members = cluster.members
µ_pec = np.sqrt(cluster_members['pmRA']**2+cluster_members['pmRA']**2)
v_pec = µ_pec.value*cluster.distance.value*4.74*(u.km/u.s)/1000
mean_velocity = np.mean(v_pec)
std_deviation = np.std(v_pec)
if plot:
# Create a range of values for x-axis
x = np.linspace(min(v_pec), max(v_pec), 100)
# Generate the bell curve using the mean and standard deviation
y = norm.pdf(x, mean_velocity, std_deviation)
# Plot the histogram of velocities
plt.hist(v_pec, bins=10, density=True, alpha=0.6, color='g')
# Plot the bell curve
plt.plot(x, y, 'r--', label=f'Mean: {mean_velocity:.2f}, Std Dev: {std_deviation:.2f}')
# Add labels and legend
plt.xlabel('Velocity')
plt.ylabel('Probability Density')
plt.legend()
plt.title(f'{len(v_pec)} members')
# Show the plot
plt.show()
return v_pec, mean_velocity, std_deviation
def columnsplot(columnx,columny, xlabel=None, ylabel=None):
from astropy.visualization import quantity_support
"""
Plot two Astropy columns on the x and y axes.
Parameters:
columnx (astropy.table.Column): Astropy column for the x-axis.
columny (astropy.table.Column): Astropy column for the y-axis.
xlabel (str, optional): Custom label for the x-axis. Defaults to None.
ylabel (str, optional): Custom label for the y-axis. Defaults to None.
"""
with quantity_support():
fig, ax = plt.subplots()
ax.plot(columnx, columny, 'o')
ax.grid(True)
if xlabel is None:
xlabel = f"{columnx.name} ({columnx.unit})"
if ylabel is None:
ylabel = f"{columny.name} ({columny.unit})"
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
plt.show()
def theoretical_isochrone(cluster,Av=None,logage=None,FeH=None,output=None, printing = True):
metallicity = FeH
"""
Retrieves the theoretical isochrone for a given cluster from CMD interface.
By default takes the parameters of the cluster as input.
Checks if data file is already present. If not, gets it.
## Parameters:
- cluster (Cluster class object): For values of extinction, age, metallicity.
- Av (float, optional): Extinction value. Defaults to None, which uses cluster.all['Av'].
- logage (float, optional): Logarithm of the age. Defaults to None, which uses cluster.all['logage'].
- metallicity (float, optional): Metallicity value. Defaults to None, which uses cluster.all['__Fe_H_'].
- output (str, optional): 'table' gets the entire table. Default is None, which returns Bp-Rp, Gmag respectively as a tuple of astropy columns.
## Returns:
tuple or astropy.table.Table: If output is None, returns a tuple containing Bp-Rp and Gmag.
If output is 'table', returns the entire table.
## Example Usage:
BP_RPtheo, Gmagtheo = theoretical_isochrone(cluster, Av=3.1, logage=6, metallicity=-0.1)
or, for a complete table,
theo_isochrone = theoretical_isochrone(cluster, Av=3.1, logage=6, metallicity=-0.1,output='table')
Also generates a file {cluster.name}\{cluster.name}_compare_data_out_Av3.1_age6.9_FeH-0.1.dat
"""
# Accessing cluster parameters
if Av is None:
Av = cluster.all['Av']
else:
Av = float(Av)
if logage is None:
logage = cluster.all['logage']
else:
logage = float(logage)
if metallicity is None:
metallicity = cluster.all['__Fe_H_']
else:
metallicity = float(metallicity)
#Check if parameters are unchanged
if (str(Av) == str(cluster.all['Av'])) and (str(logage) == str(cluster.all['logage'])) and (str(metallicity) == str(cluster.all['__Fe_H_'])):
compate_data_out_path = os.path.join(f"./Clusters/{cluster.name}",(cluster.name+'_compare_data_out.dat'))
if printing:
print(f'{cluster.name}')
print(f'Av: {str(Av)}')
print(f'logage: {str(logage)}')
print(f'[Fe/H]: {str(metallicity)}')
#otherwise, make a new file to save the new isochrone
else:
compate_data_out_path = os.path.join(f'./Clusters/{cluster.name}',(cluster.name+f'_compare_data_out_Av{str(Av)}_age{str(logage)}_FeH{str(metallicity)}.dat'))
if printing:
print(f'{cluster.name} (Changed)')
print(f"Av: {str(cluster.all['Av'])} --> {str(Av)}")
print(f"logage: {str(cluster.all['logage'])} --> {str(logage)}")
print(f"[Fe/H]: {str(cluster.all['__Fe_H_'])} --> {str(metallicity)}")
if not os.path.exists(compate_data_out_path):
if printing:
print('Getting theoretical isochrone')
s = Service()
print(f"getting isochrone form cmd3.7 with Av:{Av:.2f}, logage:{logage:.2f}, metallicity:{metallicity:.2f}")
options = webdriver.ChromeOptions()
browser = webdriver.Chrome(service=s, options=options)
browser.get('http://stev.oapd.inaf.it/cgi-bin/cmd')
#Evolutionary Tracks #from config
browser.find_element(By.XPATH,"/html/body/form/div/fieldset[1]/table/tbody/tr[3]/td[1]/input[1]").click() #PARSEC version 2.0
browser.find_element(By.XPATH,"/html/body/form/div/fieldset[1]/table/tbody/tr[5]/td/input").click() #+ COLIBRI S_37
#Phtotometric System #from config
photometricSystem = Select(browser.find_element(By.XPATH,"//select[@name='photsys_file']")) #dropdown list for available photometric systems
photometricSystem.select_by_value("YBC_tab_mag_odfnew/tab_mag_gaiaEDR3.dat") # Gaia EDR3 bands
browser.find_element(By.XPATH,"/html/body/form/div/fieldset[2]/table/tbody/tr[5]/td[1]/input").click() # VBC +new Vega for PARSEC 2.0 #As above, but adopting revised SED for Vega from Bohlin et al. (2020) (namely CALSPEC alpha_lyr_stis_010.fits).
#Phtotometric System #from config
#Circumstellar Dust
browser.find_element(By.XPATH,"/html/body/form/div/fieldset[3]/font/table/tbody/tr[3]/td[1]/input").click() #for M stars: No dust
browser.find_element(By.XPATH,"/html/body/form/div/fieldset[3]/font/table/tbody/tr[3]/td[2]/input").click() #for C stars: No dust
#Interstellar extinction
Av_field = browser.find_element(By.XPATH,"/html/body/form/div/fieldset[4]/input")
Av_field.clear()
Av_field.send_keys(str(Av))
#Long Period Variability #from config
browser.find_element(By.XPATH,"/html/body/form/div/fieldset[5]/table/tbody/tr[4]/td[1]/input").click() #3. Periods from Trabucchi et al. (2021).
#Initial Mass Function #from config
InitialMassFunction = Select(browser.find_element(By.XPATH,"//select[@name='imf_file']"))
InitialMassFunction.select_by_value("tab_imf/imf_kroupa_orig.dat")
#ages
browser.find_element(By.XPATH,"/html/body/form/div/fieldset[7]/table/tbody/tr[4]/td[1]/input").click() #click log(age/yr)
initialAge_field = browser.find_element(By.XPATH,"/html/body/form/div/fieldset[7]/table/tbody/tr[4]/td[2]/input")
initialAge_field.clear()
initialAge_field.send_keys(str(logage))
finalAge_field = browser.find_element(By.XPATH,"/html/body/form/div/fieldset[7]/table/tbody/tr[4]/td[3]/input")
finalAge_field.clear()
finalAge_field.send_keys(str(logage))
step_field = browser.find_element(By.XPATH,"/html/body/form/div/fieldset[7]/table/tbody/tr[4]/td[4]/input")
step_field.clear()
step_field.send_keys(0)
#metallicities
selectMbyH = browser.find_element(By.XPATH,"/html/body/form/div/fieldset[7]/table/tbody/tr[7]/td[1]/input").click() #click [M/H]
initialMetallicity_field = browser.find_element(By.XPATH,"/html/body/form/div/fieldset[7]/table/tbody/tr[7]/td[2]/input")
initialMetallicity_field.clear()
initialMetallicity_field.send_keys(str(metallicity))
finalMetallicity_field = browser.find_element(By.XPATH,"/html/body/form/div/fieldset[7]/table/tbody/tr[7]/td[3]/input")
finalMetallicity_field.clear()
finalMetallicity_field.send_keys(str(metallicity))
#Submit #dosen't change
browser.find_element(By.XPATH,"/html/body/form/div/input[4]").click()
browser.find_element(By.XPATH,"/html/body/form/fieldset[1]/p[1]/a").click()
data_out = browser.find_element(By.XPATH,'/html/body/pre').text
if printing:
print('Obtained isochrone from CMD')
#save the table
with open(compate_data_out_path, 'w') as f:
f.write(data_out)
if printing:
print(f"New isochrone file: {compate_data_out_path}")
#close the browser
browser.close()
else:
if printing:
print('Theoretical isochrone exists, reading it.')
# Read the DAT file into an Astropy table
theoretical_data = Table.read(compate_data_out_path, format='ascii')
# Define the new column names
new_column_names = ['Zini', 'MH', 'logAge', 'Mini', 'int_IMF', 'Mass', 'logL', 'logTe', 'logg', 'label', 'McoreTP', 'C_O', 'period0', 'period1', 'period2', 'period3', 'period4', 'pmode', 'Mloss', 'tau1m', 'X', 'Y', 'Xc', 'Xn', 'Xo', 'Cexcess', 'Z', 'Teff0', 'omega', 'angvel', 'vtaneq', 'angmom', 'Rpol', 'Req', 'mbolmag', 'G_fSBmag', 'G_BP_fSBmag', 'G_RP_fSBmag', 'G_fSB', 'G_f0', 'G_fk', 'G_i00', 'G_i05', 'G_i10', 'G_i15', 'G_i20', 'G_i25', 'G_i30', 'G_i35', 'G_i40', 'G_i45', 'G_i50', 'G_i55', 'G_i60', 'G_i65', 'G_i70', 'G_i75', 'G_i80', 'G_i85', 'G_i90', 'G_BP_fSB', 'G_BP_f0', 'G_BP_fk', 'G_BP_i00', 'G_BP_i05', 'G_BP_i10', 'G_BP_i15', 'G_BP_i20', 'G_BP_i25', 'G_BP_i30', 'G_BP_i35', 'G_BP_i40', 'G_BP_i45', 'G_BP_i50', 'G_BP_i55', 'G_BP_i60', 'G_BP_i65', 'G_BP_i70', 'G_BP_i75', 'G_BP_i80', 'G_BP_i85', 'G_BP_i90', 'G_RP_fSB', 'G_RP_f0', 'G_RP_fk', 'G_RP_i00', 'G_RP_i05', 'G_RP_i10', 'G_RP_i15', 'G_RP_i20', 'G_RP_i25', 'G_RP_i30', 'G_RP_i35', 'G_RP_i40', 'G_RP_i45', 'G_RP_i50', 'G_RP_i55', 'G_RP_i60', 'G_RP_i65', 'G_RP_i70', 'G_RP_i75', 'G_RP_i80', 'G_RP_i85', 'G_RP_i90']
# Rename the columns
for old_name, new_name in zip(theoretical_data.colnames, new_column_names):
theoretical_data.rename_column(old_name, new_name)
# Calculate BP-RP and add column at the end
theoretical_data['BP-RP'] = theoretical_data['G_BP_fSBmag']-theoretical_data['G_RP_fSBmag']
theoretical_data['Gmag'] = theoretical_data['G_fSBmag'] + 5*np.log10(cluster.distance.value)-5
BP_RPtheo = theoretical_data['BP-RP'] #BP-RP
Gmagtheo = theoretical_data['Gmag']#-cluster['Av']
if output == None:
return BP_RPtheo,Gmagtheo
if output == 'table':
return theoretical_data
def get_runaways(cluster,fs,theoretical_data,separation_factor=2,dist_filter_factor=1):
"""
Selects stars from the input astropy table fs which are runaway star candidates by tracing them back to the cluster.
If the stars happens to be in the cluster in the past 100kyr then it is selected to be a candidate.
Then the temperature of the star is estimated by comparing its color (Bp-Rp magnitude) with the theoretical isochrone that is provided.
Parameters:
- cluster (Cluster class object): The cluster compared whose runaways are to be found.
- fs or nsfs (astropy table): A table of the fast stars (relative to the cluster mean proper motion) in the region (`search_arcmin`)
around the cluster centre. Can also take `nsfs` and trace them back with the same algorithm.
- theoretical_data (astropy table ascii format downloaded from CMD): theoretical isochrone to estimate the temperatures of the runaways (or walkaways)
- separation_factor: default 2 (since separation<cluster diameter/2 implies within cluster in the past). indicates what is the minimum separation factor to be considered inside cluster.
Returns:
- run (astropy table): A table of the selected stars out of fs which can be traced back within 100kyr or less
with their temperature estimates added in a column.
Example Usage:
cluster = Cluster('Berkeley_97')
theoretical_data = theoretical_isochrone(cluster,output="table",printing=False,logage=7)
fs = cluster.read_table('fs')
runaways = get_runaways(cluster,fs,theoretical_data)
display(runaways)
"""
if fs['v_pec'].max() < 17.6:
name = "walkaways"
print("Stars given in the fs variable have their maximum v_pec <17.6km/s. These are walkaways.")
else:
name = "runaways"
file_path = f'./Clusters/{cluster.name}/{cluster.name}_{name}_all.tsv' #`name` is either "walkaways" or "runaways"
if os.path.exists(file_path):
# If the file exists, read it and return its contents
print(f'{file_path} exists. Skipping traceback and temp estimates.')
run = Table.read(file_path, format='ascii.ecsv')
print(f'{len(fs)} stars had been checked, {len(run)} traced back to cluster')
return run
else:
#do the tracebacks
#def the proper motin of cluster
times = np.linspace(-1e5, 0, 100)*u.year #time in years
cluster_motion = cluster.coordinates.apply_space_motion(dt=times)
#read the fast stars and make new table
fast_stars = QTable()
fast_stars['coordinates'] = SkyCoord(ra=fs['RA_ICRS_1'],dec=fs['DE_ICRS_1'],pm_ra_cosdec=fs['pmRA'],pm_dec=fs['pmDE'],obstime = (Time('J2000')+1*u.Myr))
fast_stars['Source'] = fs['Source']
print(f'Tracing back {len(fast_stars)} stars...')
# Preallocate array to hold separations for all stars
sep_array = np.zeros((len(fast_stars), len(times)))
for i, star in enumerate(fast_stars):
# Apply space motion for the star for all time steps
star_motion = fast_stars['coordinates'][i].apply_space_motion(dt=times)
# Calculate separation for all time steps at once
sep = cluster_motion.separation(star_motion).to(u.arcmin)
sep_array[i, :] = sep
# Convert sep_array to list of lists
sep_list = sep_array.tolist()
# Create dictionary with star names as keys and separations as values
sep_dict = {star['Source']: sep_list[i] for i, star in enumerate(fast_stars)}
# make a list of stars that were inside the cluster
selected_stars_dict = {}
threshold_separation = cluster.diameter.value/separation_factor
# threshold_separation = cluster.diameter.value #deleteme
for source_id, separations in sep_dict.items():
# Check if any separation is less than the threshold
if any(sep < threshold_separation for sep in np.array(separations)):
# Add the source id to the list
selected_stars_dict[source_id] = sep_dict[source_id]
# Make list with temp
selected_stars_dict_with_temp = {}
for source_id, separations in selected_stars_dict.items():
new_star_bp_rp = fs[fs['Source']==source_id]['BP-RP']
new_star_gmag = fs[fs['Source']==source_id]['Gmag']
# Get the BP-RP and Temperature columns from the table
temperature_column = theoretical_data['Teff0']
bp_rp_column = theoretical_data['BP-RP']
gmag_column = theoretical_data['Gmag']
# Calculate the differences between the new star's BP-RP value and all values in the table
differences_bp_rp = abs(bp_rp_column - new_star_bp_rp)
# Calculate the differences between the new star's Gmag value and all values in the table
differences_gmag = abs(gmag_column - new_star_gmag)
# iso_dist = differences_bp_rp**2+differences_gmag**2 #method 1
iso_dist = differences_bp_rp #method 2
# Find the index of the star with the closest BP-RP value
closest_star_index = np.argmin(iso_dist)
# Get the temperature of the closest star
closest_star_temperature = temperature_column[closest_star_index]
selected_stars_dict_with_temp[source_id] = (selected_stars_dict[source_id],closest_star_temperature)
sources_to_mask = list(selected_stars_dict_with_temp.keys())