Merge pull request #50 from AlecThomson/pure

Pythonize gaussft.f to gaussft.py

Author: AlecThomson

pyproject.toml

@@ -1,6 +1,36 @@
+[tool.poetry]
+name = "racs-tools"
+version = "3.0.0"
+description = "Useful scripts for RACS."
+authors = ["Alec Thomson"]
+license = "BSD"
+readme = "README.md"
+packages = [{include = "racs_tools"}]
+
+[tool.poetry.dependencies]
+python = "^3.8"
+numpy = "<2"
+astropy = "^5"
+radio_beam = "*"
+schwimmbad = "*"
+scipy = "*"
+spectral_cube = "^0.6.3"
+tqdm = "*"
+numba = "*"
+mpi4py = {version = "*", optional = true}
+
+[tool.poetry.dev-dependencies]
+black = "*"
+isort = "*"
+
+[tool.poetry.extras]
+mpi = ["mpi4py"]
+
+[tool.poetry.scripts]
+beamcon_2D = "racs_tools.beamcon_2D:cli"
+beamcon_3D = "racs_tools.beamcon_3D:cli"
+getnoise_list = "racs_tools.getnoise_list:cli"
+
 [build-system]
-requires = [
-    "setuptools>=42",
-    "wheel",
-    "numpy"]
-build-backend = "setuptools.build_meta"
+requires = ["poetry-core"]
+build-backend = "poetry.core.masonry.api"

racs_tools/au2.py

@@ -1,12 +1,11 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 """ For getting fluxes right in Jy/beam """
 __author__ = "Tessa Vernstrom"
 
 import math
 from typing import List, Tuple
 
 import numpy as np
-from scipy import *
 
 
 def gaussianDeconvolve(smaj, smin, spa, bmaj, bmin, bpa):
@@ -26,9 +25,6 @@ def gaussianDeconvolve(smaj, smin, spa, bmaj, bmin, bpa):
     0.5. Feel a little wary about that first change.
     """
 
-    import numpy as np
-    from numpy import abs, arctan2, cos, min, sin, sqrt
-
     spa = np.radians(spa)
     bpa = np.radians(bpa)
     if smaj < bmaj:
@@ -37,29 +33,29 @@ def gaussianDeconvolve(smaj, smin, spa, bmaj, bmin, bpa):
         smin = bmin
 
     alpha = (
-        (smaj * cos(spa)) ** 2
-        + (smin * sin(spa)) ** 2
-        - (bmaj * cos(bpa)) ** 2
-        - (bmin * sin(bpa)) ** 2
+        (smaj * np.cos(spa)) ** 2
+        + (smin * np.sin(spa)) ** 2
+        - (bmaj * np.cos(bpa)) ** 2
+        - (bmin * np.sin(bpa)) ** 2
     )
     beta = (
-        (smaj * sin(spa)) ** 2
-        + (smin * cos(spa)) ** 2
-        - (bmaj * sin(bpa)) ** 2
-        - (bmin * cos(bpa)) ** 2
+        (smaj * np.sin(spa)) ** 2
+        + (smin * np.cos(spa)) ** 2
+        - (bmaj * np.sin(bpa)) ** 2
+        - (bmin * np.cos(bpa)) ** 2
     )
     gamma = 2 * (
-        (smin**2 - smaj**2) * sin(spa) * cos(spa)
-        - (bmin**2 - bmaj**2) * sin(bpa) * cos(bpa)
+        (smin**2 - smaj**2) * np.sin(spa) * np.cos(spa)
+        - (bmin**2 - bmaj**2) * np.sin(bpa) * np.cos(bpa)
     )
     #    print smaj,smin
     #    print alpha,beta,gamma
     s = alpha + beta
-    t = sqrt((alpha - beta) ** 2 + gamma**2)
+    t = np.sqrt((alpha - beta) ** 2 + gamma**2)
     #    print s,t
-    dmaj = sqrt(0.5 * (s + t))
+    dmaj = np.sqrt(0.5 * (s + t))
     if s > t:
-        dmin = sqrt(0.5 * (s - t))
+        dmin = np.sqrt(0.5 * (s - t))
     else:
         dmin = 0
     #    print dmaj,dmin
@@ -71,7 +67,7 @@ def gaussianDeconvolve(smaj, smin, spa, bmaj, bmin, bpa):
     if abs(gamma) + abs(alpha - beta) == 0:
         dpa = 0
     else:
-        dpa = 0.5 * arctan2(-gamma, alpha - beta)
+        dpa = 0.5 * np.arctan2(-gamma, alpha - beta)
     #    if((s>=t)&(bmin!=smin)):
     #        dmin=sqrt (0.5 * (s - t))
 

racs_tools/beamcon_3D.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 """ Convolve ASKAP cubes to common resolution """
 __author__ = "Alec Thomson"
 

racs_tools/convolve_uv.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 """ Fast convolution in the UV domain """
 __author__ = "Wasim Raja"
 
@@ -63,8 +63,6 @@ def convolve(
         bpa=new_beam.pa.to(units.deg).value,
         u=u_image,
         v=v_image,
-        nx=nx,
-        ny=ny,
     )
     # Perform the x-ing in the FT domain
     im_f = np.fft.fft2(image)

racs_tools/gaussft.f

@@ -0,0 +1,84 @@
+#!/usr/bin/env python3
+"""
+Code to generate FT of final 2D-Gaussian to be used
+for convolving an image. The code deconvolves the input
+psf. The intrinsic psf must be specified.
+
+Python version of gaussft.f by Wasim Raja
+"""
+from typing import Tuple
+
+import numba as nb
+import numpy as np
+
+
+@nb.njit(cache=True)
+def gaussft(
+    bmin_in: float,
+    bmaj_in: float,
+    bpa_in: float,
+    bmin: float,
+    bmaj: float,
+    bpa: float,
+    u: np.ndarray,
+    v: np.ndarray,
+) -> Tuple[np.ndarray, float]:
+    """
+    Compute the Fourier transform of a 2D Gaussian for convolution.
+
+    Parameters:
+        bmin_in (float): Intrinsic psf BMIN (degrees)
+        bmaj_in (float): Intrinsic psf BMAJ (degrees)
+        bpa_in (float): Intrinsic psf BPA (degrees)
+        bmin (float): Final psf BMIN (degrees)
+        bmaj (float): Final psf BMAJ (degrees)
+        bpa (float): Final psf BPA (degrees)
+        u (np.ndarray): Fourier coordinates corresponding to image coord x
+        v (np.ndarray): Fourier coordinates corresponding to image coord y
+
+    Returns:
+        g_final (np.ndarray): Final array to be multiplied to FT(image) for convolution
+            in the FT domain.
+        g_ratio (float): Factor for flux scaling
+    """
+    deg2rad = np.pi / 180.0
+
+    bmaj_in_rad, bmin_in_rad, bpa_in_rad = (
+        bmaj_in * deg2rad,
+        bmin_in * deg2rad,
+        bpa_in * deg2rad,
+    )
+    bmaj_rad, bmin_rad, bpa_rad = bmaj * deg2rad, bmin * deg2rad, bpa * deg2rad
+
+    sx, sy = bmaj_rad / (2 * np.sqrt(2.0 * np.log(2.0))), bmin_rad / (
+        2 * np.sqrt(2.0 * np.log(2.0))
+    )
+    sx_in, sy_in = bmaj_in_rad / (2.0 * np.sqrt(2.0 * np.log(2.0))), bmin_in_rad / (
+        2.0 * np.sqrt(2.0 * np.log(2.0))
+    )
+
+    u_cosbpa, u_sinbpa = u * np.cos(bpa_rad), u * np.sin(bpa_rad)
+    u_cosbpa_in, u_sinbpa_in = u * np.cos(bpa_in_rad), u * np.sin(bpa_in_rad)
+
+    v_cosbpa, v_sinbpa = v * np.cos(bpa_rad), v * np.sin(bpa_rad)
+    v_cosbpa_in, v_sinbpa_in = v * np.cos(bpa_in_rad), v * np.sin(bpa_in_rad)
+
+    g_amp = np.sqrt(2.0 * np.pi * sx * sy)
+
+    dg_amp = np.sqrt(2.0 * np.pi * sx_in * sy_in)
+
+    g_ratio = g_amp / dg_amp
+
+    # Vectorized calculation of ur, vr, g_arg, and dg_arg
+    ur = u_cosbpa[:, np.newaxis] - v_sinbpa[np.newaxis, :]
+    vr = u_sinbpa[:, np.newaxis] + v_cosbpa[np.newaxis, :]
+    g_arg = -2.0 * np.pi**2 * ((sx * ur) ** 2 + (sy * vr) ** 2)
+
+    ur_in = u_cosbpa_in[:, np.newaxis] - v_sinbpa_in[np.newaxis, :]
+    vr_in = u_sinbpa_in[:, np.newaxis] + v_cosbpa_in[np.newaxis, :]
+    dg_arg = -2.0 * np.pi**2 * ((sx_in * ur_in) ** 2 + (sy_in * vr_in) ** 2)
+
+    # Vectorized calculation of g_final
+    g_final = g_ratio * np.exp(g_arg - dg_arg)
+
+    return g_final, g_ratio

racs_tools/gaussft.py

@@ -1,4 +1,4 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
 """ Find bad channels by checking statistics of each channel image. """
 
 import argparse