<jupyter_start><jupyter_text>Tutorial for *diskmap* This tutorial showcases the functionalities of *diskmap*. We will use a polarized scattered light image of the LkCa 15 circumstellar disk. The data were obtained with [VLT/SPHERE](https://www.eso.org/sci/facilities/paranal/instruments/sphere.html) in the $J$ band and have been published by [Thalmann et al. (2016)](https://ui.adsabs.harvard.edu/abs/2016ApJ...828L..17T/abstract).  Getting started We start by importing the Python modules that are required for this tutorial.<jupyter_code>import diskmap
import glob
import matplotlib.pyplot as plt
import numpy as np
import urllib.request
from astropy.io import fits<jupyter_output><empty_output><jupyter_text>Next, we download the $J$ band image of LkCa 15 circumstellar disk.<jupyter_code>urllib.request.urlretrieve('https://home.strw.leidenuniv.nl/~stolker/diskmap/lkca15_irdis_qphi.fits',
                           'lkca15_irdis_qphi.fits')<jupyter_output><empty_output><jupyter_text>We can read the FITS file with `astropy`.<jupyter_code>image = fits.getdata('lkca15_irdis_qphi.fits')<jupyter_output><empty_output><jupyter_text>Let's have a look at the image. The central masked region contains NaNs.<jupyter_code>plt.imshow(image, origin='lower', vmin=np.nanmin(image), vmax=0.35*np.nanmax(image))<jupyter_output><empty_output><jupyter_text>Mapping of the disk surface We will now create an instance of `DiskMap` by providing the FITS filename, the pixel scale of the IRDIS detector (12.25 mas), the inclination (50 deg) and position angle (90 deg) of the disk, the distance (160 pc), and the image type (polarized flux)..The inclination convention is such that the near side is located on the right side of the image when using an inclination between 0 and 90 deg and a position angle of 0 deg. Therefore, with an position angle of 90 deg, the near side will be in upward direction of the image, as we will also see later in the `_radius.fits` file.<jupyter_code>mapping = diskmap.DiskMap(fitsfile='lkca15_irdis_qphi.fits',
                          pixscale=0.01225,
                          inclination=50.,
                          pos_angle=90.,
                          distance=160.,
                          image_type='polarized')<jupyter_output><empty_output><jupyter_text>The scattering surface of the disk is mapped with the `map_disk` method. Here, we provide a powerlaw funtion as the (approximate) shape of the disk surface for which we assume a constant opening angle: $h(r) = 0 + 0.05r^{1}$. The argument of `radius` specifies the sampling of the radii (100 points between 1 and 500 au).For running the deprojection later on, it is important that the outer radius of the `radius` parameter is larger than the field of view of the image. This may not be possible if the disk is strongly inclined and flaring. A mapping of the full field of view is not required for the $r^2$ scaling and phase function extraction. In that case, a smaller outer radius can be used, for example the actual outer radius of the disk. The radius and scattering angle output will contain NaNs beyond the outer radius.<jupyter_code>mapping.map_disk(power_law=(0., 0.05, 1.),
                 radius=(1., 500., 100))<jupyter_output><empty_output><jupyter_text>Radius and scattering angle The available output from the `DiskMap` methods are written by calling `write_output`. The argument of `filename` contains the prefix of the output files.<jupyter_code>mapping.write_output(filename='lkca15')<jupyter_output><empty_output><jupyter_text>Let's see which FITS files have been written.<jupyter_code>glob.glob('*.fits')<jupyter_output><empty_output><jupyter_text>For simplicity with the plots, we define half the field of view in arcseconds.<jupyter_code>size = mapping.pixscale * image.shape[0]/2<jupyter_output><empty_output><jupyter_text>The deprojected radius (in au) from the disk surface to the star is stored in the `_radius.fits` file. Let's plot the image from this FITS file.<jupyter_code>radius = fits.getdata('lkca15_radius.fits')
plt.imshow(radius, origin='lower', extent=[size, -size, -size, size])
plt.xlabel('RA offset (arcsec)', fontsize=14)
plt.ylabel('Dec offset (arcsec)', fontsize=14)
cb = plt.colorbar()
cb.set_label(label='Deprojected radius (au)', size=14)<jupyter_output><empty_output><jupyter_text>Similarly, the scattering angles on the disk surface are stored in the `_scat_angle.fits` file. The scattering angle is defined as 180 degrees minus the angle between the direction from the disk surface to the star and the direction from the disk surface to the observer.<jupyter_code>scat_angle = fits.getdata('lkca15_scat_angle.fits')
plt.imshow(scat_angle, origin='lower', extent=[size, -size, -size, size])
plt.xlabel('RA offset (arcsec)', fontsize=14)
plt.ylabel('Dec offset (arcsec)', fontsize=14)
cb = plt.colorbar()
cb.set_label(label='Scattering angle (deg)', size=14)<jupyter_output><empty_output><jupyter_text>Irradiation correction Now that we have the deprojected distance from each pixel to the star, we can compute the stellar irradiation corrected (i.e. $r^2$ scaled) image. We set a maximum radius of 100 au such that the flux at large separations, which only consists of noise, is not enhanced by the scaling.<jupyter_code>mapping.r2_scaling(r_max=100.)<jupyter_output><empty_output><jupyter_text>We run again the `write_output` method such that also the r$^2$ scaled image is stored as FITS file.<jupyter_code>mapping.write_output(filename='lkca15')<jupyter_output><empty_output><jupyter_text>Let's have a look at the r$^2$ scaled image. The dynamical range is smaller compared to the regular image which brings out the disk features more clearly.<jupyter_code>r2_scaled = fits.getdata('lkca15_r2_scaled.fits')
plt.imshow(r2_scaled, origin='lower', extent=[size, -size, -size, size])
plt.xlabel('RA offset (arcsec)', fontsize=14)
plt.ylabel('Dec offset (arcsec)', fontsize=14)
cb = plt.colorbar()
cb.set_label(label='Flux (r$^2$ ADU)', size=14)<jupyter_output><empty_output><jupyter_text>Disk deprojection Next, we will use the 3D mapping of the disk surface to deproject the image with the `deproject_disk` method. The deprojection corrects therefore both for the inclination (i.e. the disk midplane) and height of the disk surface (i.e. the powerlaw profile).<jupyter_code>mapping.deproject_disk()<jupyter_output><empty_output><jupyter_text>And we write again all available output files.<jupyter_code>mapping.write_output(filename='lkca15')<jupyter_output><empty_output><jupyter_text>The deprojected image is stored in the FITS file with the `_deprojected.fits` suffix. This image shows what the disk would look like at an inclination of 0 degrees. Let's have a look at the result.<jupyter_code>deprojected = fits.getdata('lkca15_deprojected.fits')
plt.imshow(deprojected, origin='lower', extent=[size, -size, -size, size],
           vmin=np.amin(deprojected), vmax=0.25*np.amax(deprojected))
plt.xlabel('RA offset (arcsec)', fontsize=14)
plt.ylabel('Dec offset (arcsec)', fontsize=14)
cb = plt.colorbar()
cb.set_label(label='Flux (ADU)', size=14)<jupyter_output><empty_output><jupyter_text>Estimated total intensity image By assuming a bell-shaped (i.e. Rayleigh-like curve with an adjustable peak value) degree of polarization and using the scattering angles from before, we use the `total_intensity` method to convert the r$^2$-scaled, polarized intensity image into an estimated total intensity image. This method should therefore only be used if the input image is a polarized light image (i.e. `image_type='polarized'`).In this example, we assume a maximum polarization of 100% at a scattering angle of 90 degrees, which is to be expected for aggregate-like dust grains with submicron-sized monomers.<jupyter_code>mapping.total_intensity(pol_max=1.)<jupyter_output><empty_output><jupyter_text>We write again all available output, which now also includes the r$^2$-scaled, total intensity image.<jupyter_code>mapping.write_output(filename='lkca15')<jupyter_output><empty_output><jupyter_text>Let's plot the total intensity image. The forward scattering by dust grains on the near/north side of the disk is visible in this image.<jupyter_code>total_intensity = fits.getdata('lkca15_total_intensity.fits')
plt.imshow(total_intensity, origin='lower', extent=[size, -size, -size, size])
plt.xlabel('RA offset (arcsec)', fontsize=14)
plt.ylabel('Dec offset (arcsec)', fontsize=14)
cb = plt.colorbar()
cb.set_label(label='Flux (r$^2$ ADU)', size=14)<jupyter_output><empty_output><jupyter_text>Scattering phase function As a last step, we extract the scattering phase function from, that is, the normalized flux as function of scattering angle. We use the `phase_function` method and select pixels between a deprojected distance of 80 and 100 au (i.e. along the bright ring in the r$^2$-scaled image) and caculate the average flux in 30 steps between 0 and 180 degrees.<jupyter_code>mapping.phase_function(radius=(80., 100.), n_phase=30)<jupyter_output><empty_output><jupyter_text>We run again the `write_output` method to store the phase function in a text file.<jupyter_code>mapping.write_output(filename='lkca15')<jupyter_output><empty_output><jupyter_text>We can read the extracted phase function with the `loadtxt` function of `numpy`. The second and third column of the data file contain the extracted phase function and error, which in this case is the polarized phase function. The fourth and fifth column contains an extimated total intensity phase function, which assumes that the degree of polarization is bell-shaped.In case the input image is a total intensity image (i.e. `image_type='total'`), the data file contains only the regular / total intensity phase function.<jupyter_code>angle, pol_flux, pol_error, total_flux, total_error = np.loadtxt('lkca15_phase_function.dat', unpack=True)<jupyter_output><empty_output><jupyter_text>Let's plot the polarized phase function that is extracted from the r$^2$-scaled image.<jupyter_code>plt.errorbar(angle, pol_flux, yerr=pol_error)
plt.xlabel('Scattering angle (deg)', fontsize=14)
plt.ylabel('Normalized polarized flux', fontsize=14)<jupyter_output><empty_output><jupyter_text>Finally, we plot the total intensity phase function, which shows the onset of a strong forward scattering peak and a more shallow backward scattering peak.<jupyter_code>plt.errorbar(angle, total_flux, yerr=total_error)
plt.xlabel('Scattering angle (deg)', fontsize=14)
plt.ylabel('Normalized total flux', fontsize=14)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09448295508498718, 'fraction_numerical': 0.018611717785585416, 'mean_word_length': 5.834523036988968, 'pattern_counts': {'":': 0, '<': 101, '<?xml version=': 0, '>': 101, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '25', 'dir': 'jupyter-structured-clean-dedup', 'id': '498367', 'license': 'MIT', 'n_tokens_mistral': 3384, 'n_tokens_neox': 3135, 'n_words': 1120, 'path': 'docs/tutorial.ipynb', 'repo_name': 'tomasstolker/diskmap', 'size': '10531'}