Synth. Obs.: 3D Phantom

We create synthetic observations for the Magritte model of the 3D Phantom snapshot that was created in the this example.

Setup

Import the required functionalty.

from magrittetorch.model.model import Model
import magrittetorch.tools.plot as plot
import magrittetorch.algorithms.solvers as solvers
from astropy import units, constants
import os
import matplotlib.pyplot as plt
import torch

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/Phantom_3D/"

Define file names.

model_file = os.path.join(wdir, 'model_Phantom_3D.hdf5')   # 3D Phantom Magrittetorch model

Load the Magritte model.

model = Model(model_file)
model.read()

Reading model from:  /lhome/thomasc/Magrittetorch-examples/Phantom_3D/model_Phantom_3D.hdf5
Reading Magrittetorch model

print(model.parameters.npoints.get())

1034671

Model the medium

Initialize the model

model.dataCollection.infer_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)

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(device)[1]
vpix = 1500   # velocity pixel size [m/s]
nfreqs = 31   # number of frequency bins
dd   = vpix * (nfreqs-1)/2 / constants.c.value # frequency bin size [Hz]
fmin = fcen - fcen*dd
fmax = fcen + fcen*dd

image_freqs = torch.linspace(fmin, fmax, nfreqs, device=device, dtype=torch.float32)

# Compute image in the specified direction at the given frequencies
solvers.image_model(model, torch.Tensor([0,0,1]).type(torch.float32).to(device), image_freqs, device)

Plot observations

Plot the resulting channel maps with the Magrittetorch matplotlib back end.

plot.image_mpl(
    model,
    image_nr =  -1,
    zoom     = 1.3,
    npix_x   = 256,
    npix_y   = 256,
    x_unit   = units.au,
    v_unit   = units.km / units.s,
    method = 'nearest'
)
plt.show()

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 31/31 [00:22<00:00,  1.41it/s]

(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/Phantom_3D/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()#directions of the x-and y-vectors of the image
ydir = model.images[-1].image_direction_y.get()
zdir = model.images[-1].image_direction_z.get()#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:  tensor([1., 0., 0.]) tensor([0.0000e+00, 1.0000e+00, 4.3711e-08]) tensor([ 0.0000e+00, -4.3711e-08,  1.0000e+00])
# of frequencies:  31  frequencies : [2.3051970e+11 2.3052085e+11 2.3052201e+11 2.3052316e+11 2.3052432e+11
 2.3052547e+11 2.3052662e+11 2.3052778e+11 2.3052893e+11 2.3053009e+11
 2.3053124e+11 2.3053238e+11 2.3053355e+11 2.3053469e+11 2.3053586e+11
 2.3053700e+11 2.3053815e+11 2.3053931e+11 2.3054046e+11 2.3054162e+11
 2.3054277e+11 2.3054392e+11 2.3054508e+11 2.3054623e+11 2.3054739e+11
 2.3054854e+11 2.3054968e+11 2.3055085e+11 2.3055199e+11 2.3055316e+11
 2.3055430e+11] Hz