.gitlab-ci.yml

@@ -17,14 +17,15 @@ test:
   - python3 -m gwinc aLIGO -s aLIGO.png
   - python3 -m gwinc A+ -s A+.png
   - python3 -m gwinc Voyager -s Voyager.png
-  - python3 -m gwinc.test aLIGO -p -s aLIGO_test.png || true
-  - python3 -m gwinc.test A+ -p -s A+_test.png || true
+  - python3 -m gwinc.test aLIGO -t 10e-6 -k Seismic -k "Substrate Thermo-Elastic" -p -s aLIGO_test.png
+  - python3 -m gwinc.test -t 100e-6 -k Seismic -k "Substrate Thermo-Elastic" A+ -p -s A+_test.png
   after_script:
   - rm gitID.txt
   cache:
     key: "$CI_PROJECT_NAMESPACE:$CI_PROJECT_NAME:$CI_JOB_NAME"
     untracked: true
   artifacts:
+    when: always
     expire_in: 4w
     paths:
     - aLIGO.png

gwinc/noise/substratethermal.py

@@ -9,7 +9,31 @@ from ..const import BESSEL_ZEROS as zeta
 from ..const import J0M as j0m
 
 
-def carrierdensity(f, ifo):
+def carrierdensity_adiabatic(f, ifo):
+    """strain noise psd arising from charge carrier density
+    fluctuations in ITM substrate (for semiconductor substrates)."""
+
+    Omega = 2*pi*f
+    H = ifo.Materials.MassThickness
+    gammaElec = ifo.Materials.Substrate.ElectronIndexGamma
+    gammaHole = ifo.Materials.Substrate.HoleIndexGamma
+    diffElec = ifo.Materials.Substrate.ElectronDiffusion
+    diffHole = ifo.Materials.Substrate.HoleDiffusion
+    cdDens = ifo.Materials.Substrate.CarrierDensity
+    r0 = ifo.Optics.ITM.BeamRadius/np.sqrt(2)
+    T = ifo.Optics.ITM.Transmittance
+    L = ifo.Infrastructure.Length
+    Finesse = 2*pi/T
+    gPhase = 2*Finesse/pi
+
+    psdElec = 4*H*gammaElec**2*cdDens*diffElec/(pi*r0**4*Omega**2) # units are meters
+    psdHole = 4*H*gammaHole**2*cdDens*diffHole/(pi*r0**4*Omega**2) # units are meters
+    psdMeters = 2 * (psdElec + psdHole) # electrons and holes for two ITMs
+    n = psdMeters / (gPhase*L)**2
+    return n
+
+
+def carrierdensity_exact(f, ifo):
     """Strain noise arising from charge carrier density fluctuations in ITM substrate
 
     For semiconductor substrates
@@ -37,7 +61,7 @@ def carrierdensity(f, ifo):
     gPhase = 2*Finesse/pi
 
     omega = 2*pi*f
-    
+
     def integrand(k, om, D):
         return D * k**3 * exp(-k**2 * w**2/4) / (D**2 * k**4 + om**2)
     
@@ -58,8 +82,34 @@ def carrierdensity(f, ifo):
 
     return n
 
+carrierdensity = carrierdensity_adiabatic
 
-def thermorefractiveITM(f, ifo):
+
+def thermorefractiveITM_adiabatic(f, ifo):
+    """strain noise psd arising from thermorefractive
+    fluctuations in ITM substrate (for semiconductor substrates)."""
+
+    Omega = 2*pi*f
+    H = ifo.Materials.MassThickness
+    beta = ifo.Materials.Substrate.dndT
+    kappa = ifo.Materials.Substrate.MassKappa
+    rho = ifo.Materials.Substrate.MassDensity
+    C = ifo.Materials.Substrate.MassCM
+    Temp = ifo.Materials.Substrate.Temp
+    kBT = scipy.constants.k * Temp
+    r0 = ifo.Optics.ITM.BeamRadius/np.sqrt(2)
+    L = ifo.Infrastructure.Length
+    T = ifo.Optics.ITM.Transmittance
+    Finesse = 2*pi/T
+    gPhase = 2*Finesse/pi
+
+    psd = 4*H*beta**2*kappa*kBT*Temp/(pi*r0**4*Omega**2*(rho*C)**2) # units are meters
+    psdMeters = 2*psd # two ITMs
+    n = psdMeters / (gPhase*L)**2
+    return n
+
+
+def thermorefractiveITM_exact(f, ifo):
     """Strain noise from thermorefractive fluctuations in ITM substrate
 
     For semiconductor substrates.
@@ -105,6 +155,8 @@ def thermorefractiveITM(f, ifo):
 
     return n
 
+thermorefractiveITM = thermorefractiveITM_adiabatic
+
 
 def subbrownian(f, ifo):
     """Strain noise from the Brownian thermal noise due to substrate mechanical loss

gwinc/precomp.py

@@ -261,7 +261,7 @@ def dhdl(f, armlen):
 
     """
     c = scipy.constants.c
-    nu_small = 0.05
+    nu_small = 15*pi/180
     omega_arm = pi * f * armlen / c
     omega_arm_f = (1 - sin(nu_small)) * pi * f * armlen / c
     omega_arm_b = (1 + sin(nu_small)) * pi * f * armlen / c

gwinc/test/__main__.py

@@ -23,13 +23,6 @@ except ImportError:
 FLO = 5
 FHI = 6000
 NPOINTS = 3000
-FRACTIONAL_TOLERANCE = 0.01
-# comparisons to skip
-SKIP = [
-    # 'Seismic',
-    # 'Suspension Thermal',
-    # 'Newtonian Gravity',
-    ]
 
 
 def path_hash(path):
@@ -58,6 +51,8 @@ def main():
     parser.add_argument('--plot', '-p', action='store_true', help='plot differences')
     parser.add_argument('--save', '-s', help='save plot to file')
     parser.add_argument('--recalc', '-r', action='store_true', help='recalculate all traces')
+    parser.add_argument('--tolerance', '-t', help='fractional tolerance', type=float, default=1e-6)
+    parser.add_argument('--skip', '-k', action='append', help='traces to skip comparing')
     parser.add_argument('IFO', help='IFO name or description file')
     args = parser.parse_args()
 
@@ -145,11 +140,14 @@ pygwinc: {fom:.2f} Mpc""".format(
     ##############################
     # find differences
 
+    skip = args.skip
+    fractional_tolerance = args.tolerance
+
     diffs = {}
     for name, noise in noises.items():
         if name in ['Freq']:
             continue
-        if name in SKIP:
+        if skip and name in skip:
             logging.warning("SKIPPING TEST: '{}'".format(name))
             continue
 
@@ -167,11 +165,11 @@ pygwinc: {fom:.2f} Mpc""".format(
         diff = np.sqrt(mn) - np.sqrt(pn)
         frac = abs(diff / np.sqrt(pn))
 
-        if max(frac) < FRACTIONAL_TOLERANCE:
+        if max(frac) < fractional_tolerance:
             continue
 
-        logging.warning("EXCESSIVE DIFFERENCE: {:{w}} {:6.1f}%".format(
-            name, max(frac)*100, w=max([len(n) for n in noises])))
+        logging.warning("EXCESSIVE DIFFERENCE: {:{w}} {:6.1f} ppm".format(
+            name, max(frac)*1e6, w=max([len(n) for n in noises])))
         # logging.warning("  max: {:e}, min: {:e}".format(max(frac), min(frac)))
 
         diffs[name] = (mn, pn, frac)
@@ -198,22 +196,26 @@ pygwinc: {fom:.2f} Mpc""".format(
             axr.loglog(freq, frac)
             axr.grid()
             axr.axhline(y=max(frac), color='r', linestyle='--')
-            axr.text(max(freq)+4000, max(frac), '{:.1f}%'.format(max(frac)*100),
+            axr.text(max(freq)+4000, max(frac), '{:.1f} ppm'.format(max(frac)*1e6),
                      horizontalalignment='left', verticalalignment='center',
                      color='red')
 
-        axl.set_xlabel("frequency [Hz]")
-        axr.set_xlabel("frequency [Hz]")
+        if diffs:
+            axl.set_xlabel("frequency [Hz]")
+            axr.set_xlabel("frequency [Hz]")
     
-        plt.suptitle("""{} mat/py gwinc noise comparison
-noises that differ by more than {}% [(mat-py)/py]
-{}""".format(args.IFO, FRACTIONAL_TOLERANCE*100, fom_title))
+            plt.suptitle("""{} mat/py gwinc noise comparison
+noises that differ by more than {} ppm [(mat-py)/py]
+{}""".format(args.IFO, fractional_tolerance*1e6, fom_title))
+
+            if args.save:
+                plt.gcf().set_size_inches(11, (len(diffs)+1)*4)
+                plt.savefig(args.save)
+            else:
+                plt.show()
 
-        if args.save:
-            plt.gcf().set_size_inches(11, (len(diffs)+1)*4)
-            plt.savefig(args.save)
         else:
-            plt.show()
+            logging.warning("All tests passed, so no plot was generated")
 
     ##############################