# Module refract.py
import math

def refdry(nu, T, Pdry, Pvap):
    # From Miriad: Determine the complex refractivity of the dry components
    # of the atmosphere.
    #
    # Input:
    #  nu = observing frequency (Hz)
    #  T = temperature (K)
    #  Pdry = partial pressure of dry components (Pa)
    #  Pvap = partial pressure of water vapour (Pa)

    # Table of microwave oxygen lines and their parameters.
    nu0 = [ 49.452379, 49.962257, 50.474238, 50.987748, 51.503350, 52.021409,
            52.542393, 53.066906, 53.595748, 54.129999, 54.671157, 55.221365,
            55.783800, 56.264777, 56.363387, 56.968180, 57.612481, 58.323874,
            58.446589, 59.164204, 59.590982, 60.306057, 60.434775, 61.150558,
            61.800152, 62.411212, 62.486253, 62.997974, 63.568515, 64.127764,
            64.678900, 65.224067, 65.764769, 66.302088, 66.836827, 67.369595,
            67.900862, 68.431001, 68.960306, 69.489021, 70.017342, 18.750341,
            68.498350, 24.763120, 87.249370, 15.393150, 73.838730, 34.145330 ]
    a1 = [    0.12E-6,    0.34E-6,    0.94E-6,    2.46E-6,    6.08E-6,   14.14E-6,
              31.02E-6,   64.10E-6,  124.70E-6,  228.00E-6,  391.80E-6,  631.60E-6,
              953.50E-6,  548.90E-6, 1344.00E-6, 1763.00E-6, 2141.00E-6, 2386.00E-6,
              1457.00E-6, 2404.00E-6, 2112.00E-6, 2124.00E-6, 2461.00E-6, 2504.00E-6,
              2298.00E-6, 1933.00E-6, 1517.00E-6, 1503.00E-6, 1087.00E-6,  733.50E-6,
              463.50E-6,  274.80E-6,  153.00E-6,   80.09E-6,   39.46E-6,   18.32E-6,
              8.01E-6,    3.30E-6,    1.28E-6,    0.47E-6,    0.16E-6,  945.00E-6,
              67.90E-6,  638.00E-6,  235.00E-6,   99.60E-6,  671.00E-6,  180.00E-6 ]
    a2 = [  11.830, 10.720,  9.690,  8.690,  7.740,  6.840,
            6.000,  5.220,  4.480,  3.810,  3.190,  2.620,
            2.115,  0.010,  1.655,  1.255,  0.910,  0.621,
            0.079,  0.386,  0.207,  0.207,  0.386,  0.621,
            0.910,  1.255,  0.078,  1.660,  2.110,  2.620,
            3.190,  3.810,  4.480,  5.220,  6.000,  6.840,
            7.740,  8.690,  9.690, 10.720, 11.830,  0.000,
            0.020,  0.011,  0.011,  0.089,  0.079,  0.079 ]
    a3 = [   8.40E-3,  8.50E-3,  8.60E-3,  8.70E-3,  8.90E-3,  9.20E-3,
             9.40E-3,  9.70E-3, 10.00E-3, 10.20E-3, 10.50E-3, 10.79E-3,
             11.10E-3, 16.46E-3, 11.44E-3, 11.81E-3, 12.21E-3, 12.66E-3,
             14.49E-3, 13.19E-3, 13.60E-3, 13.82E-3, 12.97E-3, 12.48E-3,
             12.07E-3, 11.71E-3, 14.68E-3, 11.39E-3, 11.08E-3, 10.78E-3,
             10.50E-3, 10.20E-3, 10.00E-3,  9.70E-3,  9.40E-3,  9.20E-3,
             8.90E-3,  8.70E-3,  8.60E-3,  8.50E-3,  8.40E-3, 15.92E-3,
             19.20E-3, 19.16E-3, 19.20E-3, 18.10E-3, 18.10E-3, 18.10E-3 ]
    a4 = [  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.6, 0.6, 0.6, 0.6, 0.6, 0.6 ]
    a5 = [  5.60E-3,  5.60E-3,  5.60E-3,  5.50E-3,  5.60E-3,  5.50E-3,
            5.70E-3,  5.30E-3,  5.40E-3,  4.80E-3,  4.80E-3,  4.17E-3,
            3.75E-3,  7.74E-3,  2.97E-3,  2.12E-3,  0.94E-3, -0.55E-3,
            5.97E-3, -2.44E-3,  3.44E-3, -4.13E-3,  1.32E-3, -0.36E-3,
            -1.59E-3, -2.66E-3, -4.77E-3, -3.34E-3, -4.17E-3, -4.48E-3,
            -5.10E-3, -5.10E-3, -5.70E-3, -5.50E-3, -5.90E-3, -5.60E-3,
            -5.80E-3, -5.70E-3, -5.60E-3, -5.60E-3, -5.60E-3, -0.44E-3,
            0.00E00,  0.00E00,  0.00E00,  0.00E00,  0.00E00,  0.00E00 ]
    a6 = [  1.7,  1.7,  1.7,  1.7,  1.8,  1.8,
            1.8,  1.9,  1.8,  2.0,  1.9,  2.1,
            2.1,  0.9,  2.3,  2.5,  3.7, -3.1,
            0.8,  0.1,  0.5,  0.7, -1.0,  5.8,
            2.9,  2.3,  0.9,  2.2,  2.0,  2.0,
            1.8,  1.9,  1.8,  1.8,  1.7,  1.8,
            1.7,  1.7,  1.7,  1.7,  1.7,  0.9,
            1.0,  1.0,  1.0,  1.0,  1.0,  1.0 ]

    # Convert to the units of Liebe.
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-9

    ap = 1.4e-10 * (1.0 - 1.2e-5 * f ** 1.5)
    gamma0 = 5.6e-3 * (p + 1.1 * e) * theta ** 0.8
    nr = 2.588 * p * theta + 3.07e-4 * (1.0 / (1.0 + (f / gamma0) ** 2) - 1.0) * p * theta * theta
    ni = (2.0 * 3.07e-4 / (gamma0 * (1.0 + (f / gamma0) ** 2) * (1.0 + (f / 60.0) ** 2)) +
          ap * p * theta ** 2.5) * f * p * theta * theta

    # Sum the contributions of the lines.
    for i in xrange(0, len(nu0)):
        S = a1[i] * p * theta ** 3 * math.exp(a2[i] * (1.0 - theta))
        gamma = a3[i] * (p * theta ** (0.8 - a4[i])) + 1.1 * e * theta
        delta = a5[i] * p * theta ** a6[i]
        x = (nu0[i] - f) * (nu0[i] - f) + gamma * gamma
        y = (nu0[i] + f) * (nu0[i] + f) + gamma * gamma
        z = (nu0[i] + gamma * gamma / nu0[i])
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / nu0[i] + delta *
                       (1.0 / x - 1.0 / y) * gamma * f / nu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / nu0[i] - delta *
                       ((nu0[i] - f) / x + (nu0[i] + f) / y) * f / nu0[i])

    # Return the result.
    return complex(nr, ni)

def refvap(nu, T, Pdry, Pvap):
    # From Miriad; Determine the complex refractivity of the water vapour monomers.
    #
    # Inputs:
    #  nu = observating frequency (Hz)
    #  T = temperature (K)
    #  Pdry = partial pressure of dry components (Pa)
    #  Pvap = partial pressure of water vapour (Pa)

    # Table of the microwave water lines.
    mnu0 = [  22.235080,  67.813960, 119.995940, 183.310117, 321.225644, 325.152919,
              336.187000, 380.197372, 390.134508, 437.346667, 439.150812, 443.018295,
              448.001075, 470.888947, 474.689127, 488.491133, 503.568532, 504.482692,
              556.936002, 620.700807, 658.006500, 752.033227, 841.073593, 859.865000,
              899.407000, 902.555000, 906.205524, 916.171582, 970.315022, 987.926764 ]
    b1 = [   0.1090,  0.0011,  0.0007,   2.3000,  0.0464,  1.5400,
             0.0010, 11.9000,  0.0044,   0.0637,  0.9210,  0.1940,
             10.6000,  0.3300,  1.2800,   0.2530,  0.0374,  0.0125,
             510.0000,  5.0900,  0.2740, 250.0000,  0.0130,  0.1330,
             0.0550,  0.0380,  0.1830,   8.5600,  9.1600, 138.000 ]
    b2 = [ 2.143, 8.730, 8.347, 0.653, 6.156, 1.515,
           9.802, 1.018, 7.318, 5.015, 3.561, 5.015,
           1.370, 3.561, 2.342, 2.814, 6.693, 6.693,
           0.114, 2.150, 7.767, 0.336, 8.113, 7.989,
           7.845, 8.360, 5.039, 1.369, 1.842, 0.178 ]
    b3 = [ 27.84E-3, 27.60E-3, 27.00E-3, 28.35E-3, 21.40E-3, 27.00E-3,
           26.50E-3, 27.60E-3, 19.00E-3, 13.70E-3, 16.40E-3, 14.40E-3,
           23.80E-3, 18.20E-3, 19.80E-3, 24.90E-3, 11.50E-3, 11.90E-3,
           30.00E-3, 22.30E-3, 30.00E-3, 28.60E-3, 14.10E-3, 28.60E-3,
           28.60E-3, 26.40E-3, 23.40E-3, 25.30E-3, 24.00E-3, 28.60E-3 ]

    # Convert to the units of Liebe.
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-9

    nr = 2.39 * e * theta + 41.6 * e * theta * theta + 6.47e-6 * f ** 2.05 * e * theta ** 2.4
    ni = (0.915 * 1.40e-6 * p + 5.41e-5 * e * theta * theta * theta) * f * e * theta ** 2.5

    # Sum the contributions of the lines.
    for i in xrange(0, len(mnu0)):
        S = b1[i] * e * theta ** 3.5 * math.exp(b2[i] * (1.0 - theta))
        gamma = b3[i] * (p * theta ** 0.8 + 4.80 * e * theta)
        x = (mnu0[i] - f) * (mnu0[i] - f) + gamma * gamma
        y = (mnu0[i] + f) * (mnu0[i] + f) + gamma * gamma
        z = (mnu0[i] + gamma * gamma / mnu0[i])
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / mnu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / mnu0[i])

    # Return the result.
    return complex(nr, ni)
    
def pvapsat(T):
    # From Miriad; Determine the saturation pressure of water vapour.
    # Input:
    #  T = temperature (K)
    #
    # Output:
    #  vapour saturation pressure (Pa)

    if (T > 215):
        theta = 300.0 / T
        return 1e5 / (41.51 * (theta ** -5) * (10 **(9.384 * theta - 10.0)))
    else:
        return 0.0

def refract(t, pdry, pvap, z, n, nu, T0, el):
    # From Miriad; Compute refractive index for an atmosphere.
    # Determine the sky brightness and excess path lengths for a parallel
    # slab atmosphere. Liebe's model (1985) is used to determine the complex
    # refractive index of air.
    #
    # Input:
    #  n = the number of atmospheric layers.
    #  t = temperature of the layers. T[0] is the temperature at the lowest
    #      layer (K)
    #  Pdry = partial pressure of the dry components (Pa)
    #  Pvap = partial pressure of the water vapour components (Pa)
    #  z = height of the layer.
    #  nu = frequency of interest (Hz)
    #  T0 = astronomical brightness temperature (K)
    #  el = elevation angle of the source above the atmosphere (rad)
    #
    # Output:
    #  { 'Tb' = brightness temperature (K),
    #    'tau' = opacity (nepers)
    #    'Ldry' = excess path, dry component (m)
    #    'Lvap' = excess path, water vapour component (m) }

    # Some constants.
    HMKS = 6.6260755e-34 # Planck constant, J.s
    KMKS = 1.380658e-23 # Boltzmann constant, J/K
    CMKS = 299792458 # Speed of light, m/s

    tau = 0.0
    Tb = HMKS * nu / (KMKS * (math.exp(HMKS * nu / (KMKS * T0)) - 1))
    Ldry = 0.0
    Lvap = 0.0

    snell = math.sin(el)
    for i in xrange(n, 0, -1):
        if (i == 1):
            dz = 0.5 * (z[1] - z[0])
        elif (i == n):
            dz = 0.5 * (z[n] - z[n - 1])
        else:
            dz = 0.5 * (z[i + 1] - z[i - 1])
        Ndry = refdry(nu, t[i], pdry[i], pvap[i])
        Nvap = refvap(nu, t[i], pdry[i], pvap[i])
        nr = 1 + (Ndry.real + Nvap.real) * 1e-6
        ni = (Ndry.imag + Nvap.imag) * 1e-6
        l = dz * nr / math.sqrt(nr * nr + (snell * snell) - 1.0)
        dtau = l * 4.0 * math.pi * nu / CMKS * ni
        Tb = (Tb - t[i]) * math.exp(-dtau) + t[i]
        tau = tau + dtau
        Ldry = Ldry + l * Ndry.real * 1e-6
        Lvap = Lvap + l * Nvap.real * 1e-6

    return { 'Tb': Tb, 'tau': tau, 'Ldry': Ldry, 'Lvap': Lvap }

def calcOpacity(freq, el, t0, p0, h0):
    # From Miriad; Compute sky brightness and opacity of a model atmosphere.
    # Returns the transmissivity of the atmosphere given frequency, elevation
    # angle and meteorological data. This uses a simple model of the atmosphere
    # and Liebe's model (1985) of the complex refractive index of air.
    # Input:
    #  freq = frequency (Hz)
    #  el = elevation angle (radians)
    #  t0,p0,h0 = Met data; observatory temperature, pressure and humidity
    #             (K, Pa, fraction)
    # Output:
    #  { 'fac' = transmissivity (fraction between 0 and 1)
    #    'Tb' = sky brightness temperature (K) }

    # Atmospheric parameters.
    M = 28.96e-3
    R = 8.314
    Mv = 18e-3
    rho0 = 1e3
    g = 9.81

    # d = temperature lapse rate (K/m)
    d = 0.0065
    # z0 = water vapour scale height (m)
    z0 = 1540.0
    # zmax = max altitude of model atmosphere (m)
    zmax = 10000.0

    # Generate a model of the atmosphere - T is temperature, Pdry is
    # partial pressure of "dry" constituents, Pvap is the partial
    # pressure of the water vapour.
    N = 50 # Number of iterations through the atmosphere in height
    z = []
    T = []
    Pvap = []
    Pdry = []
    fac = []
    Tb = []
    tau = []
    ofreq = []
    for i in xrange(0, (N + 1)):
        zd = float(i) * zmax / float(N)
        z.append(zd)
        T.append(t0 / (1 + d / t0 * zd))
        P = p0 * math.exp(-1.0 * M * g / (R * t0) * (zd + 0.5 * d * zd * zd / t0))
        Pvap.append(min(h0 * math.exp(-1.0 * zd / z0) * pvapsat(T[0]),
                        pvapsat(T[i])))
        Pdry.append(P - Pvap[i])

    # Determine the transmissivity and sky brightness.
    #job_server = pp.Server(secret="hello")
    #jobs = [(inpf, job_server.submit(refract, (T, Pdry, Pvap, z, N, inpf, 2.7, el),
    #                                 (refdry, refvap, pvapsat), ("math",))) for inpf in freq]
    for i in xrange(0, len(freq)):
        #print "calculating atmosphere for frequency %.1f" % freq[i]
        resref = refract(T, Pdry, Pvap, z, N, freq[i], 2.7, el)
        fac.append(math.exp(-1.0 * resref['tau']))
        Tb.append(resref['Tb'])
        tau.append(resref['tau'])
        ofreq.append(freq[i])
    #for inpf, job in jobs:
    #    resref = job()
    #    fac.append(math.exp(-1.0 * resref['tau']))
    #    Tb.append(resref['Tb'])
    #    tau.append(resref['tau'])
    #    ofreq.append(inpf)

    return { 'fac': fac, 'Tb': Tb, 'tau': tau, 'freq': ofreq }

67dd6583829909bc34af509c3a8c3484dbbf8ca7
<mask token>


def refdry(nu, T, Pdry, Pvap):
    nu0 = [49.452379, 49.962257, 50.474238, 50.987748, 51.50335, 52.021409,
        52.542393, 53.066906, 53.595748, 54.129999, 54.671157, 55.221365, 
        55.7838, 56.264777, 56.363387, 56.96818, 57.612481, 58.323874, 
        58.446589, 59.164204, 59.590982, 60.306057, 60.434775, 61.150558, 
        61.800152, 62.411212, 62.486253, 62.997974, 63.568515, 64.127764, 
        64.6789, 65.224067, 65.764769, 66.302088, 66.836827, 67.369595, 
        67.900862, 68.431001, 68.960306, 69.489021, 70.017342, 18.750341, 
        68.49835, 24.76312, 87.24937, 15.39315, 73.83873, 34.14533]
    a1 = [1.2e-07, 3.4e-07, 9.4e-07, 2.46e-06, 6.08e-06, 1.414e-05, 
        3.102e-05, 6.41e-05, 0.0001247, 0.000228, 0.0003918, 0.0006316, 
        0.0009535, 0.0005489, 0.001344, 0.001763, 0.002141, 0.002386, 
        0.001457, 0.002404, 0.002112, 0.002124, 0.002461, 0.002504, 
        0.002298, 0.001933, 0.001517, 0.001503, 0.001087, 0.0007335, 
        0.0004635, 0.0002748, 0.000153, 8.009e-05, 3.946e-05, 1.832e-05, 
        8.01e-06, 3.3e-06, 1.28e-06, 4.7e-07, 1.6e-07, 0.000945, 6.79e-05, 
        0.000638, 0.000235, 9.96e-05, 0.000671, 0.00018]
    a2 = [11.83, 10.72, 9.69, 8.69, 7.74, 6.84, 6.0, 5.22, 4.48, 3.81, 3.19,
        2.62, 2.115, 0.01, 1.655, 1.255, 0.91, 0.621, 0.079, 0.386, 0.207, 
        0.207, 0.386, 0.621, 0.91, 1.255, 0.078, 1.66, 2.11, 2.62, 3.19, 
        3.81, 4.48, 5.22, 6.0, 6.84, 7.74, 8.69, 9.69, 10.72, 11.83, 0.0, 
        0.02, 0.011, 0.011, 0.089, 0.079, 0.079]
    a3 = [0.0084, 0.0085, 0.0086, 0.0087, 0.0089, 0.0092, 0.0094, 0.0097, 
        0.01, 0.0102, 0.0105, 0.01079, 0.0111, 0.01646, 0.01144, 0.01181, 
        0.01221, 0.01266, 0.01449, 0.01319, 0.0136, 0.01382, 0.01297, 
        0.01248, 0.01207, 0.01171, 0.01468, 0.01139, 0.01108, 0.01078, 
        0.0105, 0.0102, 0.01, 0.0097, 0.0094, 0.0092, 0.0089, 0.0087, 
        0.0086, 0.0085, 0.0084, 0.01592, 0.0192, 0.01916, 0.0192, 0.0181, 
        0.0181, 0.0181]
    a4 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6]
    a5 = [0.0056, 0.0056, 0.0056, 0.0055, 0.0056, 0.0055, 0.0057, 0.0053, 
        0.0054, 0.0048, 0.0048, 0.00417, 0.00375, 0.00774, 0.00297, 0.00212,
        0.00094, -0.00055, 0.00597, -0.00244, 0.00344, -0.00413, 0.00132, -
        0.00036, -0.00159, -0.00266, -0.00477, -0.00334, -0.00417, -0.00448,
        -0.0051, -0.0051, -0.0057, -0.0055, -0.0059, -0.0056, -0.0058, -
        0.0057, -0.0056, -0.0056, -0.0056, -0.00044, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0]
    a6 = [1.7, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.9, 1.8, 2.0, 1.9, 2.1, 2.1, 
        0.9, 2.3, 2.5, 3.7, -3.1, 0.8, 0.1, 0.5, 0.7, -1.0, 5.8, 2.9, 2.3, 
        0.9, 2.2, 2.0, 2.0, 1.8, 1.9, 1.8, 1.8, 1.7, 1.8, 1.7, 1.7, 1.7, 
        1.7, 1.7, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-09
    ap = 1.4e-10 * (1.0 - 1.2e-05 * f ** 1.5)
    gamma0 = 0.0056 * (p + 1.1 * e) * theta ** 0.8
    nr = 2.588 * p * theta + 0.000307 * (1.0 / (1.0 + (f / gamma0) ** 2) - 1.0
        ) * p * theta * theta
    ni = (2.0 * 0.000307 / (gamma0 * (1.0 + (f / gamma0) ** 2) * (1.0 + (f /
        60.0) ** 2)) + ap * p * theta ** 2.5) * f * p * theta * theta
    for i in xrange(0, len(nu0)):
        S = a1[i] * p * theta ** 3 * math.exp(a2[i] * (1.0 - theta))
        gamma = a3[i] * (p * theta ** (0.8 - a4[i])) + 1.1 * e * theta
        delta = a5[i] * p * theta ** a6[i]
        x = (nu0[i] - f) * (nu0[i] - f) + gamma * gamma
        y = (nu0[i] + f) * (nu0[i] + f) + gamma * gamma
        z = nu0[i] + gamma * gamma / nu0[i]
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / nu0[i] + delta * (
            1.0 / x - 1.0 / y) * gamma * f / nu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / nu0[i] - delta * (
            (nu0[i] - f) / x + (nu0[i] + f) / y) * f / nu0[i])
    return complex(nr, ni)


<mask token>

<mask token>


def refdry(nu, T, Pdry, Pvap):
    nu0 = [49.452379, 49.962257, 50.474238, 50.987748, 51.50335, 52.021409,
        52.542393, 53.066906, 53.595748, 54.129999, 54.671157, 55.221365, 
        55.7838, 56.264777, 56.363387, 56.96818, 57.612481, 58.323874, 
        58.446589, 59.164204, 59.590982, 60.306057, 60.434775, 61.150558, 
        61.800152, 62.411212, 62.486253, 62.997974, 63.568515, 64.127764, 
        64.6789, 65.224067, 65.764769, 66.302088, 66.836827, 67.369595, 
        67.900862, 68.431001, 68.960306, 69.489021, 70.017342, 18.750341, 
        68.49835, 24.76312, 87.24937, 15.39315, 73.83873, 34.14533]
    a1 = [1.2e-07, 3.4e-07, 9.4e-07, 2.46e-06, 6.08e-06, 1.414e-05, 
        3.102e-05, 6.41e-05, 0.0001247, 0.000228, 0.0003918, 0.0006316, 
        0.0009535, 0.0005489, 0.001344, 0.001763, 0.002141, 0.002386, 
        0.001457, 0.002404, 0.002112, 0.002124, 0.002461, 0.002504, 
        0.002298, 0.001933, 0.001517, 0.001503, 0.001087, 0.0007335, 
        0.0004635, 0.0002748, 0.000153, 8.009e-05, 3.946e-05, 1.832e-05, 
        8.01e-06, 3.3e-06, 1.28e-06, 4.7e-07, 1.6e-07, 0.000945, 6.79e-05, 
        0.000638, 0.000235, 9.96e-05, 0.000671, 0.00018]
    a2 = [11.83, 10.72, 9.69, 8.69, 7.74, 6.84, 6.0, 5.22, 4.48, 3.81, 3.19,
        2.62, 2.115, 0.01, 1.655, 1.255, 0.91, 0.621, 0.079, 0.386, 0.207, 
        0.207, 0.386, 0.621, 0.91, 1.255, 0.078, 1.66, 2.11, 2.62, 3.19, 
        3.81, 4.48, 5.22, 6.0, 6.84, 7.74, 8.69, 9.69, 10.72, 11.83, 0.0, 
        0.02, 0.011, 0.011, 0.089, 0.079, 0.079]
    a3 = [0.0084, 0.0085, 0.0086, 0.0087, 0.0089, 0.0092, 0.0094, 0.0097, 
        0.01, 0.0102, 0.0105, 0.01079, 0.0111, 0.01646, 0.01144, 0.01181, 
        0.01221, 0.01266, 0.01449, 0.01319, 0.0136, 0.01382, 0.01297, 
        0.01248, 0.01207, 0.01171, 0.01468, 0.01139, 0.01108, 0.01078, 
        0.0105, 0.0102, 0.01, 0.0097, 0.0094, 0.0092, 0.0089, 0.0087, 
        0.0086, 0.0085, 0.0084, 0.01592, 0.0192, 0.01916, 0.0192, 0.0181, 
        0.0181, 0.0181]
    a4 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6]
    a5 = [0.0056, 0.0056, 0.0056, 0.0055, 0.0056, 0.0055, 0.0057, 0.0053, 
        0.0054, 0.0048, 0.0048, 0.00417, 0.00375, 0.00774, 0.00297, 0.00212,
        0.00094, -0.00055, 0.00597, -0.00244, 0.00344, -0.00413, 0.00132, -
        0.00036, -0.00159, -0.00266, -0.00477, -0.00334, -0.00417, -0.00448,
        -0.0051, -0.0051, -0.0057, -0.0055, -0.0059, -0.0056, -0.0058, -
        0.0057, -0.0056, -0.0056, -0.0056, -0.00044, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0]
    a6 = [1.7, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.9, 1.8, 2.0, 1.9, 2.1, 2.1, 
        0.9, 2.3, 2.5, 3.7, -3.1, 0.8, 0.1, 0.5, 0.7, -1.0, 5.8, 2.9, 2.3, 
        0.9, 2.2, 2.0, 2.0, 1.8, 1.9, 1.8, 1.8, 1.7, 1.8, 1.7, 1.7, 1.7, 
        1.7, 1.7, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-09
    ap = 1.4e-10 * (1.0 - 1.2e-05 * f ** 1.5)
    gamma0 = 0.0056 * (p + 1.1 * e) * theta ** 0.8
    nr = 2.588 * p * theta + 0.000307 * (1.0 / (1.0 + (f / gamma0) ** 2) - 1.0
        ) * p * theta * theta
    ni = (2.0 * 0.000307 / (gamma0 * (1.0 + (f / gamma0) ** 2) * (1.0 + (f /
        60.0) ** 2)) + ap * p * theta ** 2.5) * f * p * theta * theta
    for i in xrange(0, len(nu0)):
        S = a1[i] * p * theta ** 3 * math.exp(a2[i] * (1.0 - theta))
        gamma = a3[i] * (p * theta ** (0.8 - a4[i])) + 1.1 * e * theta
        delta = a5[i] * p * theta ** a6[i]
        x = (nu0[i] - f) * (nu0[i] - f) + gamma * gamma
        y = (nu0[i] + f) * (nu0[i] + f) + gamma * gamma
        z = nu0[i] + gamma * gamma / nu0[i]
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / nu0[i] + delta * (
            1.0 / x - 1.0 / y) * gamma * f / nu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / nu0[i] - delta * (
            (nu0[i] - f) / x + (nu0[i] + f) / y) * f / nu0[i])
    return complex(nr, ni)


<mask token>


def pvapsat(T):
    if T > 215:
        theta = 300.0 / T
        return 100000.0 / (41.51 * theta ** -5 * 10 ** (9.384 * theta - 10.0))
    else:
        return 0.0


def refract(t, pdry, pvap, z, n, nu, T0, el):
    HMKS = 6.6260755e-34
    KMKS = 1.380658e-23
    CMKS = 299792458
    tau = 0.0
    Tb = HMKS * nu / (KMKS * (math.exp(HMKS * nu / (KMKS * T0)) - 1))
    Ldry = 0.0
    Lvap = 0.0
    snell = math.sin(el)
    for i in xrange(n, 0, -1):
        if i == 1:
            dz = 0.5 * (z[1] - z[0])
        elif i == n:
            dz = 0.5 * (z[n] - z[n - 1])
        else:
            dz = 0.5 * (z[i + 1] - z[i - 1])
        Ndry = refdry(nu, t[i], pdry[i], pvap[i])
        Nvap = refvap(nu, t[i], pdry[i], pvap[i])
        nr = 1 + (Ndry.real + Nvap.real) * 1e-06
        ni = (Ndry.imag + Nvap.imag) * 1e-06
        l = dz * nr / math.sqrt(nr * nr + snell * snell - 1.0)
        dtau = l * 4.0 * math.pi * nu / CMKS * ni
        Tb = (Tb - t[i]) * math.exp(-dtau) + t[i]
        tau = tau + dtau
        Ldry = Ldry + l * Ndry.real * 1e-06
        Lvap = Lvap + l * Nvap.real * 1e-06
    return {'Tb': Tb, 'tau': tau, 'Ldry': Ldry, 'Lvap': Lvap}


<mask token>

<mask token>


def refdry(nu, T, Pdry, Pvap):
    nu0 = [49.452379, 49.962257, 50.474238, 50.987748, 51.50335, 52.021409,
        52.542393, 53.066906, 53.595748, 54.129999, 54.671157, 55.221365, 
        55.7838, 56.264777, 56.363387, 56.96818, 57.612481, 58.323874, 
        58.446589, 59.164204, 59.590982, 60.306057, 60.434775, 61.150558, 
        61.800152, 62.411212, 62.486253, 62.997974, 63.568515, 64.127764, 
        64.6789, 65.224067, 65.764769, 66.302088, 66.836827, 67.369595, 
        67.900862, 68.431001, 68.960306, 69.489021, 70.017342, 18.750341, 
        68.49835, 24.76312, 87.24937, 15.39315, 73.83873, 34.14533]
    a1 = [1.2e-07, 3.4e-07, 9.4e-07, 2.46e-06, 6.08e-06, 1.414e-05, 
        3.102e-05, 6.41e-05, 0.0001247, 0.000228, 0.0003918, 0.0006316, 
        0.0009535, 0.0005489, 0.001344, 0.001763, 0.002141, 0.002386, 
        0.001457, 0.002404, 0.002112, 0.002124, 0.002461, 0.002504, 
        0.002298, 0.001933, 0.001517, 0.001503, 0.001087, 0.0007335, 
        0.0004635, 0.0002748, 0.000153, 8.009e-05, 3.946e-05, 1.832e-05, 
        8.01e-06, 3.3e-06, 1.28e-06, 4.7e-07, 1.6e-07, 0.000945, 6.79e-05, 
        0.000638, 0.000235, 9.96e-05, 0.000671, 0.00018]
    a2 = [11.83, 10.72, 9.69, 8.69, 7.74, 6.84, 6.0, 5.22, 4.48, 3.81, 3.19,
        2.62, 2.115, 0.01, 1.655, 1.255, 0.91, 0.621, 0.079, 0.386, 0.207, 
        0.207, 0.386, 0.621, 0.91, 1.255, 0.078, 1.66, 2.11, 2.62, 3.19, 
        3.81, 4.48, 5.22, 6.0, 6.84, 7.74, 8.69, 9.69, 10.72, 11.83, 0.0, 
        0.02, 0.011, 0.011, 0.089, 0.079, 0.079]
    a3 = [0.0084, 0.0085, 0.0086, 0.0087, 0.0089, 0.0092, 0.0094, 0.0097, 
        0.01, 0.0102, 0.0105, 0.01079, 0.0111, 0.01646, 0.01144, 0.01181, 
        0.01221, 0.01266, 0.01449, 0.01319, 0.0136, 0.01382, 0.01297, 
        0.01248, 0.01207, 0.01171, 0.01468, 0.01139, 0.01108, 0.01078, 
        0.0105, 0.0102, 0.01, 0.0097, 0.0094, 0.0092, 0.0089, 0.0087, 
        0.0086, 0.0085, 0.0084, 0.01592, 0.0192, 0.01916, 0.0192, 0.0181, 
        0.0181, 0.0181]
    a4 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6]
    a5 = [0.0056, 0.0056, 0.0056, 0.0055, 0.0056, 0.0055, 0.0057, 0.0053, 
        0.0054, 0.0048, 0.0048, 0.00417, 0.00375, 0.00774, 0.00297, 0.00212,
        0.00094, -0.00055, 0.00597, -0.00244, 0.00344, -0.00413, 0.00132, -
        0.00036, -0.00159, -0.00266, -0.00477, -0.00334, -0.00417, -0.00448,
        -0.0051, -0.0051, -0.0057, -0.0055, -0.0059, -0.0056, -0.0058, -
        0.0057, -0.0056, -0.0056, -0.0056, -0.00044, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0]
    a6 = [1.7, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.9, 1.8, 2.0, 1.9, 2.1, 2.1, 
        0.9, 2.3, 2.5, 3.7, -3.1, 0.8, 0.1, 0.5, 0.7, -1.0, 5.8, 2.9, 2.3, 
        0.9, 2.2, 2.0, 2.0, 1.8, 1.9, 1.8, 1.8, 1.7, 1.8, 1.7, 1.7, 1.7, 
        1.7, 1.7, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-09
    ap = 1.4e-10 * (1.0 - 1.2e-05 * f ** 1.5)
    gamma0 = 0.0056 * (p + 1.1 * e) * theta ** 0.8
    nr = 2.588 * p * theta + 0.000307 * (1.0 / (1.0 + (f / gamma0) ** 2) - 1.0
        ) * p * theta * theta
    ni = (2.0 * 0.000307 / (gamma0 * (1.0 + (f / gamma0) ** 2) * (1.0 + (f /
        60.0) ** 2)) + ap * p * theta ** 2.5) * f * p * theta * theta
    for i in xrange(0, len(nu0)):
        S = a1[i] * p * theta ** 3 * math.exp(a2[i] * (1.0 - theta))
        gamma = a3[i] * (p * theta ** (0.8 - a4[i])) + 1.1 * e * theta
        delta = a5[i] * p * theta ** a6[i]
        x = (nu0[i] - f) * (nu0[i] - f) + gamma * gamma
        y = (nu0[i] + f) * (nu0[i] + f) + gamma * gamma
        z = nu0[i] + gamma * gamma / nu0[i]
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / nu0[i] + delta * (
            1.0 / x - 1.0 / y) * gamma * f / nu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / nu0[i] - delta * (
            (nu0[i] - f) / x + (nu0[i] + f) / y) * f / nu0[i])
    return complex(nr, ni)


def refvap(nu, T, Pdry, Pvap):
    mnu0 = [22.23508, 67.81396, 119.99594, 183.310117, 321.225644, 
        325.152919, 336.187, 380.197372, 390.134508, 437.346667, 439.150812,
        443.018295, 448.001075, 470.888947, 474.689127, 488.491133, 
        503.568532, 504.482692, 556.936002, 620.700807, 658.0065, 
        752.033227, 841.073593, 859.865, 899.407, 902.555, 906.205524, 
        916.171582, 970.315022, 987.926764]
    b1 = [0.109, 0.0011, 0.0007, 2.3, 0.0464, 1.54, 0.001, 11.9, 0.0044, 
        0.0637, 0.921, 0.194, 10.6, 0.33, 1.28, 0.253, 0.0374, 0.0125, 
        510.0, 5.09, 0.274, 250.0, 0.013, 0.133, 0.055, 0.038, 0.183, 8.56,
        9.16, 138.0]
    b2 = [2.143, 8.73, 8.347, 0.653, 6.156, 1.515, 9.802, 1.018, 7.318, 
        5.015, 3.561, 5.015, 1.37, 3.561, 2.342, 2.814, 6.693, 6.693, 0.114,
        2.15, 7.767, 0.336, 8.113, 7.989, 7.845, 8.36, 5.039, 1.369, 1.842,
        0.178]
    b3 = [0.02784, 0.0276, 0.027, 0.02835, 0.0214, 0.027, 0.0265, 0.0276, 
        0.019, 0.0137, 0.0164, 0.0144, 0.0238, 0.0182, 0.0198, 0.0249, 
        0.0115, 0.0119, 0.03, 0.0223, 0.03, 0.0286, 0.0141, 0.0286, 0.0286,
        0.0264, 0.0234, 0.0253, 0.024, 0.0286]
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-09
    nr = (2.39 * e * theta + 41.6 * e * theta * theta + 6.47e-06 * f ** 
        2.05 * e * theta ** 2.4)
    ni = (0.915 * 1.4e-06 * p + 5.41e-05 * e * theta * theta * theta
        ) * f * e * theta ** 2.5
    for i in xrange(0, len(mnu0)):
        S = b1[i] * e * theta ** 3.5 * math.exp(b2[i] * (1.0 - theta))
        gamma = b3[i] * (p * theta ** 0.8 + 4.8 * e * theta)
        x = (mnu0[i] - f) * (mnu0[i] - f) + gamma * gamma
        y = (mnu0[i] + f) * (mnu0[i] + f) + gamma * gamma
        z = mnu0[i] + gamma * gamma / mnu0[i]
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / mnu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / mnu0[i])
    return complex(nr, ni)


def pvapsat(T):
    if T > 215:
        theta = 300.0 / T
        return 100000.0 / (41.51 * theta ** -5 * 10 ** (9.384 * theta - 10.0))
    else:
        return 0.0


def refract(t, pdry, pvap, z, n, nu, T0, el):
    HMKS = 6.6260755e-34
    KMKS = 1.380658e-23
    CMKS = 299792458
    tau = 0.0
    Tb = HMKS * nu / (KMKS * (math.exp(HMKS * nu / (KMKS * T0)) - 1))
    Ldry = 0.0
    Lvap = 0.0
    snell = math.sin(el)
    for i in xrange(n, 0, -1):
        if i == 1:
            dz = 0.5 * (z[1] - z[0])
        elif i == n:
            dz = 0.5 * (z[n] - z[n - 1])
        else:
            dz = 0.5 * (z[i + 1] - z[i - 1])
        Ndry = refdry(nu, t[i], pdry[i], pvap[i])
        Nvap = refvap(nu, t[i], pdry[i], pvap[i])
        nr = 1 + (Ndry.real + Nvap.real) * 1e-06
        ni = (Ndry.imag + Nvap.imag) * 1e-06
        l = dz * nr / math.sqrt(nr * nr + snell * snell - 1.0)
        dtau = l * 4.0 * math.pi * nu / CMKS * ni
        Tb = (Tb - t[i]) * math.exp(-dtau) + t[i]
        tau = tau + dtau
        Ldry = Ldry + l * Ndry.real * 1e-06
        Lvap = Lvap + l * Nvap.real * 1e-06
    return {'Tb': Tb, 'tau': tau, 'Ldry': Ldry, 'Lvap': Lvap}


def calcOpacity(freq, el, t0, p0, h0):
    M = 0.02896
    R = 8.314
    Mv = 0.018
    rho0 = 1000.0
    g = 9.81
    d = 0.0065
    z0 = 1540.0
    zmax = 10000.0
    N = 50
    z = []
    T = []
    Pvap = []
    Pdry = []
    fac = []
    Tb = []
    tau = []
    ofreq = []
    for i in xrange(0, N + 1):
        zd = float(i) * zmax / float(N)
        z.append(zd)
        T.append(t0 / (1 + d / t0 * zd))
        P = p0 * math.exp(-1.0 * M * g / (R * t0) * (zd + 0.5 * d * zd * zd /
            t0))
        Pvap.append(min(h0 * math.exp(-1.0 * zd / z0) * pvapsat(T[0]),
            pvapsat(T[i])))
        Pdry.append(P - Pvap[i])
    for i in xrange(0, len(freq)):
        resref = refract(T, Pdry, Pvap, z, N, freq[i], 2.7, el)
        fac.append(math.exp(-1.0 * resref['tau']))
        Tb.append(resref['Tb'])
        tau.append(resref['tau'])
        ofreq.append(freq[i])
    return {'fac': fac, 'Tb': Tb, 'tau': tau, 'freq': ofreq}

import math


def refdry(nu, T, Pdry, Pvap):
    nu0 = [49.452379, 49.962257, 50.474238, 50.987748, 51.50335, 52.021409,
        52.542393, 53.066906, 53.595748, 54.129999, 54.671157, 55.221365, 
        55.7838, 56.264777, 56.363387, 56.96818, 57.612481, 58.323874, 
        58.446589, 59.164204, 59.590982, 60.306057, 60.434775, 61.150558, 
        61.800152, 62.411212, 62.486253, 62.997974, 63.568515, 64.127764, 
        64.6789, 65.224067, 65.764769, 66.302088, 66.836827, 67.369595, 
        67.900862, 68.431001, 68.960306, 69.489021, 70.017342, 18.750341, 
        68.49835, 24.76312, 87.24937, 15.39315, 73.83873, 34.14533]
    a1 = [1.2e-07, 3.4e-07, 9.4e-07, 2.46e-06, 6.08e-06, 1.414e-05, 
        3.102e-05, 6.41e-05, 0.0001247, 0.000228, 0.0003918, 0.0006316, 
        0.0009535, 0.0005489, 0.001344, 0.001763, 0.002141, 0.002386, 
        0.001457, 0.002404, 0.002112, 0.002124, 0.002461, 0.002504, 
        0.002298, 0.001933, 0.001517, 0.001503, 0.001087, 0.0007335, 
        0.0004635, 0.0002748, 0.000153, 8.009e-05, 3.946e-05, 1.832e-05, 
        8.01e-06, 3.3e-06, 1.28e-06, 4.7e-07, 1.6e-07, 0.000945, 6.79e-05, 
        0.000638, 0.000235, 9.96e-05, 0.000671, 0.00018]
    a2 = [11.83, 10.72, 9.69, 8.69, 7.74, 6.84, 6.0, 5.22, 4.48, 3.81, 3.19,
        2.62, 2.115, 0.01, 1.655, 1.255, 0.91, 0.621, 0.079, 0.386, 0.207, 
        0.207, 0.386, 0.621, 0.91, 1.255, 0.078, 1.66, 2.11, 2.62, 3.19, 
        3.81, 4.48, 5.22, 6.0, 6.84, 7.74, 8.69, 9.69, 10.72, 11.83, 0.0, 
        0.02, 0.011, 0.011, 0.089, 0.079, 0.079]
    a3 = [0.0084, 0.0085, 0.0086, 0.0087, 0.0089, 0.0092, 0.0094, 0.0097, 
        0.01, 0.0102, 0.0105, 0.01079, 0.0111, 0.01646, 0.01144, 0.01181, 
        0.01221, 0.01266, 0.01449, 0.01319, 0.0136, 0.01382, 0.01297, 
        0.01248, 0.01207, 0.01171, 0.01468, 0.01139, 0.01108, 0.01078, 
        0.0105, 0.0102, 0.01, 0.0097, 0.0094, 0.0092, 0.0089, 0.0087, 
        0.0086, 0.0085, 0.0084, 0.01592, 0.0192, 0.01916, 0.0192, 0.0181, 
        0.0181, 0.0181]
    a4 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0, 0.0, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6]
    a5 = [0.0056, 0.0056, 0.0056, 0.0055, 0.0056, 0.0055, 0.0057, 0.0053, 
        0.0054, 0.0048, 0.0048, 0.00417, 0.00375, 0.00774, 0.00297, 0.00212,
        0.00094, -0.00055, 0.00597, -0.00244, 0.00344, -0.00413, 0.00132, -
        0.00036, -0.00159, -0.00266, -0.00477, -0.00334, -0.00417, -0.00448,
        -0.0051, -0.0051, -0.0057, -0.0055, -0.0059, -0.0056, -0.0058, -
        0.0057, -0.0056, -0.0056, -0.0056, -0.00044, 0.0, 0.0, 0.0, 0.0, 
        0.0, 0.0]
    a6 = [1.7, 1.7, 1.7, 1.7, 1.8, 1.8, 1.8, 1.9, 1.8, 2.0, 1.9, 2.1, 2.1, 
        0.9, 2.3, 2.5, 3.7, -3.1, 0.8, 0.1, 0.5, 0.7, -1.0, 5.8, 2.9, 2.3, 
        0.9, 2.2, 2.0, 2.0, 1.8, 1.9, 1.8, 1.8, 1.7, 1.8, 1.7, 1.7, 1.7, 
        1.7, 1.7, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-09
    ap = 1.4e-10 * (1.0 - 1.2e-05 * f ** 1.5)
    gamma0 = 0.0056 * (p + 1.1 * e) * theta ** 0.8
    nr = 2.588 * p * theta + 0.000307 * (1.0 / (1.0 + (f / gamma0) ** 2) - 1.0
        ) * p * theta * theta
    ni = (2.0 * 0.000307 / (gamma0 * (1.0 + (f / gamma0) ** 2) * (1.0 + (f /
        60.0) ** 2)) + ap * p * theta ** 2.5) * f * p * theta * theta
    for i in xrange(0, len(nu0)):
        S = a1[i] * p * theta ** 3 * math.exp(a2[i] * (1.0 - theta))
        gamma = a3[i] * (p * theta ** (0.8 - a4[i])) + 1.1 * e * theta
        delta = a5[i] * p * theta ** a6[i]
        x = (nu0[i] - f) * (nu0[i] - f) + gamma * gamma
        y = (nu0[i] + f) * (nu0[i] + f) + gamma * gamma
        z = nu0[i] + gamma * gamma / nu0[i]
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / nu0[i] + delta * (
            1.0 / x - 1.0 / y) * gamma * f / nu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / nu0[i] - delta * (
            (nu0[i] - f) / x + (nu0[i] + f) / y) * f / nu0[i])
    return complex(nr, ni)


def refvap(nu, T, Pdry, Pvap):
    mnu0 = [22.23508, 67.81396, 119.99594, 183.310117, 321.225644, 
        325.152919, 336.187, 380.197372, 390.134508, 437.346667, 439.150812,
        443.018295, 448.001075, 470.888947, 474.689127, 488.491133, 
        503.568532, 504.482692, 556.936002, 620.700807, 658.0065, 
        752.033227, 841.073593, 859.865, 899.407, 902.555, 906.205524, 
        916.171582, 970.315022, 987.926764]
    b1 = [0.109, 0.0011, 0.0007, 2.3, 0.0464, 1.54, 0.001, 11.9, 0.0044, 
        0.0637, 0.921, 0.194, 10.6, 0.33, 1.28, 0.253, 0.0374, 0.0125, 
        510.0, 5.09, 0.274, 250.0, 0.013, 0.133, 0.055, 0.038, 0.183, 8.56,
        9.16, 138.0]
    b2 = [2.143, 8.73, 8.347, 0.653, 6.156, 1.515, 9.802, 1.018, 7.318, 
        5.015, 3.561, 5.015, 1.37, 3.561, 2.342, 2.814, 6.693, 6.693, 0.114,
        2.15, 7.767, 0.336, 8.113, 7.989, 7.845, 8.36, 5.039, 1.369, 1.842,
        0.178]
    b3 = [0.02784, 0.0276, 0.027, 0.02835, 0.0214, 0.027, 0.0265, 0.0276, 
        0.019, 0.0137, 0.0164, 0.0144, 0.0238, 0.0182, 0.0198, 0.0249, 
        0.0115, 0.0119, 0.03, 0.0223, 0.03, 0.0286, 0.0141, 0.0286, 0.0286,
        0.0264, 0.0234, 0.0253, 0.024, 0.0286]
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-09
    nr = (2.39 * e * theta + 41.6 * e * theta * theta + 6.47e-06 * f ** 
        2.05 * e * theta ** 2.4)
    ni = (0.915 * 1.4e-06 * p + 5.41e-05 * e * theta * theta * theta
        ) * f * e * theta ** 2.5
    for i in xrange(0, len(mnu0)):
        S = b1[i] * e * theta ** 3.5 * math.exp(b2[i] * (1.0 - theta))
        gamma = b3[i] * (p * theta ** 0.8 + 4.8 * e * theta)
        x = (mnu0[i] - f) * (mnu0[i] - f) + gamma * gamma
        y = (mnu0[i] + f) * (mnu0[i] + f) + gamma * gamma
        z = mnu0[i] + gamma * gamma / mnu0[i]
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / mnu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / mnu0[i])
    return complex(nr, ni)


def pvapsat(T):
    if T > 215:
        theta = 300.0 / T
        return 100000.0 / (41.51 * theta ** -5 * 10 ** (9.384 * theta - 10.0))
    else:
        return 0.0


def refract(t, pdry, pvap, z, n, nu, T0, el):
    HMKS = 6.6260755e-34
    KMKS = 1.380658e-23
    CMKS = 299792458
    tau = 0.0
    Tb = HMKS * nu / (KMKS * (math.exp(HMKS * nu / (KMKS * T0)) - 1))
    Ldry = 0.0
    Lvap = 0.0
    snell = math.sin(el)
    for i in xrange(n, 0, -1):
        if i == 1:
            dz = 0.5 * (z[1] - z[0])
        elif i == n:
            dz = 0.5 * (z[n] - z[n - 1])
        else:
            dz = 0.5 * (z[i + 1] - z[i - 1])
        Ndry = refdry(nu, t[i], pdry[i], pvap[i])
        Nvap = refvap(nu, t[i], pdry[i], pvap[i])
        nr = 1 + (Ndry.real + Nvap.real) * 1e-06
        ni = (Ndry.imag + Nvap.imag) * 1e-06
        l = dz * nr / math.sqrt(nr * nr + snell * snell - 1.0)
        dtau = l * 4.0 * math.pi * nu / CMKS * ni
        Tb = (Tb - t[i]) * math.exp(-dtau) + t[i]
        tau = tau + dtau
        Ldry = Ldry + l * Ndry.real * 1e-06
        Lvap = Lvap + l * Nvap.real * 1e-06
    return {'Tb': Tb, 'tau': tau, 'Ldry': Ldry, 'Lvap': Lvap}


def calcOpacity(freq, el, t0, p0, h0):
    M = 0.02896
    R = 8.314
    Mv = 0.018
    rho0 = 1000.0
    g = 9.81
    d = 0.0065
    z0 = 1540.0
    zmax = 10000.0
    N = 50
    z = []
    T = []
    Pvap = []
    Pdry = []
    fac = []
    Tb = []
    tau = []
    ofreq = []
    for i in xrange(0, N + 1):
        zd = float(i) * zmax / float(N)
        z.append(zd)
        T.append(t0 / (1 + d / t0 * zd))
        P = p0 * math.exp(-1.0 * M * g / (R * t0) * (zd + 0.5 * d * zd * zd /
            t0))
        Pvap.append(min(h0 * math.exp(-1.0 * zd / z0) * pvapsat(T[0]),
            pvapsat(T[i])))
        Pdry.append(P - Pvap[i])
    for i in xrange(0, len(freq)):
        resref = refract(T, Pdry, Pvap, z, N, freq[i], 2.7, el)
        fac.append(math.exp(-1.0 * resref['tau']))
        Tb.append(resref['Tb'])
        tau.append(resref['tau'])
        ofreq.append(freq[i])
    return {'fac': fac, 'Tb': Tb, 'tau': tau, 'freq': ofreq}

# Module refract.py
import math

def refdry(nu, T, Pdry, Pvap):
    # From Miriad: Determine the complex refractivity of the dry components
    # of the atmosphere.
    #
    # Input:
    #  nu = observing frequency (Hz)
    #  T = temperature (K)
    #  Pdry = partial pressure of dry components (Pa)
    #  Pvap = partial pressure of water vapour (Pa)

    # Table of microwave oxygen lines and their parameters.
    nu0 = [ 49.452379, 49.962257, 50.474238, 50.987748, 51.503350, 52.021409,
            52.542393, 53.066906, 53.595748, 54.129999, 54.671157, 55.221365,
            55.783800, 56.264777, 56.363387, 56.968180, 57.612481, 58.323874,
            58.446589, 59.164204, 59.590982, 60.306057, 60.434775, 61.150558,
            61.800152, 62.411212, 62.486253, 62.997974, 63.568515, 64.127764,
            64.678900, 65.224067, 65.764769, 66.302088, 66.836827, 67.369595,
            67.900862, 68.431001, 68.960306, 69.489021, 70.017342, 18.750341,
            68.498350, 24.763120, 87.249370, 15.393150, 73.838730, 34.145330 ]
    a1 = [    0.12E-6,    0.34E-6,    0.94E-6,    2.46E-6,    6.08E-6,   14.14E-6,
              31.02E-6,   64.10E-6,  124.70E-6,  228.00E-6,  391.80E-6,  631.60E-6,
              953.50E-6,  548.90E-6, 1344.00E-6, 1763.00E-6, 2141.00E-6, 2386.00E-6,
              1457.00E-6, 2404.00E-6, 2112.00E-6, 2124.00E-6, 2461.00E-6, 2504.00E-6,
              2298.00E-6, 1933.00E-6, 1517.00E-6, 1503.00E-6, 1087.00E-6,  733.50E-6,
              463.50E-6,  274.80E-6,  153.00E-6,   80.09E-6,   39.46E-6,   18.32E-6,
              8.01E-6,    3.30E-6,    1.28E-6,    0.47E-6,    0.16E-6,  945.00E-6,
              67.90E-6,  638.00E-6,  235.00E-6,   99.60E-6,  671.00E-6,  180.00E-6 ]
    a2 = [  11.830, 10.720,  9.690,  8.690,  7.740,  6.840,
            6.000,  5.220,  4.480,  3.810,  3.190,  2.620,
            2.115,  0.010,  1.655,  1.255,  0.910,  0.621,
            0.079,  0.386,  0.207,  0.207,  0.386,  0.621,
            0.910,  1.255,  0.078,  1.660,  2.110,  2.620,
            3.190,  3.810,  4.480,  5.220,  6.000,  6.840,
            7.740,  8.690,  9.690, 10.720, 11.830,  0.000,
            0.020,  0.011,  0.011,  0.089,  0.079,  0.079 ]
    a3 = [   8.40E-3,  8.50E-3,  8.60E-3,  8.70E-3,  8.90E-3,  9.20E-3,
             9.40E-3,  9.70E-3, 10.00E-3, 10.20E-3, 10.50E-3, 10.79E-3,
             11.10E-3, 16.46E-3, 11.44E-3, 11.81E-3, 12.21E-3, 12.66E-3,
             14.49E-3, 13.19E-3, 13.60E-3, 13.82E-3, 12.97E-3, 12.48E-3,
             12.07E-3, 11.71E-3, 14.68E-3, 11.39E-3, 11.08E-3, 10.78E-3,
             10.50E-3, 10.20E-3, 10.00E-3,  9.70E-3,  9.40E-3,  9.20E-3,
             8.90E-3,  8.70E-3,  8.60E-3,  8.50E-3,  8.40E-3, 15.92E-3,
             19.20E-3, 19.16E-3, 19.20E-3, 18.10E-3, 18.10E-3, 18.10E-3 ]
    a4 = [  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
            0.6, 0.6, 0.6, 0.6, 0.6, 0.6 ]
    a5 = [  5.60E-3,  5.60E-3,  5.60E-3,  5.50E-3,  5.60E-3,  5.50E-3,
            5.70E-3,  5.30E-3,  5.40E-3,  4.80E-3,  4.80E-3,  4.17E-3,
            3.75E-3,  7.74E-3,  2.97E-3,  2.12E-3,  0.94E-3, -0.55E-3,
            5.97E-3, -2.44E-3,  3.44E-3, -4.13E-3,  1.32E-3, -0.36E-3,
            -1.59E-3, -2.66E-3, -4.77E-3, -3.34E-3, -4.17E-3, -4.48E-3,
            -5.10E-3, -5.10E-3, -5.70E-3, -5.50E-3, -5.90E-3, -5.60E-3,
            -5.80E-3, -5.70E-3, -5.60E-3, -5.60E-3, -5.60E-3, -0.44E-3,
            0.00E00,  0.00E00,  0.00E00,  0.00E00,  0.00E00,  0.00E00 ]
    a6 = [  1.7,  1.7,  1.7,  1.7,  1.8,  1.8,
            1.8,  1.9,  1.8,  2.0,  1.9,  2.1,
            2.1,  0.9,  2.3,  2.5,  3.7, -3.1,
            0.8,  0.1,  0.5,  0.7, -1.0,  5.8,
            2.9,  2.3,  0.9,  2.2,  2.0,  2.0,
            1.8,  1.9,  1.8,  1.8,  1.7,  1.8,
            1.7,  1.7,  1.7,  1.7,  1.7,  0.9,
            1.0,  1.0,  1.0,  1.0,  1.0,  1.0 ]

    # Convert to the units of Liebe.
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-9

    ap = 1.4e-10 * (1.0 - 1.2e-5 * f ** 1.5)
    gamma0 = 5.6e-3 * (p + 1.1 * e) * theta ** 0.8
    nr = 2.588 * p * theta + 3.07e-4 * (1.0 / (1.0 + (f / gamma0) ** 2) - 1.0) * p * theta * theta
    ni = (2.0 * 3.07e-4 / (gamma0 * (1.0 + (f / gamma0) ** 2) * (1.0 + (f / 60.0) ** 2)) +
          ap * p * theta ** 2.5) * f * p * theta * theta

    # Sum the contributions of the lines.
    for i in xrange(0, len(nu0)):
        S = a1[i] * p * theta ** 3 * math.exp(a2[i] * (1.0 - theta))
        gamma = a3[i] * (p * theta ** (0.8 - a4[i])) + 1.1 * e * theta
        delta = a5[i] * p * theta ** a6[i]
        x = (nu0[i] - f) * (nu0[i] - f) + gamma * gamma
        y = (nu0[i] + f) * (nu0[i] + f) + gamma * gamma
        z = (nu0[i] + gamma * gamma / nu0[i])
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / nu0[i] + delta *
                       (1.0 / x - 1.0 / y) * gamma * f / nu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / nu0[i] - delta *
                       ((nu0[i] - f) / x + (nu0[i] + f) / y) * f / nu0[i])

    # Return the result.
    return complex(nr, ni)

def refvap(nu, T, Pdry, Pvap):
    # From Miriad; Determine the complex refractivity of the water vapour monomers.
    #
    # Inputs:
    #  nu = observating frequency (Hz)
    #  T = temperature (K)
    #  Pdry = partial pressure of dry components (Pa)
    #  Pvap = partial pressure of water vapour (Pa)

    # Table of the microwave water lines.
    mnu0 = [  22.235080,  67.813960, 119.995940, 183.310117, 321.225644, 325.152919,
              336.187000, 380.197372, 390.134508, 437.346667, 439.150812, 443.018295,
              448.001075, 470.888947, 474.689127, 488.491133, 503.568532, 504.482692,
              556.936002, 620.700807, 658.006500, 752.033227, 841.073593, 859.865000,
              899.407000, 902.555000, 906.205524, 916.171582, 970.315022, 987.926764 ]
    b1 = [   0.1090,  0.0011,  0.0007,   2.3000,  0.0464,  1.5400,
             0.0010, 11.9000,  0.0044,   0.0637,  0.9210,  0.1940,
             10.6000,  0.3300,  1.2800,   0.2530,  0.0374,  0.0125,
             510.0000,  5.0900,  0.2740, 250.0000,  0.0130,  0.1330,
             0.0550,  0.0380,  0.1830,   8.5600,  9.1600, 138.000 ]
    b2 = [ 2.143, 8.730, 8.347, 0.653, 6.156, 1.515,
           9.802, 1.018, 7.318, 5.015, 3.561, 5.015,
           1.370, 3.561, 2.342, 2.814, 6.693, 6.693,
           0.114, 2.150, 7.767, 0.336, 8.113, 7.989,
           7.845, 8.360, 5.039, 1.369, 1.842, 0.178 ]
    b3 = [ 27.84E-3, 27.60E-3, 27.00E-3, 28.35E-3, 21.40E-3, 27.00E-3,
           26.50E-3, 27.60E-3, 19.00E-3, 13.70E-3, 16.40E-3, 14.40E-3,
           23.80E-3, 18.20E-3, 19.80E-3, 24.90E-3, 11.50E-3, 11.90E-3,
           30.00E-3, 22.30E-3, 30.00E-3, 28.60E-3, 14.10E-3, 28.60E-3,
           28.60E-3, 26.40E-3, 23.40E-3, 25.30E-3, 24.00E-3, 28.60E-3 ]

    # Convert to the units of Liebe.
    theta = 300.0 / T
    e = 0.001 * Pvap
    p = 0.001 * Pdry
    f = nu * 1e-9

    nr = 2.39 * e * theta + 41.6 * e * theta * theta + 6.47e-6 * f ** 2.05 * e * theta ** 2.4
    ni = (0.915 * 1.40e-6 * p + 5.41e-5 * e * theta * theta * theta) * f * e * theta ** 2.5

    # Sum the contributions of the lines.
    for i in xrange(0, len(mnu0)):
        S = b1[i] * e * theta ** 3.5 * math.exp(b2[i] * (1.0 - theta))
        gamma = b3[i] * (p * theta ** 0.8 + 4.80 * e * theta)
        x = (mnu0[i] - f) * (mnu0[i] - f) + gamma * gamma
        y = (mnu0[i] + f) * (mnu0[i] + f) + gamma * gamma
        z = (mnu0[i] + gamma * gamma / mnu0[i])
        nr = nr + S * ((z - f) / x + (z + f) / y - 2 / mnu0[i])
        ni = ni + S * ((1.0 / x + 1.0 / y) * gamma * f / mnu0[i])

    # Return the result.
    return complex(nr, ni)
    
def pvapsat(T):
    # From Miriad; Determine the saturation pressure of water vapour.
    # Input:
    #  T = temperature (K)
    #
    # Output:
    #  vapour saturation pressure (Pa)

    if (T > 215):
        theta = 300.0 / T
        return 1e5 / (41.51 * (theta ** -5) * (10 **(9.384 * theta - 10.0)))
    else:
        return 0.0

def refract(t, pdry, pvap, z, n, nu, T0, el):
    # From Miriad; Compute refractive index for an atmosphere.
    # Determine the sky brightness and excess path lengths for a parallel
    # slab atmosphere. Liebe's model (1985) is used to determine the complex
    # refractive index of air.
    #
    # Input:
    #  n = the number of atmospheric layers.
    #  t = temperature of the layers. T[0] is the temperature at the lowest
    #      layer (K)
    #  Pdry = partial pressure of the dry components (Pa)
    #  Pvap = partial pressure of the water vapour components (Pa)
    #  z = height of the layer.
    #  nu = frequency of interest (Hz)
    #  T0 = astronomical brightness temperature (K)
    #  el = elevation angle of the source above the atmosphere (rad)
    #
    # Output:
    #  { 'Tb' = brightness temperature (K),
    #    'tau' = opacity (nepers)
    #    'Ldry' = excess path, dry component (m)
    #    'Lvap' = excess path, water vapour component (m) }

    # Some constants.
    HMKS = 6.6260755e-34 # Planck constant, J.s
    KMKS = 1.380658e-23 # Boltzmann constant, J/K
    CMKS = 299792458 # Speed of light, m/s

    tau = 0.0
    Tb = HMKS * nu / (KMKS * (math.exp(HMKS * nu / (KMKS * T0)) - 1))
    Ldry = 0.0
    Lvap = 0.0

    snell = math.sin(el)
    for i in xrange(n, 0, -1):
        if (i == 1):
            dz = 0.5 * (z[1] - z[0])
        elif (i == n):
            dz = 0.5 * (z[n] - z[n - 1])
        else:
            dz = 0.5 * (z[i + 1] - z[i - 1])
        Ndry = refdry(nu, t[i], pdry[i], pvap[i])
        Nvap = refvap(nu, t[i], pdry[i], pvap[i])
        nr = 1 + (Ndry.real + Nvap.real) * 1e-6
        ni = (Ndry.imag + Nvap.imag) * 1e-6
        l = dz * nr / math.sqrt(nr * nr + (snell * snell) - 1.0)
        dtau = l * 4.0 * math.pi * nu / CMKS * ni
        Tb = (Tb - t[i]) * math.exp(-dtau) + t[i]
        tau = tau + dtau
        Ldry = Ldry + l * Ndry.real * 1e-6
        Lvap = Lvap + l * Nvap.real * 1e-6

    return { 'Tb': Tb, 'tau': tau, 'Ldry': Ldry, 'Lvap': Lvap }

def calcOpacity(freq, el, t0, p0, h0):
    # From Miriad; Compute sky brightness and opacity of a model atmosphere.
    # Returns the transmissivity of the atmosphere given frequency, elevation
    # angle and meteorological data. This uses a simple model of the atmosphere
    # and Liebe's model (1985) of the complex refractive index of air.
    # Input:
    #  freq = frequency (Hz)
    #  el = elevation angle (radians)
    #  t0,p0,h0 = Met data; observatory temperature, pressure and humidity
    #             (K, Pa, fraction)
    # Output:
    #  { 'fac' = transmissivity (fraction between 0 and 1)
    #    'Tb' = sky brightness temperature (K) }

    # Atmospheric parameters.
    M = 28.96e-3
    R = 8.314
    Mv = 18e-3
    rho0 = 1e3
    g = 9.81

    # d = temperature lapse rate (K/m)
    d = 0.0065
    # z0 = water vapour scale height (m)
    z0 = 1540.0
    # zmax = max altitude of model atmosphere (m)
    zmax = 10000.0

    # Generate a model of the atmosphere - T is temperature, Pdry is
    # partial pressure of "dry" constituents, Pvap is the partial
    # pressure of the water vapour.
    N = 50 # Number of iterations through the atmosphere in height
    z = []
    T = []
    Pvap = []
    Pdry = []
    fac = []
    Tb = []
    tau = []
    ofreq = []
    for i in xrange(0, (N + 1)):
        zd = float(i) * zmax / float(N)
        z.append(zd)
        T.append(t0 / (1 + d / t0 * zd))
        P = p0 * math.exp(-1.0 * M * g / (R * t0) * (zd + 0.5 * d * zd * zd / t0))
        Pvap.append(min(h0 * math.exp(-1.0 * zd / z0) * pvapsat(T[0]),
                        pvapsat(T[i])))
        Pdry.append(P - Pvap[i])

    # Determine the transmissivity and sky brightness.
    #job_server = pp.Server(secret="hello")
    #jobs = [(inpf, job_server.submit(refract, (T, Pdry, Pvap, z, N, inpf, 2.7, el),
    #                                 (refdry, refvap, pvapsat), ("math",))) for inpf in freq]
    for i in xrange(0, len(freq)):
        #print "calculating atmosphere for frequency %.1f" % freq[i]
        resref = refract(T, Pdry, Pvap, z, N, freq[i], 2.7, el)
        fac.append(math.exp(-1.0 * resref['tau']))
        Tb.append(resref['Tb'])
        tau.append(resref['tau'])
        ofreq.append(freq[i])
    #for inpf, job in jobs:
    #    resref = job()
    #    fac.append(math.exp(-1.0 * resref['tau']))
    #    Tb.append(resref['Tb'])
    #    tau.append(resref['tau'])
    #    ofreq.append(inpf)

    return { 'fac': fac, 'Tb': Tb, 'tau': tau, 'freq': ofreq }
