Image-based analysis (experimental)

Although FRAP primarily operates on the visibility plane, there are situations where visibility analysis is not suitable. For instance, if a disk exhibits significant azimuthal asymmetry, azimuthal averaging is not adequate. In such cases, FRAP can also process radial intensity profiles. While this function has not been extensively tested, we demonstrate its usage in this notebook.

The basic usage is the same as the normal procedure.

import frap
import numpyro
import pickle
import dsharp_opac
import jax.numpy as jnp

numpyro.set_host_device_count(10)

In this example, we use a mock data set, and assume a cut along the major axis.

bands = [ 'band3', 'band6' ]
data = {
    'band3' : 'data_band3.pkl',
    'band6' : 'data_band6.pkl'
}

obs = {}

for band in bands:
    with open(data[band], 'rb') as f:
        obs[band] = pickle.load(f)

he structure of the data should be like this.

obs['band6'].keys()

dict_keys(['r', 'Tb', 'Cov', 'nu'])

The unit of Tb should be brightness temprature. Cov is the covariance matrix of the noise.

import matplotlib.pyplot as plt
plt.plot( obs['band6']['r'], obs['band6']['Tb'] )

[<matplotlib.lines.Line2D at 0x113a5d450>]

Those settings are identical to the visibility-based analysis

dust_opacity = jnp.load(dsharp_opac.get_datafile('default_opacities_smooth.npz'))

disk = frap.model( incl= 6.0,
                     r_out = 1.3,
                     N_GP = 130,
                     flux_uncert = True )

D = 114.9 # pc

disk.set_parameter('log10_T',
                   free = False,
                   profile = lambda r: jnp.log10(77.0*(r*D/10)**(-0.5)) )


disk.set_parameter('q',
                   free = False,
                   profile = lambda r: 3.5 )


lengthscale = 0.03

disk.set_parameter('log10_Sigma_d', free = True,
                                    bounds = ( -6.0, 3.0 ),
                                    lengthscale = lengthscale)

disk.set_parameter('log10_a_max', free = True,
                                    bounds =  ( -2.0, 2.0 ),
                                    lengthscale = lengthscale)

Then you can use set_radialprofile() to set your radial profile. Here you also have to give the covariance matrix.

f_s = { 'band3': 0.025, 'band6': 0.05 }
f_mean = { 'band3': 1.0, 'band6': 1.0 }

for band in bands:
    disk.set_radialprofile( band = band,
                            r = [obs[band]['r']],
                            Tb = [obs[band]['Tb']],
                            nu = [obs[band]['nu']],
                            Nch = 1,
                            f_s = f_s[band],
                            f_mean = f_mean[band],
                            FWHM = 0.05,
                            Cov = [obs[band]['Cov']] )

Remaining part is the same as before!

disk.set_opacity( opac_dict = dust_opacity )

inference = frap.inference(disk)

num_warmup = 500
num_samples = 500
num_chains = 10

posterior = inference.MCMC( num_warmup = num_warmup,
                            num_samples = num_samples,
                            num_chains = num_chains )

                            mean       std    median      5.0%     95.0%     n_eff     r_hat
       g_f_scale_band3     -0.72      0.75     -0.72     -1.97      0.53    435.13      1.02
       g_f_scale_band6     -1.23      0.27     -1.23     -1.70     -0.79    538.59      1.02
  g_log10_Sigma_d_z[0]     -1.06      0.57     -0.95     -1.91     -0.23    353.63      1.02
... [Output Truncated]

plotter = frap.plot( posterior )

saved keys:  dict_keys(['warmup_f', 'posterior_f', 'warmup_g', 'posterior_g'])

plotter.sample_paths( 'posterior_f', nskip = 10, plot_kwargs = { 'alpha' : 0.01, 'color' : 'b' } )