diff --git a/py/desispec/fluxcalibration.py b/py/desispec/fluxcalibration.py
index bdd474ab7..5989361d1 100644
--- a/py/desispec/fluxcalibration.py
+++ b/py/desispec/fluxcalibration.py
@@ -18,6 +18,7 @@
 from desitarget.targets import main_cmx_or_sv
 from desispec.fiberfluxcorr import flat_to_psf_flux_correction,psf_to_fiber_flux_correction
 from desispec.gpu import is_gpu_available, NoGPU
+from desispec.maskbits import fibermask
 import scipy, scipy.sparse, scipy.ndimage
 import sys
 import time
diff --git a/py/desispec/test/test_fiberbitmask.py b/py/desispec/test/test_fiberbitmask.py
index c0cede27d..5874efa1a 100644
--- a/py/desispec/test/test_fiberbitmask.py
+++ b/py/desispec/test/test_fiberbitmask.py
@@ -73,11 +73,10 @@ def test_ambiguous_maskbits(self):
         #- POORPOSITION is ok for flats, sky, and fluxcalib; but bad for stdstars
         self.check_mask('POORPOSITION', ok_steps=['flat', 'sky', 'fluxcalib'], bad_steps=['stdstar'])
 
-        #- NEARCHARGETRAP is informative; treated as ok for everyone including sky
+        #- NEARCHARGETRAP and VARIABLETHRU are informative for fiberbitmasking;
+        #- i.e. they don't trigger masking fibers
         #- TODO: it's actually bad for faint targets and sky for a single amp, but we structurally
         #- don't have a way to encode that in FIBERSTATUS (fiber not CCD or amp)
         self.check_mask('NEARCHARGETRAP', ok_steps=['flat', 'sky', 'stdstar', 'fluxcalib'], bad_steps=[])
+        self.check_mask('VARIABLETHRU', ok_steps=['flat', 'sky', 'stdstar', 'fluxcalib'], bad_steps=[])
 
-        #- VARIABLETHRU is ok for flats because otherwise we'd block the entire fiber,
-        #- and ok to at least attempt to flux calibrate it, but it shouldn't be used for sky or stdstars
-        self.check_mask('VARIABLETHRU', ok_steps=['flat', 'fluxcalib'], bad_steps=['sky', 'stdstar'])
diff --git a/py/desispec/test/test_flux_calibration.py b/py/desispec/test/test_flux_calibration.py
index 259a3f34c..7b19644fe 100644
--- a/py/desispec/test/test_flux_calibration.py
+++ b/py/desispec/test/test_flux_calibration.py
@@ -10,7 +10,7 @@
 import numpy as np
 #import scipy.sparse
 
-#from desispec.maskbits import specmask
+from desispec.maskbits import fibermask
 from desispec.frame import Frame
 from desispec.fluxcalibration import normalize_templates
 from desispec.fluxcalibration import FluxCalib
@@ -204,6 +204,7 @@ def test_isStdStar(self):
         fm['CMX_TARGET'] = np.zeros(10, dtype=int)
         fm['CMX_TARGET'][0:2] = cmx_mask.STD_FAINT
         fm['CMX_TARGET'][2:4] = cmx_mask.SV0_STD_FAINT
+        fm['FIBERSTATUS'] = 0
         self.assertEqual(main_cmx_or_sv(fm)[2], 'cmx')
         self.assertEqual(np.count_nonzero(isStdStar(fm)), 4)
 
@@ -213,6 +214,7 @@ def test_isStdStar(self):
         fm['SV1_MWS_TARGET'] = np.zeros(10, dtype=int)
         fm['SV1_DESI_TARGET'][0:2] = sv1_desi_mask.STD_FAINT
         fm['SV1_MWS_TARGET'][2:4] = sv1_mws_mask.GAIA_STD_FAINT
+        fm['FIBERSTATUS'] = 0
         self.assertEqual(main_cmx_or_sv(fm)[2], 'sv1')
         self.assertEqual(np.count_nonzero(isStdStar(fm)), 4)
 
@@ -224,11 +226,21 @@ def test_isStdStar(self):
         fm['DESI_TARGET'][2:4] = desi_mask.STD_BRIGHT
         fm['DESI_TARGET'][4:6] |= desi_mask.MWS_ANY
         fm['MWS_TARGET'][4:6] = sv1_mws_mask.GAIA_STD_FAINT
+        fm['FIBERSTATUS'] = 0
         self.assertEqual(main_cmx_or_sv(fm)[2], 'main')
         self.assertEqual(np.count_nonzero(isStdStar(fm)), 6)
         self.assertEqual(np.count_nonzero(isStdStar(fm, bright=False)), 4)
         self.assertEqual(np.count_nonzero(isStdStar(fm, bright=True)), 2)
 
+        #- VARIABLETHRU should be excluded
+        fm['FIBERSTATUS'] = 0
+        ii = isStdStar(fm)
+        fm['FIBERSTATUS'][0] = fibermask.VARIABLETHRU
+        jj = isStdStar(fm)
+        self.assertTrue(ii[0])
+        self.assertFalse(jj[0])
+
+
 
     def test_main(self):
         pass
diff --git a/py/desispec/test/test_sky.py b/py/desispec/test/test_sky.py
index 08953641b..9ec1ebfdd 100644
--- a/py/desispec/test/test_sky.py
+++ b/py/desispec/test/test_sky.py
@@ -5,9 +5,11 @@
 import unittest
 
 import numpy as np
-from desispec.sky import compute_sky, subtract_sky, SkyModel
+from astropy.table import Table
+from desispec.sky import compute_sky, subtract_sky, SkyModel, get_sky_fibers
 from desispec.resolution import Resolution
 from desispec.frame import Frame
+from desispec.maskbits import fibermask
 import desispec.io
 
 import desispec.scripts.sky as skyscript
@@ -155,6 +157,29 @@ def test_sky_slice(self):
         self.assertEqual(sky2.throughput_corrections, None)
         self.assertEqual(sky2.nrej, sky1.nrej)
 
+    def test_get_sky_fibers(self):
+        """Test desispec.sky.get_sky_fibers"""
+        fibermap = Table()
+        fibermap['TARGETID'] = np.arange(5)
+        fibermap['OBJTYPE'] = ['TGT', 'SKY', 'TGT', 'SKY', 'TGT']
+        fibermap['FIBERSTATUS'] = 0
+
+        #- Test default behavior
+        skyfibers = list(get_sky_fibers(fibermap))
+        self.assertEqual(skyfibers, [1,3])
+
+        #- Test overrides
+        skyfibers = list(get_sky_fibers(fibermap, override_sky_targetids=[0,1]))
+        self.assertEqual(skyfibers, [0,1])
+
+        skyfibers = list(get_sky_fibers(fibermap, exclude_sky_targetids=[0,1]))
+        self.assertEqual(skyfibers, [3,])
+
+        #- SKY fibers with VARIABLETHRU should be excluded
+        fibermap['FIBERSTATUS'][3] = fibermask.VARIABLETHRU
+        skyfibers = list(get_sky_fibers(fibermap))
+        self.assertEqual(skyfibers, [1,])
+
     def test_main(self):
         pass