precomp.py 7.32 KB
Newer Older
1
from __future__ import division
2
from numpy import pi, sqrt, sin, exp
Christopher Wipf's avatar
Christopher Wipf committed
3 4
from scipy.io import loadmat
import scipy.special
5
import logging
6
import copy
Christopher Wipf's avatar
Christopher Wipf committed
7

8
from . import suspension
9
from .const import CONSTANTS
10
from .struct import Struct
11
from .noise.coatingthermal import getCoatDopt
Christopher Wipf's avatar
Christopher Wipf committed
12 13


Jameson Graef Rollins's avatar
Jameson Graef Rollins committed
14
def precompIFO(f, ifoin, PRfixed=True):
15
    """Add precomputed data to the IFO model.
16

Christopher Wipf's avatar
Christopher Wipf committed
17 18 19 20
    To prevent recomputation of these precomputed data, if the
    ifo argument contains ifo.gwinc.PRfixed, and this matches
    the argument PRfixed, no changes are made.

Jameson Graef Rollins's avatar
Jameson Graef Rollins committed
21 22 23
    This function DOES NOT modify IFO. Its return value is a copy of
    IFO with precomputations filled out.

24
    """
25
    ifo = copy.deepcopy(ifoin)
26

27 28
    if 'gwinc' not in ifo:
        ifo.gwinc = Struct()
29

Christopher Wipf's avatar
Christopher Wipf committed
30 31
    ifo.gwinc.PRfixed = PRfixed

32 33
    ##############################
    # derived temp
34

35
    if 'Temp' not in ifo.Materials.Substrate:
36
        ifo.Materials.Substrate.Temp = ifo.Infrastructure.Temp
37

38 39 40 41 42
    ##############################
    # suspension vertical-horizontal coupling

    if 'VHCoupling' not in ifo.Suspension:
        ifo.Suspension.VHCoupling = Struct()
43
        ifo.Suspension.VHCoupling.theta = ifo.Infrastructure.Length / CONSTANTS['R_earth']
44

45 46 47 48
    ##############################
    # optics values

    # calculate optics' parameters
Christopher Wipf's avatar
Christopher Wipf committed
49 50
    ifo.Materials.MirrorVolume = pi*ifo.Materials.MassRadius**2 * \
                                 ifo.Materials.MassThickness
51
    ifo.Materials.MirrorMass = ifo.Materials.MirrorVolume * \
52
                               ifo.Materials.Substrate.MassDensity
Christopher Wipf's avatar
Christopher Wipf committed
53 54 55
    ifo.Optics.ITM.Thickness = ifo.Materials.MassThickness

    # coating layer optical thicknesses - mevans 2 May 2008
56
    if 'CoatLayerOpticalThickness' not in ifo.Optics.ITM:
57 58 59 60 61 62 63 64
        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)
65

66 67 68 69 70 71 72 73 74 75 76 77 78 79
    ##############################
    # beam parameters

    armlen = ifo.Infrastructure.Length

    # g-factors
    g1 = 1 - armlen / ifo.Optics.Curvature.ITM
    g2 = 1 - armlen / ifo.Optics.Curvature.ETM
    gcav = sqrt(g1 * g2 * (1 - g1 * g2))
    gden = g1 - 2 * g1 * g2 + g2

    if (g1 * g2 * (1 - g1 * g2)) <= 0:
        raise Exception('Unstable arm cavity g-factors.  Change ifo.Optics.Curvature')
    elif gcav < 1e-3:
80
        logging.warning('Nearly unstable arm cavity g-factors.  Reconsider ifo.Optics.Curvature')
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

    ws = sqrt(armlen * ifo.Laser.Wavelength / pi)
    w1 = ws * sqrt(abs(g2) / gcav)
    w2 = ws * sqrt(abs(g1) / gcav)

    # waist size
    w0 = ws * sqrt(gcav / abs(gden))
    zr = pi * w0**2 / ifo.Laser.Wavelength
    z1 = armlen * g2 * (1 - g1) / gden
    z2 = armlen * g1 * (1 - g2) / gden

    ifo.Optics.ITM.BeamRadius = w1
    ifo.Optics.ETM.BeamRadius = w2

    ##############################
    # calc power and IFO parameters

Christopher Wipf's avatar
Christopher Wipf committed
98
    pbs, parm, finesse, prfactor, Tpr = precompPower(ifo, PRfixed)
99

Christopher Wipf's avatar
Christopher Wipf committed
100 101 102 103
    ifo.gwinc.pbs = pbs
    ifo.gwinc.parm = parm
    ifo.gwinc.finesse = finesse
    ifo.gwinc.prfactor = prfactor
104 105 106 107 108 109
    ifo.gwinc.gITM = g1
    ifo.gwinc.gETM = g2
    ifo.gwinc.BeamWaist = w0
    ifo.gwinc.BeamRayleighRange = zr
    ifo.gwinc.BeamWaistToITM = z1
    ifo.gwinc.BeamWaistToETM = z2
Christopher Wipf's avatar
Christopher Wipf committed
110 111
    ifo.Optics.PRM.Transmittance = Tpr

112 113 114
    ##############################
    # saved seismic spectrum

115
    if 'darmSeiSusFile' in ifo.Seismic and ifo.Seismic.darmSeiSusFile:
116 117 118
        darmsei = loadmat(ifo.Seismic.darmSeiSusFile)
        ifo.Seismic.darmseis_f = darmsei['darmseis_f'][0]
        ifo.Seismic.darmseis_x = darmsei['darmseis_x'][0]
Christopher Wipf's avatar
Christopher Wipf committed
119

120 121 122 123
    # --------------------------------------------------------
    # Suspension
    # if the suspension code supports different temps for the stages
    fname = eval('suspension.susp{}'.format(ifo.Suspension.Type))
Jameson Graef Rollins's avatar
Jameson Graef Rollins committed
124
    hForce, vForce, hTable, vTable = fname(f, ifo)
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144

    try:
        # full TF (conventional)
        ifo.Suspension.hForce = hForce.fullylossy
        ifo.Suspension.vForce = vForce.fullylossy
        # TFs with each stage lossy
        ifo.Suspension.hForce_singlylossy = hForce.singlylossy
        ifo.Suspension.vForce_singlylossy = vForce.singlylossy
    except:
        ifo.Suspension.hForce = hForce
        ifo.Suspension.vForce = vForce

    ifo.Suspension.hTable = hTable
    ifo.Suspension.vTable = vTable

    ##############################
    # strain to length conversion
    #
    # used for converting displacement to strain in all noise calculations

Jameson Graef Rollins's avatar
Jameson Graef Rollins committed
145
    ifo.gwinc.dhdl_sqr, ifo.gwinc.sinc_sqr = dhdl(f, ifo.Infrastructure.Length)
146

Christopher Wipf's avatar
Christopher Wipf committed
147 148 149
    return ifo


150
def precompPower(ifo, PRfixed=True):
151
    """Compute power on beamsplitter, finesse, and power recycling factor.
Christopher Wipf's avatar
Christopher Wipf committed
152

153
    """
Christopher Wipf's avatar
Christopher Wipf committed
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
    c       = scipy.constants.c
    pin     = ifo.Laser.Power
    t1      = sqrt(ifo.Optics.ITM.Transmittance)
    r1      = sqrt(1 - ifo.Optics.ITM.Transmittance)
    r2      = sqrt(1 - ifo.Optics.ETM.Transmittance)
    t5      = sqrt(ifo.Optics.PRM.Transmittance)
    r5      = sqrt(1 - ifo.Optics.PRM.Transmittance)
    loss    = ifo.Optics.Loss                          # single TM loss
    bsloss  = ifo.Optics.BSLoss
    acoat   = ifo.Optics.ITM.CoatingAbsorption
    pcrit   = ifo.Optics.pcrit

    # Finesse, effective number of bounces in cavity, power recycling factor
    finesse = 2*pi / (t1**2 + 2*loss)        # arm cavity finesse
    neff    = 2 * finesse / pi

    # Arm cavity reflectivity with finite loss
    garm = t1 / (1 - r1*r2*sqrt(1-2*loss))  # amplitude gain wrt input field
    rarm = r1 - t1 * r2 * sqrt(1-2*loss) * garm

174
    if PRfixed:
Christopher Wipf's avatar
Christopher Wipf committed
175 176 177 178 179 180 181 182 183 184 185 186 187 188
        Tpr = ifo.Optics.PRM.Transmittance  # use given value
    else:
        Tpr = 1-(rarm*sqrt(1-bsloss))**2 # optimal recycling mirror transmission
        t5 = sqrt(Tpr)
        r5 = sqrt(1 - Tpr)
    prfactor = t5**2 / (1 + r5 * rarm * sqrt(1-bsloss))**2

    pbs  = pin * prfactor          # BS power from input power
    parm = pbs * garm**2 / 2       # arm power from BS power

    asub = 1.3*2*ifo.Optics.ITM.Thickness*ifo.Optics.SubstrateAbsorption
    pbsl = 2*pcrit/(asub+acoat*neff) # bs power limited by thermal lensing

    if pbs > pbsl:
189
        logging.warning('P_BS exceeds BS Thermal limit!')
Christopher Wipf's avatar
Christopher Wipf committed
190 191 192 193

    return pbs, parm, finesse, prfactor, Tpr


194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
def dhdl(f, armlen):
    """Strain to length conversion for noise power spetra

    This takes into account the GW wavelength and is only important
    when this is comparable to the detector arm length.

    From R. Schilling, CQG 14 (1997) 1513-1519, equation 5,
    with n = 1, nu = 0.05, ignoring overall phase and cos(nu)^2.
    A small value of nu is used instead of zero to avoid infinities.

    Returns the square of the dh/dL function, and the same divided by
    the arm length squared.

    """
    c = scipy.constants.c
209
    nu_small = 15*pi/180
210 211 212 213 214 215 216 217 218
    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
    sinc_sqr = 4 / abs(sin(omega_arm_f) * exp(-1j * omega_arm) / omega_arm_f
                       + sin(omega_arm_b) * exp(1j * omega_arm) / omega_arm_b)**2
    # keep DC value equal to 1
    sinc_sqr /= sinc_sqr[0]
    dhdl_sqr = sinc_sqr / armlen**2
    return dhdl_sqr, sinc_sqr