Synth. Obs.: 3D Phantom - Reduction
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.
[1]:
import magritte.core as magritte # Core functionality
import magritte.plot as plot # Plotting
import magritte.tools as tools # Save fits
import os
from astropy import units # 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.
[2]:
wdir = "/lhome/thomasc/Magritte-examples/Phantom_3D/"
Define file names.
[3]:
model_file = os.path.join(wdir, 'model_Phantom_3D_red.hdf5') # reduced 3D Phantom Magritte model
Load the Magritte model.
[4]:
model = magritte.Model(model_file)
-------------------------------------------
Reading Model...
-------------------------------------------
model file = /lhome/thomasc/Magritte-examples/Phantom_3D/model_Phantom_3D_red.hdf5
-------------------------------------------
Reading parameters...
Reading points...
Reading rays...
Reading boundary...
Reading chemistry...
Reading species...
Reading thermodynamics...
Reading temperature...
Reading turbulence...
Reading lines...
Reading lineProducingSpecies...
Reading linedata...
read num 0
read sym CO
nlev = 41
nrad = 1
Reading collisionPartner...
Reading collisionPartner...
Reading quadrature...
Reading radiation...
Reading frequencies...
Not using scattering!
-------------------------------------------
Model read, parameters:
-------------------------------------------
npoints = 100159
nrays = 2
nboundary = 386
nfreqs = 31
nspecs = 3
nlspecs = 1
nlines = 1
nquads = 31
-------------------------------------------
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.
[5]:
model.compute_spectral_discretisation ()
model.compute_inverse_line_widths ()
model.compute_LTE_level_populations ()
Computing spectral discretisation...
Computing inverse line widths...
Computing LTE level populations...
[5]:
0
In this example we will work with the LTE level populations and do not demand statistical equilibrium.
[6]:
# Iterate level populations until statistical equilibrium
# model.compute_level_populations_sparse (True, 20)
Make synthetic observations
Now we can make synthetic observations of the model.
[7]:
fcen = model.lines.lineProducingSpecies[0].linedata.frequency[0]
vpix = 1500 # velocity pixel size [m/s]
dd = vpix * (model.parameters.nfreqs()-1)/2 / magritte.CC
fmin = fcen - fcen*dd
fmax = fcen + fcen*dd
model.compute_spectral_discretisation (fmin, fmax, 31)#bins using the specified amount of frequency bins (31). Can be any integer >=1
model.compute_image_new (model.parameters.hnrays()-1, 512, 512)#using a resolution of 512x512 for the image.
#Instead of definining a ray index [0, nrays-1], you can also define a ray direction for the imager
#model.compute_image_new (rx, ry, rz, 512, 512)#in which (rx, ry, rz) is the (normalized) ray direction
Computing spectral discretisation...
Computing image new...
[7]:
0
Plot observations
Plot the resulting channel maps with the Magritte matplotlib back end.
[8]:
plot.image_mpl(
model,
image_nr = -1,
zoom = 1.3,
npix_x = 240,
npix_y = 240,
x_unit = units.au,
v_unit = units.km / units.s
)
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 31/31 [00:14<00:00, 2.15it/s]
[8]:
<function magritte.plot.image_mpl.<locals>.<lambda>(v)>
(The plot is only interactive in a live notebook.)
Save the image cube in a fits file.
[9]:
tools.save_fits(model)
Written file to: /lhome/thomasc/Magritte-examples/Phantom_3D/images/image.fits
(Optional: To create your own plots) Overview of data stored in the Image object
[10]:
import numpy as np
xdir = np.array(model.images[-1].image_direction_x)#directions of the x-and y-vectors of the image
ydir = np.array(model.images[-1].image_direction_y)
zdir = np.array(model.images[-1].image_direction_z)#this is direction in which we observe the object
print("image directions: ", xdir, ydir, zdir)
nfreqs = np.array(model.images[-1].nfreqs) #number of frequency bins
freqs = np.array(model.images[-1].freqs) #frequency bins [Hz]
print("# of frequencies: ", nfreqs, " frequencies :", freqs)
ImX = np.array(model.images[-1].ImX)#X position in image [m]
ImY = np.array(model.images[-1].ImY)#Y position in image [m]
I = np.array(model.images[-1].I)#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: [1. 0. 0.] [0. 1. 0.] [0. 0. 1.]
# of frequencies: 31 frequencies : [1.15262550e+11 1.15263127e+11 1.15263704e+11 1.15264281e+11
1.15264857e+11 1.15265434e+11 1.15266011e+11 1.15266588e+11
1.15267165e+11 1.15267741e+11 1.15268318e+11 1.15268895e+11
1.15269472e+11 1.15270048e+11 1.15270625e+11 1.15271202e+11
1.15271779e+11 1.15272355e+11 1.15272932e+11 1.15273509e+11
1.15274086e+11 1.15274662e+11 1.15275239e+11 1.15275816e+11
1.15276393e+11 1.15276969e+11 1.15277546e+11 1.15278123e+11
1.15278700e+11 1.15279276e+11 1.15279853e+11]
[ ]: