From 2159f428c9e26f133b3dafee8e7ce7aae4fd78ce Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 17 Aug 2017 09:59:47 -0600
Subject: [PATCH 01/16] Add initial code for full VCS modelling

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx  | 301 ++++++++++++++++++
 turbustat/statistics/vca_vcs/vcs_models.py    | 167 ++++++++++
 2 files changed, 468 insertions(+)
 create mode 100644 turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
 create mode 100644 turbustat/statistics/vca_vcs/vcs_models.py

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
new file mode 100644
index 00000000..f4729b6e
--- /dev/null
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -0,0 +1,301 @@
+
+import numpy as np
+from scipy.integrate import quad
+from scipy.special import gamma, erf, erfc
+cimport cython
+
+from libc.math cimport sqrt, exp, sin, cos, atan
+
+cdef double pi = np.pi
+
+
+def C_eps(double r, double k_cut, double alphae, double norm_factor):
+    '''
+
+    k_cut is k0 for alphae > 3 and k1 for alphae<3
+    '''
+
+    cdef double Int
+
+    # Steep
+    if alphae > 3:
+        # Eq. B2
+        Int = Int3(r, k_cut, alphae)[0]
+
+        return 1 - 4 * pi * Int / norm_factor
+    elif alphae == 3:
+        # No density dependence.
+        return 1. / norm_factor
+    # Shallow
+    elif alphae > 1:
+        # For alphae between 1 and 3
+
+        Int = Int4(r, k_cut, alphae)[0]
+
+        return 1 - 4 * pi * Int / norm_factor
+
+    else:
+        raise ValueError("Solution not defined for alphae <= 1.")
+
+
+def Dz(double r, double theta, double V0, double k0, double alphav):
+    '''
+    Eq. A3 to A5.
+
+    I've only included the solutions for a solenoidal velocity
+    field here. The solution for a potential field is defined in Eqs. A8 & A9
+    and can be easily included here. Do they need to be?
+    '''
+
+    cdef double intone, inttwo, I_C, I_S
+
+    intone = Int1(r, k0, alphav)[0]
+    inttwo = Int2(r, k0, alphav)[0]
+
+    I_C = (4 / 3.) * intone
+    I_S = 2 * (inttwo - intone / 3.)
+
+    return 4 * pi * V0**2 * r**(alphav - 3) * \
+        (I_C * cos(theta)**2 + I_S * sin(theta)**2)
+
+
+def F_eps_norm(double alphae, double k_cut):
+    '''
+    Normalization for F_eps.
+
+    If alphae > 3, k_cut is k0 (Eq. B1). If alphae < 3, k_cut is k1 (Eq. B5).
+
+    '''
+
+    cdef double prefactor
+
+    prefactor = 2 * pi * k_cut**(3 - alphae)
+
+    if alphae > 3:
+        return prefactor * gamma(0.5 * (alphae - 3))
+    elif alphae < 3:
+        return prefactor * gamma(0.5 * (3 - alphae))
+    else:
+        # If alphae == 3, the density spectrum has no contribution
+        return 1.
+
+# Integrals in Dz and C_eps
+
+def Int1(double r, double k0, double alphav):
+    '''
+    Integral in Eq. A4
+    '''
+    def integrand(double q):
+        cdef double out
+        out = q**(2 - alphav) * exp(-(k0 * r / q)**2) * \
+            (1 - (3 / q**2) * (sin(q) / q - cos(q)))
+        return out
+
+    cdef double value, err
+
+    value, err = quad(integrand, 0, np.inf)
+
+    return value, err
+
+
+def Int2(double r, double k0, double alphav):
+    '''
+    First integral in Eq. A5
+    '''
+
+    def integrand(double q):
+        cdef double out
+        out = q**(2 - alphav) * exp(-(k0 * r / q)**2) * \
+            (1 - sin(q) / q)
+        return out
+
+    cdef double value, err
+
+    value, err = quad(integrand, 0, np.inf)
+
+    return value, err
+
+
+def Int3(double r, double k0, double alphae):
+    '''
+    Eq. B2 (w/o 4pi constant)
+    '''
+
+    def integrand(double k):
+        cdef double out
+        out = k**(1 - alphae) * exp(-(k0 / k)**2) * \
+            (sin(k * r) / r)
+        return out
+
+    cdef double value, err
+
+    value, err = quad(integrand, 0, np.inf)
+
+    return value, err
+
+
+def Int4(double r, double k1, double alphae):
+    '''
+    Eq. B7
+    '''
+
+    def integrand(double q):
+        cdef double out
+        out = q**(1 - alphae) * exp(-(q / (k1 * r))**2) * sin(q)
+        return out
+
+    cdef double value, err
+
+    value, err = quad(integrand, 0, np.inf)
+
+    return value, err
+
+
+# Window functions
+
+def gaussian_beam(double theta, double theta_0):
+    '''
+    Normalized Gaussian
+
+    Note that Eq. 33 and 34 in CL06 are missing the sqrt in the normalization.
+    '''
+    cdef double out
+    out = exp(-(theta / theta_0)**2) / sqrt(pi * theta_0**2)
+    return out
+
+
+def gaussian_autocorr(double R, double z_0, double theta_0):
+    '''
+    Gaussian autocorrelation function for a circular Gaussian beam defined
+    in the projected frame (:math:`\vec{R}` in CL06). This is the solution for
+    Eq. 32.
+
+    For W_b:
+    :math:`-\frac{e^{-\frac{R^2}{2 \theta ^2 z^2}} \left(\sqrt{\pi } R e^{\frac{R^2}{4 \theta ^2 z^2}} \text{erfc}\left(\frac{R}{2 \theta  z}\right)-2 \theta  z\right)}{4 \theta  z^3}`
+
+    w_b was solved for using Mathematica. My by-hand solution had the same
+    asymptotic behaviour, but looked so much worse than the compacter form
+    above.
+
+    '''
+
+    cdef double ratio_term, exp_term
+
+    ratio_term = R / (theta_0 * z_0)
+    exp_term = ratio_term**2
+
+    term1 = - 1. / (4 * theta_0 * z_0**3)
+    term2 = sqrt(pi) * R * exp(- 0.25 * exp_term) * \
+        erfc(0.5 * ratio_term)
+    term3 = 2 * theta_0 * z_0 * exp(- 0.5 * exp_term)
+
+    return term1 * (term2 - term3)
+
+
+def slab_autocorr(double z, double z_0, double z_1):
+    '''
+    Slab model. The function is 1 for z within z_0 and z_1. See Eq. 40.
+    This is a normalized version, so the values are 1 / (z_1 - z_0).
+
+    The autocorrelation function is just the square, so
+    :math:`(z_1 - z_0)^{-2}`.
+
+    '''
+
+    cdef double out
+
+    if z >= z_0 and z <= z_1:
+        out = (z_1 - z_0)**2
+        return out
+
+
+def pencil_beam_gaussian_z(double z, double sigma_z):
+    '''
+    Pencil beam with a Gaussian w_eps_a.
+    '''
+
+    return gaussian_beam(z, 2 * sigma_z)
+
+
+def pencil_beam_slab_z(double z, double z0, double z1):
+    '''
+    Pencil beam with a slab w_eps_a.
+    '''
+
+    return slab_autocorr(z, z0, z1)
+
+
+def gaussian_beam_slab_z_crossing(double R, double z, double z0, double z1,
+                                  double theta0):
+    '''
+    Eq. 41 -- Gaussian beam with w_eps as tophat between z0 and z1.
+
+    Same as Eq.14 in Chepurnov+10 (called g(r) there).
+
+    Appropriate for limited resolution with nearby emission.
+
+    '''
+
+    cdef double term1, p, term_a, term_b, term2, term3, term4
+
+    term1 = - 1 / (2 * pi**0.5 * theta0 * R * z)
+
+    p = (z1 + z0) / (z1 - z0)
+
+    term_a = sqrt(2 * z0**2 + p * z**2)
+    term_b = sqrt(2 * z1**2 - p * z**2)
+
+    term2 = atan(1 + 2 * z0 / z) + atan(1 - 2 * z1 / z)
+    term3 = 1 / term_a - 1 / term_b
+
+    term4 = erf(R / (theta0 * term_a) - erf(R / (theta0 * term_b)))
+
+    return term1 * (term2 / term3) * term4
+
+
+def gaussian_beam_slab_z_parallel(double R, double z, double z0, double z1,
+                                  double theta0):
+    '''
+    Gaussian beam with w_eps_a as a slab (or tophat) with parallel LOS.
+
+    B/c this is in the parallel limit, the two component for the beam and
+    the structure along the LOS are independent. This allows us to easily
+    write out the components without needing to solve for the special case,
+    as was done for gaussian_beam_slab_z_crossing.
+
+    '''
+
+    return slab_autocorr(z, z0, z1) * gaussian_autocorr(R, z0, theta0)
+
+
+def gaussian_beam_gaussian_z_parallel(double R, double z, double z_0,
+                                      double sigma_z, double theta0):
+    '''
+    Solution to Eq. 31 assuming a Gaussian beam and a Gaussian for w_eps.
+
+    * Gaussian beam has width theta0 (not FWHM, 2 sigma)
+    * Gaussian for structure along LOS at a distance of z0 and a width of
+      sigma_z
+
+    The auto-correlation function is presented for both of these.
+
+    For w_eps:
+    :math:`e^{-z^2 / 4\sigma_z^2} / \sqrt{4 \pi \sigma_z^2}`
+
+    This is independent of z_0 as it is the distance to the object for this
+    definition. The autocorrelation function is for the object thickness along
+    the LOS, and so should be independent of distance.
+
+    For W_b:
+    :math:`-\frac{e^{-\frac{R^2}{2 \theta ^2 z^2}} \left(\sqrt{\pi } R e^{\frac{R^2}{4 \theta ^2 z^2}} \text{erfc}\left(\frac{R}{2 \theta  z}\right)-2 \theta  z\right)}{4 \theta  z^3}`
+
+    '''
+
+    cdef double w_eps_a, w_b_a
+
+    # See form above. Equal to 2 * sigma_z for "theta_0"
+    w_eps_a = gaussian_beam(z, 2 * sigma_z)
+
+    w_b_a = gaussian_autocorr(R, z_0, theta0)
+
+    return w_eps_a * w_b_a
diff --git a/turbustat/statistics/vca_vcs/vcs_models.py b/turbustat/statistics/vca_vcs/vcs_models.py
new file mode 100644
index 00000000..51612db4
--- /dev/null
+++ b/turbustat/statistics/vca_vcs/vcs_models.py
@@ -0,0 +1,167 @@
+
+import numpy as np
+from scipy.integrate import tplquad
+from astropy import constants as const
+import astropy.units as u
+from functools import partial
+import math
+
+try:
+    import vegas
+    VEGAS_IMPORT_FLAG = True
+except ImportError:
+    VEGAS_IMPORT_FLAG = False
+
+
+from vcs_model_helpers import (Dz, C_eps, F_eps_norm, pencil_beam_slab_z,
+                               pencil_beam_gaussian_z,
+                               gaussian_beam_slab_z_parallel,
+                               gaussian_beam_slab_z_crossing,
+                               gaussian_beam_gaussian_z_parallel)
+
+
+def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
+       beam_type='gaussian', object_type='gaussian', los_type='parallel',
+       z0=None, z1=None, sigma_z=None, theta0=None, k1=None,
+       integration_type='mc', **integration_kwargs):
+    '''
+    VCS model.
+    '''
+
+    fk2 = f_k(kv, T)**2
+
+    window = partial(window_function, beam_type=beam_type,
+                     object_type=object_type, los_type=los_type,
+                     sigma_z=sigma_z, z0=z0, z1=z1, theta0=theta0)
+
+    # Set bounds based on the inputs
+    if object_type == "gaussian":
+        z_bounds = [z0 - 5 * sigma_z, z0 + 5 * sigma_z]
+    else:
+        # Only consider values within the slab
+        z_bounds = [z0, z1]
+
+    if beam_type == "gaussian":
+        # becomes highly attenuated near the beam size
+        r_bounds = [0, 1.5 * theta0 * z0]
+    else:
+        # Pencil beam
+        raise NotImplementedError("A reasonable bound for R still needs to be"
+                                  " implemented!")
+        r_bounds = [0, ]
+
+    if alphae == 3:
+
+        def integrand(z, theta, r):
+
+            # Extra factor of r on each for the cylindrical Jacobian
+            # No density contribution
+            return window(r, z) * \
+                math.exp(-0.5 * kv**2 * Dz(r, theta, V0, k0, alphav))  # -
+                       # 1.0j * kv * b * z) * r
+
+    elif alphae > 3:
+
+        norm_factor = F_eps_norm(alphae, k0)
+
+        def integrand(z, theta, r):
+
+            # Extra factor of r on each for the cylindrical Jacobian
+            # Steep density
+            return window(r, z) * \
+                C_eps(r, k0, alphae, norm_factor) * \
+                math.exp(-0.5 * kv**2 * Dz(r, theta, V0, k0, alphav))  # -
+                       # 1.0j * kv * b * z) * r
+    else:
+
+        norm_factor = F_eps_norm(alphae, k1)
+
+        def integrand(z, theta, r):
+
+            # Extra factor of r on each for the cylindrical Jacobian
+            # Shallow density
+            return window(r, z) * \
+                C_eps(r, k1, alphae, norm_factor) * \
+                math.exp(-0.5 * kv**2 * Dz(r, theta, V0, k0, alphav))  # -
+                       # 1.0j * kv * b * z) * r
+
+    if integration_type == "quad":
+        # Integration order is from the last to first arguments
+        # So z, theta, then r
+        value, error = tplquad(integrand, 0, r_bounds[1],
+                               lambda r: 0.0, lambda r: 2 * np.pi,
+                               lambda r, theta: z_bounds[0],
+                               lambda r, theta: z_bounds[1],
+                               **integration_kwargs)
+    elif integration_type == 'mc':
+        if not VEGAS_IMPORT_FLAG:
+            raise ImportError("Monte Carlo integration require the vegas "
+                              "package.")
+        integ = vegas.Integrator([z_bounds,
+                                  [0, 2 * np.pi],
+                                  r_bounds])
+
+        def wrap_integrand(vals):
+            z, r, theta = vals
+            return integrand(z, r, theta)
+
+        result = integ(wrap_integrand, **integration_kwargs)
+        value = result.mean
+        error = np.sqrt(result.var)
+
+    return P0 * fk2 * value, P0 * fk2 * error
+    # return P0 * value, P0 * error
+
+
+k_B = const.k_B.to(u.J / u.K).value
+m_p = const.m_p.to(u.kg).value
+
+
+def f_k(kv, T):
+    return np.exp(- 0.5 * k_B * T * kv**2 / m_p)
+
+
+# Window functions and related
+
+
+def window_function(r, z, beam_type, object_type, los_type, sigma_z=None,
+                    z0=None, z1=None, theta0=None):
+    '''
+    Return the window auto-correlation function for a variety of cases.
+    '''
+
+    beam_types = ['pencil', 'gaussian']
+    object_types = ['slab', 'gaussian']
+    los_types = ['parallel', 'crossing']
+
+    if beam_type not in beam_types:
+        raise ValueError("beam_type must be one of {}".format(beam_types))
+    if object_type not in object_types:
+        raise ValueError("object_type must be one of "
+                         "{}".format(object_types))
+    if los_type not in los_types:
+        raise ValueError("los_type must be one of {}".format(los_types))
+
+    if beam_type == 'pencil':
+        # los_type has no effect hasa
+        if object_type == 'slab':
+            return pencil_beam_slab_z(z, z0, z1)
+        else:
+            # Gaussian case
+            return pencil_beam_gaussian_z(z, sigma_z)
+    else:
+        if object_type == 'slab':
+            if los_type == 'parallel':
+                return gaussian_beam_slab_z_parallel(r, z, z0, z1, theta0)
+            else:
+                # Crossing
+                return gaussian_beam_slab_z_crossing(r, z, z0, z1, theta0)
+        else:
+            # Gaussian
+            if los_type == 'parallel':
+                return gaussian_beam_gaussian_z_parallel(r, z, z0, sigma_z,
+                                                         theta0)
+            else:
+                raise NotImplementedError("The gaussian beam, gaussian LOS "
+                                          "structure in the parallel limit is"
+                                          " not implemented.")

From 74133f006b11ace28def0a00d5d9ff00cf6a0ff5 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 31 Aug 2017 17:51:46 -0600
Subject: [PATCH 02/16] Add integral of F_eps needed for normalization

---
 turbustat/statistics/vca_vcs/vcs_model_helpers.pyx | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index f4729b6e..45cc26ab 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -151,6 +151,14 @@ def Int4(double r, double k1, double alphae):
     return value, err
 
 
+def Int5(double k, double k0, double alphae):
+    '''
+    Integral of F_eps^2. Needed for the normalization of C_eps
+    '''
+
+    return = 0.5 * k0**(0.5 * (1 - alphae)) * gamma(0.5 * (alphae - 1))
+
+
 # Window functions
 
 def gaussian_beam(double theta, double theta_0):

From 68b9e93c8a6b73ce2f4bde3e8410699a4c763672 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 1 Sep 2017 11:17:50 -0600
Subject: [PATCH 03/16] Actually already wrote this function in F_eps_norm

---
 turbustat/statistics/vca_vcs/vcs_model_helpers.pyx | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 45cc26ab..521c1604 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -150,15 +150,6 @@ def Int4(double r, double k1, double alphae):
 
     return value, err
 
-
-def Int5(double k, double k0, double alphae):
-    '''
-    Integral of F_eps^2. Needed for the normalization of C_eps
-    '''
-
-    return = 0.5 * k0**(0.5 * (1 - alphae)) * gamma(0.5 * (alphae - 1))
-
-
 # Window functions
 
 def gaussian_beam(double theta, double theta_0):

From fd52391a0d576d31e3f4f4afe7f66dfdc2e77c94 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 1 Sep 2017 18:08:28 -0600
Subject: [PATCH 04/16] Fixes for C_eps

---
 turbustat/statistics/vca_vcs/vcs_model_helpers.pyx | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 521c1604..2fd2bcf0 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -22,17 +22,17 @@ def C_eps(double r, double k_cut, double alphae, double norm_factor):
         # Eq. B2
         Int = Int3(r, k_cut, alphae)[0]
 
-        return 1 - 4 * pi * Int / norm_factor
+        return 1 + 4 * pi * r**(alphae - 3) * Int / norm_factor
     elif alphae == 3:
         # No density dependence.
-        return 1. / norm_factor
+        return 1.
     # Shallow
     elif alphae > 1:
         # For alphae between 1 and 3
 
         Int = Int4(r, k_cut, alphae)[0]
 
-        return 1 - 4 * pi * Int / norm_factor
+        return 1 + 4 * pi * r**(alphae - 3) * Int / norm_factor
 
     else:
         raise ValueError("Solution not defined for alphae <= 1.")

From 16836f54844a7f8a183adcd84a1c0f9fe00ec79b Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 7 Sep 2017 19:49:20 -0600
Subject: [PATCH 05/16] Several critical VCS fixes

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx  | 21 +++++---
 turbustat/statistics/vca_vcs/vcs_models.py    | 49 ++++++++++---------
 2 files changed, 41 insertions(+), 29 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 2fd2bcf0..19c6a54f 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -38,7 +38,7 @@ def C_eps(double r, double k_cut, double alphae, double norm_factor):
         raise ValueError("Solution not defined for alphae <= 1.")
 
 
-def Dz(double r, double theta, double V0, double k0, double alphav):
+def Dz(double R, double z, double V0, double k0, double alphav):
     '''
     Eq. A3 to A5.
 
@@ -47,7 +47,9 @@ def Dz(double r, double theta, double V0, double k0, double alphav):
     and can be easily included here. Do they need to be?
     '''
 
-    cdef double intone, inttwo, I_C, I_S
+    cdef double intone, inttwo, I_C, I_S, costheta, sintheta, r
+
+    r = np.sqrt(R**2 + z**2)
 
     intone = Int1(r, k0, alphav)[0]
     inttwo = Int2(r, k0, alphav)[0]
@@ -55,8 +57,11 @@ def Dz(double r, double theta, double V0, double k0, double alphav):
     I_C = (4 / 3.) * intone
     I_S = 2 * (inttwo - intone / 3.)
 
+    costheta = z / r
+    sintheta = R / r
+
     return 4 * pi * V0**2 * r**(alphav - 3) * \
-        (I_C * cos(theta)**2 + I_S * sin(theta)**2)
+        (I_C * costheta**2 + I_S * sintheta**2)
 
 
 def F_eps_norm(double alphae, double k_cut):
@@ -121,10 +126,9 @@ def Int3(double r, double k0, double alphae):
     Eq. B2 (w/o 4pi constant)
     '''
 
-    def integrand(double k):
+    def integrand(double q):
         cdef double out
-        out = k**(1 - alphae) * exp(-(k0 / k)**2) * \
-            (sin(k * r) / r)
+        out = q**(1 - alphae) * exp(-(k0 * r / q)**2) * sin(q)
         return out
 
     cdef double value, err
@@ -137,6 +141,7 @@ def Int3(double r, double k0, double alphae):
 def Int4(double r, double k1, double alphae):
     '''
     Eq. B7
+
     '''
 
     def integrand(double q):
@@ -206,6 +211,8 @@ def slab_autocorr(double z, double z_0, double z_1):
     if z >= z_0 and z <= z_1:
         out = (z_1 - z_0)**2
         return out
+    else:
+        return 0.0
 
 
 def pencil_beam_gaussian_z(double z, double sigma_z):
@@ -293,7 +300,7 @@ def gaussian_beam_gaussian_z_parallel(double R, double z, double z_0,
     cdef double w_eps_a, w_b_a
 
     # See form above. Equal to 2 * sigma_z for "theta_0"
-    w_eps_a = gaussian_beam(z, 2 * sigma_z)
+    w_eps_a = exp(-(0.25 * (z - z_0) / sigma_z)**2) / sqrt(4 * pi * sigma_z**2)
 
     w_b_a = gaussian_autocorr(R, z_0, theta0)
 
diff --git a/turbustat/statistics/vca_vcs/vcs_models.py b/turbustat/statistics/vca_vcs/vcs_models.py
index 51612db4..432b4c05 100644
--- a/turbustat/statistics/vca_vcs/vcs_models.py
+++ b/turbustat/statistics/vca_vcs/vcs_models.py
@@ -1,6 +1,6 @@
 
 import numpy as np
-from scipy.integrate import tplquad
+from scipy.integrate import dblquad
 from astropy import constants as const
 import astropy.units as u
 from functools import partial
@@ -36,7 +36,7 @@ def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
 
     # Set bounds based on the inputs
     if object_type == "gaussian":
-        z_bounds = [z0 - 5 * sigma_z, z0 + 5 * sigma_z]
+        z_bounds = [max(z0 - 5 * sigma_z, 0), z0 + 5 * sigma_z]
     else:
         # Only consider values within the slab
         z_bounds = [z0, z1]
@@ -52,58 +52,63 @@ def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
 
     if alphae == 3:
 
-        def integrand(z, theta, r):
+        def integrand(z, R):
 
             # Extra factor of r on each for the cylindrical Jacobian
             # No density contribution
-            return window(r, z) * \
-                math.exp(-0.5 * kv**2 * Dz(r, theta, V0, k0, alphav))  # -
+
+            return 2 * np.pi * R * window(R, z) * \
+                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
 
     elif alphae > 3:
 
         norm_factor = F_eps_norm(alphae, k0)
 
-        def integrand(z, theta, r):
+        def integrand(z, R):
 
             # Extra factor of r on each for the cylindrical Jacobian
             # Steep density
-            return window(r, z) * \
+
+            r = np.sqrt(R**2 + z**2)
+
+            return 2 * np.pi * R * window(R, z) * \
                 C_eps(r, k0, alphae, norm_factor) * \
-                math.exp(-0.5 * kv**2 * Dz(r, theta, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
     else:
 
         norm_factor = F_eps_norm(alphae, k1)
 
-        def integrand(z, theta, r):
+        def integrand(z, R):
 
             # Extra factor of r on each for the cylindrical Jacobian
             # Shallow density
-            return window(r, z) * \
+
+            r = np.sqrt(R**2 + z**2)
+
+            return 2 * np.pi * R * window(R, z) * \
                 C_eps(r, k1, alphae, norm_factor) * \
-                math.exp(-0.5 * kv**2 * Dz(r, theta, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
 
     if integration_type == "quad":
         # Integration order is from the last to first arguments
         # So z, theta, then r
-        value, error = tplquad(integrand, 0, r_bounds[1],
-                               lambda r: 0.0, lambda r: 2 * np.pi,
-                               lambda r, theta: z_bounds[0],
-                               lambda r, theta: z_bounds[1],
+        value, error = dblquad(integrand, 0, r_bounds[1],
+                               lambda r: z_bounds[0],
+                               lambda r: z_bounds[1],
                                **integration_kwargs)
     elif integration_type == 'mc':
         if not VEGAS_IMPORT_FLAG:
             raise ImportError("Monte Carlo integration require the vegas "
                               "package.")
         integ = vegas.Integrator([z_bounds,
-                                  [0, 2 * np.pi],
                                   r_bounds])
 
         def wrap_integrand(vals):
-            z, r, theta = vals
-            return integrand(z, r, theta)
+            z, R = vals
+            return integrand(z, R)
 
         result = integ(wrap_integrand, **integration_kwargs)
         value = result.mean
@@ -124,7 +129,7 @@ def f_k(kv, T):
 # Window functions and related
 
 
-def window_function(r, z, beam_type, object_type, los_type, sigma_z=None,
+def window_function(R, z, beam_type, object_type, los_type, sigma_z=None,
                     z0=None, z1=None, theta0=None):
     '''
     Return the window auto-correlation function for a variety of cases.
@@ -152,14 +157,14 @@ def window_function(r, z, beam_type, object_type, los_type, sigma_z=None,
     else:
         if object_type == 'slab':
             if los_type == 'parallel':
-                return gaussian_beam_slab_z_parallel(r, z, z0, z1, theta0)
+                return gaussian_beam_slab_z_parallel(R, z, z0, z1, theta0)
             else:
                 # Crossing
-                return gaussian_beam_slab_z_crossing(r, z, z0, z1, theta0)
+                return gaussian_beam_slab_z_crossing(R, z, z0, z1, theta0)
         else:
             # Gaussian
             if los_type == 'parallel':
-                return gaussian_beam_gaussian_z_parallel(r, z, z0, sigma_z,
+                return gaussian_beam_gaussian_z_parallel(R, z, z0, sigma_z,
                                                          theta0)
             else:
                 raise NotImplementedError("The gaussian beam, gaussian LOS "

From 7d0ce9ea56f351bdf87e445925a039ee88929739 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 14 Sep 2017 08:42:58 -0600
Subject: [PATCH 06/16] Various fixes for Dz, Ceps, and window functions

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx  | 51 +++++++++++++------
 1 file changed, 35 insertions(+), 16 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 19c6a54f..a40b8048 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -30,6 +30,12 @@ def C_eps(double r, double k_cut, double alphae, double norm_factor):
     elif alphae > 1:
         # For alphae between 1 and 3
 
+        # Don't use b/c there's a small jump vs. the numerical integration
+        # (But the numerical integration's error contains the actual value)
+        # Use the Inf approximation on large scales.
+        # if r >= 1 / k_cut:
+        #     Int = Int4_inf(k_cut, alphae)
+        # else:
         Int = Int4(r, k_cut, alphae)[0]
 
         return 1 + 4 * pi * r**(alphae - 3) * Int / norm_factor
@@ -64,6 +70,15 @@ def Dz(double R, double z, double V0, double k0, double alphav):
         (I_C * costheta**2 + I_S * sintheta**2)
 
 
+def Dz_simp(double R, double z, double V0, double k0, double alphav):
+
+    cdef double costheta, sintheta, r
+
+    r = np.sqrt(R**2 + z**2)
+
+    return 4 * pi * V0**2 * r**(alphav - 3) / (r**(alphav - 3) + k0**(3 - alphav))
+
+
 def F_eps_norm(double alphae, double k_cut):
     '''
     Normalization for F_eps.
@@ -93,7 +108,7 @@ def Int1(double r, double k0, double alphav):
     def integrand(double q):
         cdef double out
         out = q**(2 - alphav) * exp(-(k0 * r / q)**2) * \
-            (1 - (3 / q**2) * (sin(q) / q - cos(q)))
+            (1 - (3 / q**2) * ((sin(q) / q) - cos(q)))
         return out
 
     cdef double value, err
@@ -111,7 +126,7 @@ def Int2(double r, double k0, double alphav):
     def integrand(double q):
         cdef double out
         out = q**(2 - alphav) * exp(-(k0 * r / q)**2) * \
-            (1 - sin(q) / q)
+            (1 - (sin(q) / q))
         return out
 
     cdef double value, err
@@ -155,6 +170,18 @@ def Int4(double r, double k1, double alphae):
 
     return value, err
 
+def Int4_inf(double k1, double alphae):
+    '''
+    Analytical form for Eq. B9
+
+    Use for r >> 1 / k1.
+
+    Valid for alphae > 1, alphae < 3.
+    '''
+
+    return gamma(2 - alphae) * sin(alphae * pi * 0.5)
+
+
 # Window functions
 
 def gaussian_beam(double theta, double theta_0):
@@ -175,7 +202,7 @@ def gaussian_autocorr(double R, double z_0, double theta_0):
     Eq. 32.
 
     For W_b:
-    :math:`-\frac{e^{-\frac{R^2}{2 \theta ^2 z^2}} \left(\sqrt{\pi } R e^{\frac{R^2}{4 \theta ^2 z^2}} \text{erfc}\left(\frac{R}{2 \theta  z}\right)-2 \theta  z\right)}{4 \theta  z^3}`
+    :math:`\frac{1}{2 \pi z_0^2 \theta_0^2} \exp^{-R^2 / 2 \theta_0^2 z_0^2}`
 
     w_b was solved for using Mathematica. My by-hand solution had the same
     asymptotic behaviour, but looked so much worse than the compacter form
@@ -183,17 +210,12 @@ def gaussian_autocorr(double R, double z_0, double theta_0):
 
     '''
 
-    cdef double ratio_term, exp_term
-
-    ratio_term = R / (theta_0 * z_0)
-    exp_term = ratio_term**2
+    cdef double ratio_term
 
-    term1 = - 1. / (4 * theta_0 * z_0**3)
-    term2 = sqrt(pi) * R * exp(- 0.25 * exp_term) * \
-        erfc(0.5 * ratio_term)
-    term3 = 2 * theta_0 * z_0 * exp(- 0.5 * exp_term)
+    ratio_term = R / (2 * theta_0 * z_0)
 
-    return term1 * (term2 - term3)
+    return 0.5 * np.exp(- ratio_term**2) / (z_0 * theta_0)**2
+    # return 0.5 * np.exp(-0.25 * ratio_term**2)
 
 
 def slab_autocorr(double z, double z_0, double z_1):
@@ -292,15 +314,12 @@ def gaussian_beam_gaussian_z_parallel(double R, double z, double z_0,
     definition. The autocorrelation function is for the object thickness along
     the LOS, and so should be independent of distance.
 
-    For W_b:
-    :math:`-\frac{e^{-\frac{R^2}{2 \theta ^2 z^2}} \left(\sqrt{\pi } R e^{\frac{R^2}{4 \theta ^2 z^2}} \text{erfc}\left(\frac{R}{2 \theta  z}\right)-2 \theta  z\right)}{4 \theta  z^3}`
-
     '''
 
     cdef double w_eps_a, w_b_a
 
     # See form above. Equal to 2 * sigma_z for "theta_0"
-    w_eps_a = exp(-(0.25 * (z - z_0) / sigma_z)**2) / sqrt(4 * pi * sigma_z**2)
+    w_eps_a = gaussian_beam(z, 2 * sigma_z)
 
     w_b_a = gaussian_autocorr(R, z_0, theta0)
 

From 33b8431d4b1e019449485cc042434ef43fc453a0 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 14 Sep 2017 11:44:39 -0600
Subject: [PATCH 07/16] Use analytic solution for integral in C_eps for shallow
 spectra

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx  | 44 +++++++++++--------
 1 file changed, 25 insertions(+), 19 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index a40b8048..c7ac7b93 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -1,7 +1,7 @@
 
 import numpy as np
 from scipy.integrate import quad
-from scipy.special import gamma, erf, erfc
+from scipy.special import gamma, erf, erfc, hyp1f1
 cimport cython
 
 from libc.math cimport sqrt, exp, sin, cos, atan
@@ -13,6 +13,14 @@ def C_eps(double r, double k_cut, double alphae, double norm_factor):
     '''
 
     k_cut is k0 for alphae > 3 and k1 for alphae<3
+
+
+    When alphae<3, this uses the analytic solution to Eq. B7,
+
+    :math:`\int_0^\inf q^{1 - \alpha_e} \sin(q) \exp[-q^2/(k_1 r)^2]`
+
+    :math:`0.5 (k_1 r)^{3 - \alpha_e} \Gamma(0.5 * (3 - alpha_e) {\rm F1}(0.5 (3 - alphae), 3/2, -0.25 (k_1 r)^2))`
+
     '''
 
     cdef double Int
@@ -28,18 +36,16 @@ def C_eps(double r, double k_cut, double alphae, double norm_factor):
         return 1.
     # Shallow
     elif alphae > 1:
-        # For alphae between 1 and 3
+        # For alphae between 1 and 3. See full solution above and in B7.
+        # The form is simplified here for efficiency
 
-        # Don't use b/c there's a small jump vs. the numerical integration
-        # (But the numerical integration's error contains the actual value)
-        # Use the Inf approximation on large scales.
-        # if r >= 1 / k_cut:
-        #     Int = Int4_inf(k_cut, alphae)
-        # else:
-        Int = Int4(r, k_cut, alphae)[0]
+        cdef double kr, nu
 
-        return 1 + 4 * pi * r**(alphae - 3) * Int / norm_factor
+        nu = 3 - alphae
+        kr = k1 * r
 
+        return 1 + 2 * pi * k1**nu * gamma(0.5 * nu) * \
+            hyp1f1(0.5 * nu, 1.5, -0.25 * kr**2) / norm_factor
     else:
         raise ValueError("Solution not defined for alphae <= 1.")
 
@@ -155,20 +161,21 @@ def Int3(double r, double k0, double alphae):
 
 def Int4(double r, double k1, double alphae):
     '''
-    Eq. B7
+    Analytic solution to Eq. B7,
+
+    :math:`\int_0^\inf q^{1 - \alpha_e} \sin(q) \exp[-q^2/(k_1 r)^2]`
+
+    :math:`0.5 (k_1 r)^{3 - \alpha_e} \Gamma(0.5 * (3 - alpha_e) {\rm F1}(0.5 (3 - alphae), 3/2, -0.25 (k_1 r)^2))`
 
     '''
 
-    def integrand(double q):
-        cdef double out
-        out = q**(1 - alphae) * exp(-(q / (k1 * r))**2) * sin(q)
-        return out
+    cdef double nu, kr
 
-    cdef double value, err
+    nu = 3 - alphae
+    kr = k1 * r
 
-    value, err = quad(integrand, 0, np.inf)
+    return 0.5 * kr**nu * gamma(0.5 * nu) * hyp1f1(0.5 * nu, 1.5, -0.25 * kr**2)
 
-    return value, err
 
 def Int4_inf(double k1, double alphae):
     '''
@@ -181,7 +188,6 @@ def Int4_inf(double k1, double alphae):
 
     return gamma(2 - alphae) * sin(alphae * pi * 0.5)
 
-
 # Window functions
 
 def gaussian_beam(double theta, double theta_0):

From df36e2f46a3c9f9dd2ffc4669dfcd32c1bd458be Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 14 Sep 2017 14:08:39 -0600
Subject: [PATCH 08/16] Add analytic solution for steep emissivity spectrum

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx  | 32 +++++++++++++++----
 1 file changed, 26 insertions(+), 6 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index c7ac7b93..4959f346 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -2,6 +2,8 @@
 import numpy as np
 from scipy.integrate import quad
 from scipy.special import gamma, erf, erfc, hyp1f1
+from mpmath import hyper
+
 cimport cython
 
 from libc.math cimport sqrt, exp, sin, cos, atan
@@ -23,12 +25,13 @@ def C_eps(double r, double k_cut, double alphae, double norm_factor):
 
     '''
 
-    cdef double Int
+    cdef double Int, kr, nu
 
     # Steep
     if alphae > 3:
         # Eq. B2
-        Int = Int3(r, k_cut, alphae)[0]
+        # Int = Int3(r, k_cut, alphae)[0]
+        Int = Int3_analytic(r, k_cut, alphae)
 
         return 1 + 4 * pi * r**(alphae - 3) * Int / norm_factor
     elif alphae == 3:
@@ -39,12 +42,10 @@ def C_eps(double r, double k_cut, double alphae, double norm_factor):
         # For alphae between 1 and 3. See full solution above and in B7.
         # The form is simplified here for efficiency
 
-        cdef double kr, nu
-
         nu = 3 - alphae
-        kr = k1 * r
+        kr = k_cut * r
 
-        return 1 + 2 * pi * k1**nu * gamma(0.5 * nu) * \
+        return 1 + 2 * pi * k_cut**nu * gamma(0.5 * nu) * \
             hyp1f1(0.5 * nu, 1.5, -0.25 * kr**2) / norm_factor
     else:
         raise ValueError("Solution not defined for alphae <= 1.")
@@ -159,6 +160,23 @@ def Int3(double r, double k0, double alphae):
     return value, err
 
 
+def Int3_analytic(double r, double k0, double alphae):
+    '''
+
+    Analytic solution for Eq. B2
+
+    '''
+
+    cdef double term1, term2, nu, kr2
+
+    nu = 3 - alphae
+    kr2 = (k0 * r)**2
+
+    term1 = 0.5 * (k0 * r)**nu * gamma(- 0.5 * nu) * float(hyper([], [1.5, 2.5 - 0.5 * alphae], 0.25 * kr2))
+    term2 = gamma(2 - alphae) * float(hyper([], [0.5 * (- 1 + alphae), 0.5 * alphae], 0.25 * kr2)) * sin(0.5 * alphae * pi)
+
+    return term1 + term2
+
 def Int4(double r, double k1, double alphae):
     '''
     Analytic solution to Eq. B7,
@@ -167,6 +185,8 @@ def Int4(double r, double k1, double alphae):
 
     :math:`0.5 (k_1 r)^{3 - \alpha_e} \Gamma(0.5 * (3 - alpha_e) {\rm F1}(0.5 (3 - alphae), 3/2, -0.25 (k_1 r)^2))`
 
+    Verified against mathematica numerical integration.
+
     '''
 
     cdef double nu, kr

From dbb5ed5b6e71beb12e35ab2f48b8cfc75008d901 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 14 Sep 2017 22:16:36 -0600
Subject: [PATCH 09/16] Add analytical version of Int1 and Int2; use flint for
 a fast general hypergeometric function (need to make optional!)

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx  | 62 ++++++++++++++++++-
 1 file changed, 59 insertions(+), 3 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 4959f346..79f83c38 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -2,7 +2,18 @@
 import numpy as np
 from scipy.integrate import quad
 from scipy.special import gamma, erf, erfc, hyp1f1
-from mpmath import hyper
+
+
+import flint
+hyper = flint.arb.hypgeom
+
+# ImportError = Exception
+
+# try:
+#     import flint
+#     hyper = flint.arb.hypergeom
+# except ImportError:
+#     from mpmath import hyper
 
 cimport cython
 
@@ -64,8 +75,11 @@ def Dz(double R, double z, double V0, double k0, double alphav):
 
     r = np.sqrt(R**2 + z**2)
 
-    intone = Int1(r, k0, alphav)[0]
-    inttwo = Int2(r, k0, alphav)[0]
+    # intone = Int1(r, k0, alphav)[0]
+    # inttwo = Int2(r, k0, alphav)[0]
+
+    intone = Int1_analytic(r, k0, alphav)
+    inttwo = Int2_analytic(r, k0, alphav)
 
     I_C = (4 / 3.) * intone
     I_S = 2 * (inttwo - intone / 3.)
@@ -125,6 +139,28 @@ def Int1(double r, double k0, double alphav):
     return value, err
 
 
+def Int1_analytic(double r, double k0, double alphav):
+    '''
+    Analytical solution to Eq. A4
+    '''
+
+    cdef double kr, kr2, nu, term1, term2, term3, term4, terma, termb
+
+    kr2 = (k0 * r)**2
+    kr = k0 * r
+    nu = alphav - 3
+
+    term1 = 3 * nu * float(hyper([], [0.5, 1.5 - 0.5 * alphav], 0.25 * kr2))
+    term2 = 3 * nu * float(hyper([], [1.5, 1.5 - 0.5 * alphav], 0.25 * kr2))
+    terma = 0.25 * kr**(1 - alphav) * gamma(0.5 * nu) * (2 * kr2 + term1 - term2)
+
+    term3 = float(hyper([], [0.5 * (1 + alphav), 1 + 0.5 * alphav], 0.25 * kr2))
+    term4 = alphav * float(hyper([], [0.5 * (1 + alphav), 0.5 * alphav], 0.25 * kr2))
+    termb = 3 * gamma(-alphav) * sin(alphav * pi * 0.5) * (term3 - term4)
+
+    return terma + termb
+
+
 def Int2(double r, double k0, double alphav):
     '''
     First integral in Eq. A5
@@ -143,6 +179,26 @@ def Int2(double r, double k0, double alphav):
     return value, err
 
 
+def Int2_analytic(double r, double k0, double alphav):
+    '''
+    Analytical solution to first integral in Eq. A5
+    '''
+
+    cdef double kr, kr3, kr2, term1, term2, term3, nu
+
+    kr3 = (k0 * r)**3
+    kr2 = (k0 * r)**2
+    kr = k0 * r
+    nu = alphav - 3
+
+
+    term1 = - kr3 * gamma(0.5 * nu)
+    term2 = kr3 * gamma(0.5 * nu) * float(hyper([], [1.5, 2.5 - 0.5 * alphav], 0.25 * kr2))
+    term3 = 2 * kr**alphav * gamma(2 - alphav) * float(hyper([], [0.5 * (alphav - 1), 0.5 * alphav], 0.25 * kr2)) * sin(alphav * pi * 0.5)
+
+    return -0.5 * kr**-alphav * (term1 + term2 + term3)
+
+
 def Int3(double r, double k0, double alphae):
     '''
     Eq. B2 (w/o 4pi constant)

From b58d023ae63347cb7152b005815d9e4b9cb05c23 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Thu, 14 Sep 2017 22:18:08 -0600
Subject: [PATCH 10/16] Various updates and changes to the VCS main model
 functions

---
 turbustat/statistics/vca_vcs/vcs_models.py | 32 +++++++++++++++-------
 1 file changed, 22 insertions(+), 10 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_models.py b/turbustat/statistics/vca_vcs/vcs_models.py
index 432b4c05..3c8ed2c7 100644
--- a/turbustat/statistics/vca_vcs/vcs_models.py
+++ b/turbustat/statistics/vca_vcs/vcs_models.py
@@ -17,7 +17,8 @@
                                pencil_beam_gaussian_z,
                                gaussian_beam_slab_z_parallel,
                                gaussian_beam_slab_z_crossing,
-                               gaussian_beam_gaussian_z_parallel)
+                               gaussian_beam_gaussian_z_parallel,
+                               Dz_simp)
 
 
 def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
@@ -28,7 +29,10 @@ def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
     VCS model.
     '''
 
-    fk2 = f_k(kv, T)**2
+    if T > 0:
+        fk2 = f_k(kv, T)**2
+    else:
+        fk2 = 1.
 
     window = partial(window_function, beam_type=beam_type,
                      object_type=object_type, los_type=los_type,
@@ -36,19 +40,24 @@ def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
 
     # Set bounds based on the inputs
     if object_type == "gaussian":
-        z_bounds = [max(z0 - 5 * sigma_z, 0), z0 + 5 * sigma_z]
+        z_bounds = [0, 1 / k0]
     else:
         # Only consider values within the slab
         z_bounds = [z0, z1]
 
-    if beam_type == "gaussian":
+    if beam_type == "gaussian" or beam_type == 'none':
         # becomes highly attenuated near the beam size
-        r_bounds = [0, 1.5 * theta0 * z0]
+        R_bounds = [0, 1.5 * theta0 * z0]
+
+    elif beam_type == 'none':
+        R_bounds = [0, 1 / k0]
     else:
         # Pencil beam
+        # TODO: In this case, you only need to integrate over z. R's contribution
+        # is a delta-fcn
         raise NotImplementedError("A reasonable bound for R still needs to be"
                                   " implemented!")
-        r_bounds = [0, ]
+        R_bounds = [0, ]
 
     if alphae == 3:
 
@@ -95,7 +104,7 @@ def integrand(z, R):
     if integration_type == "quad":
         # Integration order is from the last to first arguments
         # So z, theta, then r
-        value, error = dblquad(integrand, 0, r_bounds[1],
+        value, error = dblquad(integrand, 0, R_bounds[1],
                                lambda r: z_bounds[0],
                                lambda r: z_bounds[1],
                                **integration_kwargs)
@@ -104,7 +113,7 @@ def integrand(z, R):
             raise ImportError("Monte Carlo integration require the vegas "
                               "package.")
         integ = vegas.Integrator([z_bounds,
-                                  r_bounds])
+                                  R_bounds])
 
         def wrap_integrand(vals):
             z, R = vals
@@ -115,7 +124,6 @@ def wrap_integrand(vals):
         error = np.sqrt(result.var)
 
     return P0 * fk2 * value, P0 * fk2 * error
-    # return P0 * value, P0 * error
 
 
 k_B = const.k_B.to(u.J / u.K).value
@@ -127,6 +135,8 @@ def f_k(kv, T):
 
 
 # Window functions and related
+def no_window(*args, **kwargs):
+    return 1.
 
 
 def window_function(R, z, beam_type, object_type, los_type, sigma_z=None,
@@ -135,7 +145,7 @@ def window_function(R, z, beam_type, object_type, los_type, sigma_z=None,
     Return the window auto-correlation function for a variety of cases.
     '''
 
-    beam_types = ['pencil', 'gaussian']
+    beam_types = ['pencil', 'gaussian', 'none']
     object_types = ['slab', 'gaussian']
     los_types = ['parallel', 'crossing']
 
@@ -147,6 +157,8 @@ def window_function(R, z, beam_type, object_type, los_type, sigma_z=None,
     if los_type not in los_types:
         raise ValueError("los_type must be one of {}".format(los_types))
 
+    if beam_type == 'none':
+        return no_window(R, z)
     if beam_type == 'pencil':
         # los_type has no effect hasa
         if object_type == 'slab':

From e43955aef32054ca5246a1f2c40a2805c135e6c4 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 15 Sep 2017 17:00:17 -0600
Subject: [PATCH 11/16] Simplify Dz a bit

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx  | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 79f83c38..27d7d59f 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -71,7 +71,7 @@ def Dz(double R, double z, double V0, double k0, double alphav):
     and can be easily included here. Do they need to be?
     '''
 
-    cdef double intone, inttwo, I_C, I_S, costheta, sintheta, r
+    cdef double intone, inttwo, I_C, I_S, r
 
     r = np.sqrt(R**2 + z**2)
 
@@ -84,11 +84,13 @@ def Dz(double R, double z, double V0, double k0, double alphav):
     I_C = (4 / 3.) * intone
     I_S = 2 * (inttwo - intone / 3.)
 
-    costheta = z / r
-    sintheta = R / r
-
-    return 4 * pi * V0**2 * r**(alphav - 3) * \
-        (I_C * costheta**2 + I_S * sintheta**2)
+    # Equation looks like:
+    # 4 * pi * V0**2 * r**(alphav - 3) * \
+    #     (I_C * costheta**2 + I_S * sintheta**2)
+    # Since costheta = z /r and sintheta = R/r, can
+    # pull 1/r^2 out of brackets to get:
+    return 4 * pi * V0**2 * r**(alphav - 5) * \
+        (I_C * z**2 + I_S * R**2)
 
 
 def Dz_simp(double R, double z, double V0, double k0, double alphav):
@@ -286,10 +288,6 @@ def gaussian_autocorr(double R, double z_0, double theta_0):
     For W_b:
     :math:`\frac{1}{2 \pi z_0^2 \theta_0^2} \exp^{-R^2 / 2 \theta_0^2 z_0^2}`
 
-    w_b was solved for using Mathematica. My by-hand solution had the same
-    asymptotic behaviour, but looked so much worse than the compacter form
-    above.
-
     '''
 
     cdef double ratio_term
@@ -297,7 +295,6 @@ def gaussian_autocorr(double R, double z_0, double theta_0):
     ratio_term = R / (2 * theta_0 * z_0)
 
     return 0.5 * np.exp(- ratio_term**2) / (z_0 * theta_0)**2
-    # return 0.5 * np.exp(-0.25 * ratio_term**2)
 
 
 def slab_autocorr(double z, double z_0, double z_1):

From 98d46ca0e69c8b25d87e869d6c79c15f292a0dde Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 15 Sep 2017 20:18:49 -0600
Subject: [PATCH 12/16] Use hyp1f1 from flint; another Dz function to test

---
 .../statistics/vca_vcs/vcs_model_helpers.pyx     | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 27d7d59f..79859d00 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -1,11 +1,12 @@
 
 import numpy as np
 from scipy.integrate import quad
-from scipy.special import gamma, erf, erfc, hyp1f1
+from scipy.special import gamma, erf, erfc  # , hyp1f1
 
 
 import flint
 hyper = flint.arb.hypgeom
+hyp1f1 = flint.arb.hypgeom_1f1
 
 # ImportError = Exception
 
@@ -95,11 +96,20 @@ def Dz(double R, double z, double V0, double k0, double alphav):
 
 def Dz_simp(double R, double z, double V0, double k0, double alphav):
 
-    cdef double costheta, sintheta, r
+    cdef double r
 
     r = np.sqrt(R**2 + z**2)
 
-    return 4 * pi * V0**2 * r**(alphav - 3) / (r**(alphav - 3) + k0**(3 - alphav))
+    return V0**2 * r**(alphav - 3) / (r**(alphav - 3) + k0**(3 - alphav))
+
+
+def Dz_simp2(double R, double z, double V0, double k0, double alphav):
+
+    cdef double m
+
+    m = alphav - 3
+
+    return V0**2  * (R**2 + z**2)**(0.5 * m - 1) * ((1 + 0.5 * m) * R**2 + z**2)
 
 
 def F_eps_norm(double alphae, double k_cut):

From afa323914e7afa496ec5e892d0709013d1538f33 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 15 Sep 2017 20:21:41 -0600
Subject: [PATCH 13/16] Missing a factor of pi in gaussian autocorrelation
 function

---
 turbustat/statistics/vca_vcs/vcs_model_helpers.pyx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
index 79859d00..af22e2de 100644
--- a/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
+++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx
@@ -304,7 +304,7 @@ def gaussian_autocorr(double R, double z_0, double theta_0):
 
     ratio_term = R / (2 * theta_0 * z_0)
 
-    return 0.5 * np.exp(- ratio_term**2) / (z_0 * theta_0)**2
+    return 0.5 * np.exp(- ratio_term**2) / (pi * (z_0 * theta_0)**2)
 
 
 def slab_autocorr(double z, double z_0, double z_1):

From 94dffc5311949536b2b6298e36c8f746013288de Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Fri, 15 Sep 2017 20:22:31 -0600
Subject: [PATCH 14/16] Use Dz_simp for now (not completely sure it's right,
 though); normalization for the no window function case

---
 turbustat/statistics/vca_vcs/vcs_models.py | 26 +++++++++++++---------
 1 file changed, 16 insertions(+), 10 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_models.py b/turbustat/statistics/vca_vcs/vcs_models.py
index 3c8ed2c7..9f45f02d 100644
--- a/turbustat/statistics/vca_vcs/vcs_models.py
+++ b/turbustat/statistics/vca_vcs/vcs_models.py
@@ -18,7 +18,7 @@
                                gaussian_beam_slab_z_parallel,
                                gaussian_beam_slab_z_crossing,
                                gaussian_beam_gaussian_z_parallel,
-                               Dz_simp)
+                               Dz_simp, Dz_simp2)
 
 
 def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
@@ -40,14 +40,14 @@ def P1(kv, alphav, alphae, P0, k0, V0, T, b=0,
 
     # Set bounds based on the inputs
     if object_type == "gaussian":
-        z_bounds = [0, 1 / k0]
+        z_bounds = [0, max(2 / float(k0), 3 * sigma_z)]
     else:
         # Only consider values within the slab
         z_bounds = [z0, z1]
 
-    if beam_type == "gaussian" or beam_type == 'none':
+    if beam_type == "gaussian":
         # becomes highly attenuated near the beam size
-        R_bounds = [0, 1.5 * theta0 * z0]
+        R_bounds = [0, max(5. * theta0 * z0, 1 / k0)]
 
     elif beam_type == 'none':
         R_bounds = [0, 1 / k0]
@@ -67,7 +67,10 @@ def integrand(z, R):
             # No density contribution
 
             return 2 * np.pi * R * window(R, z) * \
-                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz_simp(R, z, V0, k0, alphav))  # -
+
+            # return 2 * np.pi * R * window(R, z) * \
+            #     math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
 
     elif alphae > 3:
@@ -83,7 +86,7 @@ def integrand(z, R):
 
             return 2 * np.pi * R * window(R, z) * \
                 C_eps(r, k0, alphae, norm_factor) * \
-                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz_simp(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
     else:
 
@@ -98,7 +101,7 @@ def integrand(z, R):
 
             return 2 * np.pi * R * window(R, z) * \
                 C_eps(r, k1, alphae, norm_factor) * \
-                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz_simp(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
 
     if integration_type == "quad":
@@ -135,8 +138,11 @@ def f_k(kv, T):
 
 
 # Window functions and related
-def no_window(*args, **kwargs):
-    return 1.
+def no_window(R, z, sigma_z, z0, theta0):
+    '''
+    Only normalized for gaussian windows in R and z!
+    '''
+    return 1. / (4 * np.sqrt(np.pi) * sigma_z * z0**2 * theta0**2)
 
 
 def window_function(R, z, beam_type, object_type, los_type, sigma_z=None,
@@ -158,7 +164,7 @@ def window_function(R, z, beam_type, object_type, los_type, sigma_z=None,
         raise ValueError("los_type must be one of {}".format(los_types))
 
     if beam_type == 'none':
-        return no_window(R, z)
+        return no_window(R, z, sigma_z, z0, theta0)
     if beam_type == 'pencil':
         # los_type has no effect hasa
         if object_type == 'slab':

From c45060f476c723ba87b91f35b58dd2df4673482e Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Mon, 18 Sep 2017 15:51:06 -0600
Subject: [PATCH 15/16] Add in a mu parameter for the thermal line width

---
 turbustat/statistics/vca_vcs/vcs_models.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/turbustat/statistics/vca_vcs/vcs_models.py b/turbustat/statistics/vca_vcs/vcs_models.py
index 9f45f02d..9ffa7577 100644
--- a/turbustat/statistics/vca_vcs/vcs_models.py
+++ b/turbustat/statistics/vca_vcs/vcs_models.py
@@ -67,7 +67,7 @@ def integrand(z, R):
             # No density contribution
 
             return 2 * np.pi * R * window(R, z) * \
-                math.exp(-0.5 * kv**2 * Dz_simp(R, z, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
 
             # return 2 * np.pi * R * window(R, z) * \
             #     math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
@@ -86,7 +86,7 @@ def integrand(z, R):
 
             return 2 * np.pi * R * window(R, z) * \
                 C_eps(r, k0, alphae, norm_factor) * \
-                math.exp(-0.5 * kv**2 * Dz_simp(R, z, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
     else:
 
@@ -101,7 +101,7 @@ def integrand(z, R):
 
             return 2 * np.pi * R * window(R, z) * \
                 C_eps(r, k1, alphae, norm_factor) * \
-                math.exp(-0.5 * kv**2 * Dz_simp(R, z, V0, k0, alphav))  # -
+                math.exp(-0.5 * kv**2 * Dz(R, z, V0, k0, alphav))  # -
                        # 1.0j * kv * b * z) * r
 
     if integration_type == "quad":
@@ -133,8 +133,8 @@ def wrap_integrand(vals):
 m_p = const.m_p.to(u.kg).value
 
 
-def f_k(kv, T):
-    return np.exp(- 0.5 * k_B * T * kv**2 / m_p)
+def f_k(kv, T, mu=1.3):
+    return np.exp(- 0.5 * k_B * T * kv**2 / (mu * m_p))
 
 
 # Window functions and related

From ef4f3637f0deac61ae4c7ad2de2fd0f6f6c8b738 Mon Sep 17 00:00:00 2001
From: e-koch <koch.eric.w@gmail.com>
Date: Wed, 20 Sep 2017 11:27:50 -0600
Subject: [PATCH 16/16] Try imposing bounds on theta in the Gaussian fit

---
 turbustat/statistics/pca/width_estimate.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/turbustat/statistics/pca/width_estimate.py b/turbustat/statistics/pca/width_estimate.py
index f1d7dd14..6b0fe48e 100644
--- a/turbustat/statistics/pca/width_estimate.py
+++ b/turbustat/statistics/pca/width_estimate.py
@@ -400,6 +400,9 @@ def fit_2D_gaussian(xmat, ymat, z):
                                          'y_mean': True}) + \
         astropy_models.Const2D(amplitude=[np.percentile(z, 10)])
 
+    # Impose restrictions on theta
+    g.theta_0.bounds = (-np.pi, np.pi)
+
     fit_g = fitting.LevMarLSQFitter()
     output = fit_g(g, xmat, ymat, z)
     cov = fit_g.fit_info['param_cov']