`PPVCube`

analysis module to project fields along a given line of sight traveling at different line-of-sight velocities, to "mock-up" what would be seen in observations.

In [1]:

```
from yt.config import ytcfg
import yt
import numpy as np
from yt.analysis_modules.ppv_cube.api import PPVCube
import yt.units as u
```

Density: $\rho(r) \propto r^{\alpha}$

Velocity: $v_{\theta}(r) \propto \frac{r}{1+(r/r_0)^{\beta}}$

where for simplicity we won't worry about the normalizations of these profiles.

First, we'll set up the grid and the parameters of the profiles:

In [2]:

```
# increasing the resolution will make the images in this notebook more visually appealing
nx,ny,nz = (64, 64, 64) # domain dimensions
R = 10. # outer radius of disk, kpc
r_0 = 3. # scale radius, kpc
beta = 1.4 # for the tangential velocity profile
alpha = -1. # for the radial density profile
x, y = np.mgrid[-R:R:nx*1j,-R:R:ny*1j] # cartesian coordinates of x-y plane of disk
r = np.sqrt(x*x+y*y) # polar coordinates
theta = np.arctan2(y, x) # polar coordinates
```

`velx`

and `vely`

. Everywhere outside the disk, all fields are set to zero.

In [3]:

```
dens = np.zeros((nx,ny,nz))
dens[:,:,nz//2-3:nz//2+3] = (r**alpha).reshape(nx,ny,1) # the density profile of the disk
temp = np.zeros((nx,ny,nz))
temp[:,:,nz//2-3:nz//2+3] = 1.0e5 # Isothermal
vel_theta = 100.*r/(1.+(r/r_0)**beta) # the azimuthal velocity profile of the disk
velx = np.zeros((nx,ny,nz))
vely = np.zeros((nx,ny,nz))
velx[:,:,nz//2-3:nz//2+3] = (-vel_theta*np.sin(theta)).reshape(nx,ny,1) # convert polar to cartesian
vely[:,:,nz//2-3:nz//2+3] = (vel_theta*np.cos(theta)).reshape(nx,ny,1) # convert polar to cartesian
dens[r > R] = 0.0
temp[r > R] = 0.0
velx[r > R] = 0.0
vely[r > R] = 0.0
```

`load_uniform_grid`

. We'll define the width of the grid to be `2*R`

kpc, which will be equal to 1 `code_length`

.

In [4]:

```
data = {}
data["density"] = (dens,"g/cm**3")
data["temperature"] = (temp, "K")
data["velocity_x"] = (velx, "km/s")
data["velocity_y"] = (vely, "km/s")
data["velocity_z"] = (np.zeros((nx,ny,nz)), "km/s") # zero velocity in the z-direction
bbox = np.array([[-0.5,0.5],[-0.5,0.5],[-0.5,0.5]]) # bbox of width 1 on a side with center (0,0,0)
ds = yt.load_uniform_grid(data, (nx,ny,nz), length_unit=(2*R,"kpc"), nprocs=1, bbox=bbox)
```

To get a sense of what the data looks like, we'll take a slice through the middle of the disk:

In [5]:

```
slc = yt.SlicePlot(ds, "z", ["density","velocity_x","velocity_y","velocity_magnitude"])
```

In [6]:

```
slc.set_log("velocity_x", False)
slc.set_log("velocity_y", False)
slc.set_log("velocity_magnitude", False)
slc.set_unit("velocity_magnitude", "km/s")
slc.show()
```

`PPVCube`

class. First, let's assume we rotate our viewing angle 60 degrees from face-on, from along the z-axis into the x-axis. We'll create a normal vector:

In [7]:

```
i = 60.*np.pi/180.
L = [np.sin(i),0.0,np.cos(i)]
```

`(vmin, vmax, nbins, units)`

, which specifies a linear range of `nbins`

velocity bins from `vmin`

to `vmax`

in units of `units`

. We may also optionally specify the dimensions of the data cube with the `dims`

argument.

In [8]:

```
cube = PPVCube(ds, L, "density", (-150.,150.,50,"km/s"), dims=200, method="sum")
```

`length_unit`

:

In [9]:

```
cube.write_fits("cube.fits", clobber=True, length_unit="kpc")
```

Or one can use the `sky_scale`

and `sky_center`

keywords to set up the coordinates in RA and Dec:

In [10]:

```
sky_scale = (1.0, "arcsec/kpc")
sky_center = (30., 45.) # RA, Dec in degrees
cube.write_fits("cube_sky.fits", clobber=True, sky_scale=sky_scale, sky_center=sky_center)
```

In [11]:

```
ds_cube = yt.load("cube.fits")
```

In [12]:

```
# Specifying no center gives us the center slice
slc = yt.SlicePlot(ds_cube, "z", ["density"])
slc.show()
```

In [13]:

```
# Picking different velocities for the slices
new_center = ds_cube.domain_center
new_center[2] = ds_cube.spec2pixel(-100.*u.km/u.s)
slc = yt.SlicePlot(ds_cube, "z", ["density"], center=new_center)
slc.show()
```

In [14]:

```
new_center[2] = ds_cube.spec2pixel(70.0*u.km/u.s)
slc = yt.SlicePlot(ds_cube, "z", ["density"], center=new_center)
slc.show()
```

In [15]:

```
new_center[2] = ds_cube.spec2pixel(-30.0*u.km/u.s)
slc = yt.SlicePlot(ds_cube, "z", ["density"], center=new_center)
slc.show()
```

In [16]:

```
prj = yt.ProjectionPlot(ds_cube, "z", ["density"], method="sum")
prj.set_log("density", True)
prj.set_zlim("density", 1.0e-3, 0.2)
prj.show()
```

`thermal_broad`

keyword allows one to simulate thermal line broadening based on the temperature, and the `atomic_weight`

argument is used to specify the atomic weight of the particle that is doing the emitting.

In [17]:

```
cube2 = PPVCube(ds, L, "density", (-150.,150.,50,"km/s"), dims=200, thermal_broad=True,
atomic_weight=12.0, method="sum")
cube2.write_fits("cube2.fits", clobber=True, length_unit="kpc")
```

Taking a slice of this cube shows:

In [18]:

```
ds_cube2 = yt.load("cube2.fits")
new_center = ds_cube2.domain_center
new_center[2] = ds_cube2.spec2pixel(70.0*u.km/u.s)
slc = yt.SlicePlot(ds_cube2, "z", ["density"], center=new_center)
slc.show()
```

In [19]:

```
new_center[2] = ds_cube2.spec2pixel(-100.*u.km/u.s)
slc = yt.SlicePlot(ds_cube2, "z", ["density"], center=new_center)
slc.show()
```

where we can see the emission has been smeared into this velocity slice from neighboring slices due to the thermal broadening.

Finally, the "velocity" or "spectral" axis of the cube can be changed to a different unit, such as wavelength, frequency, or energy:

In [20]:

```
print (cube2.vbins[0], cube2.vbins[-1])
cube2.transform_spectral_axis(400.0,"nm")
print (cube2.vbins[0], cube2.vbins[-1])
```

In [21]:

```
cube2.reset_spectral_axis()
print (cube2.vbins[0], cube2.vbins[-1])
```