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'] 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..af22e2de --- /dev/null +++ b/turbustat/statistics/vca_vcs/vcs_model_helpers.pyx @@ -0,0 +1,415 @@ + +import numpy as np +from scipy.integrate import quad +from scipy.special import gamma, erf, erfc # , hyp1f1 + + +import flint +hyper = flint.arb.hypgeom +hyp1f1 = flint.arb.hypgeom_1f1 + +# ImportError = Exception + +# try: +# import flint +# hyper = flint.arb.hypergeom +# except ImportError: +# from mpmath import hyper + +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 + + + 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, kr, nu + + # Steep + if alphae > 3: + # Eq. B2 + # 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: + # No density dependence. + return 1. + # Shallow + elif alphae > 1: + # For alphae between 1 and 3. See full solution above and in B7. + # The form is simplified here for efficiency + + nu = 3 - alphae + kr = k_cut * r + + 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.") + + +def Dz(double R, double z, 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, r + + r = np.sqrt(R**2 + z**2) + + # 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.) + + # 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): + + cdef double r + + r = np.sqrt(R**2 + z**2) + + 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): + ''' + 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 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 + ''' + + 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 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) + ''' + + def integrand(double q): + cdef double out + out = q**(1 - alphae) * exp(-(k0 * r / q)**2) * sin(q) + return out + + cdef double value, err + + value, err = quad(integrand, 0, np.inf) + + 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, + + :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))` + + Verified against mathematica numerical integration. + + ''' + + cdef double nu, kr + + nu = 3 - alphae + kr = k1 * r + + return 0.5 * kr**nu * gamma(0.5 * nu) * hyp1f1(0.5 * nu, 1.5, -0.25 * kr**2) + + +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): + ''' + 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{1}{2 \pi z_0^2 \theta_0^2} \exp^{-R^2 / 2 \theta_0^2 z_0^2}` + + ''' + + cdef double ratio_term + + ratio_term = R / (2 * theta_0 * z_0) + + 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): + ''' + 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 + else: + return 0.0 + + +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. + + ''' + + 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..9ffa7577 --- /dev/null +++ b/turbustat/statistics/vca_vcs/vcs_models.py @@ -0,0 +1,190 @@ + +import numpy as np +from scipy.integrate import dblquad +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, + Dz_simp, Dz_simp2) + + +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. + ''' + + 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, + sigma_z=sigma_z, z0=z0, z1=z1, theta0=theta0) + + # Set bounds based on the inputs + if object_type == "gaussian": + 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": + # becomes highly attenuated near the beam size + R_bounds = [0, max(5. * theta0 * z0, 1 / k0)] + + 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, ] + + if alphae == 3: + + def integrand(z, R): + + # Extra factor of r on each for the cylindrical Jacobian + # No density contribution + + return 2 * np.pi * R * window(R, z) * \ + 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)) # - + # 1.0j * kv * b * z) * r + + elif alphae > 3: + + norm_factor = F_eps_norm(alphae, k0) + + def integrand(z, R): + + # Extra factor of r on each for the cylindrical Jacobian + # Steep density + + 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, z, V0, k0, alphav)) # - + # 1.0j * kv * b * z) * r + else: + + norm_factor = F_eps_norm(alphae, k1) + + def integrand(z, R): + + # Extra factor of r on each for the cylindrical Jacobian + # Shallow density + + 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, 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 = 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, + R_bounds]) + + def wrap_integrand(vals): + z, R = vals + return integrand(z, R) + + result = integ(wrap_integrand, **integration_kwargs) + value = result.mean + error = np.sqrt(result.var) + + return P0 * fk2 * value, P0 * fk2 * 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, mu=1.3): + return np.exp(- 0.5 * k_B * T * kv**2 / (mu * m_p)) + + +# Window functions and related +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, + z0=None, z1=None, theta0=None): + ''' + Return the window auto-correlation function for a variety of cases. + ''' + + beam_types = ['pencil', 'gaussian', 'none'] + 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 == 'none': + return no_window(R, z, sigma_z, z0, theta0) + 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.")