Synth. Obs.: Analytic spiral

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

Define file names.

model_file = os.path.join(wdir, 'model_analytic_spiral.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_spiral/model_analytic_spiral.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, 20, True)

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]
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)
# Ray orthogonal to image plane
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.

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%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 51/51 [00:33<00:00,  1.55it/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_spiral/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.15258712e+11 1.15259212e+11 1.15259711e+11 1.15260211e+11
 1.15260711e+11 1.15261211e+11 1.15261710e+11 1.15262210e+11
 1.15262710e+11 1.15263209e+11 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 1.15279200e+11 1.15279700e+11 1.15280200e+11
 1.15280699e+11 1.15281199e+11 1.15281699e+11 1.15282199e+11
 1.15282698e+11 1.15283198e+11 1.15283698e+11] Hz