diff --git a/gwinc/noise/coatingthermal.py b/gwinc/noise/coatingthermal.py
index 29e641e531b742df4ebee1bc4b18f9857fa48c69..41a642a43efd3e42fa1476ab819cea37a2b772d4 100644
--- a/gwinc/noise/coatingthermal.py
+++ b/gwinc/noise/coatingthermal.py
@@ -3,7 +3,7 @@
 '''
 from __future__ import division, print_function
 import numpy as np
-from numpy import pi, sum, zeros, exp, real, imag, sqrt, sin, cos, sinh, cosh, polyfit, roots, min, ceil, log
+from numpy import pi, exp, real, imag, sqrt, sin, cos, sinh, cosh, ceil, log
 
 from .. import const
 from ..const import BESSEL_ZEROS as zeta
@@ -53,10 +53,10 @@ def coating_brownian(f, materials, wavelength, wBeam, dOpt):
     nL = coat.Indexlown
 
     # make vectors of material properties
-    nN = zeros(len(dOpt))
-    yN = zeros(len(dOpt))
-    pratN = zeros(len(dOpt))
-    phiN = zeros(len(dOpt))
+    nN = np.zeros(len(dOpt))
+    yN = np.zeros(len(dOpt))
+    pratN = np.zeros(len(dOpt))
+    phiN = np.zeros(len(dOpt))
 
     # make simple alternating structure (low-index, high-index doublets)
     # (adapted from the more general calculation in
@@ -77,7 +77,7 @@ def coating_brownian(f, materials, wavelength, wBeam, dOpt):
 
     # geometrical thickness of each layer and total
     dGeo = wavelength * np.asarray(dOpt) / nN
-    #dCoat = sum(dGeo)
+    #dCoat = np.sum(dGeo)
 
     ###################################
     # Brownian
@@ -107,12 +107,12 @@ def coating_brownian(f, materials, wavelength, wBeam, dOpt):
     else:
         SbrZ = (4 * kBT / (pi * wBeam**2 * w)) * \
                (1 - pratsub - 2 * pratsub**2) / Ysub
-        SbrZ *= (sum(dGeo[::2] * brLayer[::2] * phiN[::2]) * (f / 100)**(coat.Philown_slope) +
-                 sum(dGeo[1::2] * brLayer[1::2] * phiN[1::2]) * (f / 100)**(coat.Phihighn_slope))
+        SbrZ *= (np.sum(dGeo[::2] * brLayer[::2] * phiN[::2]) * (f / 100)**(coat.Philown_slope) +
+                 np.sum(dGeo[1::2] * brLayer[1::2] * phiN[1::2]) * (f / 100)**(coat.Phihighn_slope))
 
     # for the record: evaluate summation in eq 1 of PhysRevD.91.042002
     # normalized by total coating thickness to make it unitless
-    #cCoat = sum(dGeo * brLayer * phiN) / dCoat
+    #cCoat = np.sum(dGeo * brLayer * phiN) / dCoat
 
     return SbrZ
 
@@ -339,17 +339,17 @@ def getCoatLayers(materials, wavelength, dOpt):
              * ( ((1 + sigC) / (1 + sigS)) + (1 - 2 * sigS) * Y_C / Y_S )
         return ce
 
-    aLayer = zeros(Nlayer)
+    aLayer = np.zeros(Nlayer)
     aLayer[::2] = alphaL * getExpansionRatio(Y_L, sigL, Y_S, sigS)
     aLayer[1::2] = alphaH * getExpansionRatio(Y_H, sigH, Y_S, sigS)
 
     # and beta
-    bLayer = zeros(Nlayer)
+    bLayer = np.zeros(Nlayer)
     bLayer[::2] = betaL
     bLayer[1::2] = betaH
 
     # and refractive index
-    nLayer = zeros(Nlayer)
+    nLayer = np.zeros(Nlayer)
     nLayer[::2] = nL
     nLayer[1::2] = nH
 
@@ -357,7 +357,7 @@ def getCoatLayers(materials, wavelength, dOpt):
     dLayer = wavelength * np.asarray(dOpt) / nLayer
 
     # and sigma correction
-    sLayer = zeros(Nlayer)
+    sLayer = np.zeros(Nlayer)
     sLayer[::2] = alphaL * (1 + sigL) / (1 - sigL)
     sLayer[1::2] = alphaH * (1 + sigH) / (1 - sigH)
 
@@ -396,9 +396,9 @@ def getCoatAvg(materials, wavelength, dOpt):
     junk1, junk2, junk3, dLayer, junk4 = getCoatLayers(materials, wavelength, dOpt)
 
     # heat capacity
-    dc = sum(dLayer)
-    dL = sum(dLayer[::2])
-    dH = sum(dLayer[1::2])
+    dc = np.sum(dLayer)
+    dL = np.sum(dLayer[::2])
+    dH = np.sum(dLayer[1::2])
     Cc = (C_L * dL + C_H * dH) / dc
 
     # thermal diffusivity
@@ -452,10 +452,10 @@ def getCoatTOPhase(nIn, nOut, nLayer, dOpt, aLayer, bLayer, sLayer):
     dphi_dd = 4 * pi * dcdp
 
     # thermo-refractive coupling
-    dphi_TR = sum(dphi_dd * (bLayer + sLayer * nLayer) * dGeo)
+    dphi_TR = np.sum(dphi_dd * (bLayer + sLayer * nLayer) * dGeo)
 
     # thermo-elastic
-    dphi_TE = 4 * pi * sum(aLayer * dGeo)
+    dphi_TE = 4 * pi * np.sum(aLayer * dGeo)
 
     # total
     dphi_dT = dphi_TR + dphi_TE
@@ -482,8 +482,8 @@ def getCoatFiniteCorr(materials, wavelength, wBeam, dOpt):
 
     """
     # parameter extraction
-    R = materials.MassRadius      #substrate radius
-    H = materials.MassThickness   #substrate thickness
+    R = materials.MassRadius
+    H = materials.MassThickness
 
     alphaS = materials.Substrate.MassAlpha
     C_S = materials.Substrate.MassCM * materials.Substrate.MassDensity
@@ -503,8 +503,8 @@ def getCoatFiniteCorr(materials, wavelength, wBeam, dOpt):
     nH = materials.Coating.Indexhighn
 
     # coating sums
-    dL = wavelength * sum(dOpt[::2]) / nL
-    dH = wavelength * sum(dOpt[1::2]) / nH
+    dL = wavelength * np.sum(dOpt[::2]) / nL
+    dH = wavelength * np.sum(dOpt[1::2]) / nH
     dc = dH + dL
 
     # AVERAGE SPECIFIC HEAT (simple volume average for coating)
@@ -543,8 +543,8 @@ def getCoatFiniteCorr(materials, wavelength, wBeam, dOpt):
          (1 + sigS) * (1 - Qm)**2 / ((1 - Qm)**2 - 4 * km**2 * H**2 * Qm)
 
     # eq 90 and 91
-    S1 = (12 * R**2 / H**2) * sum(pm / zeta**2)
-    S2 = sum(pm**2 * Lm**2)
+    S1 = (12 * R**2 / H**2) * np.sum(pm / zeta**2)
+    S2 = np.sum(pm**2 * Lm**2)
     P = (Xr - 2 * sigS * Yr - Cr + S1 * (Cr - Yr * (1 - sigS)))**2 + S2
 
     # eq 60 and 70
@@ -571,13 +571,13 @@ def getCoatDopt(materials, T, dL, dCap=0.5):
     def getTrans(materials, Ndblt, dL, dH, dCap, dTweak):
 
         # the optical thickness vector
-        dOpt = zeros(2 * Ndblt)
+        dOpt = np.zeros(2 * Ndblt)
         dOpt[0] = dCap
         dOpt[1::2] = dH
         dOpt[2::2] = dL
 
         N = dTweak.size
-        T = zeros(N)
+        T = np.zeros(N)
         for n in range(N):
             dOpt[-1] = dTweak[n]
             r = getCoatRefl(materials, dOpt)[0]
@@ -590,13 +590,13 @@ def getCoatDopt(materials, T, dL, dCap=0.5):
 
         # tweak bottom layer
         Tn = getTrans(materials, Ndblt, dL, dH, dCap, dScan)
-        pf = polyfit(dScan, Tn - T, Nfit)
-        rts = roots(pf)
+        pf = np.polyfit(dScan, Tn - T, Nfit)
+        rts = np.roots(pf)
         if not any((imag(rts) == 0) & (rts > 0)):
             dTweak = None
             Td = 0
             return dTweak, Td
-        dTweak = real(min(rts[(imag(rts) == 0) & (rts > 0)]))
+        dTweak = real(np.min(rts[(imag(rts) == 0) & (rts > 0)]))
 
         # compute T for this dTweak
         Td = getTrans(materials, Ndblt, dL, dH, dCap, np.array([dTweak]))
@@ -664,7 +664,7 @@ def getCoatDopt(materials, T, dL, dCap=0.5):
 
     ########################
     # return dOpt vector
-    dOpt = zeros(2 * Ndblt)
+    dOpt = np.zeros(2 * Ndblt)
     dOpt[0] = dCap
     dOpt[1::2] = dH
     dOpt[2::2] = dL
@@ -692,7 +692,7 @@ def getCoatRefl(materials, dOpt):
     Nlayer = len(dOpt)
 
     # refractive index of input, coating, and output materials
-    nAll = zeros(Nlayer + 2)
+    nAll = np.zeros(Nlayer + 2)
     nAll[0] = 1  # vacuum input
     nAll[1::2] = nL
     nAll[2::2] = nH
@@ -730,8 +730,8 @@ def getCoatRefl2(nIn, nOut, nLayer, dOpt):
     r = (nAll[:-1] - nAll[1:]) / (nAll[:-1] + nAll[1:])
 
     # combine reflectivities
-    rbar = zeros(r.size, dtype=complex)
-    ephi = zeros(r.size, dtype=complex)
+    rbar = np.zeros(r.size, dtype=complex)
+    ephi = np.zeros(r.size, dtype=complex)
 
     # round-trip phase in each layer
     ephi[0] = 1
diff --git a/gwinc/noise/quantum.py b/gwinc/noise/quantum.py
index 65651c1e9869440b4326cf902d1befad35ffc2a1..c019e2a006e25f7b69e19ed5a08ea58bbaa01a9a 100644
--- a/gwinc/noise/quantum.py
+++ b/gwinc/noise/quantum.py
@@ -3,7 +3,7 @@
 '''
 from __future__ import division
 import numpy as np
-from numpy import pi, sqrt, arctan, sin, cos, exp, size, ones, zeros, log10, conj, sum
+from numpy import pi, sqrt, arctan, sin, cos, exp, log10, conj
 import logging
 
 from .. import const
@@ -225,8 +225,8 @@ def shotrad(f, ifo):
     # and compute the final noise
     def HDnoise(eta):
         vHD = np.array([[sin(eta), cos(eta)]])
-        n = coeff * np.squeeze(sum(abs(getProdTF(vHD, Mnoise))**2, axis=1)) / \
-            np.squeeze(sum(abs(getProdTF(vHD, Msig))**2, axis=1))
+        n = coeff * np.squeeze(np.sum(abs(getProdTF(vHD, Mnoise))**2, axis=1)) / \
+            np.squeeze(np.sum(abs(getProdTF(vHD, Msig))**2, axis=1))
         return n
 
     if etaRMS <= 0:
@@ -243,10 +243,10 @@ def shotrad(f, ifo):
     #    sym(M) = real(M * M')
     #
     # it is also the same as taking the sum of the squared directly
-    #   n = zeros(1, numel(f));
+    #   n = np.zeros(1, numel(f));
     #   for k = 1:numel(f)
-    #     n(k) = coeff(k) * sum(abs((vHD * Msqz(:,:,k))).^2) ./ ...
-    #       sum(abs((vHD * Msig(:,:,k))).^2);
+    #     n(k) = coeff(k) * np.sum(abs((vHD * Msqz(:,:,k))).^2) ./ ...
+    #       np.sum(abs((vHD * Msig(:,:,k))).^2);
     #   end
 
     return n * ifo.Infrastructure.Length**2
@@ -335,12 +335,12 @@ def shotradSignalRecycled(f, ifo):
     ID = np.array([[np.ones(nf), np.zeros(nf)], [np.zeros(nf), np.ones(nf)]])
 
     # transfer matrices for dark port input and signal field
-    Mifo = zeros((2,2,nf), dtype=complex)
-    Msig = zeros((2,1,nf), dtype=complex)
+    Mifo = np.zeros((2, 2, nf), dtype=complex)
+    Msig = np.zeros((2, 1, nf), dtype=complex)
 
     # transfer matrices for SEC and arm loss fields
-    Mp = zeros((2,2,nf), dtype=complex)
-    Mn = zeros((2,2,nf), dtype=complex)
+    Mp = np.zeros((2, 2, nf), dtype=complex)
+    Mn = np.zeros((2, 2, nf), dtype=complex)
 
     # SRC rotation matrix
     SEr = np.array([[np.tile(cos(phi), nf), np.tile(sin(phi), nf)],
@@ -368,7 +368,7 @@ def shotradSignalRecycled(f, ifo):
 
     # the following code is generated by compile_ARM_RES_TF()
     # and is equivalent to the following:
-    # RES = zeros((2,2,nf), dtype=complex)
+    # RES = np.zeros((2,2,nf), dtype=complex)
     # RES = np.linalg.pinv(ID.transpose((2,0,1)) - rITM * ARM.transpose((2,0,1))).transpose((1,2,0))
     x0 = exp_2jOmegaL_c*rITM*tArm
     x1 = -x0 + 1
@@ -390,7 +390,7 @@ def shotradSignalRecycled(f, ifo):
 
     # the following code is generated by compile_SEC_RES_TF()
     # and is equivalent to the following:
-    # RES = zeros((2,2,nf), dtype=complex)
+    # RES = np.zeros((2,2,nf), dtype=complex)
     # RES = np.linalg.pinv(ID.transpose((2,0,1)) + rho * SEC.transpose((2,0,1))).transpose((1,2,0))
     x0 = cos(phi)
     x1 = exp_2jOmegaL_c*tArm*(R + T)
@@ -433,7 +433,7 @@ def shotradSignalRecycled(f, ifo):
     Msig[1, 0, :] = -Msig[1, 0, :]
 
     def adapt_to_gwinc(Mx):
-        My = zeros(Mx.shape, dtype=complex)
+        My = np.zeros(Mx.shape, dtype=complex)
         My[0, 0, :] = Mx[1, 1, :]
         My[1, 1, :] = Mx[0, 0, :]
         My[0, 1, :] = -Mx[1, 0, :]
@@ -686,7 +686,7 @@ def sqzFilterCavityChain(f, params, Mn):
 
     """
     # make an identity TF
-    Mc = make2x2TF(ones(f.shape), 0, 0, 1)
+    Mc = make2x2TF(np.ones(f.shape), 0, 0, 1)
 
     # loop through the filter cavites
     for k in range(params.size):
@@ -800,7 +800,7 @@ def sqzFilterCavity(f, Lcav, Ti, Te, Lrt, fdetune, MinR, MinT=1):
 
     # reflected fields
     if MinR == []:
-        Mnoise = zeros((2, 0, f.size))
+        Mnoise = np.zeros((2, 0, f.size))
     else:
         if np.isscalar(MinR):
             Mnoise = Mr * MinR
diff --git a/gwinc/noise/suspensionthermal.py b/gwinc/noise/suspensionthermal.py
index c5b4910ac6e6d6c783e0ebdf695cb8bf7b32be4a..34c7cbf9c1fa53cef60df2a56f63d10163aa9d18 100644
--- a/gwinc/noise/suspensionthermal.py
+++ b/gwinc/noise/suspensionthermal.py
@@ -3,7 +3,7 @@
 '''
 from __future__ import division
 import numpy as np
-from numpy import pi, imag, zeros
+from numpy import pi, imag
 
 from .. import const
 
@@ -30,7 +30,7 @@ def suspension_thermal(f, sus):
     # and vertical to beamline coupling angle
     theta = sus.VHCoupling.theta
 
-    noise = zeros((1, f.size))
+    noise = np.zeros((1, f.size))
 
     # if the temperature is uniform along the suspension
     if 'Temp' in sus:
@@ -67,7 +67,7 @@ def suspension_thermal(f, sus):
         # Thermal Noise Calculation
         ##########################################################
 
-        dxdF = zeros(hForce.shape, dtype=complex)
+        dxdF = np.zeros(hForce.shape, dtype=complex)
         for n, stage in enumerate(reversed(sus.Stage)):
             # add up the contribution from each stage