Synth. Obs.: 3D AMRVAC

We create synthetic observations for the Magritte model of the 3D AMRVAC 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 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/AMRVAC_3D/"

Define file names.

model_file = os.path.join(wdir, 'model_AMRVAC_3D.hdf5')   # 3D AMRVAC Magritte model

Load the Magritte model.

model = Model(model_file)
model.read()

Reading model from:  /lhome/thomasc/Magrittetorch-examples/AMRVAC_3D/model_AMRVAC_3D.hdf5
Reading Magrittetorch model

Model the medium

Initialize the model by setting up a spectral discretisation, computing the inverse line widths and initializing the level populations with their LTE values.

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.

fcen = model.sources.lines.lineProducingSpecies[0].linedata.frequency.get()[0]
vpix = 1300   # velocity pixel size [m/s]
nfreqs = 31   # number of frequency channels
dd   = vpix * (nfreqs-1)/2 / constants.c.value
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.52704628, -0.32686023, 0.78446454]).type(torch.float32).to(device), image_freqs, device)

Plot observations

Plot the resulting channel maps with the Magritte 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
)

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 31/31 [00:24<00:00,  1.26it/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/AMRVAC_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([-0.5270, -0.8498,  0.0000]) tensor([ 0.6667, -0.4134, -0.6202]) tensor([ 0.5270, -0.3269,  0.7845])
# of frequencies:  31  frequencies : [1.15263709e+11 1.15264209e+11 1.15264709e+11 1.15265208e+11
 1.15265708e+11 1.15266208e+11 1.15266707e+11 1.15267207e+11
 1.15267707e+11 1.15268207e+11 1.15268706e+11 1.15269206e+11
 1.15269706e+11 1.15270205e+11 1.15270705e+11 1.15271205e+11
 1.15271705e+11 1.15272204e+11 1.15272704e+11 1.15273204e+11
 1.15273703e+11 1.15274203e+11 1.15274703e+11 1.15275203e+11
 1.15275702e+11 1.15276202e+11 1.15276702e+11 1.15277201e+11
 1.15277701e+11 1.15278201e+11 1.15278701e+11] Hz