Synth. Obs.: 3D AMRVAC - Reduction
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.
[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/AMRVAC_3D/"
Define file names.
[3]:
model_file = os.path.join(wdir, 'model_AMRVAC_3D_red.hdf5') # Reduced 3D AMRVAC Magritte model
Load the Magritte model.
[4]:
model = magritte.Model(model_file)
-------------------------------------------
Reading Model...
-------------------------------------------
model file = /lhome/thomasc/Magritte-examples/AMRVAC_3D/model_AMRVAC_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 = 32138
nrays = 12
nboundary = 1538
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 = 1300 # 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)#bins the frequency spectrum [fmin, fmax] into model.parameters.nfreqs bins.
# 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 = 256,
npix_y = 256,
x_unit = units.au,
v_unit = units.km / units.s
)
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 31/31 [00:15<00:00, 1.95it/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/AMRVAC_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: [-0.52704628 -0.84983659 0. ] [ 0.66666667 -0.41344912 -0.62017367] [ 0.52704628 -0.32686023 0.78446454]
# of frequencies: 31 frequencies : [1.15263704e+11 1.15264204e+11 1.15264704e+11 1.15265204e+11
1.15265703e+11 1.15266203e+11 1.15266703e+11 1.15267203e+11
1.15267703e+11 1.15268203e+11 1.15268703e+11 1.15269202e+11
1.15269702e+11 1.15270202e+11 1.15270702e+11 1.15271202e+11
1.15271702e+11 1.15272202e+11 1.15272701e+11 1.15273201e+11
1.15273701e+11 1.15274201e+11 1.15274701e+11 1.15275201e+11
1.15275700e+11 1.15276200e+11 1.15276700e+11 1.15277200e+11
1.15277700e+11 1.15278200e+11 1.15278700e+11]
[ ]: