Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org on the morning of Tuesday 11th August 2020, starting at approximately 9am PDT. It is expected to take around 20 minutes and there will be a short period of downtime (less than five minutes) towards the end of the maintenance window. Please direct any comments, questions, or concerns to computing-help@ligo.org.

Commit e17950bd authored by Sean Leavey's avatar Sean Leavey

Initial commit

parent 4acc3a49
# dielectric-thermal-noise # Ditherm: Dielectric Thermal Noise Calculator
Dielectric mirror thermal noise This is a tool to calculate Brownian thermal noise in dielectric mirror stacks. In theory, it supports noise calculations for n-layer stacks of any material, but it has only been tested with two-material stacks using data from the [Advanced LIGO](https://www.advancedligo.mit.edu/) noise calculator, [GWINC](https://awiki.ligo-wa.caltech.edu/aLIGO/GWINC) (that link requires albert.einstein style credentials, available to members of the LIGO Scientific Community only).
## Notes ##
I provide this software without any kind of guarantee that it provides the correct answer. This software is in no way affiliated with or endorsed by any university or organisation - it is purely the work of a private individual.
In particular, the extension from the two-material Brownian noise calculation (mostly based on Harry et al., 2002) to an n-material calculation was performed simply by inspecting the formulae. I don't know if this is reasonable or not. I suspect that for coatings with more than a few (~3) materials and with extremely small or large (outwith the range wavelength * 0.15 to wavelength * 0.35) optical thicknesses per layer, the results will diverge from reality. Use at your own risk!
## Testing ##
To test the software against precomputed values from GWINC, run:
`python ditherm test`
from the root Ditherm directory (the same directory as this readme).
## Future Work ##
- Implement phase correction to Brownian calculation as per the functions on [this page](https://awiki.ligo-wa.caltech.edu/aLIGO/GWINC)
- Look at proper way to combine multiple materials by reading these papers:
- http://journals.aps.org/prd/abstract/10.1103/PhysRevD.91.042002
- http://journals.aps.org/prd/abstract/10.1103/PhysRevD.87.082001
- http://journals.aps.org/prd/abstract/10.1103/PhysRevD.91.042001
## Credits ##
This implementation was written by Sean Leavey, but the calculations are based on a number of sources:
- http://dx.doi.org/10.1088/0264-9381/19/5/305
- http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=851432 (for the correction to perpendicular component of the loss angle from Harry et al. 2002)
- http://dx.doi.org/10.1088/0264-9381/24/2/008
- GWINC (for the simplified version of the parallel component of Poisson's ratio)
- Personal discussion with colleagues
Sean Leavey
https://github.com/SeanDS/
\ No newline at end of file
import sys
if len(sys.argv) > 1:
if sys.argv[1] == 'test':
import unittest
loader = unittest.TestLoader()
suite = loader.discover('.')
runner = unittest.runner.TextTestRunner()
runner.run(suite)
\ No newline at end of file
from __future__ import division
import materials
class Layer(object):
def __init__(self, material, d):
self.material = material
self.d = d
@property
def material(self):
return self.__material
@material.setter
def material(self, material):
self.__material = material
@property
def d(self):
"""
Layer physical thickness.
"""
return self.__d
@d.setter
def d(self, d):
self.__d = d
class Bilayer(object):
def __init__(self, layerA, layerB):
self.layerA = layerA
self.layerB = layerB
@property
def layerA(self):
return self.__layerA
@layerA.setter
def layerA(self, layerA):
self.__layerA = layerA
@property
def layerB(self):
return self.__layerB
@layerB.setter
def layerB(self, layerB):
self.__layerB = layerB
def d(self):
"""
Total bilayer physical thickness.
"""
return self.layerA.d + self.layerB.d
def yPara(self):
"""
Parallel component of Young's modulus.
"""
return 1 / self.d() * (self.layerA.d * self.layerA.material.Y + self.layerB.d * self.layerB.material.Y)
def yPerp(self):
"""
Perpendicular component of Young's modulus.
"""
return self.d() / (self.layerA.d / self.layerA.material.Y + self.layerB.d / self.layerB.material.Y)
def phiPara(self):
"""
Parallel component of loss angle.
"""
return 1 / (self.d() * self.yPara()) * (self.layerA.material.Y * self.layerA.material.phi * self.layerA.d + self.layerB.material.Y * self.layerB.material.phi * self.layerB.d)
def phiPerp(self):
"""
Perpendicular component of loss angle.
"""
return self.yPerp() / self.d() * (self.layerA.d * self.layerA.material.phi / self.layerA.material.Y + self.layerB.d * self.layerB.material.phi / self.layerB.material.Y)
def sigmaPara(self):
"""
Parallel component of Poisson ratio.
"""
return (self.layerA.d * self.layerA.material.sigma * self.layerA.material.Y * (1 - self.layerA.material.sigma ** 2) + self.layerB.d * self.layerB.material.sigma * self.layerB.material.Y * (1 - self.layerB.material.sigma ** 2)) / (self.layerA.d * self.layerA.material.Y * (1 - self.layerB.material.sigma ** 2) + self.layerB.d * self.layerB.material.Y * (1 - self.layerA.material.sigma ** 2))
def sigmaPerp(self):
"""
Perpendicular component of Poisson ratio.
"""
return (self.layerA.material.sigma * self.layerA.material.Y * self.layerA.d + self.layerB.material.sigma * self.layerB.material.Y * self.layerB.d) / (self.layerA.material.Y * self.layerA.d + self.layerB.material.Y * self.layerB.d)
\ No newline at end of file
from __future__ import division
class Material(object):
def __init__(self, name, Y, sigma, phi, n):
self.name = name
self.Y = Y
self.sigma = sigma
self.phi = phi
self.n = n
@property
def name(self):
return self.__name
@name.setter
def name(self, name):
self.__name = name
@property
def Y(self):
return self.__Y
@Y.setter
def Y(self, Y):
self.__Y = Y
@property
def sigma(self):
return self.__sigma
@sigma.setter
def sigma(self, sigma):
self.__sigma = sigma
@property
def phi(self):
return self.__phi
@phi.setter
def phi(self, phi):
self.__phi = phi
@property
def n(self):
return self.__n
@n.setter
def n(self, n):
self.__n = n
class SilicaCoating(Material):
def __init__(self, *args, **kwargs):
super(SilicaCoating, self).__init__("Silica", 72e9, 0.17, 4e-5, 1.45)
class TantalaCoating(Material):
def __init__(self, *args, **kwargs):
super(TantalaCoating, self).__init__("Tantala", 140e9, 0.23, 3.8e-4, 2.06)
class TitaniumTantalaCoating(Material):
def __init__(self, *args, **kwargs):
super(TitaniumTantalaCoating, self).__init__("Titanium doped Tantala", 140e9, 0.23, 3e-4, 2.06)
class SilicaSubstrate(Material):
def __init__(self, *args, **kwargs):
super(SilicaSubstrate, self).__init__("Silica", 72e9, 0.17, 5e-9, 1.45)
\ No newline at end of file
from __future__ import division
import numpy as np
import layers
import materials
class Stack(object):
def __init__(self, layers, substrate):
self.layers = layers
self.substrate = substrate
@property
def layers(self):
return self.__layers
@layers.setter
def layers(self, layers):
self.__layers = layers
@property
def substrate(self):
return self.__substrate
@substrate.setter
def substrate(self, substrate):
self.__substrate = substrate
def d(self):
"""
Total thickness.
"""
return sum(layer.d for layer in self.layers)
def yPara(self):
"""
Total parallel Young's modulus.
"""
return 1 / self.d() * sum([layer.d * layer.material.Y for layer in self.layers])
def yPerp(self):
"""
Total perpendicular Young's modulus.
"""
return self.d() / sum([layer.d / layer.material.Y for layer in self.layers])
def phiPara(self):
"""
Total parallel loss angle.
"""
return 1 / (self.d() * self.yPara()) * sum([layer.material.Y * layer.material.phi * layer.d for layer in self.layers])
def phiPerp(self):
"""
Total perpendicular loss angle.
"""
return self.yPerp() / self.d() * sum([layer.d * layer.material.phi / layer.material.Y for layer in self.layers])
def sigmaPara(self):
"""
Total stack parallel Poisson ratio.
"""
return np.mean([layer.material.sigma for layer in self.layers])
def sigmaPerp(self):
"""
Total perpendicular Poisson ratio.
"""
return sum([layer.material.sigma * layer.material.Y * layer.d for layer in self.layers]) / sum([layer.material.Y * layer.d for layer in self.layers])
def phi(self, beamSize):
"""
Effective loss angle.
"""
return (self.d() / (np.sqrt(np.pi) * beamSize * self.yPerp()) *
(self.phiPerp() *
(self.substrate.Y / (1 - self.substrate.sigma ** 2) -
2 * self.sigmaPerp() ** 2 * self.substrate.Y * self.yPara() /
(self.yPerp() * (1 - self.substrate.sigma ** 2) * (1 - self.sigmaPara()))) +
self.yPara() * self.sigmaPerp() * (1 - 2 * self.substrate.sigma) /
((1 - self.sigmaPara()) * (1 - self.substrate.sigma)) *
(self.phiPara() - self.phiPerp()) +
self.yPara() * self.yPerp() * (1 + self.substrate.sigma) *
(self.phiPara() * (1 - 2 * self.substrate.sigma) ** 2) /
(self.substrate.Y * (1 - self.sigmaPara() ** 2) * (1 - self.substrate.sigma))))
def brownianNoise(self, freq, beamSize, temperature):
k = 1.3806503e-23;
return 2 * k * temperature / (np.sqrt(np.pi ** 3) * freq * beamSize * self.substrate.Y) * (1 - self.substrate.sigma ** 2) * self.phi(beamSize)
# GWINC version, which does the same thing
#c = self.d()*(1-self.substrate.sigma**2)/(np.pi*beamSize**2)*((1/(self.yPerp()*(1-self.substrate.sigma**2))-2*self.sigmaPerp()**2*self.yPara()/(self.yPerp()**2*(1-self.substrate.sigma**2)*(1-self.sigmaPara())))*self.phiPerp() + self.yPara()*self.sigmaPerp()*(1-2*self.substrate.sigma)/(self.yPerp()*self.substrate.Y*(1-self.sigmaPara())*(1-self.substrate.sigma))*(self.phiPara()-self.phiPerp())+self.yPara()*(1+self.substrate.sigma)*(1-2*self.substrate.sigma)**2/(self.substrate.Y**2*(1-self.sigmaPara()**2)*(1-self.substrate.sigma))*self.phiPara())
#return 4 * k * temperature * c / (2 * np.pi * freq)
\ No newline at end of file
This diff is collapsed.
% Script to extract calculated Brownian noise values from GWINC.
% Run from within the GWINC directory.
% load precomputed IFO model
load('ifo.mat');
% frequencies to calculate
f = logspace(-2, 5, 1000);
% calculate Brownian noise for mirrors
SbrITM = getCoatBrownian(f, ifo, 'ITM');
SbrETM = getCoatBrownian(f, ifo, 'ETM');
% save CSV file
csvwrite('aligo.csv', [f', SbrITM', SbrETM']);
\ No newline at end of file
from __future__ import division
import sys
import os
sys.path.append('..')
from unittest import TestCase
import materials
import layers
import stacks
import numpy as np
class AdvancedLigoEtmStack(stacks.Stack):
def __init__(self, wavelength):
substrate = materials.Material("Silica Substrate", 7.27e10, 0.167, 5e-9, 1.45)
coatingA = materials.Material("Silica Coating", 7.2e10, 0.17, 4e-5, 1.45)
coatingB = materials.Material("Titanium Tantala Coating", 1.4e11, 0.23, 2.3e-4, 2.06539)
layerA = layers.Layer(coatingA, 0.27 * wavelength / coatingA.n)
layerB = layers.Layer(coatingB, 0.23 * wavelength / coatingB.n)
topLayer = layers.Layer(coatingA, 0.5 * wavelength / coatingA.n)
baseLayer = layers.Layer(coatingA, 0.163870186147445 * wavelength / coatingB.n)
theseLayers = np.array([topLayer, layerB] + [layerA, layerB] * 17 + [layerA, baseLayer])
super(AdvancedLigoEtmStack, self).__init__(theseLayers, substrate)
class AdvancedLigoItmStack(stacks.Stack):
def __init__(self, wavelength):
substrate = materials.Material("Silica Substrate", 7.27e10, 0.167, 5e-9, 1.45)
coatingA = materials.Material("Silica Coating", 7.2e10, 0.17, 4e-5, 1.45)
coatingB = materials.Material("Titanium Tantala Coating", 1.4e11, 0.23, 2.3e-4, 2.06539)
layerA = layers.Layer(coatingA, 0.308 * wavelength / coatingA.n)
layerB = layers.Layer(coatingB, 0.192 * wavelength / coatingB.n)
topLayer = layers.Layer(coatingA, 0.5 * wavelength / coatingA.n)
baseLayer = layers.Layer(coatingA, 0.186957128029190 * wavelength / coatingB.n)
theseLayers = np.array([topLayer, layerB] + [layerA, layerB] * 7 + [layerA, baseLayer])
super(AdvancedLigoItmStack, self).__init__(theseLayers, substrate)
class TestBrownianNoise(TestCase):
def setUp(self):
self.etmStack = AdvancedLigoEtmStack(1064e-9)
self.itmStack = AdvancedLigoItmStack(1064e-9)
def test_d(self):
self.assertAlmostEqual(self.etmStack.d(), 6.150299645767251e-6, delta=1e-20)
def test_y_para(self):
gwincVal = 9.651384066646355e10
self.assertAlmostEqual(self.etmStack.yPara(), gwincVal, delta=0.3*gwincVal)
def test_y_perp(self):
gwincVal = 8.728318664479849e10
self.assertAlmostEqual(self.etmStack.yPerp(), gwincVal, delta=0.3*gwincVal)
def test_phi_para(self):
gwincVal = 1.393560882693337e-4
self.assertAlmostEqual(self.etmStack.phiPara(), gwincVal, delta=0.05*gwincVal)
def test_phi_perp(self):
gwincVal = 8.270302150752522e-5
self.assertAlmostEqual(self.etmStack.phiPerp(), gwincVal, delta=0.05*gwincVal)
def test_sigma_para(self):
gwincVal = 0.2
self.assertAlmostEqual(self.etmStack.sigmaPara(), gwincVal, delta=0.01*gwincVal)
def test_sigma_perp(self):
gwincVal = 0.201375606821895
self.assertAlmostEqual(self.etmStack.sigmaPerp(), gwincVal, delta=0.01*gwincVal)
def test_brownian_noise_itm(self):
# load GWINC exported data
gwincData = np.genfromtxt(os.path.join(os.path.dirname(__file__), 'gwinc', 'aligo.csv'), delimiter=',')
wITM = 55e-3
T = 290
# calculate noise
dithermITM = self.itmStack.brownianNoise(gwincData[:, 0], wITM, T)
self.assertTrue(np.allclose(dithermITM, gwincData[:, 1], rtol=0.05, atol=0))
def test_brownian_noise_etm(self):
# load GWINC exported data
gwincData = np.genfromtxt(os.path.join(os.path.dirname(__file__), 'gwinc', 'aligo.csv'), delimiter=',')
wETM = 62e-3
T = 290
# calculate noise
dithermETM = self.etmStack.brownianNoise(gwincData[:, 0], wETM, T)
self.assertTrue(np.allclose(dithermETM, gwincData[:, 2], rtol=0.05, atol=0))
\ No newline at end of file
"""
Brownian noise in Advanced LIGO ETM
"""
from __future__ import division
import sys
sys.path.append('..')
import stacks
import matplotlib.pyplot as plt
import numpy as np
###
# Parameters
# wavelength
wavelength = 1064e-9 # [m]
# frequency range
f = np.logspace(1, 4, 1000) # [Hz]
# temperature
T = 290 # [K]
# beam sizes (distance from centre where power drops to 1/e)
wITM = 55e-3 # [m]
wETM = 62e-3 # [m]
###
# Calculation
stackITM = stacks.AdvancedLigoItmStack(wavelength)
stackETM = stacks.AdvancedLigoEtmStack(wavelength)
# calculate Brownian noise amplitude for frequency range
brownianNoiseITM = stackITM.brownianNoise(f, wITM, T)
brownianNoiseETM = stackETM.brownianNoise(f, wETM, T)
# combined noise for interferometer
brownianNoiseCombined = 2 * (brownianNoiseITM + brownianNoiseETM);
###
# Plot
# log-log plot
plt.loglog(f, np.sqrt(brownianNoiseITM), f, np.sqrt(brownianNoiseETM), f, np.sqrt(brownianNoiseCombined))
# axis labels, etc.
plt.legend(['ITM', 'ETM', 'Advanced LIGO Total'])
plt.xlabel('Frequency [Hz]')
plt.ylabel('Displacement-equivalent noise [m / sqrt(Hz)]')
plt.title('Brownian noise in Advanced LIGO ETM')
# add grid
plt.grid(True)
# display
plt.show()
\ No newline at end of file
from __future__ import division
import sys
sys.path.append('../..')
import ditherm.materials as materials
import ditherm.layers as layers
import ditherm.stacks as stacks
import numpy as np
class AdvancedLigoItmStack(stacks.Stack):
def __init__(self, wavelength):
substrate = materials.Material("Silica Substrate", 7.27e10, 0.167, 5e-9, 1.45)
coatingA = materials.Material("Silica Coating", 7.2e10, 0.17, 4e-5, 1.45)
coatingB = materials.Material("Titanium Tantala Coating", 1.4e11, 0.23, 2.3e-4, 2.06539)
layerA = layers.Layer(coatingA, 0.308 * wavelength / coatingA.n)
layerB = layers.Layer(coatingB, 0.192 * wavelength / coatingB.n)
topLayer = layers.Layer(coatingA, 0.5 * wavelength / coatingA.n)
baseLayer = layers.Layer(coatingA, 0.186957128029190 * wavelength / coatingB.n)
theseLayers = np.array([topLayer, layerB] + [layerA, layerB] * 7 + [layerA, baseLayer])
super(AdvancedLigoItmStack, self).__init__(theseLayers, substrate)
class AdvancedLigoEtmStack(stacks.Stack):
def __init__(self, wavelength):
substrate = materials.Material("Silica Substrate", 7.27e10, 0.167, 5e-9, 1.45)
coatingA = materials.Material("Silica Coating", 7.2e10, 0.17, 4e-5, 1.45)
coatingB = materials.Material("Titanium Tantala Coating", 1.4e11, 0.23, 2.3e-4, 2.06539)
layerA = layers.Layer(coatingA, 0.27 * wavelength / coatingA.n)
layerB = layers.Layer(coatingB, 0.23 * wavelength / coatingB.n)
topLayer = layers.Layer(coatingA, 0.5 * wavelength / coatingA.n)
baseLayer = layers.Layer(coatingA, 0.163870186147445 * wavelength / coatingB.n)
theseLayers = np.array([topLayer, layerB] + [layerA, layerB] * 17 + [layerA, baseLayer])
super(AdvancedLigoEtmStack, self).__init__(theseLayers, substrate)
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment