diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..ef8d369716351225ab87491064b8cf70fd5074b2
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,8 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# Backups
+*~
+
diff --git a/IFOModel.py b/IFOModel.py
new file mode 100644
index 0000000000000000000000000000000000000000..98b0d5e447277b5fcc3743b8adafee649e1c987c
--- /dev/null
+++ b/IFOModel.py
@@ -0,0 +1,364 @@
+from __future__ import division, print_function
+from numpy import pi, NaN
+from scipy.io.matlab.mio5_params import mat_struct
+
+
+def IFOModel():
+ """IFOMODEL returns a structure describing an IFO for use in
+ benchmark programs and noise simulator. Part of the gwinc
+ package, which provides science-grounded figures of merit for
+ comparing interferometric gravitational wave detector designs.
+
+
+
+ parameters for quad pendulum suspension updated 3rd May 2006, NAR
+ References:
+ LIGO-T000012-00-D
+ * Differentiate between silica and sapphire substrate absorption
+ * Change ribbon suspension aspect ratio
+ * Change pendulum frequency
+ * References:
+ * 1. Electro-Optic Handbook, Waynant & Ediger (McGraw-Hill: 1993)
+ * 2. LIGO/GEO data/experience
+ * 3. Suspension reference design, LIGO-T000012-00
+ * 4. Quartz Glass for Optics Data and Properties, Heraeus data sheet,
+ * numbers for suprasil
+ * 5. Y.S. Touloukian (ed), Thermophysical Properties of Matter
+ * (IFI/Plenum,1970)
+ * 6. Marvin J. Weber (ed) CRC Handbook of laser science and technology,
+ * Vol 4, Pt 2
+ * 7. R.S. Krishnan et al.,Thermal Expansion of Crystals, Pergamon Press
+ * 8. P. Klocek, Handbook of infrared and optical materials, Marcel Decker,
+ * 1991
+ * 9. Rai Weiss, electronic log from 5/10/2006
+ * 10. Wikipedia online encyclopedia, 2006
+ * 11. D.K. Davies, The Generation and Dissipation of Static Charge on
+ * dielectrics in a Vacuum, page 29
+ * 12. Gretarsson & Harry, Gretarsson thesis
+ * 13. Fejer
+ * 14. Braginsky"""
+
+ ifo = mat_struct()
+
+ ## Infrastructure----------------------------------------------------------
+
+ ifo.Infrastructure = mat_struct()
+ ifo.Infrastructure.ResidualGas = mat_struct()
+ ifo.Infrastructure.Length = 3995 # m;
+ ifo.Infrastructure.ResidualGas.pressure = 4.0e-7 # Pa;
+ ifo.Infrastructure.ResidualGas.mass = 3.35e-27 # kg; Mass of H_2 (ref. 10)
+ ifo.Infrastructure.ResidualGas.polarizability = 7.8e-31 # m^3; Gas polarizability
+
+ ## Physical and other constantMaterialss; All in SI units------------------
+
+ ifo.Constants = mat_struct()
+ #ifo.Constants.E0 = 8.8541878176e-12; % F / m; Permittivity of Free Space
+ #ifo.Constants.hbar = 1.054572e-34; % J - s; (Plancks constant) / (2 * pi)
+ #ifo.Constants.kB = 1.380658e-23; % J / K; Boltzman Constant
+ #ifo.Constants.h = ifo.Constants.hbar * 2 * pi; % J - s; Planks constant
+ #ifo.Constants.R = 8.31447215; % J / (K * mol); Gas Constant
+ #ifo.Constants.m_e = 9.10938291e-31; % kg; electron mass
+ #ifo.Constants.c = 2.99792458e8; % m / s; speed of light in vacuum
+ ifo.Constants.Temp = 290 # K; Temperature of the Vacuum
+ #ifo.Constants.yr = 365.2422 * 86400; % sec; Seconds in a year
+ #ifo.Constants.M_earth = 5.972e24; % mass of Earth [kg]
+ #ifo.Constants.R_earth = 6.3781e6; % radius of Earth [m]
+ #ifo.Constants.fs = 16384; % Sampling frequency (Hz)
+ #ifo.Constants.AU = 149597870700; % m; Astronomical unit, IAU 2012 Resolution B2
+ #ifo.Constants.parsec = ifo.Constants.AU * (648000 / pi); % m, IAU 2015 Resolution B2
+ #ifo.Constants.Mpc = ifo.Constants.parsec * 1e6; % m, IAU 2015 Resolution B2
+ #ifo.Constants.SolarMassParameter = 1.3271244e20; % m^3 / s^2; G * MSol, IAU 2015 Resolution B3
+ #ifo.Constants.G = 6.67408e-11; % m^3 / (kg s^2); Grav. const
+ # % http://arxiv.org/abs/1507.07956
+ #ifo.Constants.MSol = ifo.Constants.SolarMassParameter / ifo.Constants.G; % kg; Solar mass
+ #ifo.Constants.g = 9.806; % m / s^2; grav. acceleration
+ # % http://physics.nist.gov/cuu/Constants/
+ #ifo.Constants.H0 = 67110; % ms^( - 1); Hubble const.
+ # % http://arxiv.org/pdf/1303.5076v3.pdf
+ #ifo.Constants.omegaM = 0.3175; % Mass density parameter
+ # % http://arxiv.org/pdf/1303.5076v3.pdf
+ #ifo.Constants.omegaLambda = 1 - ifo.Constants.omegaM; % Cosmological constant density parameter
+ # % omegaK = 0 (flat universe) is assumed
+
+ ## Parameter describing thermal lensing --------------------------------------
+ # The presumably dominant effect of a thermal lens in the ITMs is an increased
+ # mode mismatch into the SRC, and thus an increased effective loss of the SRC.
+ # This increase is estimated by calculating the round-trip loss S in the SRC as
+ # 1-S = |<Psi|exp(i*phi)|Psi>|^2, where
+ # |Psi> is the beam hitting the ITM and
+ # phi = P_coat*phi_coat + P_subs*phi_subs
+ # with phi_coat & phi__subs the specific lensing profiles
+ # and P_coat & P_subst the power absorbed in coating and substrate
+ #
+ # This expression can be expanded to 2nd order and is given by
+ # S= s_cc P_coat^2 + 2*s_cs*P_coat*P_subst + s_ss*P_subst^2
+ # s_cc, s_cs and s_ss where calculated analytically by Phil Wilems (4/2007)
+ ifo.TCS = mat_struct()
+ ifo.TCS.s_cc=7.024 # Watt^-2
+ ifo.TCS.s_cs=7.321 # Watt^-2
+ ifo.TCS.s_ss=7.631 # Watt^-2
+
+ # The hardest part to model is how efficient the TCS system is in
+ # compensating this loss. Thus as a simple Ansatz we define the
+ # TCS efficiency TCSeff as the reduction in effective power that produces
+ # a phase distortion. E.g. TCSeff=0.99 means that the compensated distortion
+ # of 1 Watt absorbed is eqivalent to the uncompensated distortion of 10mWatt.
+ # The above formula thus becomes:
+ # S= s_cc P_coat^2 + 2*s_cs*P_coat*P_subst + s_ss*P_subst^2 * (1-TCSeff)^2
+ #
+ # To avoid iterative calculation we define TCS.SCRloss = S as an input
+ # and calculate TCSeff as an output.
+ # TCS.SRCloss is incorporated as an additional loss in the SRC
+ ifo.TCS.SRCloss = 0.00
+
+
+ ## Seismic and Gravity Gradient Parameters---------------------------------
+ ifo.Seismic = mat_struct()
+ ifo.Seismic.Site = 'LHO' # LHO or LLO (only used for Newtonian noise)
+ ifo.Seismic.darmSeiSusFile = 'seismic.mat' # .mat file containing predictions for darm displacement
+ ifo.Seismic.KneeFrequency = 10 # Hz; freq where 'flat' noise rolls off
+ ifo.Seismic.LowFrequencyLevel = 1e-9 # m/rtHz; seismic noise level below f_knee
+ ifo.Seismic.Gamma = .8 # abruptness of change at f_knee
+ ifo.Seismic.Rho = 1.8e3 # kg/m^3; density of the ground nearby
+ ifo.Seismic.Beta = 0.5 # quiet times beta = 0.35-0.60
+ # noisy times beta = 0.15-1.4
+ ifo.Seismic.Omicron = 1 # Feedforward cancellation factor
+
+ ifo.Seismic.darmSeiSusFile = 'CryogenicLIGO/Sensitivity/GWINC/seismic.mat'
+
+ ## Suspension: SI Units----------------------------------------------------
+ ifo.Suspension = mat_struct()
+ ifo.Suspension.BreakStress = 750e6 # Pa; ref. K. Strain
+ ifo.Suspension.Temp = 300
+ ifo.Suspension.VHCoupling = mat_struct()
+ ifo.Suspension.VHCoupling.theta = 1e-3 # vertical-horizontal x-coupling
+
+ # new Silicon parameters added for gwincDev RA 10-04-20 ~~~~~~~~~~~~~~~~~~~
+ # ref ---- http://design.caltech.edu/Research/MEMS/siliconprop.html
+ # all properties should be for T ~ 20 K
+ ifo.Suspension.Silicon = mat_struct()
+ ifo.Suspension.Silicon.Rho = 2330 # Kg/m^3; density
+ ifo.Suspension.Silicon.C = 772 # J/kg/K heat capacity
+ ifo.Suspension.Silicon.K = 4980 # W/m/K thermal conductivity
+ ifo.Suspension.Silicon.Alpha = 1e-9 # 1/K thermal expansion coeff
+
+ # from Gysin, et. al. PRB (2004) E(T) = E0 - B*T*exp(-T0/T)
+ # E0 = 167.5e9 Pa T0 = 317 K B = 15.8e6 Pa/K
+ ifo.Suspension.Silicon.dlnEdT = 2.5e-10 # (1/K) dlnE/dT T=20K
+
+ ifo.Suspension.Silicon.Phi = 2e-9 # Nawrodt (2010) loss angle 1/Q
+ ifo.Suspension.Silicon.Y = 150e9 # Pa Youngs Modulus
+ ifo.Suspension.Silicon.Dissdepth = 1.5e-3 # 10x smaller surface loss depth (Nawrodt (2010))
+ ifo.Suspension.FiberType = 0 # 0 = round, 1 = ribbons
+ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ ifo.Suspension.Silica = mat_struct()
+ ifo.Suspension.Silica.Rho = 2200 # Kg/m^3;
+ ifo.Suspension.Silica.C = 772 # J/Kg/K;
+ ifo.Suspension.Silica.K = 1.38 # W/m/kg;
+ ifo.Suspension.Silica.Alpha = 3.9e-7 # 1/K;
+ ifo.Suspension.Silica.dlnEdT = 1.52e-4 # (1/K), dlnE/dT
+ ifo.Suspension.Silica.Phi = 4.1e-10 # from G Harry e-mail to NAR 27April06 dimensionless units
+ ifo.Suspension.Silica.Y = 72e9 # Pa; Youngs Modulus
+ ifo.Suspension.Silica.Dissdepth = 1.5e-2 # from G Harry e-mail to NAR 27April06
+
+ ifo.Suspension.C70Steel = mat_struct()
+ ifo.Suspension.C70Steel.Rho = 7800
+ ifo.Suspension.C70Steel.C = 486
+ ifo.Suspension.C70Steel.K = 49
+ ifo.Suspension.C70Steel.Alpha = 12e-6
+ ifo.Suspension.C70Steel.dlnEdT = -2.5e-4
+ ifo.Suspension.C70Steel.Phi = 2e-4
+ ifo.Suspension.C70Steel.Y = 212e9 # measured by MB for one set of wires
+
+ ifo.Suspension.MaragingSteel = mat_struct()
+ ifo.Suspension.MaragingSteel.Rho = 7800
+ ifo.Suspension.MaragingSteel.C = 460
+ ifo.Suspension.MaragingSteel.K = 20
+ ifo.Suspension.MaragingSteel.Alpha = 11e-6
+ ifo.Suspension.MaragingSteel.dlnEdT = 0
+ ifo.Suspension.MaragingSteel.Phi = 1e-4
+ ifo.Suspension.MaragingSteel.Y = 187e9
+ # consistent with measured blade spring constants NAR
+
+ ifo.Suspension.Type = 'Quad' # 0 for cylindrical suspension
+
+ # Note stage numbering: mirror is at beginning of stack, not end
+ ifo.Suspension.Stage1 = mat_struct()
+ ifo.Suspension.Stage2 = mat_struct()
+ ifo.Suspension.Stage3 = mat_struct()
+ ifo.Suspension.Stage4 = mat_struct()
+
+ ifo.Suspension.Stage1.Mass = 39.6 # kg; current numbers May 2006 NAR
+ ifo.Suspension.Stage2.Mass = 39.6
+ ifo.Suspension.Stage3.Mass = 21.8
+ ifo.Suspension.Stage4.Mass = 22.1
+
+ ifo.Suspension.Stage1.Length = 0.602 # m; current numbers May 2006 NAR
+ ifo.Suspension.Stage2.Length = 0.341 # m;
+ ifo.Suspension.Stage3.Length = 0.277 # m;
+ ifo.Suspension.Stage4.Length = 0.416 # m;
+
+ ifo.Suspension.Stage1.Dilution = NaN
+ ifo.Suspension.Stage2.Dilution = 106 # updated May06 NAR
+ ifo.Suspension.Stage3.Dilution = 80
+ ifo.Suspension.Stage4.Dilution = 87
+
+ ifo.Suspension.Stage1.K = NaN #
+ ifo.Suspension.Stage2.K = 5200 # N/m; vertical spring constant
+ ifo.Suspension.Stage3.K = 3900 # N/m; vertical spring constant
+ ifo.Suspension.Stage4.K = 3400 # N/m; vertical spring constant
+
+ ifo.Suspension.Stage1.WireRadius = NaN
+ ifo.Suspension.Stage2.WireRadius = 310e-6 # current numbers May 2006 NAR
+ ifo.Suspension.Stage3.WireRadius = 350e-6
+ ifo.Suspension.Stage4.WireRadius = 520e-6
+
+ # For Ribbon suspension
+ ifo.Suspension.Ribbon = mat_struct()
+ ifo.Suspension.Fiber = mat_struct()
+ ifo.Suspension.Ribbon.Thickness = 115e-6 # m;
+ ifo.Suspension.Ribbon.Width = 1150e-6 # m;
+ ifo.Suspension.Fiber.Radius = 205e-6 # m;
+
+ ifo.Suspension.Stage1.Blade = NaN # blade thickness
+ ifo.Suspension.Stage2.Blade = 4200e-6 # current numbers May 2006 NAR
+ ifo.Suspension.Stage3.Blade = 4600e-6
+ ifo.Suspension.Stage4.Blade = 4300e-6
+
+ ifo.Suspension.Stage1.NWires = 4
+ ifo.Suspension.Stage2.NWires = 4
+ ifo.Suspension.Stage3.NWires = 4
+ ifo.Suspension.Stage4.NWires = 2
+
+ ## Dielectric coating material parameters----------------------------------
+
+ ## high index material: tantala
+ ifo.Materials = mat_struct()
+ ifo.Materials.Coating = mat_struct()
+ ifo.Materials.Coating.Yhighn = 140e9
+ ifo.Materials.Coating.Sigmahighn = 0.23
+ ifo.Materials.Coating.CVhighn = 2.1e6 # Crooks et al, Fejer et al
+ ifo.Materials.Coating.Alphahighn = 3.6e-6 # 3.6e-6 Fejer et al, 5e-6 from Braginsky
+ ifo.Materials.Coating.Betahighn = 1.4e-5 # dn/dT, value Gretarrson (G070161)
+ ifo.Materials.Coating.ThermalDiffusivityhighn = 33 # Fejer et al
+ ifo.Materials.Coating.Phihighn = 2.3e-4
+ ifo.Materials.Coating.Indexhighn = 2.06539
+
+ ## low index material: silica
+ ifo.Materials.Coating.Ylown = 72e9
+ ifo.Materials.Coating.Sigmalown = 0.17
+ ifo.Materials.Coating.CVlown = 1.6412e6 # Crooks et al, Fejer et al
+ ifo.Materials.Coating.Alphalown = 5.1e-7 # Fejer et al
+ ifo.Materials.Coating.Betalown = 8e-6 # dn/dT, (ref. 14)
+ ifo.Materials.Coating.ThermalDiffusivitylown = 1.38 # Fejer et al
+ ifo.Materials.Coating.Philown = 4.0e-5
+ ifo.Materials.Coating.Indexlown = 1.45
+
+ ## Substrate Material parameters--------------------------------------------
+
+ ifo.Materials.Substrate = mat_struct()
+ ifo.Materials.Substrate.c2 = 7.6e-12 # Coeff of freq depend. term for bulk mechanical loss, 7.15e-12 for Sup2
+ ifo.Materials.Substrate.MechanicalLossExponent=0.77 # Exponent for freq dependence of silica loss, 0.822 for Sup2
+ ifo.Materials.Substrate.Alphas = 5.2e-12 # Surface loss limit (ref. 12)
+ ifo.Materials.Substrate.MirrorY = 7.27e10 # N/m^2; Youngs modulus (ref. 4)
+ ifo.Materials.Substrate.MirrorSigma = 0.167 # Kg/m^3; Poisson ratio (ref. 4)
+ ifo.Materials.Substrate.MassDensity = 2.2e3 # Kg/m^3; (ref. 4)
+ ifo.Materials.Substrate.MassAlpha = 3.9e-7 # 1/K; thermal expansion coeff. (ref. 4)
+ ifo.Materials.Substrate.MassCM = 739 # J/Kg/K; specific heat (ref. 4)
+ ifo.Materials.Substrate.MassKappa = 1.38 # J/m/s/K; thermal conductivity (ref. 4)
+ ifo.Materials.Substrate.RefractiveIndex = 1.45 # mevans 25 Apr 2008
+
+ ifo.Materials.MassRadius = 0.17 # m;
+ ifo.Materials.MassThickness = 0.200 # m; Peter F 8/11/2005
+
+ ## Laser-------------------------------------------------------------------
+ ifo.Laser = mat_struct()
+ ifo.Laser.Wavelength = 1.064e-6 # m;
+ ifo.Laser.Power = 125 # W;
+
+ ## Optics------------------------------------------------------------------
+ ifo.Optics = mat_struct()
+ ifo.Optics.Type = 'SignalRecycled'
+
+ ifo.Optics.ITM = mat_struct()
+ ifo.Optics.ETM = mat_struct()
+ ifo.Optics.PRM = mat_struct()
+ ifo.Optics.SRM = mat_struct()
+ ifo.Optics.SRM.CavityLength = 55 # m; ITM to SRM distance
+ ifo.Optics.PhotoDetectorEfficiency = 0.95 # photo-detector quantum efficiency
+ ifo.Optics.Loss = 37.5e-6 # average per mirror power loss
+ ifo.Optics.BSLoss = 0.5e-3 # power loss near beamsplitter
+ ifo.Optics.coupling = 1.0 # mismatch btwn arms & SRC modes; used to
+ # calculate an effective r_srm
+ ifo.Optics.Curvature = mat_struct()
+ ifo.Optics.Curvature.ITM = 1970 # ROC of ITM
+ ifo.Optics.Curvature.ETM = 2192 # ROC of ETM
+ #ifo.Optics.SubstrateAbsorption = 0.5e-4 # 1/m; bulk absorption coef (ref. 2)
+ ifo.Optics.SubstrateAbsorption = 0.3e-4 # 1/m; 0.3 ppm/cm for Hereaus
+ ifo.Optics.pcrit = 10 # W; tolerable heating power (factor 1 ATC)
+
+ # factor of 2.5 added to simulate LNG modes - remove after new LNG code is added
+ ifo.Optics.ITM.BeamRadius = 0.055 # m; 1/e^2 power radius
+ ifo.Optics.ETM.BeamRadius = 0.062 # m; 1/e^2 power radius
+
+ ifo.Optics.ITM.CoatingAbsorption = 0.5e-6 # absorption of ITM
+ ifo.Optics.ITM.SubstrateAbsorption = 0.3e-4 # 1/m; 0.3 ppm/cm for Hereaus
+ #ifo.Optics.ITM.SubstrateAbsorption = 15e-4 # 1/m; 15 ppm/cm for CZ Si
+ ifo.Optics.ITM.Transmittance = 0.014 # Transmittance of ITM
+ ifo.Optics.ETM.Transmittance = 5e-6 # Transmittance of ETM
+ ifo.Optics.SRM.Transmittance = 0.20 # Transmittance of SRM
+ ifo.Optics.PRM.Transmittance = 0.03
+
+ # coating layer optical thicknesses - mevans June 2008
+ ifo.Optics.ITM.CoatingThicknessLown = 0.308
+ ifo.Optics.ITM.CoatingThicknessCap = 0.5
+
+ ifo.Optics.ETM.CoatingThicknessLown = 0.27
+ ifo.Optics.ETM.CoatingThicknessCap = 0.5
+
+ #ifo.Optics.SRM.Tunephase = 0.23 # SRM tuning, 795 Hz narrowband
+ ifo.Optics.SRM.Tunephase = 0.0 # SRM tuning
+ ifo.Optics.Quadrature = mat_struct()
+ ifo.Optics.Quadrature.dc = pi/2 # demod/detection/homodyne phase
+
+ ## Squeezer Parameters------------------------------------------------------
+
+ # Define the squeezing you want:
+ # None = ignore the squeezer settings
+ # Freq Independent = nothing special (no filter cavties)
+ # Freq Dependent = applies the specified filter cavites
+ # Optimal = find the best squeeze angle, assuming no output filtering
+ # OptimalOptimal = optimal squeeze angle, assuming optimal readout phase
+ ifo.Squeezer = mat_struct()
+ ifo.Squeezer.Type = 'None'
+ ifo.Squeezer.AmplitudedB = 10 # SQZ amplitude [dB]
+ ifo.Squeezer.InjectionLoss = 0.05 # power loss to sqz
+ ifo.Squeezer.SQZAngle = 0 # SQZ phase [radians]
+
+ # Parameters for frequency dependent squeezing
+ ifo.Squeezer.FilterCavity = mat_struct()
+ ifo.Squeezer.FilterCavity.fdetune = -14.5 # detuning [Hz]
+ ifo.Squeezer.FilterCavity.L = 100 # cavity length
+ ifo.Squeezer.FilterCavity.Ti = 0.12e-3 # input mirror trasmission [Power]
+ ifo.Squeezer.FilterCavity.Te = 0 # end mirror trasmission
+ ifo.Squeezer.FilterCavity.Lrt = 100e-6 # round-trip loss in the cavity
+ ifo.Squeezer.FilterCavity.Rot = 0 # phase rotation after cavity
+
+ ## Variational Output Parameters--------------------------------------------
+ # Define the output filter cavity chain
+ # None = ignore the output filter settings
+ # Chain = apply filter cavity chain
+ # Optimal = find the best readout phase
+ ifo.OutputFilter = mat_struct()
+ ifo.OutputFilter.Type = 'None'
+
+ ifo.OutputFilter.FilterCavity = mat_struct()
+ ifo.OutputFilter.FilterCavity.fdetune = -30 # detuning [Hz]
+ ifo.OutputFilter.FilterCavity.L = 4000 # cavity length
+ ifo.OutputFilter.FilterCavity.Ti = 10e-3 # input mirror trasmission [Power]
+ ifo.OutputFilter.FilterCavity.Te = 0 # end mirror trasmission
+ ifo.OutputFilter.FilterCavity.Lrt = 100e-6 # round-trip loss in the cavity
+ ifo.OutputFilter.FilterCavity.Rot = 0 # phase rotation after cavity
+
+ return ifo
diff --git a/IFOModel_123_2000.py b/IFOModel_123_2000.py
index 05393f22b6d5220d6eb78bada2815276cf13777e..df8db04bdc50a213ea03c5aa43f0720070fc7dbf 100644
--- a/IFOModel_123_2000.py
+++ b/IFOModel_123_2000.py
@@ -1,4 +1,4 @@
-from __future__ import division
+from __future__ import division, print_function
from numpy import pi, NaN
from util import SpotSizes
import scipy.constants
diff --git a/__init__.py b/__init__.py
index ac15f06e2e1fa7dd362062cba85b6c67ff96c393..48b19c23c5877bd4574c4d7b6eed40daa2597076 100644
--- a/__init__.py
+++ b/__init__.py
@@ -1,3 +1,3 @@
from gwinc import gwinc
+from IFOModel import IFOModel
from IFOModel_123_2000 import IFOModel_123_2000
-
diff --git a/gwinc.py b/gwinc.py
index 1abdb963c852d1bb56e60a5687b3da0a5c544583..fa52f794c870c6e1eab9c97f642695147777f4f6 100644
--- a/gwinc.py
+++ b/gwinc.py
@@ -82,10 +82,10 @@ def gwinc(flo, fhi, ifoin, source=None, fig=False):
# quad cases, for backward compatability
if ifo.Suspension.Type in (0, 1):
- hForce, vForce, hTable, vTable = suspQuad(f, ifo)
+ hForce, vForce, hTable, vTable = noise.suspensionthermal.suspQuad(f, ifo)
elif ifo.Suspension.Type == 'Quad':
ifo.Suspension.Type = 0
- hForce, vForce, hTable, vTable = suspQuad(f, ifo)
+ hForce, vForce, hTable, vTable = noise.suspensionthermal.suspQuad(f, ifo)
else:
fname = noise.suspensionthermal.__dict__['susp' + ifo.Suspension.Type]
hForce, vForce, hTable, vTable = fname(f, ifo)
diff --git a/noise/coatingthermal.py b/noise/coatingthermal.py
index d41aa54985f280dd0ac5f09a778cc4f593b1c4d6..06055d7aeb8c76f9fc0858b28be9a40a50ad1e52 100644
--- a/noise/coatingthermal.py
+++ b/noise/coatingthermal.py
@@ -2,7 +2,7 @@ from __future__ import division, print_function
import scipy.constants
import scipy.special
import numpy as np
-from numpy import pi, sum, zeros, exp, imag, sqrt, sin, cos, sinh, cosh
+from numpy import pi, sum, zeros, exp, real, imag, sqrt, sin, cos, sinh, cosh, polyfit, roots, max, min, ceil, log
def coatbrownian(f, ifo):
@@ -599,3 +599,173 @@ def getCoatFiniteCorr(ifo, wBeam, dOpt):
# eq 92
Cfsm = sqrt((r0**2 * P) / (2 * R**2 * (1 + sigS)**2 * LAMBDA**2))
return Cfsm
+
+
+def getCoatDopt(ifo, T, dL, dCap=0.5):
+ """get coating layer optical thicknesses to match desired transmission
+
+ ifo = gwinc IFO model, or vector of 3 refractive indexes [nS, nL, nH]
+ nS, nL, nH = refractive index of substrate, low and high-n materials
+ nL is used for the even layers, nH for odd layers
+ this algorithm may fail if nS > nL
+ opticName = name of the Optic struct to use for T, dL and dCap
+ T = power transmission of coating
+ dL = optical thickness of low-n layers
+ high-n layers have dH = 0.5 - dL
+ dCap = first layer (low-n) thickness (default 0.5)
+
+ dOpt = optical thickness vector Nlayer x 1
+
+ Examples:
+ dOpt = getCoatDopt(ifo, 'ITM');
+ dOpt = getCoatDopt(ifo, 0.014, 0.3);
+ dOpt = getCoatDopt(ifo, 0.014, 0.3, 0.5);
+ dOpt = getCoatDopt([1.45, 1.45, 2.06], 0.014, 0.3, 0.5);"""
+
+ # helper functions
+
+ ##############################################
+ def getTrans(ifo, Ndblt, dL, dH, dCap, dTweak):
+
+ # the optical thickness vector
+ dOpt = zeros(2 * Ndblt)
+ dOpt[0] = dCap
+ dOpt[1::2] = dH
+ dOpt[2::2] = dL
+
+ N = dTweak.size
+ T = zeros(N)
+ for n in range(N):
+ dOpt[-1] = dTweak[n]
+ r = getCoatRefl(ifo, dOpt)
+ T[n] = 1 - abs(r**2)
+
+ return T
+
+ ##############################################
+ def getTweak(ifo, T, Ndblt, dL, dH, dCap, dScan, Nfit):
+
+ # tweak bottom layer
+ Tn = getTrans(ifo, Ndblt, dL, dH, dCap, dScan)
+ pf = polyfit(dScan, Tn - T, Nfit)
+ rts = 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)]))
+
+ # compute T for this dTweak
+ Td = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dTweak]))
+
+ return dTweak, Td
+
+ # plot for debugging
+ # plot(dScan, [Tn - T, polyval(pf, dScan)], dTweak, Td - T, 'ro')
+ # grid on
+ # legend('T exact', 'T fit', 'T at dTweak')
+ # title(sprintf('%d doublets', Ndblt))
+ # pause(1)
+
+ # get IFO model stuff (or create it for other functions)
+ pS = ifo.Materials.Substrate
+ pC = ifo.Materials.Coating
+
+ nS = pS.RefractiveIndex
+ nL = pC.Indexlown
+ nH = pC.Indexhighn
+
+ ########################
+ # find number of quarter-wave layers required, as first guess
+ nR = nH / nL
+ a1 = (2 - T + 2 * sqrt(1 - T)) / (nR * nH * T)
+ Ndblt = int(ceil(log(a1) / (2 * log(nR))))
+
+ # search through number of doublets to find Ndblt
+ # which gives T lower than required
+ dH = 0.5 - dL
+ Tn = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dH]))
+ while Tn < T and Ndblt > 1:
+ # strange, but T is too low... remove doublets
+ Ndblt = Ndblt - 1
+ Tn = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dH]))
+ while Tn > T and Ndblt < 1e3:
+ # add doublets until T > tN
+ Ndblt = Ndblt + 1
+ Tn = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dH]))
+
+ ########################
+ # tweak bottom layer
+ delta = 0.01
+ dScan = np.arange(0, 0.25+delta, delta)
+ dTweak = getTweak(ifo, T, Ndblt, dL, dH, dCap, dScan, 5)[0]
+
+ if not dTweak:
+ if nS > nL:
+ raise Exception('Coating tweak layer not sufficient since nS > nL.')
+ else:
+ raise Exception('Coating tweak layer not found... very strange.')
+
+ # now that we are close, get a better result with a linear fit
+ delta = 0.001
+ dScan = np.linspace(dTweak - 3*delta, dTweak + 3*delta, 7)
+ dTweak, Td = getTweak(ifo, T, Ndblt, dL, dH, dCap, dScan, 3)
+
+ # negative values are bad
+ if dTweak < 0.01:
+ dTweak = 0.01
+
+ # check the result
+ if abs(log(Td / T)) > 1e-3:
+ print('Exact coating tweak layer not found... %g%% error.' % abs(log(Td / T)))
+
+ ########################
+ # return dOpt vector
+ dOpt = zeros((2 * Ndblt, 1))
+ dOpt[0] = dCap
+ dOpt[1::2] = dH
+ dOpt[2::2] = dL
+ dOpt[-1] = dTweak
+
+ return dOpt
+
+
+def getCoatRefl(ifo, dOpt):
+ """returns amplitude reflectivity, with phase, of a coating
+
+ ifo = parameter struct from IFOmodel.m
+ dOpt = optical thickness / lambda of each layer
+ = geometrical thickness * refractive index / lambda
+
+ This code is taken from ewa/multidiel1.m"""
+
+ pS = ifo.Materials.Substrate
+ pC = ifo.Materials.Coating
+
+ nS = pS.RefractiveIndex
+ nL = pC.Indexlown
+ nH = pC.Indexhighn
+
+ Nlayer = dOpt.size
+
+ # refractive index of input, coating, and output materials
+ nAll = zeros(Nlayer + 2)
+ nAll[0] = 1 # vacuum input
+ nAll[1::2] = nL
+ nAll[2::2] = nH
+ nAll[-1] = nS # substrate output
+
+ # reflectivity of each interface
+ rInterface = (nAll[:-1] - nAll[1:]) / (nAll[:-1] + nAll[1:])
+
+ # combine layers as small FP cavities, starting with last reflectivity
+ rCoat = rInterface[-1]
+ for n in reversed(range(dOpt.size)):
+ # one-way phase in this layer
+ phi = 2 * pi * dOpt[n]
+
+ # accumulate reflectivity
+ rCoat = rCoat * exp(-2j * phi)
+ rCoat = (rInterface[n] + rCoat) / (1 + rInterface[n] * rCoat)
+
+ return rCoat
diff --git a/noise/suspensionthermal.py b/noise/suspensionthermal.py
index 3b0676a282dde151b42cb6c42577c2f0974d9991..1dfdafdd64aa1b21ef60382749a54c62645b2a8e 100644
--- a/noise/suspensionthermal.py
+++ b/noise/suspensionthermal.py
@@ -18,7 +18,7 @@ def suspR(f, ifo):
theta = ifo.Suspension.VHCoupling.theta
noise = zeros((1, f.size))
- if len(Temp) == 1: # if the temperature is uniform along the suspension
+ if np.isscalar(Temp) or len(Temp) == 1: # if the temperature is uniform along the suspension
##########################################################
# Suspension TFs
##########################################################
@@ -112,6 +112,43 @@ def suspBQuad(f, ifo):
modification for the different temperatures between the stages
K.Arai Mar 1., 2012"""
+ # helper functions
+ def construct_eom_matrix(k, m, f):
+ """construct a matrix for eq of motion
+ k is the array for the spring constants
+ f is the freq vector"""
+
+ w = 2*pi * f
+
+ A = zeros((4,4,f.size), dtype=complex)
+
+ A[0,1,:] = -k[1,:]; A[1,0,:] = A[1,2,:]
+ A[1,2,:] = -k[2,:]; A[2,1,:] = A[1,2,:]
+ A[2,3,:] = -k[3,:]; A[3,2,:] = A[2,3,:]
+ A[0,0,:] = k[0,:] + k[1,:] - m[0] * w**2
+ A[1,1,:] = k[1,:] + k[2,:] - m[1] * w**2
+ A[2,2,:] = k[2,:] + k[3,:] - m[2] * w**2
+ A[3,3,:] = k[3,:] - m[3] * w**2
+ return A
+
+ def calc_transfer_functions(A, B, k, f):
+
+ X = zeros([B.size,A[0,0,:].size], dtype=complex)
+
+ for j in range(A[0,0,:].size):
+ X[:,j] = np.linalg.solve(A[:,:,j], B)
+
+ # transfer function from the force on the TM to TM motion
+ hForce = zeros(f.shape, dtype=complex)
+ hForce[:] = X[3,:]
+
+ # transfer function from the table motion to TM motion
+ hTable = zeros(f.shape, dtype=complex)
+ hTable[:] = X[0,:]
+ hTable = hTable * k[0,:]
+
+ return hForce, hTable
+
# default arguments
fiberType = 0
@@ -120,7 +157,7 @@ def suspBQuad(f, ifo):
kB = scipy.constants.k
Temp = ifo.Suspension.Temp
- if len(Temp) == 1:
+ if np.isscalar(Temp) or len(Temp) == 1:
Temp = [Temp, Temp, Temp, Temp]
# if only one temp is given, use it for all stages
@@ -364,39 +401,319 @@ def suspBQuad(f, ifo):
return hForce, vForce, hTable, vTable #, Ah, Av
-def construct_eom_matrix(k, m, f):
- """construct a matrix for eq of motion
- k is the array for the spring constants
- f is the freq vector"""
+def suspQuad(f, ifo):
+ """Suspension for quadruple pendulum
+
+ Silica ribbons used in mirror suspension stage; steel in others
+ Violin modes included
+ Adapted from code by Morag Casey (Matlab) and Geppo Cagnoli (Maple)
+
+ f = frequency vector
+ ifo = IFO model
+ fiberType = suspension sub type 0 => round fibers, otherwise ribbons
+
+ hForce, vForce = transfer functions from the force on the TM to TM motion
+ these should have the correct losses for the mechanical system such
+ that the thermal noise is
+ dxdF = force on TM along beam line to position of TM along beam line
+ = hForce + theta^2 * vForce
+ = admittance / (i * w)
+ where theta = ifo.Suspension.VHCoupling.theta.
+ Since this is just suspension thermal noise, the TM internal
+ modes and coating properties should not be included.
+
+ hTable, vTable = TFs from support motion to TM motion
+
+ Ah = horizontal equations of motion
+ Av = vertical equations of motion
+
+ modification for the different temperatures between the stages
+ K.Arai Mar 1., 2012"""
- w = 2*pi * f
+ # helper functions
+ def construct_eom_matrix(k, m, f):
+ """construct a matrix for eq of motion
+ k is the array for the spring constants
+ f is the freq vector"""
+
+ w = 2*pi * f
+
+ A = zeros((4,4,f.size), dtype=complex)
+
+ A[0,1,:] = -k[1,:]; A[1,0,:] = A[1,2,:]
+ A[1,2,:] = -k[2,:]; A[2,1,:] = A[1,2,:]
+ A[2,3,:] = -k[3,:]; A[3,2,:] = A[2,3,:]
+ A[0,0,:] = k[0,:] + k[1,:] - m[0] * w**2
+ A[1,1,:] = k[1,:] + k[2,:] - m[1] * w**2
+ A[2,2,:] = k[2,:] + k[3,:] - m[2] * w**2
+ A[3,3,:] = k[3,:] - m[3] * w**2
+ return A
+
+ def calc_transfer_functions(A, B, k, f):
- A = zeros((4,4,f.size), dtype=complex)
+ X = zeros([B.size,A[0,0,:].size], dtype=complex)
- A[0,1,:] = -k[1,:]; A[1,0,:] = A[1,2,:]
- A[1,2,:] = -k[2,:]; A[2,1,:] = A[1,2,:]
- A[2,3,:] = -k[3,:]; A[3,2,:] = A[2,3,:]
- A[0,0,:] = k[0,:] + k[1,:] - m[0] * w**2
- A[1,1,:] = k[1,:] + k[2,:] - m[1] * w**2
- A[2,2,:] = k[2,:] + k[3,:] - m[2] * w**2
- A[3,3,:] = k[3,:] - m[3] * w**2
- return A
+ for j in range(A[0,0,:].size):
+ X[:,j] = np.linalg.solve(A[:,:,j], B)
+ # transfer function from the force on the TM to TM motion
+ hForce = zeros(f.shape, dtype=complex)
+ hForce[:] = X[3,:]
-def calc_transfer_functions(A, B, k, f):
+ # transfer function from the table motion to TM motion
+ hTable = zeros(f.shape, dtype=complex)
+ hTable[:] = X[0,:]
+ hTable = hTable * k[0,:]
+
+ return hForce, hTable
+
+ # default arguments
+ fiberType = 0
+
+ # Assign Physical Constants
+ g = scipy.constants.g
+ kB = scipy.constants.k
- X = zeros([B.size,A[0,0,:].size], dtype=complex)
+ Temp = ifo.Suspension.Temp
+ Temp = ifo.Suspension.Temp
+ if np.isscalar(Temp) or len(Temp) == 1:
+ Temp = [Temp, Temp, Temp, Temp]
+ # if only one temp is given, use it for all stages
+
+ alpha_si = ifo.Suspension.Silica.Alpha # coeff. thermal expansion
+ beta_si = ifo.Suspension.Silica.dlnEdT # temp. dependence Youngs modulus
+ rho = ifo.Suspension.Silica.Rho # mass density
+ C = ifo.Suspension.Silica.C
+ K = ifo.Suspension.Silica.K # W/(m kg)
+ ds = ifo.Suspension.Silica.Dissdepth # surface loss dissipation depth
+
+ rho_st = ifo.Suspension.C70Steel.Rho
+ C_st = ifo.Suspension.C70Steel.C
+ K_st = ifo.Suspension.C70Steel.K
+ Y_st = ifo.Suspension.C70Steel.Y
+ alpha_st = ifo.Suspension.C70Steel.Alpha
+ beta_st = ifo.Suspension.C70Steel.dlnEdT
+ phi_steel = ifo.Suspension.C70Steel.Phi
+
+ rho_m = ifo.Suspension.MaragingSteel.Rho
+ C_m = ifo.Suspension.MaragingSteel.C
+ K_m = ifo.Suspension.MaragingSteel.K
+ Y_m = ifo.Suspension.MaragingSteel.Y
+ alpha_m = ifo.Suspension.MaragingSteel.Alpha
+ beta_m = ifo.Suspension.MaragingSteel.dlnEdT
+ phi_marag = ifo.Suspension.MaragingSteel.Phi
+
+
+ # Begin parameter assignment
+
+ # Note that I'm counting stages differently than Morag. Morag's
+ # counting is reflected in the variable names in this funcion; my
+ # counting is reflected in the index into Stage().
+ # Morag's count has stage "n" labeled as 1 and the mirror as stage 4.
+ # I'm counting the mirror as stage 1 and proceeding up. The reason
+ # for the change is my assumption that th eimplication of referring
+ # to stage "n" is that, once you get far enough away from the
+ # mirror, you might have additional stages but not change their
+ # characteristics. The simplest implementation of this would be to
+ # work through the stages sequenctially, starting from 1, until one
+ # reached the end, and then repeat the final stage as many times as
+ # desired. What I've done with the reordering is prepare for the
+ # day when we might do that.
+
+ theta = ifo.Suspension.VHCoupling.theta
+
+ m1 = ifo.Suspension.Stage4.Mass
+ m2 = ifo.Suspension.Stage3.Mass
+ m3 = ifo.Suspension.Stage2.Mass
+ m4 = ifo.Suspension.Stage1.Mass
+
+ M1 = m1 + m2 + m3 + m4 # mass supported by stage n
+ M2 = m2 + m3 + m4 # mass supported by stage ...
+ M3 = m3 + m4 # mass supported by stage ...
+
+ L1 = ifo.Suspension.Stage4.Length
+ L2 = ifo.Suspension.Stage3.Length
+ L3 = ifo.Suspension.Stage2.Length
+ L4 = ifo.Suspension.Stage1.Length
+
+ dil1 = ifo.Suspension.Stage4.Dilution
+ dil2 = ifo.Suspension.Stage3.Dilution
+ dil3 = ifo.Suspension.Stage2.Dilution
+
+ kv10 = ifo.Suspension.Stage4.K # N/m, vert. spring constant,
+ kv20 = ifo.Suspension.Stage3.K
+ kv30 = ifo.Suspension.Stage2.K
+
+
+ # Correction for the pendulum restoring force
+ # replaced m1->M1, m2->M2, m3->M3
+ # K. Arai Feb. 29, 2012
+ kh10 = M1*g/L1 # N/m, horiz. spring constant, stage n
+ kh20 = M2*g/L2 # N/m, horiz. spring constant, stage 1
+ kh30 = M3*g/L3 # N/m, horiz. spring constant, stage 2
+ kh40 = m4*g/L4 # N/m, horiz. spring constant, last stage
+
+ r_st1 = ifo.Suspension.Stage4.WireRadius
+ r_st2 = ifo.Suspension.Stage3.WireRadius
+ r_st3 = ifo.Suspension.Stage2.WireRadius
+
+ t_m1 = ifo.Suspension.Stage4.Blade
+ t_m2 = ifo.Suspension.Stage3.Blade
+ t_m3 = ifo.Suspension.Stage2.Blade
- for j in range(A[0,0,:].size):
- X[:,j] = np.linalg.solve(A[:,:,j], B)
+ N1 = ifo.Suspension.Stage4.NWires # number of wires in stage n
+ N2 = ifo.Suspension.Stage3.NWires # Number of wires in stage 1
+ N3 = ifo.Suspension.Stage2.NWires # Number of wires in stage 1
+ N4 = ifo.Suspension.Stage1.NWires # Number of wires in stage 1
- # transfer function from the force on the TM to TM motion
- hForce = zeros(f.shape, dtype=complex)
- hForce[:] = X[3,:]
+ Y_si = ifo.Suspension.Silica.Y
- # transfer function from the table motion to TM motion
- hTable = zeros(f.shape, dtype=complex)
- hTable[:] = X[0,:]
- hTable = hTable * k[0,:]
+ if ifo.Suspension.FiberType == 0:
+ r_fib = ifo.Suspension.Fiber.Radius
+ xsect = pi*r_fib**2 # cross-sectional area
+ II4 = r_fib**4*pi/4 # x-sectional moment of inertia
+ mu_v = 2/r_fib # mu/(V/S), vertical motion
+ mu_h = 4/r_fib # mu/(V/S), horizontal motion
+ tau_si = 7.372e-2*rho*C*(4*xsect/pi)/K # TE time constant
+ else:
+ W = ifo.Suspension.Ribbon.Width
+ t = ifo.Suspension.Ribbon.Thickness
+ xsect = W * t
+ II4 = (W * t**3)/12
+ mu_v = 2 * (W + t)/(W*t)
+ mu_h = (3 * N4 * W + t)/(N4*W + t)*2*(W+t)/(W*t)
+ tau_si = (rho*C*t**2)/(K*pi**2)
+
+ # loss factor, last stage suspension, vertical
+ phiv4 = ifo.Suspension.Silica.Phi*(1 + mu_v*ds)
+ Y_si_v = Y_si * (1 + 1j*phiv4) # Vertical Young's modulus, silica
+
+ T4 = m4*g / N4 # Tension in last stage
+
+ # TE time constant, steel wire 1-3
+ # WHAT IS THIS CONSTANT 7.37e-2?
+ tau_steel1 = 7.37e-2*(rho_st*C_st*(2*r_st1)**2)/K_st
+ tau_steel2 = 7.37e-2*(rho_st*C_st*(2*r_st2)**2)/K_st
+ tau_steel3 = 7.37e-2*(rho_st*C_st*(2*r_st3)**2)/K_st
+
+ # TE time constant, maraging blade 1
+ tau_marag1 = (rho_m*C_m*t_m1**2)/(K_m*pi**2)
+ tau_marag2 = (rho_m*C_m*t_m2**2)/(K_m*pi**2)
+ tau_marag3 = (rho_m*C_m*t_m3**2)/(K_m*pi**2)
+
+ # vertical delta, maraging
+ delta_v1 = Y_m*alpha_m**2*Temp[0]/(rho_m*C_m)
+ delta_v2 = delta_v1
+ delta_v3 = delta_v1
+
+ # horizontal delta, steel, stage n
+ delta_h1 = Y_st*(alpha_st-beta_st*g*M1/(N1*pi*r_st1**2*Y_st))**2
+ delta_h1 = delta_h1*Temp[0]/(rho_st*C_st)
+
+ delta_h2 = Y_st*(alpha_st-beta_st*g*M2/(N2*pi*r_st2**2*Y_st))**2
+ delta_h2 = delta_h2*Temp[1]/(rho_st*C_st)
+
+ delta_h3 = Y_st*(alpha_st-beta_st*g*M3/(N3*pi*r_st3**2*Y_st))**2
+ delta_h3 = delta_h3*Temp[2]/(rho_st*C_st)
+
+ # solutions to equations of motion
+ B = np.array([ 0, 0, 0, 1]).T
+
+ w = 2*pi*f
+
+ # thermoelastic correction factor, silica
+ delta_s = Y_si*(alpha_si-beta_si*T4/(xsect*Y_si))**2*Temp[3]/(rho*C)
+
+
+ # vertical loss factor, maraging
+ phiv1 = phi_marag+delta_v1*tau_marag1*w/(1+w**2*tau_marag1**2)
+ phiv2 = phi_marag+delta_v2*tau_marag2*w/(1+w**2*tau_marag2**2)
+ phiv3 = phi_marag+delta_v3*tau_marag3*w/(1+w**2*tau_marag3**2)
+
+ # horizontal loss factor, steel, stage n
+ phih1 = phi_steel+delta_h1*tau_steel1*w/(1+w**2*tau_steel1**2)
+ phih2 = phi_steel+delta_h2*tau_steel2*w/(1+w**2*tau_steel2**2)
+ phih3 = phi_steel+delta_h3*tau_steel3*w/(1+w**2*tau_steel3**2)
+
+ kv1 = kv10*(1 + 1j*phiv1) # stage n spring constant, vertical
+ kv2 = kv20*(1 + 1j*phiv2) # stage 1 spring constant, vertical
+ kv3 = kv30*(1 + 1j*phiv3) # stage 2 spring constant, vertical
+
+ kh1 = kh10*(1 + 1j*phih1/dil1) # stage n spring constant, horizontal
+ kh2 = kh20*(1 + 1j*phih2/dil2) # stage 1 spring constant, horizontal
+ kh3 = kh30*(1 + 1j*phih3/dil3) # stage 2 spring constant, horizontal
+
+ # loss factor, last stage suspension, horizontal
+ phih4 = ifo.Suspension.Silica.Phi*(1 + mu_h * ds) + \
+ delta_s * (tau_si * w/(1 + tau_si**2*w**2))
+
+ # violin mode calculations
+ Y_si_h = Y_si * (1 + 1j*phih4) # Horizontal Young's modulus
+ simp1 = sqrt(rho/Y_si_h)*w # simplification factor 1 q
+ simp2 = sqrt(rho * xsect *w**2/T4) # simplification factor 2 p
+
+ # simplification factor 3 kk
+ simp3 = sqrt(T4 * (1 + II4 * xsect * Y_si_h * w**2 / T4**2) / (Y_si_h * II4))
+
+ a = simp3 * cos(simp2 * L4) # simplification factor a
+ b = sin(simp2 * L4) # simplification factor b
+
+ # vertical spring constant, last stage
+ kv40 = abs(N4 * Y_si_v * xsect / L4) # this seems to not be used ??
+ kv4 = N4 * Y_si_v * xsect * simp1 / (tan(simp1 * L4))
+
+ # numerator, horiz spring constant, last stage
+ kh4num = N4*II4*Y_si_h*simp2*simp3*(simp2**2+simp3**2)*(a+simp2*b)
+ # denominator, horiz spring constant, last stage
+ kh4den = (2 * simp2 * a + (simp2**2 - simp3**2) * b)
+ # horizontal spring constant, last stage
+ kh4 = -kh4num / kh4den
+
+ ###############################################################
+ # Equations of motion for the system
+ ###############################################################
+
+ m_list=np.hstack((m1, m2, m3, m4)) # array of the mass
+ kh_list=np.vstack((kh1, kh2, kh3, kh4)) # array of the horiz spring constants
+ kv_list=np.vstack((kv1, kv2, kv3, kv4)) # array of the vert spring constants
+
+ # Calculate TFs turning on the loss of each stage one by one
+ hForce = mat_struct()
+ vForce = mat_struct()
+ hForce.singlylossy = zeros((4, f.size), dtype=complex)
+ vForce.singlylossy = zeros((4, f.size), dtype=complex)
+ for ii in range(4): # specify the stage to turn on the loss
+ # horizontal
+ k_list = kh_list
+ # only the imaginary part of the specified stage is used.
+ k_list = real(k_list) + 1j*imag(np.vstack((k_list[0,:]*(ii==0), k_list[1,:]*(ii==1), k_list[2,:]*(ii==2), k_list[3,:]*(ii==3))))
+ # construct Eq of motion matrix
+ Ah = construct_eom_matrix(k_list, m_list, f)
+ # calculate TFs
+ hForce.singlylossy[ii,:], hTable = calc_transfer_functions(Ah, B, k_list, f)
+
+ # vertical
+ k_list = kv_list
+ # only the imaginary part of the specified stage is used.
+ k_list = real(k_list) + 1j*imag(np.vstack((k_list[0,:]*(ii==0), k_list[1,:]*(ii==1), k_list[2,:]*(ii==2), k_list[3,:]*(ii==3))))
+ # construct Eq of motion matrix
+ Av = construct_eom_matrix(k_list, m_list, f)
+ # calculate TFs
+ vForce.singlylossy[ii,:], vTable = calc_transfer_functions(Av, B, k_list, f)
+
+ # horizontal
+ k_list=kh_list # all of the losses are on
+ # construct Eq of motion matrix
+ Ah = construct_eom_matrix(k_list, m_list, f)
+ # calculate TFs
+ hForce.fullylossy, hTable = calc_transfer_functions(Ah, B, k_list, f)
+
+ # vertical
+ k_list=kv_list # all of the losses are on
+ # construct Eq of motion matrix
+ Av = construct_eom_matrix(k_list, m_list, f)
+ # calculate TFs
+ vForce.fullylossy, vTable = calc_transfer_functions(Av, B, k_list, f)
+
+ return hForce, vForce, hTable, vTable #, Ah, Av
- return hForce, hTable
diff --git a/util.py b/util.py
index 2cd1073558cd59fd3086ac36cba56c60fcff0d19..e6ff1b613ccbdfe5cae482f73201040ffc25190a 100644
--- a/util.py
+++ b/util.py
@@ -3,6 +3,7 @@ from numpy import pi, sqrt
from scipy.io.matlab.mio5_params import mat_struct
from scipy.io import loadmat
import scipy.special
+from noise.coatingthermal import getCoatDopt
def SpotSizes(g1, g2, L, lambda_):
"""calculates spot sizes using FP cavity parameters
@@ -50,8 +51,14 @@ def precompIFO(ifo, PRfixed):
# coating layer optical thicknesses - mevans 2 May 2008
if 'CoatLayerOpticalThickness' not in ifo.Optics.ITM.__dict__:
- ifo.Optics.ITM.CoatLayerOpticalThickness = getCoatDopt(ifo, 'ITM')
- ifo.Optics.ETM.CoatLayerOpticalThickness = getCoatDopt(ifo, 'ETM')
+ T = ifo.Optics.ITM.Transmittance
+ dL = ifo.Optics.ITM.CoatingThicknessLown
+ dCap = ifo.Optics.ITM.CoatingThicknessCap
+ ifo.Optics.ITM.CoatLayerOpticalThickness = getCoatDopt(ifo, T, dL, dCap=dCap)
+ T = ifo.Optics.ETM.Transmittance
+ dL = ifo.Optics.ETM.CoatingThicknessLown
+ dCap = ifo.Optics.ETM.CoatingThicknessCap
+ ifo.Optics.ETM.CoatLayerOpticalThickness = getCoatDopt(ifo, T, dL, dCap=dCap)
# compute power on BS
pbs, parm, finesse, prfactor, Tpr = precompPower(ifo, PRfixed)