From 3237c57edc6e9226eab2233b936b0f438d373a92 Mon Sep 17 00:00:00 2001 From: Rocketpack23 Date: Wed, 11 Feb 2026 13:36:18 -0800 Subject: [PATCH 1/4] Added IRTF NSFCam L filter, updated filters.py and synthetic.py --- filt_func/nsfcam/L.dat | 390 +++++++++++++++++++++++++++++++++++ spisea/filters.py | 86 ++++---- spisea/synthetic.py | 446 ++++++++++++++++------------------------- 3 files changed, 606 insertions(+), 316 deletions(-) create mode 100755 filt_func/nsfcam/L.dat diff --git a/filt_func/nsfcam/L.dat b/filt_func/nsfcam/L.dat new file mode 100755 index 00000000..43589640 --- /dev/null +++ b/filt_func/nsfcam/L.dat @@ -0,0 +1,390 @@ +30156.8 0.000962257 +30175 0.00120592 +30193.2 0.00134921 +30211.5 0.00127125 +30229.7 0.00150251 +30248 0.00161243 +30266.3 0.00188732 +30284.7 0.00170279 +30303 0.00204206 +30321.4 0.00192571 +30339.8 0.00215411 +30358.2 0.0022676 +30376.7 0.00244951 +30395.1 0.00242066 +30413.6 0.0025301 +30432.1 0.00248027 +30450.7 0.0023284 +30469.2 0.00224257 +30487.8 0.00251317 +30506.4 0.0024693 +30525 0.00250125 +30543.7 0.00249219 +30562.3 0.00242495 +30581 0.00247026 +30599.8 0.00256109 +30618.5 0.00265122 +30637.3 0.00268054 +30656 0.0026195 +30674.8 0.00252914 +30693.7 0.0025599 +30712.5 0.00271273 +30731.4 0.00289178 +30750.3 0.0028379 +30769.2 0.0027976 +30788.2 0.00280762 +30807.1 0.00294924 +30826.1 0.00277042 +30845.2 0.00315022 +30864.2 0.00310421 +30883.3 0.0032618 +30902.3 0.00327802 +30921.5 0.00363994 +30940.6 0.00354338 +30959.8 0.00367117 +30978.9 0.00364733 +30998.1 0.00382757 +31017.4 0.00397182 +31036.6 0.00422025 +31055.9 0.00454044 +31075.2 0.00470376 +31094.5 0.00489855 +31113.9 0.00507116 +31133.3 0.00551271 +31152.6 0.00583959 +31172.1 0.00608921 +31191.5 0.00644302 +31211 0.00705314 +31230.5 0.00754762 +31250 0.00817514 +31269.5 0.0087316 +31289.1 0.00952125 +31308.7 0.0100901 +31328.3 0.0112662 +31348 0.012145 +31367.6 0.0135512 +31387.3 0.0147848 +31407 0.0162675 +31426.8 0.0178838 +31446.5 0.0200837 +31466.3 0.0222387 +31486.1 0.0248327 +31506 0.0277917 +31525.9 0.0313036 +31545.7 0.0353098 +31565.7 0.0398614 +31585.6 0.0450525 +31605.6 0.0513139 +31625.6 0.0586095 +31645.6 0.0669079 +31665.6 0.0764895 +31685.7 0.0874972 +31705.8 0.100346 +31725.9 0.115176 +31746 0.132822 +31766.2 0.152965 +31786.4 0.176419 +31806.6 0.202869 +31826.9 0.232942 +31847.1 0.26663 +31867.4 0.304328 +31887.8 0.346514 +31908.1 0.39176 +31928.5 0.440715 +31948.9 0.49065 +31969.3 0.54304 +31989.8 0.595515 +32010.2 0.646355 +32030.7 0.69369 +32051.3 0.736929 +32071.8 0.776407 +32092.4 0.810945 +32113 0.840741 +32133.7 0.865207 +32154.3 0.882942 +32175 0.895848 +32195.8 0.907165 +32216.5 0.914134 +32237.3 0.918469 +32258.1 0.918946 +32278.9 0.921683 +32299.7 0.921975 +32320.6 0.922349 +32341.5 0.920204 +32362.5 0.919714 +32383.4 0.918919 +32404.4 0.919387 +32425.4 0.917863 +32446.5 0.919708 +32467.5 0.918676 +32488.6 0.920314 +32509.8 0.921267 +32530.9 0.924305 +32552.1 0.924901 +32573.3 0.92679 +32594.5 0.927062 +32615.8 0.92867 +32637.1 0.930581 +32658.4 0.932568 +32679.7 0.933752 +32701.1 0.934407 +32722.5 0.934966 +32743.9 0.935562 +32765.4 0.935877 +32786.9 0.935577 +32808.4 0.937084 +32829.9 0.935535 +32851.5 0.935269 +32873.1 0.933562 +32894.7 0.934837 +32916.4 0.934186 +32938.1 0.93382 +32959.8 0.930758 +32981.5 0.932153 +33003.3 0.930245 +33025.1 0.931608 +33046.9 0.929196 +33068.8 0.928553 +33090.7 0.927891 +33112.6 0.927366 +33134.5 0.925955 +33156.5 0.925616 +33178.5 0.924155 +33200.5 0.9242 +33222.6 0.924963 +33244.7 0.926667 +33266.8 0.926907 +33288.9 0.927095 +33311.1 0.927768 +33333.3 0.927881 +33355.6 0.927298 +33377.8 0.928207 +33400.1 0.929544 +33422.5 0.92875 +33444.8 0.931329 +33467.2 0.931714 +33489.6 0.933643 +33512.1 0.933482 +33534.5 0.934703 +33557 0.933673 +33579.6 0.936075 +33602.2 0.935964 +33624.7 0.937671 +33647.4 0.937828 +33670 0.93953 +33692.7 0.93949 +33715.4 0.940382 +33738.2 0.939342 +33761 0.939697 +33783.8 0.939215 +33806.6 0.939933 +33829.5 0.939316 +33852.4 0.93844 +33875.3 0.936932 +33898.3 0.937429 +33921.3 0.935991 +33944.3 0.934474 +33967.4 0.93386 +33990.5 0.933919 +34013.6 0.934497 +34036.8 0.932878 +34059.9 0.932388 +34083.2 0.930515 +34106.4 0.930335 +34129.7 0.926779 +34153 0.928042 +34176.3 0.927017 +34199.7 0.927041 +34223.1 0.925245 +34246.6 0.925419 +34270 0.922434 +34293.6 0.922526 +34317.1 0.920628 +34340.7 0.922238 +34364.3 0.922019 +34387.9 0.92312 +34411.6 0.922412 +34435.3 0.922376 +34459 0.922568 +34482.8 0.92339 +34506.6 0.923071 +34530.4 0.922723 +34554.3 0.923598 +34578.1 0.924315 +34602.1 0.92512 +34626 0.925369 +34650 0.926424 +34674.1 0.926304 +34698.1 0.927742 +34722.2 0.927231 +34746.4 0.928697 +34770.5 0.928282 +34794.7 0.928515 +34818.9 0.928225 +34843.2 0.931308 +34867.5 0.93055 +34891.8 0.931965 +34916.2 0.931467 +34940.6 0.933056 +34965 0.931511 +34989.5 0.93317 +35014 0.932156 +35038.5 0.933175 +35063.1 0.93295 +35087.7 0.933402 +35112.4 0.93237 +35137 0.932794 +35161.7 0.933112 +35186.5 0.933091 +35211.3 0.932481 +35236.1 0.931611 +35260.9 0.931829 +35285.8 0.931225 +35310.7 0.931367 +35335.7 0.930309 +35360.7 0.930468 +35385.7 0.928837 +35410.8 0.929125 +35435.9 0.927485 +35461 0.928169 +35486.2 0.926087 +35511.4 0.926735 +35536.6 0.924653 +35561.9 0.925792 +35587.2 0.923333 +35612.5 0.924337 +35637.9 0.922857 +35663.3 0.923206 +35688.8 0.922176 +35714.3 0.922241 +35739.8 0.920767 +35765.4 0.921537 +35791 0.921272 +35816.6 0.921205 +35842.3 0.920734 +35868 0.921399 +35893.8 0.922207 +35919.5 0.922209 +35945.4 0.923026 +35971.2 0.923312 +35997.1 0.924458 +36023.1 0.923587 +36049 0.925046 +36075 0.925553 +36101.1 0.927519 +36127.2 0.927259 +36153.3 0.929243 +36179.5 0.929186 +36205.6 0.931319 +36231.9 0.931036 +36258.2 0.933156 +36284.5 0.932827 +36310.8 0.935527 +36337.2 0.935068 +36363.6 0.936368 +36390.1 0.93617 +36416.6 0.938059 +36443.1 0.938573 +36469.7 0.939709 +36496.4 0.93988 +36523 0.939998 +36549.7 0.94011 +36576.4 0.940587 +36603.2 0.940885 +36630 0.940585 +36656.9 0.939915 +36683.8 0.938672 +36710.7 0.93903 +36737.7 0.937852 +36764.7 0.937917 +36791.8 0.936105 +36818.9 0.936005 +36846 0.933661 +36873.2 0.933678 +36900.4 0.930873 +36927.6 0.931034 +36954.9 0.927961 +36982.2 0.927815 +37009.6 0.92528 +37037 0.925085 +37064.5 0.922289 +37092 0.92209 +37119.5 0.92012 +37147.1 0.919756 +37174.7 0.917933 +37202.4 0.917613 +37230.1 0.916242 +37257.8 0.915813 +37285.6 0.914763 +37313.4 0.914014 +37341.3 0.913402 +37369.2 0.912757 +37397.2 0.911949 +37425.1 0.910474 +37453.2 0.909624 +37481.3 0.906503 +37509.4 0.903579 +37537.5 0.898389 +37565.7 0.893606 +37594 0.884823 +37622.3 0.876159 +37650.6 0.863052 +37679 0.849638 +37707.4 0.830408 +37735.8 0.810358 +37764.4 0.784773 +37792.9 0.757942 +37821.5 0.7256 +37850.1 0.692032 +37878.8 0.654269 +37907.5 0.615479 +37936.3 0.57444 +37965.1 0.533411 +37993.9 0.491468 +38022.8 0.450241 +38051.8 0.410215 +38080.7 0.371899 +38109.8 0.335394 +38138.8 0.301191 +38167.9 0.269683 +38197.1 0.240773 +38226.3 0.214541 +38255.5 0.190484 +38284.8 0.169015 +38314.2 0.149827 +38343.6 0.132494 +38373 0.117145 +38402.5 0.103844 +38432 0.0915825 +38461.5 0.081001 +38491.1 0.071492 +38520.8 0.0633302 +38550.5 0.055959 +38580.2 0.0494764 +38610 0.043674 +38639.9 0.0386024 +38669.8 0.0340512 +38699.7 0.0300536 +38729.7 0.0265405 +38759.7 0.0234788 +38789.8 0.0208454 +38819.9 0.0184088 +38850 0.0161228 +38880.2 0.014158 +38910.5 0.0123668 +38940.8 0.0108354 +38971.2 0.00944018 +39001.6 0.00825095 +39032 0.00721598 +39062.5 0.00630403 +39093 0.00532269 +39123.6 0.00468755 +39154.3 0.00382352 +39185 0.00340962 +39215.7 0.00293064 +39246.5 0.00240946 +39277.3 0.00178313 +39308.2 0.00182629 +39339.1 0.00121212 +39370.1 0.00120258 +39401.1 0.000939846 diff --git a/spisea/filters.py b/spisea/filters.py index 23c029b3..23221c6b 100755 --- a/spisea/filters.py +++ b/spisea/filters.py @@ -29,7 +29,7 @@ def get_nirc2_filt(name): while len(idx) != 0: wavelength[idx+1] += 1.0e-8 - + diff = np.diff(wavelength) idx = np.where(diff <= 0)[0] #print( 'Duplicate entry loop' ) @@ -68,7 +68,7 @@ def get_2mass_filt(name): name='2MASS_{0}'.format(name)) return spectrum - + def get_vista_filt(name): """ @@ -79,21 +79,21 @@ def get_vista_filt(name): t = Table.read('{0}/vista/VISTA_Filters_at80K_forETC_{1}.dat'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find VISTA filter file {0}/vista/VISTA_Filters_at80K_forETC_{1}.dat'.format(filters_dir, name)) + raise ValueError('Could not find VISTA filter file {0}/vista/VISTA_Filters_at80K_forETC_{1}.dat'.format(filters_dir, name)) # Wavelength must be in angstroms, transmission in fraction wave = t['col1'] * 10 trans = t['col2'] * 0.01 - + # Change any negative numbers to 0, as well as anything shortward # of 0.4 microns or longward of 2.9 microns # (no VISTA filter transmissions beyond these boundaries) bad = np.where( (trans < 0) | (wave < 4000) | (wave > 29000) ) trans[bad] = 0 - + # Now we can define the VISTA filter bandpass objects spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='VISTA_{0}'.format(name)) - + return spectrum def get_decam_filt(name): @@ -104,18 +104,18 @@ def get_decam_filt(name): try: t = Table.read('{0}/decam/DECam_filters.txt'.format(filters_dir), format='ascii') t.rename_column('Y', 'y') - + cols = np.array(t.keys()) idx = np.where(cols == name)[0][0] trans = t[cols[idx]] except: - raise ValueError('Could not find DECAM filter {0} in {1}/decam/DECam_filters.txt'.format(name, filters_dir)) + raise ValueError('Could not find DECAM filter {0} in {1}/decam/DECam_filters.txt'.format(name, filters_dir)) # Limit to unmasked regions only mask = np.ma.getmask(trans) good = np.where(mask == False) - + # Convert wavelengths from nm to angstroms, while eliminating masked regions wave = t['wavelength'][good] * 10. trans = trans[good] @@ -126,7 +126,7 @@ def get_decam_filt(name): return spectrum -def get_PS1_filt(name): +def get_PS1_filt(name): """ Define PS1 filter as pysynphot object """ @@ -145,11 +145,11 @@ def get_PS1_filt(name): trans = t[cols[idx]] except: - raise ValueError('Could not find PS1 filter {0} in {1}/ps1'.format(name, filters_dir)) + raise ValueError('Could not find PS1 filter {0} in {1}/ps1'.format(name, filters_dir)) # Convert wavelengths from nm to angstroms wave = t['wave'] * 10. - + spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='ps1_{0}'.format(name)) return spectrum @@ -161,7 +161,7 @@ def get_jwst_filt(name): try: t = Table.read('{0}/jwst/{1}.txt'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find JWST filter {0} in {1}/jwst'.format(name, filters_dir)) + raise ValueError('Could not find JWST filter {0} in {1}/jwst'.format(name, filters_dir)) # Convert wavelengths to angstroms wave = t['microns'] * 10**4. @@ -170,10 +170,10 @@ def get_jwst_filt(name): # Change any negative numbers to 0 bad = np.where(trans < 0) trans[bad] = 0 - + spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='jwst_{0}'.format(name)) - return spectrum + return spectrum def get_Johnson_Glass_filt(name): """ @@ -182,7 +182,7 @@ def get_Johnson_Glass_filt(name): try: t = Table.read('{0}/Johnson_Glass/{1}.txt'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find Johnson-Glass filter {0} in {1}/Johnson_Glass'.format(name, filters_dir)) + raise ValueError('Could not find Johnson-Glass filter {0} in {1}/Johnson_Glass'.format(name, filters_dir)) # Convert wavelengths to angstroms wave = t['col1'] * 10. @@ -191,10 +191,10 @@ def get_Johnson_Glass_filt(name): # Change any negative numbers to 0 bad = np.where(trans < 0) trans[bad] = 0 - + spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='jg_{0}'.format(name)) - return spectrum + return spectrum def get_nirc1_filt(name): """ @@ -203,12 +203,12 @@ def get_nirc1_filt(name): try: t = Table.read('{0}/nirc1/{1}.txt'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find NIRC1 filter {0} in {1}/nirc1'.format(name, filters_dir)) + raise ValueError('Could not find NIRC1 filter {0} in {1}/nirc1'.format(name, filters_dir)) # Convert wavelengths to angstroms wave = t['col1'] * 10**4 trans = t['col2'] - + # Lets fix wavelength array for duplicate values or negative vals; # delete these entries diff = np.diff(wave) @@ -222,14 +222,14 @@ def get_nirc1_filt(name): diff = np.diff(wave) idx = np.where(diff <= 0)[0] - + # Change any negative transmission vals to 0 bad = np.where(trans < 0) trans[bad] = 0 spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='nirc1_{0}'.format(name)) - return spectrum + return spectrum def get_ctio_osiris_filt(name): """ @@ -238,7 +238,7 @@ def get_ctio_osiris_filt(name): try: t = Table.read('{0}/CTIO_OSIRIS/{1}.txt'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find CTIO/OSIRIS filter {0} in {1}/CTIO_OSIRIS'.format(name, filters_dir)) + raise ValueError('Could not find CTIO/OSIRIS filter {0} in {1}/CTIO_OSIRIS'.format(name, filters_dir)) # Convert wavelengths to angstroms wave = t['col1'] * 10**4 @@ -247,7 +247,7 @@ def get_ctio_osiris_filt(name): # Change any negative numbers to 0 bad = np.where(trans < 0) trans[bad] = 0 - + spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='ctio_osiris_{0}'.format(name)) return spectrum @@ -259,7 +259,7 @@ def get_naco_filt(name): try: t = Table.read('{0}/naco/{1}.dat'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find NACO filter {0} in {1}/naco'.format(name, filters_dir)) + raise ValueError('Could not find NACO filter {0} in {1}/naco'.format(name, filters_dir)) # Convert wavelengths to angstroms wave = t['col1'] * 10**4 @@ -268,7 +268,7 @@ def get_naco_filt(name): # Change any negative numbers to 0 bad = np.where(trans < 0) trans[bad] = 0 - + spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='naco_{0}'.format(name)) return spectrum @@ -282,7 +282,7 @@ def get_ubv_filt(name): except: raise ValueError('Could not find ubv filter {0} in {1}/ubv'.format(name, filters_dir)) - # Convert wavelength from nm to angstroms + # Convert wavelength from nm to angstroms wave = t[t.keys()[0]] * 10. # Convert transmission to ratio (from percent) trans = t[t.keys()[1]] / 100. @@ -336,8 +336,8 @@ def get_keck_osiris_filt(name): def get_gaia_filt(version, name): """ Define Gaia filters as pysynphot object. - To avoid confusion, we will only support - the revised DR2 zeropoints from + To avoid confusion, we will only support + the revised DR2 zeropoints from Evans+18. version: specify dr1, dr2, or dr2_rev @@ -357,7 +357,7 @@ def get_gaia_filt(version, name): path = '{0}/gaia/dr2_rev/'.format(filters_dir) else: raise ValueError('GAIA filter version {0} not understood. Please use dr1, dr2, or dr2_rev'.format(version)) - + # Get the filter info try: t = Table.read('{0}/Gaia_passbands.txt'.format(path), format='ascii') @@ -380,7 +380,7 @@ def get_gaia_filt(version, name): bad = np.where(trans > 90) trans[bad] = 0 except: - raise ValueError('Could not find Gaia filter {0}'.format(name)) + raise ValueError('Could not find Gaia filter {0}'.format(name)) # Convert wavelengths to angstroms (from nm) wave = t['LAMBDA'] * 10 @@ -437,7 +437,7 @@ def get_rubin_filt(name): try: t = Table.read('{0}/rubin/{1}.dat'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find Rubin LSST filter file {0}/rubin/{1}.dat'.format(filters_dir, name)) + raise ValueError('Could not find Rubin LSST filter file {0}/total_/{1}.dat'.format(filters_dir, name)) wavelength = t[t.keys()[0]] transmission = t[t.keys()[1]] @@ -451,24 +451,20 @@ def get_rubin_filt(name): return spectrum -def get_euclid_filt(name): +def get_nsfcam_filt(name): + """ - Define the Euclid filters as a pysynphot spectrum object + Define irtf nsfcam filters as pysynphot object """ - # Read in filter info try: - t = Table.read('{0}/euclid/{1}.dat'.format(filters_dir, name), format='ascii') + t = Table.read('{0}/nsfcam/{1}.dat'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find Euclid filter file {0}/euclid/{1}.dat'.format(filters_dir, name)) - - wavelength = t[t.keys()[0]] - transmission = t[t.keys()[1]] + raise ValueError('Could not find nsfcam filter {0} in {1}/nsfcam'.format(name, filters_dir)) - # Convert wavelength to Angstroms - wavelength = wavelength * 10 + # Wavelength already in angstrom and and transmission in fraction + wave = t['col1'] + trans = t['col2'] - # Make spectrum object - spectrum = pysynphot.ArrayBandpass(wavelength, transmission, waveunits='angstrom', - name='euclid_{0}'.format(name)) + spectrum = pysynphot.ArrayBandpass(wave, trans, waveunits='angstrom', name='nsfcam_{0}'.format(name)) return spectrum diff --git a/spisea/synthetic.py b/spisea/synthetic.py index fa8220f6..882a18d3 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -37,7 +37,7 @@ def Vega(): # Use Vega as our zeropoint... assume V=0.03 mag and all colors = 0.0 # These parameters are defined in Girardi+02 - vega = atm.get_kurucz_atmosphere(temperature=9550, + vega = atm.get_kurucz_atmosphere(temperature=9550, gravity=3.95, metallicity=-0.5) @@ -47,8 +47,8 @@ def Vega(): # This is (R/d)**2 as reported by Girardi et al. 2002, page 198, col 1. # and is used to convert to flux observed at Earth. - vega *= 6.247e-17 - + vega *= 6.247e-17 + return vega vega = Vega() @@ -56,13 +56,13 @@ def Vega(): class Cluster(object): """ Base class to create a cluster with user-specified isochrone, - imf, ifmr, and total mass. + imf, ifmr, and total mass. Parameters ----------- iso: isochrone object SPISEA isochrone object - + imf: imf object SPISEA IMF object @@ -90,17 +90,17 @@ def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=False, self.ifmr = ifmr self.cluster_mass = cluster_mass self.seed = seed - + return - + class ResolvedCluster(Cluster): """ Cluster sub-class that produces a *resolved* stellar cluster. - A table is output with the synthetic photometry and intrinsic - properties of the individual stars (or stellar systems, if + A table is output with the synthetic photometry and intrinsic + properties of the individual stars (or stellar systems, if mutliplicity is used in the IMF object). - If multiplicity is used, than a second table is produced that + If multiplicity is used, than a second table is produced that contains the properties of the companion stars independent of their primary stars. @@ -108,7 +108,7 @@ class ResolvedCluster(Cluster): ----------- iso: isochrone object SPISEA isochrone object - + imf: imf object SPISEA IMF object @@ -120,12 +120,6 @@ class ResolvedCluster(Cluster): produced by the cluster at the given isochrone age. Otherwise, no compact remnants are produced. - keep_low_mass_stars: boolean (default False) - If True, the cluster will not cut out stars below the isochrone grid - on initial mass. They are assigned a current mass equal to their initial - mass, a phase of 98, and no other evolutionary properties or photometry. - If False, stars below the isochrone initial mass limit are cut out. - seed: int If set to non-None, all random sampling will be seeded with the specified seed, forcing identical output. @@ -135,7 +129,7 @@ class ResolvedCluster(Cluster): True for verbose output. """ def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=True, - seed=None, keep_low_mass_stars=False): + seed=None): Cluster.__init__(self, iso, imf, cluster_mass, ifmr=ifmr, verbose=verbose, seed=seed) # Provide a user warning is random seed is set @@ -143,7 +137,7 @@ def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=True, print('WARNING: random seed set to %i' % seed) t1 = time.time() - ##### + ##### # Sample the IMF to build up our cluster mass. ##### mass, isMulti, compMass, sysMass = imf.generate_cluster(cluster_mass, @@ -161,27 +155,27 @@ def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=True, for ikey in interp_keys: self.iso_interps[ikey] = interpolate.interp1d(self.iso.points['mass'], self.iso.points[ikey], kind='linear', bounds_error=False, fill_value=np.nan) - - ##### + + ##### # Make a table to contain all the information about each stellar system. ##### star_systems = self._make_star_systems_table(mass, isMulti, sysMass) - + # Trim out bad systems; specifically, stars with masses outside those provided # by the model isochrone (except for compact objects). - star_systems, compMass = self._remove_bad_systems(star_systems, compMass, keep_low_mass_stars) + star_systems, compMass = self._remove_bad_systems(star_systems, compMass) - ##### + ##### # Make a table to contain all the information about companions. ##### if self.imf.make_multiples: companions = self._make_companions_table(star_systems, compMass) - + ##### # Save our arrays to the object ##### self.star_systems = star_systems - + if self.imf.make_multiples: self.companions = companions @@ -192,13 +186,13 @@ def set_filter_names(self): Set filter column names """ filt_names = [] - + for col_name in self.iso.points.colnames: if 'm_' in col_name: filt_names.append(col_name) return filt_names - + def _make_star_systems_table(self, mass, isMulti, sysMass): """ Make a star_systems table and get synthetic photometry for each primary star. @@ -245,7 +239,7 @@ def _make_star_systems_table(self, mass, isMulti, sysMass): for ii in range(len(bad[0])): print('WARNING: changing phase {0} to 5'.format(star_systems['phase'][bad[0][ii]])) star_systems['phase'][bad] = 5 - + for filt in self.filt_names: star_systems[filt] = self.iso_interps[filt](star_systems['mass']) @@ -253,14 +247,14 @@ def _make_star_systems_table(self, mass, isMulti, sysMass): # Make Remnants # Note: Some models already have WDs in them. If they do, then they shouldn't # be handled by this code here (because their Teff > 0). - # + # # Remnants have flux = 0 in all bands if they are generated here. - ##### + ##### if self.ifmr != None: # Identify compact objects as those with Teff = 0 or with phase > 100. highest_mass_iso = self.iso.points['mass'].max() idx_rem = np.where((np.isnan(star_systems['Teff'])) & (star_systems['mass'] > highest_mass_iso))[0] - + # Calculate remnant mass and ID for compact objects; update remnant_id and # remnant_mass arrays accordingly if 'metallicity_array' in inspect.getfullargspec(self.ifmr.generate_death_mass).args: @@ -282,13 +276,13 @@ def _make_star_systems_table(self, mass, isMulti, sysMass): star_systems[filt][idx_rem_good] = np.full(len(idx_rem_good), np.nan) return star_systems - + def _make_companions_table(self, star_systems, compMass): N_systems = len(star_systems) - + ##### - # MULTIPLICITY + # MULTIPLICITY # Make a second table containing all the companion-star masses. # This table will be much longer... here are the arrays: # sysIndex - the index of the system this star belongs too @@ -312,17 +306,17 @@ def _make_companions_table(self, star_systems, compMass): companions.add_column( Column(np.empty(N_comp_tot, dtype=float), name='metallicity') ) for filt in self.filt_names: companions.add_column( Column(np.empty(N_comp_tot, dtype=float), name=filt) ) - + if isinstance(self.imf._multi_props, multiplicity.MultiplicityResolvedDK): companions.add_column( Column(np.zeros(N_comp_tot, dtype=float), name='log_a') ) companions.add_column( Column(np.zeros(N_comp_tot, dtype=float), name='e') ) companions.add_column( Column(np.zeros(N_comp_tot, dtype=float), name='i', description = 'degrees') ) companions.add_column( Column(np.zeros(N_comp_tot, dtype=float), name='Omega') ) companions.add_column( Column(np.zeros(N_comp_tot, dtype=float), name='omega') ) - + for ii in range(len(companions)): companions['log_a'][ii] = self.imf._multi_props.log_semimajoraxis(star_systems['mass'][companions['system_idx'][ii]]) - + companions['e'] = self.imf._multi_props.random_e(np.random.rand(N_comp_tot)) companions['i'], companions['Omega'], companions['omega'] = self.imf._multi_props.random_keplarian_parameters(np.random.rand(N_comp_tot),np.random.rand(N_comp_tot),np.random.rand(N_comp_tot)) @@ -330,7 +324,7 @@ def _make_companions_table(self, star_systems, compMass): # Make an array that maps system index (ii), companion index (cc) to # the place in the 1D companions array. N_comp_max = N_companions.max() - + comp_index = np.zeros((N_systems, N_comp_max), dtype=int) kk = 0 for ii in range(N_systems): @@ -346,9 +340,9 @@ def _make_companions_table(self, star_systems, compMass): idx = np.where(N_companions >= cc)[0] # Get the location in the companions array for each system and - # the cc'th companion. + # the cc'th companion. cdx = comp_index[idx, cc-1] - + companions['mass'][cdx] = [compMass[ii][cc-1] for ii in idx] comp_mass = companions['mass'][cdx] @@ -397,13 +391,13 @@ def _make_companions_table(self, star_systems, compMass): # as np.nan. Otherwise, add fluxes together good = np.where( (f1 != 0) | (f2 != 0) ) bad = np.where( (f1 == 0) & (f2 == 0) ) - + star_systems[filt][idx[good]] = -2.5 * np.log10(f1[good] + f2[good]) star_systems[filt][idx[bad]] = np.nan ##### # Make Remnants with flux = 0 in all bands. - ##### + ##### if self.ifmr != None: # Identify compact objects as those with Teff = 0 or with masses above the max iso mass highest_mass_iso = self.iso.points['mass'].max() @@ -417,7 +411,7 @@ def _make_companions_table(self, star_systems, compMass): metallicity_array=companions['metallicity'][cdx_rem]) else: r_mass_tmp, r_id_tmp = self.ifmr.generate_death_mass(mass_array=companions['mass'][cdx_rem]) - + # Drop remnants where it is not relevant (e.g. not a compact object or # outside mass range IFMR is defined for) @@ -426,7 +420,7 @@ def _make_companions_table(self, star_systems, compMass): companions['mass_current'][cdx_rem_good] = r_mass_tmp[good] companions['phase'][cdx_rem_good] = r_id_tmp[good] - + # Give remnants a magnitude of nan, so they can be filtered out later when calculating flux. for filt in self.filt_names: companions[filt][cdx_rem_good] = np.full(len(cdx_rem_good), np.nan) @@ -439,25 +433,19 @@ def _make_companions_table(self, star_systems, compMass): if len(idx) != N_comp_tot and self.verbose: print( 'Found {0:d} companions out of stellar mass range'.format(N_comp_tot - len(idx))) - # For low-mass stars and substellar objects below isochrone, assume no mass loss and set phase to 98 - low_mass_idxs = (companions['mass']0: - assert companions['mass'][idx].min() > 0 + assert companions['mass'][idx].min() > 0 return companions - - def _remove_bad_systems(self, star_systems, compMass, keep_low_mass_stars): + + def _remove_bad_systems(self, star_systems, compMass): """ Helper function to remove stars with masses outside the isochrone - mass range from the cluster. These stars are identified by having + mass range from the cluster. These stars are identified by having a Teff = 0, as set up by _make_star_systems_table_interp. - If self.ifmr == None, then both high and low-mass bad systems are - removed. If self.ifmr != None, then we will save the high mass systems + If self.ifmr == None, then both high and low-mass bad systems are + removed. If self.ifmr != None, then we will save the high mass systems since they will be plugged into an ifmr later. """ N_systems = len(star_systems) @@ -466,41 +454,23 @@ def _remove_bad_systems(self, star_systems, compMass, keep_low_mass_stars): # Convert nan_to_num to avoid errors on greater than, less than comparisons star_systems_teff_non_nan = np.nan_to_num(star_systems['Teff'], nan=-99) star_systems_phase_non_nan = np.nan_to_num(star_systems['phase'], nan=-99) - if (self.ifmr == None) and (not keep_low_mass_stars): - print('Remove low mass stars below grid and compact objects') + if self.ifmr == None: # Keep only those stars with Teff assigned. idx = np.where(star_systems_teff_non_nan > 0)[0] - elif not keep_low_mass_stars: - print('Remove low mass stars, keep compact objects') - # Keep stars (with Teff) and any other compact objects (with phase info). - idx = np.where( (star_systems_teff_non_nan > 0) | (star_systems_phase_non_nan >= 0) )[0] - elif self.ifmr == None: - print('Remove compact objects, keep low mass stars below grid') - # Keep stars (with Teff) and objects below mass grid - idx = np.where( (star_systems_teff_non_nan > 0) | ((star_systems['mass'] 0) | (star_systems_phase_non_nan >= 0) | - ((star_systems['mass'] 0) | (star_systems_phase_non_nan >= 0) )[0] if len(idx) != N_systems and self.verbose: print( 'Found {0:d} stars out of mass range'.format(N_systems - len(idx))) - if keep_low_mass_stars: - lm_idx = np.where(star_systems['mass']0: - self.verbose = True - print('Missing photometry for',len(comp_filters),'filter - recomputing these columns:',comp_filters) - - print('Loading stellar spectra') - # Initialize output for stellar spectra - self.spec_list = [] - # For each isochrone point, extract the synthetic photometry. - for ii in range(len(self.points['Teff'])): - # Loop is currently taking about 0.11 s per iteration - gravity = float( self.points['logg'][ii] ) - L = float( (self.points['L'][ii]*units.W).cgs / (units.erg / units.s)) # in erg/s - T = float( self.points['Teff'][ii] ) # in Kelvin - R = float( (self.points['R'][ii]*units.m).to('pc') / units.pc) # in pc - phase = int(self.points['phase'][ii]) - # Get the atmosphere model now. Wavelength is in Angstroms - # This is the time-intensive call... everything else is negligable. - # If source is a star, pull from star atmospheres. If it is a WD, - # pull from WD atmospheres - if phase == 101: - star = wd_atm_func(temperature=T, gravity=gravity, metallicity=metallicity, - verbose=False) - else: - star = atm_func(temperature=T, gravity=gravity, metallicity=metallicity, - rebin=rebin) - # Trim wavelength range down to appropriate range - star = spectrum.trimSpectrum(star, wave_range[0], wave_range[1]) - # Convert into flux observed at Earth (unreddened) - star *= (R / self.points.meta["DISTANCE"])**2 # in erg s^-1 cm^-2 A^-1 - # Redden the spectrum. This doesn't take much time at all. - red = red_law.reddening(AKs).resample(star.wave) - star *= red - # Save the final spectrum to our spec_list for later use. - self.spec_list.append(star) - - self.make_photometry(rebin=rebin, vega=vega, comp_filters=comp_filters) - return - def make_photometry(self, rebin=True, vega=vega, comp_filters=None): - """ + def make_photometry(self, rebin=True, vega=vega): + """ Make synthetic photometry for the specified filters. This function udpates the self.points table to include new columns with the photometry. - + """ startTime = time.time() @@ -1110,16 +1035,12 @@ def make_photometry(self, rebin=True, vega=vega, comp_filters=None): npoints = len(self.points) verbose_fmt = 'M = {0:7.3f} Msun T = {1:5.0f} K m_{2:s} = {3:4.2f}' - #Calculate all filters, or select filters - if comp_filters is None: - comp_filters = self.filters - # Loop through the filters, get filter info, make photometry for # all stars in this filter. - for ii in comp_filters: + for ii in self.filters: prt_fmt = 'Starting filter: {0:s} Elapsed time: {1:.2f} seconds' print( prt_fmt.format(ii, time.time() - startTime)) - + filt = get_filter_info(ii, rebin=rebin, vega=vega) filt_name = get_filter_col_name(ii) @@ -1127,15 +1048,15 @@ def make_photometry(self, rebin=True, vega=vega, comp_filters=None): col_name = 'm_' + filt_name mag_col = Column(np.zeros(npoints, dtype=float), name=col_name) self.points.add_column(mag_col) - + # Loop through each star in the isochrone and do the filter integration print('Starting synthetic photometry') for ss in range(npoints): star = self.spec_list[ss] # These are already extincted, observed spectra. star_mag = mag_in_filter(star, filt) - + self.points[col_name][ss] = star_mag - + if (self.verbose and (ss % 100) == 0): print( verbose_fmt.format(self.points['mass'][ss], self.points['Teff'][ss], filt_name, star_mag)) @@ -1152,34 +1073,27 @@ def make_photometry(self, rebin=True, vega=vega, comp_filters=None): def check_save_file(self, evo_model, atm_func, red_law): """ - Check to see if save_file exists, as saved by the save_file - and save_file_legacy objects. If the filename exists, check the + Check to see if save_file exists, as saved by the save_file + and save_file_legacy objects. If the filename exists, check the meta-data as well. returns a boolean: True is file exists, false otherwise """ out_bool = False - + if os.path.exists(self.save_file) | os.path.exists(self.save_file_legacy): try: tmp = Table.read(self.save_file) except: tmp = Table.read(self.save_file_legacy) - - + + # See if the meta-data matches: evo model, atm_func, redlaw if ( (tmp.meta['EVOMODEL'] == type(evo_model).__name__) & (tmp.meta['ATMFUNC'] == atm_func.__name__) & (tmp.meta['REDLAW'] == red_law.name) ): out_bool = True - # Check model version if it was logged - if 'EVOMODELVERSION' in tmp.meta: - if tmp.meta['EVOMODELVERSION']!=evo_model.model_version_name: - out_bool=False - else: - print(f"No version information found for evolution model {type(evo_model).__name__}.") - return out_bool def plot_CMD(self, mag1, mag2, savefile=None): @@ -1193,7 +1107,7 @@ def plot_CMD(self, mag1, mag2, savefile=None): mag2 : string The name of the second magnitude column to be plotted. savefile : string (default None) - If a savefile is specified, then the plot will be saved to that file. + If a savefile is specified, then the plot will be saved to that file. """ plt.clf() plt.plot(self.points[mag1] - self.points[mag2], self.points[mag1], @@ -1201,7 +1115,7 @@ def plot_CMD(self, mag1, mag2, savefile=None): plt.gca().invert_yaxis() plt.xlabel(mag1 + ' - ' + mag2 + ' (mag)') plt.ylabel(mag1 + ' (mag)') - + fmt_title = 'logAge={0:.2f}, d={1:.2f} kpc, AKs={2:.2f}' plt.title(fmt_title.format(self.points.meta['LOGAGE'], self.points.meta['DISTANCE']/1e3, @@ -1209,34 +1123,34 @@ def plot_CMD(self, mag1, mag2, savefile=None): if savefile != None: plt.savefig(savefile) - + return def plot_mass_magnitude(self, mag, savefile=None): """ Make a standard mass-luminosity relation plot for this isochrone. - + Parameters ---------- mag : string The name of the magnitude column to be plotted. savefile : string (default None) - If a savefile is specified, then the plot will be saved to that file. + If a savefile is specified, then the plot will be saved to that file. """ plt.clf() plt.semilogx(self.points['mass'], self.points[mag], 'k.') plt.gca().invert_yaxis() plt.xlabel(r'Mass (M$_\odot$)') plt.ylabel(mag + ' (mag)') - + fmt_title = 'logAge={0:.2f}, d={1:.2f} kpc, AKs={2:.2f}' plt.title(fmt_title.format(self.points.meta['LOGAGE'], self.points.meta['DISTANCE']/1e3, self.points.meta['AKS'])) - + if savefile != None: plt.savefig(savefile) - + return #===================================================# @@ -1258,7 +1172,7 @@ def __init__(self, logAge, distance, evo_model=default_evo_model, Functions on this object: apply_reddening (Apply reddening with defined redlaw and AKs) make_photometry (make synthetic photometry for spectra) - + Parameters ---------- logAge : float @@ -1296,13 +1210,13 @@ def __init__(self, logAge, distance, evo_model=default_evo_model, spectrum. This is very useful to save computation time down the road. """ - t1 = time.time() + t1 = time.time() c = constants # Get solar metallicity models for a population at a specific age. # Takes about 0.1 seconds. - evol = evo_model.isochrone(age=10**logAge) # solar metallicity - + evol = evo_model.isochrone(age=10**logAge) # solar metallicity + # Eliminate cases where log g is less than 0 idx = np.where(evol['logg'] > 0) evol = evol[idx] @@ -1313,8 +1227,8 @@ def __init__(self, logAge, distance, evo_model=default_evo_model, evol = evol[idx] if max_mass != None: idx = np.where(evol['mass'] <= max_mass) - evol = evol[idx] - + evol = evol[idx] + # Trim down the table by selecting every Nth point where # N = mass sampling factor. evol = evol[::mass_sampling] @@ -1349,31 +1263,30 @@ def __init__(self, logAge, distance, evo_model=default_evo_model, # Get the atmosphere model now. Wavelength is in Angstroms # This is the time-intensive call... everything else is negligable. star = atm_func(temperature=T, gravity=gravity) - + # Trim wavelength range down to JHKL range (0.5 - 5.2 microns) star = spectrum.trimSpectrum(star, wave_range[0], wave_range[1]) # Convert into flux observed at Earth (unreddened) star *= (R / distance)**2 # in erg s^-1 cm^-2 A^-1 - - # Save the final spectrum to our spec_list for later use. + + # Save the final spectrum to our spec_list for later use. self.spec_list.append(star) # Append all the meta data to the summary table. - + tab.meta['ATMFUNC'] = atm_func.__name__ tab.meta['EVOMODEL'] = type(evo_model).__name__ - tab.meta['EVOMODELVERSION'] = evo_model.model_version_name tab.meta['LOGAGE'] = logAge tab.meta['DISTANCE'] = distance tab.meta['WAVEMIN'] = wave_range[0] tab.meta['WAVEMAX'] = wave_range[1] self.points = tab - + t2 = time.time() print('Isochrone generation took {0:f} s.'.format(t2-t1)) - + return def apply_reddening(self, AKs, extinction_law, dAKs=0, dist='uniform', dAKs_max=None): @@ -1385,7 +1298,7 @@ def apply_reddening(self, AKs, extinction_law, dAKs=0, dist='uniform', dAKs_max= ---------- AKs: float Total extinction in AKs - + extinction_law: SPISEA extinction object Extinction law to be used on the spectra @@ -1395,13 +1308,13 @@ def apply_reddening(self, AKs, extinction_law, dAKs=0, dist='uniform', dAKs_max= dAKs_max: float or None If not none, defines the maximum |dAKs| a star can - have in gaussian distribution case + have in gaussian distribution case dist: string, 'uniform' or 'gaussian' Distribution to draw differential reddening from. If uniform, dAKs will cut off at Aks +/- dAKs. Otherwise, will draw from Gaussian of width AKs +/- dAks - + """ self.AKs = np.ones(len(self.spec_list)) # Apply reddening to each object in the spec list @@ -1430,7 +1343,7 @@ def apply_reddening(self, AKs, extinction_law, dAKs=0, dist='uniform', dAKs_max= else: AKs_act = AKs - red = extinction_law.reddening(AKs_act).resample(star.wave) + red = extinction_law.reddening(AKs_act).resample(star.wave) star *= red # Update the spectrum in spec list @@ -1443,7 +1356,7 @@ def apply_reddening(self, AKs, extinction_law, dAKs=0, dist='uniform', dAKs_max= return def make_photometry(self, filters, rebin=True): - """ + """ Make synthetic photometry for the specified filters. This function udpates the self.points table to include new columns with the photometry. @@ -1453,11 +1366,11 @@ def make_photometry(self, filters, rebin=True): filters : dictionary A dictionary containing the filter name (for the output columns) and the filter specification string that can be processed by pysynphot. - + rebin: boolean - True to rebin filter function (only used if non-zero transmission points are + True to rebin filter function (only used if non-zero transmission points are larger than 1500 points) - + """ npoints = len(self.points) @@ -1474,35 +1387,35 @@ def make_photometry(self, filters, rebin=True): col_name = 'mag_' + filt_name mag_col = Column(np.zeros(npoints, dtype=float), name=col_name) self.points.add_column(mag_col) - + # Loop through each star in the isochrone and do the filter integration for ss in range(npoints): star = self.spec_list[ss] # These are already extincted, observed spectra. star_mag = mag_in_filter(star, filt) - + self.points[col_name][ss] = star_mag - - + + endTime = time.time() print( ' Time taken: {0:.2f} seconds'.format(endTime - ts)) return def get_filter_info(name, vega=vega, rebin=True): - """ + """ Define filter functions, setting ZP according to Vega spectrum. Input name is the SPISEA obs_string """ tmp = name.split(',') filterName = tmp[-1] - + if name.startswith('nirc2'): filt = filters.get_nirc2_filt(filterName) elif name.startswith('2mass'): filt = filters.get_2mass_filt(filterName) - + elif name.startswith('vista'): filt = filters.get_vista_filt(filterName) @@ -1517,10 +1430,10 @@ def get_filter_info(name, vega=vega, rebin=True): elif name.startswith('jg'): filt = filters.get_Johnson_Glass_filt(filterName) - + elif name.startswith('nirc1'): filt = filters.get_nirc1_filt(filterName) - + elif name.startswith('ctio_osiris'): filt = filters.get_ctio_osiris_filt(filterName) @@ -1532,7 +1445,7 @@ def get_filter_info(name, vega=vega, rebin=True): elif name.startswith('ukirt'): filt = filters.get_ukirt_filt(filterName) - + elif name.startswith('keck_osiris'): filt = filters.get_keck_osiris_filt(filterName) @@ -1545,20 +1458,20 @@ def get_filter_info(name, vega=vega, rebin=True): elif name.startswith('hawki'): filt = filters.get_hawki_filt(filterName) - + elif name.startswith('rubin'): filt = filters.get_rubin_filt(filterName) - elif name.startswith('euclid'): - filt = filters.get_euclid_filt(filterName) - + elif name.startswith('nsfcam'): + filt = filters.get_nsfcam_filt(filterName) + else: # Otherwise, look for the filter info in the cdbs/mtab and cdbs/comp files try: filt = ObsBandpass(name) except: raise Exception('Filter {0} not understood. Check spelling and make sure cdbs/mtab and cdbs/comp files are up to date'.format(name)) - + # Convert to ArraySpectralElement for resampling. filt = spectrum.ArraySpectralElement(filt.wave, filt.throughput, waveunits=filt.waveunits, @@ -1576,14 +1489,14 @@ def get_filter_info(name, vega=vega, rebin=True): # Otherwise, throw an error idx = np.where(filt.throughput > 0.001)[0] if (min(filt.wave[idx]) < min(vega.wave)) | (max(filt.wave[idx]) > max(vega.wave)): - raise ValueError('Vega spectrum doesnt cover filter wavelength range!') + raise ValueError('Vega spectrum doesnt cover filter wavelength range!') vega_obs = obs.Observation(vega, filt, binset=filt.wave, force='taper') #vega_flux = vega_obs.binflux.sum() diff = np.diff(vega_obs.binwave) diff = np.append(diff, diff[-1]) vega_flux = np.sum(vega_obs.binflux * diff) - + vega_mag = 0.03 filt.flux0 = vega_flux @@ -1594,7 +1507,7 @@ def get_filter_info(name, vega=vega, rebin=True): def get_filter_col_name(obs_str): """ - Get standard column name for synthetic photometry based on + Get standard column name for synthetic photometry based on the input string. The input string is expected to be an appropriate SPISEA obs_string """ @@ -1613,7 +1526,7 @@ def get_filter_col_name(obs_str): filt_name = 'hst_{0}'.format(tmp[-1]) else: filt_name = '{0}_{1}'.format(tmp[0], tmp[1]) - + return filt_name def get_obs_str(col): @@ -1623,7 +1536,7 @@ def get_obs_str(col): """ # Remove the trailing m_ name = col[2:] - + # Define dictionary for filters filt_list = {'hst_f127m': 'wfc3,ir,f127m', 'hst_f139m': 'wfc3,ir,f139m', 'hst_f153m': 'wfc3,ir,f153m', 'hst_f814w': 'acs,wfc1,f814w', 'hst_f125w': 'wfc3,ir,f125w', 'hst_f160w': 'wfc3,ir,f160w', @@ -1645,7 +1558,7 @@ def get_obs_str(col): 'jwst_F187N': 'jwst,F187N', 'jwst_F200W': 'jwst,F200W', 'jwst_F210M': 'jwst,F210M', - 'jwst_F250M': 'jwst,F250M', + 'jwst_F250M': 'jwst,F250M', 'jwst_F277W': 'jwst,F277W', 'jwst_F300M': 'jwst,F300M', 'jwst_F322W2': 'jwst,F322W2', @@ -1665,14 +1578,10 @@ def get_obs_str(col): 'nirc2_FeII': 'nirc2,FeII', 'nirc2_Brgamma': 'nirc2,Brgamma', '2mass_J': '2mass,J', '2mass_H': '2mass,H', '2mass_Ks': '2mass,Ks', 'ubv_U':'ubv,U', 'ubv_B':'ubv,B', 'ubv_V':'ubv,V', 'ubv_R':'ubv,R', - 'ubv_I':'ubv,I', + 'ubv_I':'ubv,I', 'jg_J': 'jg,J', 'jg_H': 'jg,H', 'jg_K': 'jg,K', 'nirc1_K':'nirc1,K', 'nirc1_H':'nirc1,H', 'naco_J':'naco,J', 'naco_H':'naco,H', 'naco_Ks':'naco,Ks', - 'naco_IB_2.00': 'naco,IB_2.00', 'naco_IB_2.03':'naco,IB_2.03', 'naco_IB_2.06':'naco,IB_2.06', - 'naco_IB_2.24':'naco,IB_2.24', 'naco_IB_2.27':'naco,IB_2.27', - 'naco_IB_2.30':'naco,IB_2.30', 'naco_IB_2.33':'naco,IB_2.33', - 'naco_IB_2.36':'naco,IB_2.36', 'ukirt_J':'ukirt,J', 'ukirt_H':'ukirt,H', 'ukirt_K':'ukirt,K', 'ctio_osiris_H': 'ctio_osiris,H', 'ctio_osiris_K': 'ctio_osiris,K', 'ztf_g':'ztf,g', 'ztf_r':'ztf,r', 'ztf_i':'ztf,i', @@ -1686,7 +1595,7 @@ def get_obs_str(col): 'roman_f106': 'roman,wfi,f106', 'roman_f129': 'roman,wfi,f129', 'roman_f158': 'roman,wfi,f158', - 'roman_f146': 'roman,wfi,f146', + 'roman_w146': 'roman,wfi,w146', 'roman_f213': 'roman,wfi,f213', 'roman_f184': 'roman,wfi,f184', 'rubin_g':'rubin,g', @@ -1694,14 +1603,10 @@ def get_obs_str(col): 'rubin_r':'rubin,r', 'rubin_u':'rubin,u', 'rubin_z':'rubin,z', - 'rubin_y':'rubin,y', - 'euclid_VIS':'euclid,VIS', - 'euclid_Y':'euclid,Y', - 'euclid_J':'euclid,J', - 'euclid_H':'euclid,H'} + 'rubin_y':'rubin,y'} obs_str = filt_list[name] - + return obs_str def rebin_spec(wave, specin, wavnew): @@ -1714,7 +1619,7 @@ def rebin_spec(wave, specin, wavnew): f = np.ones(len(wave)) filt = spectrum.ArraySpectralElement(wave, f, waveunits='angstrom') obs_f = obs.Observation(spec, filt, binset=wavnew, force='taper') - + return obs_f.binflux def make_isochrone_grid(age_arr, AKs_arr, dist_arr, evo_model=default_evo_model, @@ -1725,7 +1630,7 @@ def make_isochrone_grid(age_arr, AKs_arr, dist_arr, evo_model=default_evo_model, 'wfc3,ir,f153m']): """ Wrapper routine to generate a grid of isochrones of different ages, - extinctions, and distances. + extinctions, and distances. Parameters: ---------- @@ -1737,7 +1642,7 @@ def make_isochrone_grid(age_arr, AKs_arr, dist_arr, evo_model=default_evo_model, dist_arr: array Array of distances to loop over (pc) - + evo_models: SPISEA evolution object Which evolution models to use @@ -1754,7 +1659,7 @@ def make_isochrone_grid(age_arr, AKs_arr, dist_arr, evo_model=default_evo_model, Mass sampling of isochrone, relative to original mass sampling filters: dictionary - Which filters to do the synthetic photometry on + Which filters to do the synthetic photometry on """ print( '**************************************') print( 'Start generating isochrones') @@ -1799,7 +1704,7 @@ def mag_in_filter(star, filt): diff = np.diff(star_in_filter.binwave) diff = np.append(diff, diff[-1]) star_flux = np.sum(star_in_filter.binflux * diff) - + star_mag = -2.5 * math.log10(star_flux / filt.flux0) + filt.mag0 return star_mag @@ -1822,10 +1727,10 @@ def match_model_masses(isoMasses, starMasses): idx = np.where(dm_frac > 0.1)[0] indices[idx] = -1 - + return indices - + def get_evo_model_by_string(evo_model_string): return getattr(evolution, evo_model_string) @@ -1836,7 +1741,7 @@ def calc_ab_vega_filter_conversion(filt_str): AB and Vega magnitudes for a given filter: m_AB - m_vega - Note: this conversion is just the vega magnitude in + Note: this conversion is just the vega magnitude in AB system Parameters: @@ -1856,14 +1761,14 @@ def calc_ab_vega_filter_conversion(filt_str): filt_wave = filt.wave filt_mu = c / filt_wave s_filt = filt.throughput - + # Interpolate the filter function, determine what the # filter function is at the exact sampling of the # vega spectrum (in freq space) filt_interp = scipy.interpolate.interp1d(filt_mu, s_filt, kind='linear', bounds_error=False, fill_value=0) s_interp = filt_interp(vega_mu) - + # Now for the m_ab calculation mu_diff = np.diff(vega_mu) numerator = np.sum(vega_flux_mu[:-1] * s_interp[:-1] * mu_diff) @@ -1881,10 +1786,10 @@ def calc_ab_vega_filter_conversion(filt_str): # fill_value=0) #s_interp = filt_interp(vega.wave) - # Calculate the numerator + # Calculate the numerator #diff = np.diff(vega.wave) #numerator2 = np.sum((vega.wave[:-1]**2. / c) * vega.flux[:-1] * s_interp[:-1] * diff) - + # Now we need to intergrate the filter response for the denominator #denominator2 = np.sum(s_interp[:-1] * diff) @@ -1892,4 +1797,3 @@ def calc_ab_vega_filter_conversion(filt_str): #vega_mag_ab2 = -2.5 * np.log10(numerator2 / denominator2) - 48.6 return vega_mag_ab - From c100459d619b000d6866ce77a37c1d3c10d3c059 Mon Sep 17 00:00:00 2001 From: Rocketpack23 Date: Wed, 11 Feb 2026 13:45:25 -0800 Subject: [PATCH 2/4] Edited filters.py to include irtf/nsfcam --- spisea/filters.py | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/spisea/filters.py b/spisea/filters.py index 23221c6b..e2a2aee9 100755 --- a/spisea/filters.py +++ b/spisea/filters.py @@ -437,7 +437,7 @@ def get_rubin_filt(name): try: t = Table.read('{0}/rubin/{1}.dat'.format(filters_dir, name), format='ascii') except: - raise ValueError('Could not find Rubin LSST filter file {0}/total_/{1}.dat'.format(filters_dir, name)) + raise ValueError('Could not find Rubin LSST filter file {0}/rubin/{1}.dat'.format(filters_dir, name)) wavelength = t[t.keys()[0]] transmission = t[t.keys()[1]] @@ -451,6 +451,28 @@ def get_rubin_filt(name): return spectrum +def get_euclid_filt(name): + """ + Define the Euclid filters as a pysynphot spectrum object + """ + # Read in filter info + try: + t = Table.read('{0}/euclid/{1}.dat'.format(filters_dir, name), format='ascii') + except: + raise ValueError('Could not find Euclid filter file {0}/euclid/{1}.dat'.format(filters_dir, name)) + + wavelength = t[t.keys()[0]] + transmission = t[t.keys()[1]] + + # Convert wavelength to Angstroms + wavelength = wavelength * 10 + + # Make spectrum object + spectrum = pysynphot.ArrayBandpass(wavelength, transmission, waveunits='angstrom', + name='euclid_{0}'.format(name)) + + return spectrum + def get_nsfcam_filt(name): """ From e3f520c558f5a419c31ad62c3dc65fb38fac12ce Mon Sep 17 00:00:00 2001 From: Rocketpack23 Date: Wed, 11 Feb 2026 13:47:11 -0800 Subject: [PATCH 3/4] Edited synthetic.py to include irtf/nsfcam --- spisea/synthetic.py | 124 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 111 insertions(+), 13 deletions(-) diff --git a/spisea/synthetic.py b/spisea/synthetic.py index 882a18d3..4c98dc57 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -120,6 +120,12 @@ class ResolvedCluster(Cluster): produced by the cluster at the given isochrone age. Otherwise, no compact remnants are produced. + keep_low_mass_stars: boolean (default False) + If True, the cluster will not cut out stars below the isochrone grid + on initial mass. They are assigned a current mass equal to their initial + mass, a phase of 98, and no other evolutionary properties or photometry. + If False, stars below the isochrone initial mass limit are cut out. + seed: int If set to non-None, all random sampling will be seeded with the specified seed, forcing identical output. @@ -129,7 +135,7 @@ class ResolvedCluster(Cluster): True for verbose output. """ def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=True, - seed=None): + seed=None, keep_low_mass_stars=False): Cluster.__init__(self, iso, imf, cluster_mass, ifmr=ifmr, verbose=verbose, seed=seed) # Provide a user warning is random seed is set @@ -163,7 +169,7 @@ def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=True, # Trim out bad systems; specifically, stars with masses outside those provided # by the model isochrone (except for compact objects). - star_systems, compMass = self._remove_bad_systems(star_systems, compMass) + star_systems, compMass = self._remove_bad_systems(star_systems, compMass, keep_low_mass_stars) ##### # Make a table to contain all the information about companions. @@ -433,13 +439,19 @@ def _make_companions_table(self, star_systems, compMass): if len(idx) != N_comp_tot and self.verbose: print( 'Found {0:d} companions out of stellar mass range'.format(N_comp_tot - len(idx))) + # For low-mass stars and substellar objects below isochrone, assume no mass loss and set phase to 98 + low_mass_idxs = (companions['mass'] 0 + if len(idx)>0: + assert companions['mass'][idx].min() > 0 return companions - def _remove_bad_systems(self, star_systems, compMass): + def _remove_bad_systems(self, star_systems, compMass, keep_low_mass_stars): """ Helper function to remove stars with masses outside the isochrone mass range from the cluster. These stars are identified by having @@ -454,16 +466,34 @@ def _remove_bad_systems(self, star_systems, compMass): # Convert nan_to_num to avoid errors on greater than, less than comparisons star_systems_teff_non_nan = np.nan_to_num(star_systems['Teff'], nan=-99) star_systems_phase_non_nan = np.nan_to_num(star_systems['phase'], nan=-99) - if self.ifmr == None: + if (self.ifmr == None) and (not keep_low_mass_stars): + print('Remove low mass stars below grid and compact objects') # Keep only those stars with Teff assigned. idx = np.where(star_systems_teff_non_nan > 0)[0] - else: + elif not keep_low_mass_stars: + print('Remove low mass stars, keep compact objects') # Keep stars (with Teff) and any other compact objects (with phase info). idx = np.where( (star_systems_teff_non_nan > 0) | (star_systems_phase_non_nan >= 0) )[0] + elif self.ifmr == None: + print('Remove compact objects, keep low mass stars below grid') + # Keep stars (with Teff) and objects below mass grid + idx = np.where( (star_systems_teff_non_nan > 0) | ((star_systems['mass'] 0) | (star_systems_phase_non_nan >= 0) | + ((star_systems['mass']0: + self.verbose = True + print('Missing photometry for',len(comp_filters),'filter - recomputing these columns:',comp_filters) + + print('Loading stellar spectra') + # Initialize output for stellar spectra + self.spec_list = [] + # For each isochrone point, extract the synthetic photometry. + for ii in range(len(self.points['Teff'])): + # Loop is currently taking about 0.11 s per iteration + gravity = float( self.points['logg'][ii] ) + L = float( (self.points['L'][ii]*units.W).cgs / (units.erg / units.s)) # in erg/s + T = float( self.points['Teff'][ii] ) # in Kelvin + R = float( (self.points['R'][ii]*units.m).to('pc') / units.pc) # in pc + phase = int(self.points['phase'][ii]) + # Get the atmosphere model now. Wavelength is in Angstroms + # This is the time-intensive call... everything else is negligable. + # If source is a star, pull from star atmospheres. If it is a WD, + # pull from WD atmospheres + if phase == 101: + star = wd_atm_func(temperature=T, gravity=gravity, metallicity=metallicity, + verbose=False) + else: + star = atm_func(temperature=T, gravity=gravity, metallicity=metallicity, + rebin=rebin) + # Trim wavelength range down to appropriate range + star = spectrum.trimSpectrum(star, wave_range[0], wave_range[1]) + # Convert into flux observed at Earth (unreddened) + star *= (R / self.points.meta["DISTANCE"])**2 # in erg s^-1 cm^-2 A^-1 + # Redden the spectrum. This doesn't take much time at all. + red = red_law.reddening(AKs).resample(star.wave) + star *= red + # Save the final spectrum to our spec_list for later use. + self.spec_list.append(star) + + self.make_photometry(rebin=rebin, vega=vega, comp_filters=comp_filters) + return - def make_photometry(self, rebin=True, vega=vega): + def make_photometry(self, rebin=True, vega=vega, comp_filters=None): """ Make synthetic photometry for the specified filters. This function udpates the self.points table to include new columns with the @@ -1035,9 +1110,13 @@ def make_photometry(self, rebin=True, vega=vega): npoints = len(self.points) verbose_fmt = 'M = {0:7.3f} Msun T = {1:5.0f} K m_{2:s} = {3:4.2f}' + #Calculate all filters, or select filters + if comp_filters is None: + comp_filters = self.filters + # Loop through the filters, get filter info, make photometry for # all stars in this filter. - for ii in self.filters: + for ii in comp_filters: prt_fmt = 'Starting filter: {0:s} Elapsed time: {1:.2f} seconds' print( prt_fmt.format(ii, time.time() - startTime)) @@ -1094,6 +1173,13 @@ def check_save_file(self, evo_model, atm_func, red_law): (tmp.meta['REDLAW'] == red_law.name) ): out_bool = True + # Check model version if it was logged + if 'EVOMODELVERSION' in tmp.meta: + if tmp.meta['EVOMODELVERSION']!=evo_model.model_version_name: + out_bool=False + else: + print(f"No version information found for evolution model {type(evo_model).__name__}.") + return out_bool def plot_CMD(self, mag1, mag2, savefile=None): @@ -1277,6 +1363,7 @@ def __init__(self, logAge, distance, evo_model=default_evo_model, tab.meta['ATMFUNC'] = atm_func.__name__ tab.meta['EVOMODEL'] = type(evo_model).__name__ + tab.meta['EVOMODELVERSION'] = evo_model.model_version_name tab.meta['LOGAGE'] = logAge tab.meta['DISTANCE'] = distance tab.meta['WAVEMIN'] = wave_range[0] @@ -1462,9 +1549,12 @@ def get_filter_info(name, vega=vega, rebin=True): elif name.startswith('rubin'): filt = filters.get_rubin_filt(filterName) + elif name.startswith('euclid'): + filt = filters.get_euclid_filt(filterName) + elif name.startswith('nsfcam'): filt = filters.get_nsfcam_filt(filterName) - + else: # Otherwise, look for the filter info in the cdbs/mtab and cdbs/comp files try: @@ -1582,6 +1672,10 @@ def get_obs_str(col): 'jg_J': 'jg,J', 'jg_H': 'jg,H', 'jg_K': 'jg,K', 'nirc1_K':'nirc1,K', 'nirc1_H':'nirc1,H', 'naco_J':'naco,J', 'naco_H':'naco,H', 'naco_Ks':'naco,Ks', + 'naco_IB_2.00': 'naco,IB_2.00', 'naco_IB_2.03':'naco,IB_2.03', 'naco_IB_2.06':'naco,IB_2.06', + 'naco_IB_2.24':'naco,IB_2.24', 'naco_IB_2.27':'naco,IB_2.27', + 'naco_IB_2.30':'naco,IB_2.30', 'naco_IB_2.33':'naco,IB_2.33', + 'naco_IB_2.36':'naco,IB_2.36', 'ukirt_J':'ukirt,J', 'ukirt_H':'ukirt,H', 'ukirt_K':'ukirt,K', 'ctio_osiris_H': 'ctio_osiris,H', 'ctio_osiris_K': 'ctio_osiris,K', 'ztf_g':'ztf,g', 'ztf_r':'ztf,r', 'ztf_i':'ztf,i', @@ -1595,7 +1689,7 @@ def get_obs_str(col): 'roman_f106': 'roman,wfi,f106', 'roman_f129': 'roman,wfi,f129', 'roman_f158': 'roman,wfi,f158', - 'roman_w146': 'roman,wfi,w146', + 'roman_f146': 'roman,wfi,f146', 'roman_f213': 'roman,wfi,f213', 'roman_f184': 'roman,wfi,f184', 'rubin_g':'rubin,g', @@ -1603,7 +1697,11 @@ def get_obs_str(col): 'rubin_r':'rubin,r', 'rubin_u':'rubin,u', 'rubin_z':'rubin,z', - 'rubin_y':'rubin,y'} + 'rubin_y':'rubin,y', + 'euclid_VIS':'euclid,VIS', + 'euclid_Y':'euclid,Y', + 'euclid_J':'euclid,J', + 'euclid_H':'euclid,H'} obs_str = filt_list[name] From 1d9cdc2ba4fb8f6660b7064e6f7d2460a5466e7a Mon Sep 17 00:00:00 2001 From: Rocketpack23 Date: Thu, 5 Mar 2026 11:18:24 -0800 Subject: [PATCH 4/4] updated filters.rst and dictionary in synthetic.py --- docs/filters.rst | 26 +++++++++++++++++--------- spisea/synthetic.py | 5 +++-- 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/docs/filters.rst b/docs/filters.rst index 918f27f4..6efb28fc 100644 --- a/docs/filters.rst +++ b/docs/filters.rst @@ -7,16 +7,16 @@ Photometric Filters The user can specify what filters are used for synthetic photometry when defining the :ref:`isochrone_objects`. Each filter is identified by a unique string, and an array of such strings -are passed into the Isochrone call. +are passed into the Isochrone call. For example:: - + # Use the HST WFC3-IR F127M and F153M filters, along with NIRC2 Kp filt_list = ['wfc3,ir,f127m', 'wfc3,ir,f153m', 'nirc2,Kp'] my_iso = synthetic.IsochronePhot(logAge, AKs, dist, metallicity=0, evo_model=evo_model, atm_func=atm_func, red_law=red_law, filters=filt_list) - + These strings follow the format ``,``. Note that there is no space after the comma, and case matters. @@ -40,7 +40,7 @@ Available filters: * JWST * Keck NIRC * Keck NIRC2 -* NACO +* NACO * PanStarrs 1 * Roman Space Telescope * UKIRT @@ -48,11 +48,11 @@ Available filters: * VISTA * ZTF - + Filter Sets ------------ - + **2MASS** `Two-Micron Sky Survey `_ @@ -83,7 +83,7 @@ Example: ``'decam,r'`` **Euclid** -Euclid space telescope `NISP filters `_ +Euclid space telescope `NISP filters `_ and `VIS single filter `_ Filters: VIS, Y, J, H @@ -146,10 +146,10 @@ Example: ``'jg,K'`` JWST NIRCam filters, downloaded from `NIRCam website `_. The filter functions in the nircam_throughputs/modAB_mean/nrc_plus_ote folder is used. -Filters: F070W, F090W, F115W, F140M, F150W, F150W2, F162M, F164N, F182M, F187N, F200W, F210M, F212N, F250M, F277W, F300M, F322W2, F323N, F335M, F356W, F360M, F405N, F410M, F430M, F444W, F460M, F466N, F470N, F480M +Filters: F070W, F090W, F115W, F140M, F150W, F150W2, F162M, F164N, F182M, F187N, F200W, F210M, F212N, F250M, F277W, F300M, F322W2, F323N, F335M, F356W, F360M, F405N, F410M, F430M, F444W, F460M, F466N, F470N, F480M Example: ``'jwst,F356W'`` - + **Keck NIRC** @@ -237,3 +237,11 @@ Example: ``'vista,Y'`` Filters: g, r, i Example: ``'ztf,g'`` + +**IRTF** + +`IRTF NSFCam `_ + +Filters: L + +Example: ``'nsfcam,L'`` diff --git a/spisea/synthetic.py b/spisea/synthetic.py index 4c98dc57..171d0bda 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -1554,7 +1554,7 @@ def get_filter_info(name, vega=vega, rebin=True): elif name.startswith('nsfcam'): filt = filters.get_nsfcam_filt(filterName) - + else: # Otherwise, look for the filter info in the cdbs/mtab and cdbs/comp files try: @@ -1701,7 +1701,8 @@ def get_obs_str(col): 'euclid_VIS':'euclid,VIS', 'euclid_Y':'euclid,Y', 'euclid_J':'euclid,J', - 'euclid_H':'euclid,H'} + 'euclid_H':'euclid,H', + 'nsfcam_L':'nsfcam,L'} obs_str = filt_list[name]