Synth. Obs.: Analytic disk
We create synthetic observations for the Magritte model of the analytic spiral 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/Analytic_disk/"
Define file names.
[3]:
model_file = os.path.join(wdir, 'model_analytic_disk.hdf5') # Analytic spiral Magritte model
Load the Magritte model.
[4]:
model = magritte.Model(model_file)
-------------------------------------------
Reading Model...
-------------------------------------------
model file = /lhome/thomasc/Magritte-examples/Analytic_disk/model_analytic_disk.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 = 114051
nrays = 2
nboundary = 3000
nfreqs = 51
nspecs = 3
nlspecs = 1
nlines = 1
nquads = 51
-------------------------------------------
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 = 1.0e+3 # velocity pixel size [m/s]
dd = vpix * (model.parameters.nfreqs()-1)/2 / magritte.CC
fmin = fcen - fcen*dd
fmax = fcen + fcen*dd
# Ray orthogonal to image plane
ray_nr = 0
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 (ray_nr, 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 matplotlib. Note that the resolutions do not match, as we here plot a zoomed in version.
[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%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 51/51 [00:23<00:00, 2.21it/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/Analytic_disk/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. -1. 0.] [ 0. 0. -1.] [1. 0. 0.]
# of frequencies: 51 frequencies : [1.15261589e+11 1.15261974e+11 1.15262358e+11 1.15262743e+11
1.15263127e+11 1.15263512e+11 1.15263896e+11 1.15264281e+11
1.15264665e+11 1.15265050e+11 1.15265434e+11 1.15265819e+11
1.15266203e+11 1.15266588e+11 1.15266972e+11 1.15267357e+11
1.15267741e+11 1.15268126e+11 1.15268510e+11 1.15268895e+11
1.15269279e+11 1.15269664e+11 1.15270048e+11 1.15270433e+11
1.15270817e+11 1.15271202e+11 1.15271586e+11 1.15271971e+11
1.15272355e+11 1.15272740e+11 1.15273124e+11 1.15273509e+11
1.15273893e+11 1.15274278e+11 1.15274662e+11 1.15275047e+11
1.15275431e+11 1.15275816e+11 1.15276200e+11 1.15276585e+11
1.15276969e+11 1.15277354e+11 1.15277738e+11 1.15278123e+11
1.15278507e+11 1.15278892e+11 1.15279276e+11 1.15279661e+11
1.15280045e+11 1.15280430e+11 1.15280814e+11]
[ ]: