Synth. Obs.: Analytic disk

We create synthetic observations for the Magritte model of the analytic spiral that was created in the this example.

Setup

Import the required functionalty.

import magrittetorch.model.model     as magritte   # Core functionality
import magrittetorch.algorithms.solvers as solvers       # Plotting
import magrittetorch.tools.plot    as plot      # Save fits
import os
import torch

from astropy import units, constants              # Unit conversions

Define a working directory (you will have to change this). We assume here that the scripts of the this example have already been executed and go back to that working directory.

wdir = "/lhome/thomasc/Magrittetorch-examples/Analytic_disk/"

Define file names.

model_file = os.path.join(wdir, 'model_analytic_disk.hdf5')   # Analytic spiral Magritte model

Load the Magritte model.

model = magritte.Model(model_file)
model.read(legacy_mode=False) # Load Magrittetorch model

Reading model from:  /lhome/thomasc/Magrittetorch-examples/Analytic_disk/model_analytic_disk.hdf5
Reading Magrittetorch model

Model the medium

Initialize the model.

model.dataCollection.infer_data() # Correctly initialize all data

In this example we will work with the LTE level populations and do not demand statistical equilibrium.

# If you have a GPU, we can use it to speed up the computations
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

# Getting the full model on GPU might be a bit too much for its memory, so you might want to use your CPU instead for the computations
# device = torch.device("cpu")

# Iterate level populations until statistical equilibrium
# solvers.compute_level_populations(model, device=device, 20, True)

Make synthetic observations

Now we can make synthetic observations of the model.

# Get the line frequency corresponding to the CO 2-1 line
fcen = model.sources.lines.lineProducingSpecies[0].linedata.frequency.get()[0]
vpix = 1.0e+3 # velocity pixel size [m/s]
dd   = vpix * 25 / constants.c.si.value
fmin = fcen - fcen*dd
fmax = fcen + fcen*dd

# Make sure that all input tensors lie on the same device (CPU or GPU)
device = torch.device("cpu")
freqs = torch.linspace(fmin, fmax, 51, dtype=torch.float32, device=device) # 25 frequency bins
raydir = torch.Tensor([0.0, 1.0, 0.0], device=device).type(torch.float32)
ray_nr = 3
raydir = model.geometry.rays.direction.get(device)[ray_nr]
solvers.image_model(model, raydir, freqs, torch.device("cpu"), Nxpix=256, Nypix=256)#using a resolution of 256x256 for the image.

Plot observations

Plot the resulting channel maps with matplotlib. Note that the resolutions do not match, as we here plot a zoomed in version.

plot.image_mpl(
    model,
    image_nr =  -1,
    zoom     = 1.3,
    npix_x   = 256,
    npix_y   = 256,cd
    x_unit   = units.au,
    v_unit   = units.km / units.s)

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 51/51 [00:25<00:00,  2.00it/s]

<function magrittetorch.tools.plot.image_mpl.<locals>.<lambda>(v)>

(The plot is only interactive in a live notebook.)

Save the image cube in a fits file.

plot.save_fits(model)

WARNING: No regularly spaced frequency bins!
Written file to: /lhome/thomasc/Magrittetorch-examples/Analytic_disk/images/image.fits

(Optional: To create your own plots) Overview of data stored in the Image object

xdir = model.images[-1].image_direction_x.get_astropy()#directions of the x-and y-vectors of the image
ydir = model.images[-1].image_direction_y.get_astropy()
zdir = model.images[-1].image_direction_z.get_astropy()#this is direction in which we observe the object
print("image directions: ", xdir, ydir, zdir)
nfreqs = model.images[-1].nfreqs.get() #number of frequency bins
freqs = model.images[-1].freqs.get_astropy() #frequency bins [Hz]
print("# of frequencies: ", nfreqs, " frequencies :", freqs)
ImX = model.images[-1].imX.get_astropy()#X position in image [m]
ImY = model.images[-1].imY.get_astropy()#Y position in image [m]
I = model.images[-1].I.get_astropy()#Intensity at the corresponding ImX, ImY position [W/(m^2*Hz*sr)], at a given frequency bin
# print("Intensities :", I, " ImX:", ImX, "ImY:", ImY) #prints a lot of output

image directions:  [ 0.6126789  -0.79033196  0.        ] [-0.6535241 -0.5066231 -0.5623516] [ 0.44444442  0.34454095 -0.8268982 ]
# of frequencies:  51  frequencies : [1.15261596e+11 1.15261981e+11 1.15262366e+11 1.15262751e+11
 1.15263136e+11 1.15263521e+11 1.15263906e+11 1.15264283e+11
 1.15264668e+11 1.15265053e+11 1.15265438e+11 1.15265823e+11
 1.15266208e+11 1.15266593e+11 1.15266978e+11 1.15267363e+11
 1.15267748e+11 1.15268133e+11 1.15268518e+11 1.15268895e+11
 1.15269280e+11 1.15269665e+11 1.15270050e+11 1.15270435e+11
 1.15270820e+11 1.15271205e+11 1.15271590e+11 1.15271975e+11
 1.15272360e+11 1.15272745e+11 1.15273130e+11 1.15273515e+11
 1.15273892e+11 1.15274277e+11 1.15274662e+11 1.15275047e+11
 1.15275432e+11 1.15275817e+11 1.15276202e+11 1.15276587e+11
 1.15276972e+11 1.15277357e+11 1.15277742e+11 1.15278127e+11
 1.15278504e+11 1.15278889e+11 1.15279274e+11 1.15279659e+11
 1.15280044e+11 1.15280429e+11 1.15280814e+11] Hz