Merge pull request #18 from radionets-project/mojave_simulation

Mojave simulation

README.md

@@ -1,2 +1,30 @@
 # radiosim
 Simulation of radio skies to create astrophysical data sets.
+This repository is part of the [`radionets-project`](https://github.com/radionets-project).
+
+## Installation
+
+This repository is built as a python package. We recommend creating a mamba environment to handle the dependencies of all packages.
+You can create one by running the following command in this repository:
+```
+$ mamba env create -f environment.yml
+```
+You can start the environment with
+```
+$ mamba activate radiosim
+```
+after the installation.
+
+## Usage
+
+There are currently three supported simulation types: 
+1. `survey` full sky simulation
+2. `jet` extended source
+3. `mojave` MOJAVE like extended source
+
+In the `radiosim` environment you can start the simulation with
+```
+$ radiosim path/to/rc/file.toml
+```
+You can find an exemplary file in `rc/default_simulation.toml`.
+The simulations will be saved as `.h5` files.

pyproject.toml

@@ -30,6 +30,7 @@ classifiers=[
 dependencies = [
   "astropy",
   "click",
+  "scikit-image",
   "h5py",
   "matplotlib",
   "numpy",

radiosim/gauss.py

@@ -1,6 +1,8 @@
 import numpy as np
+from numpy.typing import ArrayLike
 from astropy.convolution import Gaussian2DKernel
 from astropy.nddata.utils import Cutout2D
+from scipy.stats import norm, skewnorm
 
 
 def gauss(params: list, size: int):
@@ -82,3 +84,60 @@ def rotgauss(x, y):
         return g
 
     return rotgauss(*np.indices((size, size)))
+
+
+def skewed_gauss(
+    size: int,
+    x: float,
+    y: float,
+    amp: float,
+    width: float,
+    length: float,
+    a: float,
+) -> ArrayLike:
+    """
+    Generate a skewed 2d normal distribution.
+
+    Parameters
+    ----------
+    size: int
+        length of the square image
+    x: float
+        x coordinate of the center
+    y: float
+        y coordinate of the center
+    amp: float
+        maximal amplitude of the distribution
+    width: float
+        width of the distribution (perpendicular to the skewed function)
+    length: float
+        length of the distribution
+    a: float
+        skewness parameter
+
+    Returns
+    -------
+    skew_gauss: ArrayLike
+        two dimensional skewed gaussian distribution
+    """
+
+    def skewfunc(x: float, **args) -> ArrayLike:
+        func = skewnorm.pdf(x, **args)
+        func /= func.max()
+        return func
+
+    vals = np.arange(size)
+    normal = norm.pdf(vals, loc=size / 2, scale=width).reshape(-1, 1)
+    normal /= normal.max()
+
+    skew = skewfunc(vals, a=a, loc=size / 2, scale=length)
+
+    skew_gauss = normal * skew
+    skew_gauss *= amp
+    skew_gauss[np.isclose(skew_gauss, 0)] = 0
+
+    skew_gauss = np.roll(
+        skew_gauss, (y - int(size / 2), x - int(size / 2)), axis=(0, 1)
+    )
+
+    return skew_gauss