Oftentimes it is necessary to get a symmetric coordinate array that spans $N$ elements from -width to +width, but makes sure that the middle point lands on $0$. The indices go from $0$ to $N -1.$ linspace returns the end points inclusive, wheras we want to leave out the right endpoint, because we are sampling the function in a cyclic manner.

A dictionary of parameter lists for conversion.

Define a grid object which stores all of these variables This will not change for the duration of the run.

Create a grid object using a logarithmic then linear then logarithmic radial spacing

Read from a dictionary, then choose how to make the grid based upon the arguments.

Hold the ray-tracing grid.

Hold the ray-tracing grid.

Parameters for the NUKER model.

Parameters for the standard model.

Parameters for the vertical model.

Parameters for the vertical model with variable slope.

Sigma(r::Float64, pars::ParametersNuker)

Calculate the gas surface density using the Nuker profile.

Calculate the gas surface density

Used to turn a dictionary of parameter values (from config.yaml) directly into a parameter type. Generally used for synthesis and plotting command line scripts.

Unroll a vector of parameter values into a parameter type.

Generate a boolean mask corresponding to the velocities which fall inside the user-specified mask ranges. arr is an array of [start, end] pairs, like arr = [[1.0, 3.0], [4.0, 5.2]] vels are the velocities corresponding to the dataset, like vels = [1.2, 2.4, 3.6], etc.

The common sense priors that apply to all parameter values

This function only needs to be run once, upon setup.

Write the wavelength sampling file. Only run on setup

The image module contains various data types for reading and holding images produced by the radiative transfer programs (via RADMC-3D), as well as routines for processing these images.

RawImage(data, pixsize_x, pixsize_y, lams)

Hold the raw output from RADMC-3D in a 3D array (npixy, npixx, nlam). RawImage reflects the RADMC convention that both x and y are increasing with array index. This means that to display the image as RADMC intends it, you must set the first array element to the lower left corner.

SkyImage is a holder that has both RA and DEC increasing with array index This convention is necessary for the FFT step However, to display this image in the traditional sky convention (North up, East to the left), you must set the first array element to the lower left corner and flip the array along the RA axis: fliplr(data) or flipdim(data, 2)

SkyImage constructor for just a single frame

TausurfImage is designed to hold the results of the RADMC-3D tausurf operation. This is The distance in cm above or below the plane tangent to the observer, which intersects the origin of the model.

From the RADMC3D manual: The image output file image.out will now contain, for each pixel, the position along the ray in centimeters where τ = τs. The zero point is the surface perpendicular to the direction of observation, going through the pointing position (which is, by default the origin (0, 0, 0)). Positive values mean that the surface is closer to the observer than the plane, while negative values mean that the surface is behind the plane. So for this datastructure, it's the same thing as RawImage, just instead of intensity, we have distance above/below plane.

Encapsulates the 3D position of the pixels representing the tau=1 surface, in the same datashape as the image. For each pixel, this is the x, y, or z position.

Storage for zeroth moment map

Subtraction for SkyImages

convert(::Type{ZerothMoment}, img::SkyImage)

Convert a SkyImage to a ZerothMoment map.

imToSky(img::RawImage, dpc::Float64)

Convert a RawImage to a SkyImage. Assumes dpc is parsecs.

Take an image and integrate all the frames to create a spatially-integrated spectrum

imread(file="image.out")

Read the image file (default=image.out) and return it as an Image object, which contains the fluxes in Jy/pixel, the sizes and locations of the pixels in arcseconds, and the wavelengths (in microns) corresponding to the images

Calculate the integrated line flux

taureadImg(file="image_tausurf.out")

Like imread, but for tausurf. Pixels that have no $\tau$ surface are set to NaN.

Read the (x,y,z) positions of the $\tau=1$ pixels.

DataVis(fname::AbstractString, index::Int, flagged::Bool=false)

Read just one channel of visibilities from the HDF5 file.

DataVis(fname::AbstractString, mask::BitArray, flagged::Bool=false)

Read just a subset of channels from the HDF5 file and return an array of DataVis.

DataVis(fname::AbstractString, indices::Vector{Int}, flagged::Bool=false)

Read just a subset of channels from the HDF5 file and return an array of DataVis.

DataVis(fname::AbstractString, flagged::Bool=false)

Read all the visibilities from an HDF5 string and return them as an array of DataVis objects flagged means whether we should be loading the flagged points or not.

DataVisReal(dvis::DataVis)

Convert a DataVis object, with visibilities stored as Comlpex128, into an object with visibilities stored with real and complex values stored separately as Float64.

ModelVis(dvis::DataVis, fmvis::FullModelVis)

Given a DataSet and a FullModelVis, go through and interpolate at the u,v locations of the DataVis.

ModelVisReal(mvis::ModelVis)

Convert a ModelVis object, with visibilities stored as ComplexF64, into a ModelVisReal, with real and imaginary visibility values stored separately as Float64.

RawModelVis

Type to store the product of FFT'ing a single channel of a SkyImage.

-(vis1::FullModelVis, vis2::FullModelVis)

Subtract two FullModelVis's.

conj!(dv::DataVis)

Import the complex conjugate function from Base, and extend it to work on a DataVis. I think this is necessary because the SMA and ALMA baseline conventions are swapped from what I'm using in the NRAO Synthesis Summer School textbook.

conj!(dvarr::Array{DataVis, 1})

Apply $conj!$ to an entire DataVis array.

ModelVis2DataVis(mvis::ModelVis)

Collapse a ModelVis into a DataVis in order to write to disk.

ModelVisRotate(dvis::DataVis, fmvis::FullModelVis, PA::Real)

Given a DataSet, FullModelVis, and a desired position angle rotation, use the Fourier rotation theorem to sample the image at the rotated baselines and produce a rotated (sampled) model. A positive PA value (in degrees) means that the new model is rotated $PA$ number of degrees counter-clockwise in the image plane (towards East, from North).

ResidVis(dvarr::Array{DataVis, 1}, mvarr)

Return a model visibility file that contains the residuals.

chi2(dvis::DataVis, mvis::ModelVis)

Calculate the $\chi^2$ value of the current model. This function is not used in inference (see lnprob), but it is sometimes a useful number to report.

copy_flags(source::AbstractString, dest::AbstractString)

Copy the flags from one dataset to another.

fftfreq(n::Int, d::Float64)

After numpy.fft.fftfreq

f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (d*n)   if n is even
f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (d*n)   if n is odd.

fillModelVis(vis::RawModelVis)

This function is designed to copy the partial arrays in RawModelVis into a full image for easy plotting. This means that the u axis can remain the same but we'll need to make the complex conjugate of the v axis.

get_nyquist_pixel(max_base::Float64, angular_width::Float64)

Determine how many pixels the image needs at a given distance in order to satisfy the Nyquist sampling theorem.

max_base is in kilolambda. angular_width is in radians.

get_qq(dv::DataVis)

Given a DataVis, convert the Cartesian spatial frequency coordinates $(u,v)$ into a radial spatial coordinate $q$.

interpolate_uv(u::Float64, v::Float64, vis::FullModelVis)

Interpolates a dense grid of visibilities (e.g., from FFT of an image) to a specfic (u,v) point using spheroidal functions in a band-limited manner designed to reduce aliasing.

lnprob(dvis::DataVis, mvis::DataVis)

Lnprob for computing the difference between two datasets. Does not check whether the baselines of the two datasets are the same. Use this function at your own discretion.

lnprob(dvis::DataVis, mvis::ModelVis, mu_RA::Real, mu_DEC::Real)

Calculate the lnlikelihood using a single loop to do the phase shifting and summation. No extra memory needed for re-copying the model.

lnprob(dvis::DataVis, mvis::ModelVis)

Likelihood function used to calculate the likelihood of the data visibilities given a current model of them. Asserts that the model visibilities were sampled at the same baselines as the data.

lnprob_real(dvis::DataVis, mvis::ModelVis, mu_RA::Float64, mu_DEC::Float64)

Calculate the lnlikelihood using a single loop to do the phase shifting and summation. No extra memory needed for re-copying the model. Uses Euler formula to do the calculation, rather than exponentiating a complex number.

max_baseline(dvarr)

Determine the maximum uu or vv baseline contained in the dataset, so we know at what resolution we will need to synthesize the images.

returned in kilolambda.

phase_shift!(fvis::FullModelVis, mu_RA, mu_DEC)

Given a new model centroid in the image plane (in arcseconds), shift the model visibilities by corresponding amount.

phase_shift!(mvis::ModelVis, mu_RA, mu_DEC)

Given a new model centroid in the image plane (in arcseconds), shift the model visibilities by corresponding amount.

plan_interpolate(dvis::DataVis, uu::Vector{Float64}, vv::Vector{Float64})

Return a function that is used to interpolate the visibilities, in the spirit of interpolate_uv but much faster. Closures save time and money!

The uu and vv vectors are the visibility coordinates for the dense visibility model.

The point is that if the distance to the source and the size of the sythesized image are not changing, then we will always be interpolating from a dense Visibility grid that has the exact same U and V spacings, meaning that the weighting terms used to evaluate the interpolation for a specific visibility can be cached.

So, using a closure, those weighting terms are calculated once and then cached for further use.

rfftfreq(n::Int, d::Float64)

Return the frequencies corresponding to the output of the real FFT. After numpy.fft.rfftfreq f = [0, 1, ..., n/2-1, n/2] / (dn) if n is even f = [0, 1, ..., (n-1)/2-1, (n-1)/2] / (dn) if n is odd.

transform(img::SkyImage, index::Int=1)

Transform the SkyImage produced by RADMC using FFT.

write(dv::DataVis, fname::AbstractString)

Take in a visibility data set and then write it to the HDF5 file The HDF5 file actually expects multi-channel data, so instead we will need to store all of this information with arrays of shape (..., 1) [an extra trailing] dimension of 1.

write(dvarr::Array{DataVis, 1}, fname::AbstractString)

Write an array of DataVis to a single HDF5 file.

corrfun!(img::SkyImage, mu_RA, mu_DEC)

Apply the correction function to the SkyImage, with an offset in RA and DEC.

corrfun!(img::SkyImage)

Apply the correction function to (and mutate) a SkyImage in place.

corrfun(img::SkyImage)

Apply the correction function to a SkyImage, but return a copy of the image, leaving the original unchanged.

corrfun(eta::T) where {T}

Gridding correction function, but able to be passed either floating point numbers or vectors of Float64.

corrfun(eta::T, alpha::Float64) where {T}

Gridding correction function, used to pre-divide the image to correct for the effect of the gcffun. This function is also the Fourier transform of gcffun.

gcffun(eta::T) where {T}

The gridding convolution function, used to do the convolution and interpolation of the visibilities in the Fourier domain. This is also the Fourier transform of corrfun.

gcffun(eta::T, alpha::Float64) where {T}

The gridding convolution function, with variable $\alpha$.

spheroid(eta, alpha)

Prolate spheroidal wavefunction, assuming that $m = 6$. This allows arguments for $\eta < 1.0 + 10^{-7}$, but returns 0.0 (i.e., the spheroid window is truncated).

spheroid(eta)

spheroid function which assumes $\alpha$ = 1.0, $m=6$, built for speed.

Provide routines for making fake Gaussian models which have analytic Fourier transforms, for the purposes of testing.

FTGauss(uu::AbstractVector{Float64}, vv::AbstractVector{Float64}, p::Vector{Float64}, k::Int)

p is a length 5 vector of mu_RA, mu_DEC, sigma_RA, sigma_DEC, rho.

Given two arrays of u and v coordinates in [kλ], fill an array with the analytic FT of imageGauss evaluated at every pairwise (u,v) pair.

FTGauss(uu::Float64, vv::Float64, p::Vector{Float64}, k::Real)

Given u and v coordinates in [kλ], evaluate the analytic FT of the aforementioned Gaussian.

p is a length 5 vector of [mu_RA, mu_DEC, sigma_x, sigma_y, rho] in units of arcseconds, corresponding to the image plane.

N.B. Here Sigma refers to the (same) covariance matrix in the image domain.

This function always returns a complex value.

imageGauss(ll::AbstractVector{Float64}, mm::AbstractVector{Float64}, p::Vector{Float64}, k::Real; theta=0)

Given two arrays of $l$ and $m$ coordinates, corresponding to $x$ and $y$, fill an array of the Gaussian image following the MATLAB convention for images (each row corresponds to a different y value).

p0 is a vector of [mu_RA, mu_DEC, sigma_x, sigma_y, rho] in units of arcseconds. $\rho$ is the correlation of the Gaussian, ranging from 0 to 1.

\[\rho = \frac{\sigma_{xy}}{\sigma_x \sigma_y}\]

$\mu_\alpha$ and $\mu_\delta$ are the locations of the centroid emission relative to the image origin (RA=0, DEC=0).

$k$ is a scaling pre-factor to adjust the amplitude of the Gaussian.

The image intensity as a function of $l$ and $m$ (sky plane coordinates) is

\[ I(l,m) = \frac{k}{2 \pi \sqrt{|\boldsymbol{\Sigma}|}} \exp \left \{ -\frac{1}{2} \boldsymbol{R}^\mathrm{T} \boldsymbol{\Sigma}^{-1} \boldsymbol{R} \right \}\]

Where $\boldsymbol{\Sigma}$ is

\[\boldsymbol{\Sigma} = \left [ \begin{array}{cc} \sigma_x & \sigma_{xy} \\ \sigma_{xy} & \sigma_y \\ \end{array} \right ]\]

$\theta$ is a desired angle to rotate the Gaussian about the origin. To make sense of the order of operations, the rotation must occur about the origin. A positive value means a counter-clockwise rotation of the Gaussian from North towards East.