Demo of the use of satellite spots

Satellite spots are patterns put on the deformable mirror (DM1) and are used in post-processing to find the center of the star. In this tutorial, we will show how to apply them, the possible parameters and the different modes

#import necessary packages
from corgisim import scene
from corgisim import instrument
import matplotlib.pyplot as plt
import numpy as np
import proper
import roman_preflight_proper
roman_preflight_proper.copy_here()

WARNING: AstropyDeprecationWarning: The TestRunner class is deprecated and may be removed in a future version.
        Use pytest instead. [astropy.tests.runner]
WARNING: AstropyDeprecationWarning: The TestRunnerBase class is deprecated and may be removed in a future version.
        Use pytest instead. [astropy.utils.decorators]
Matplotlib is building the font cache; this may take a moment.

corgidrp: Configuration file written to /home/docs/.corgidrp/corgidrp.cfg. Please edit if you want things stored in different locations.

/home/docs/checkouts/readthedocs.org/user_builds/corgisim/envs/latest/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Let’s start by simulating an image without satellite spots

#define host star properties
Vmag = 8
sptype = 'G0V'
cgi_mode = 'excam'
host_star_properties = {'Vmag': Vmag, 'spectral_type': sptype, 'magtype':'vegamag'}

#initialize scene
base_scene = scene.Scene(host_star_properties)

#define coronagraph properties
cgi_mode = 'excam'
cor_type = 'hlc'
bandpass = '1'
cases = ['2e-9']
rootname = 'hlc_ni_' + cases[0]
dm1 = proper.prop_fits_read( roman_preflight_proper.lib_dir + '/examples/'+rootname+'_dm1_v.fits' )
dm2 = proper.prop_fits_read( roman_preflight_proper.lib_dir + '/examples/'+rootname+'_dm2_v.fits' )

##  Define the polaxis parameter. Use 10 for non-polaxis cases only, as other options are not yet implemented.
polaxis = 10
# output_dim define the size of the output image
output_dim = 101

optics_keywords = {'cor_type':cor_type, 'use_errors':1, 'polaxis':10, 'output_dim':output_dim,\
                    'use_dm1':1, 'dm1_v':dm1, 'use_dm2':1, 'dm2_v':dm2,'use_fpm':1, 'use_lyot_stop':1,  'use_field_stop':1 }

optics = instrument.CorgiOptics(cgi_mode, bandpass, optics_keywords=optics_keywords, if_quiet=True)

#generate host star psf
sim_scene = optics.get_host_star_psf(base_scene)
image_star_corgi = sim_scene.host_star_image.data

#plot psf
plt.imshow(image_star_corgi, origin='lower')
plt.title('Host star without satellite spots')
co = plt.colorbar()

CorgiOptics initialized with proper keywords.

Now let’s define our satellite spots. To do this, you need to choose:

• The number of pairs of satellite spots (1 or 2)

• The separation from the center in lambda/D

• The angle position (0 and 90 degrees or 45 and 135 degrees)

• The desired contrast

• The wavelength (the default value is the central wavelength of the bandpass of your optics)

• The sign (positive or negative, default is positive)

You can input any combination of the values for the parameters, but this image lists the expected combinations.

Let’s put satellite spots on our optics. We can either create a new object or add them to the existing optics.

satspot_keywords = {'num_pairs':2, 'sep_lamD': 6.25, 'angle_deg': [0,90], 'contrast': 1e-6}

#Adding them to the existing object
optics.add_satspot(satspot_keywords=satspot_keywords)

fig = plt.figure(figsize=(12,3))
plt.subplot(133)
plt.imshow(optics.optics_keywords['dm1_v'] - dm1)
plt.colorbar()
plt.title('Added Pattern')

plt.subplot(132)
plt.imshow(optics.optics_keywords['dm1_v'])
plt.colorbar()
plt.title('Satellite spots added')
plt.subplot(131)
plt.imshow(dm1)
plt.colorbar()
plt.title('Initial state')

fig.suptitle('Deformable Mirrors')
fig.subplots_adjust(top=0.8)

Let’s see what an on-axis star looks like after going through the optics with the satellite spots.

sim_scene_with_spots = optics.get_host_star_psf(base_scene)
image_star_with_spots = sim_scene_with_spots.host_star_image.data

#plot psf
plt.imshow(image_star_with_spots, origin='lower')
plt.title('Star with satellite spots')
co = plt.colorbar()

Let’s see what happens if we introduce an offset.

wavelength = 0.575e-6    #band1
lam_D = np.degrees(wavelength/2.3)*3600*1000 # in mas
shift = [2*lam_D,0] # shift in [x,y]
optics_keywords ={'cor_type':cor_type, 'use_errors':1, 'polaxis':polaxis, 'output_dim':output_dim,\
                'use_dm1':1, 'dm1_v':dm1, 'use_dm2':1, 'dm2_v':dm2,'use_fpm':1, 'use_lyot_stop':1,  'use_field_stop':1,\
                'source_x_offset_mas': shift[0], 'source_y_offset_mas': shift[1]}

optics_with_shift = instrument.CorgiOptics(cgi_mode, bandpass, optics_keywords=optics_keywords,if_quiet=True)

# Let's take one image with the satellite spots and one without
sim_scene = optics_with_shift.get_host_star_psf(base_scene)
image_star = sim_scene.host_star_image.data

# Add positive satspots
optics_with_shift.add_satspot(satspot_keywords=satspot_keywords)
sim_scene_with_pos_spots = optics_with_shift.get_host_star_psf(base_scene)
image_star_with_pos_spots = sim_scene_with_pos_spots.host_star_image.data
fig = plt.figure(figsize=(12,4))
plt.subplot(121)
plt.imshow(image_star)
plt.colorbar()
plt.title('Image Star Without Satellite Spots')
plt.subplot(122)
plt.imshow(image_star_with_pos_spots)
plt.colorbar()
plt.title('Image Star With Satellite Spots')

CorgiOptics initialized with proper keywords.

Text(0.5, 1.0, 'Image Star With Satellite Spots')

As you can see, the light of the star is too bright to see the satellite spots if the star is misaligned. We can substract the star image without satellite spots…

plt.imshow(image_star_with_pos_spots - image_star)
plt.colorbar()
plt.title('Image Satellite Spots minus Image Star')

Text(0.5, 1.0, 'Image Satellite Spots minus Image Star')

… But the resulting image isn’t quite good enough to get the position of the satellite spots. We need to take another image with negative satellite spots, then take the median image of our satellite spots images, and substract the star image from that (This is a simplified version of what the DRP does). Let’s do just that!

# First, we remove the positive satellilte spots
optics_with_shift.remove_satspot(satspot_keywords=satspot_keywords)

# Then, add the negative satellite spots
satspot_keywords['sign'] = "negative"
optics_with_shift.add_satspot(satspot_keywords=satspot_keywords)

sim_scene_with_neg_spots = optics_with_shift.get_host_star_psf(base_scene)
image_star_with_neg_spots = sim_scene_with_neg_spots.host_star_image.data

#Take the median and remove the star image
image_med = np.median(np.stack([image_star_with_pos_spots,image_star_with_neg_spots]), axis = 0)
image_spots = image_med - image_star

plt.imshow(image_spots)
plt.colorbar()
plt.title('Image Satellite Spots')

Text(0.5, 1.0, 'Image Satellite Spots')

Now we can check that the center of the satellite spots is at the location of the star.

from pyklip.fakes import gaussfit2d

xcen = ycen =(output_dim-1)/2

pix_scale = 0.0218*1000 #mas/pix
star_x = xcen+shift[0]/pix_scale
star_y = ycen+shift[1]/pix_scale

fitresult=[]
for angle in satspot_keywords['angle_deg']:
    sep = satspot_keywords['sep_lamD']*lam_D # mas

    guess_x1 = star_x+sep*np.cos(np.radians(angle))/pix_scale
    guess_x2 = star_x-sep*np.cos(np.radians(angle))/pix_scale

    guess_y1 = star_y+sep*np.sin(np.radians(angle))/pix_scale
    guess_y2 = star_y-sep*np.sin(np.radians(angle))/pix_scale

    for guess_x,guess_y in zip([guess_x1,guess_x2],[guess_y1,guess_y2]):
        peak,fwhm,x,y=gaussfit2d(image_spots,guess_x,guess_y)
        fitresult.append([peak,fwhm,x,y])

#### check location, averaged coordinates the four satellite spots are located at the FoV center
guessx_star = np.mean(np.array(fitresult)[:,2])
guessy_star = np.mean(np.array(fitresult)[:,3])

#Distance between the actual star center and the measured star center
dist = np.sqrt((guessx_star-star_x)**2+(guessy_star-star_y)**2)

print('The center of the satellite spots is measured at', dist, 'pixels from the star')

The center of the satellite spots is measured at 0.24109886561858315 pixels from the star

Wide Field of View

You can also use add satellite spots wen using a shaped-pupil coronograph

bandpass= '4'
cor_type = 'spc-wide'
cases = ['2e-8']
rootname = 'spc-wide_ni_' + cases[0]
dm1 = proper.prop_fits_read( roman_preflight_proper.lib_dir + '/examples/'+rootname+'_dm1_v.fits' )
dm2 = proper.prop_fits_read( roman_preflight_proper.lib_dir + '/examples/'+rootname+'_dm2_v.fits' )

optics_keywords = {'cor_type':cor_type, 'use_errors':1, 'polaxis':10, 'output_dim':201,\
                    'use_dm1':1, 'dm1_v':dm1, 'use_dm2':1, 'dm2_v':dm2,'use_fpm':1, 'use_lyot_stop':1,  'use_field_stop':1 }
#For the wide field of view, the separation of the satellite spot is different. We aslo use the 45/135 angles to illustrate their use
satspot_keywords = {'num_pairs':2, 'sep_lamD': 13, 'angle_deg': [45,135], 'contrast': 1e-6}

optics = instrument.CorgiOptics(cgi_mode, bandpass, optics_keywords=optics_keywords, satspot_keywords=satspot_keywords, if_quiet=True)

#generate host star psf
sim_scene_spc = optics.get_host_star_psf(base_scene)
image_star_spc = sim_scene_spc.host_star_image.data

plt.imshow(image_star_spc)
plt.colorbar()
plt.title('Image Satellite Spots SPC')

CorgiOptics initialized with proper keywords.

Text(0.5, 1.0, 'Image Satellite Spots SPC')

Polarimetry Mode

You can also add them in polarimetry mode (beware that this utilisation has not been formally tested)

bandpass = '1'
cor_type = 'hlc'
rootname = 'hlc_ni_2e-9'
#define which wollaston prism to use
prism = 'POL0'
#initialize instrument and generate 0/90 PSF pair
optics_keywords = {'cor_type':cor_type, 'use_errors':1, 'polaxis':10, 'output_dim':output_dim, 'prism':prism,\
                    'use_dm1':1, 'dm1_v':dm1, 'use_dm2':1, 'dm2_v':dm2,'use_fpm':1, 'use_lyot_stop':1,  'use_field_stop':1 }
satspot_keywords = {'num_pairs':2, 'sep_lamD': 6.25, 'angle_deg': [0,90], 'contrast': 1e-5}

optics = instrument.CorgiOptics(cgi_mode, bandpass, optics_keywords=optics_keywords,satspot_keywords=satspot_keywords, if_quiet=True)
sim_scene = optics.get_host_star_psf(base_scene)
image_star_horizontal = sim_scene.host_star_image.data[0]
image_star_vertical = sim_scene.host_star_image.data[1]

optics.remove_satspot(satspot_keywords=satspot_keywords)
sim_scene = optics.get_host_star_psf(base_scene)
image_star_horizontal_without_satellite_spots = sim_scene.host_star_image.data[0]
image_star_vertical_without_satellite_spot = sim_scene.host_star_image.data[1]

CorgiOptics initialized with proper keywords.

fig = plt.figure(figsize=(12,8))
plt.subplot(221)
plt.imshow(image_star_horizontal)
plt.colorbar()
plt.title('0 degree polarized light intensity \n Star and satellite spots')
plt.subplot(222)
plt.imshow(image_star_vertical)
plt.colorbar()
plt.title('90 degree polarized light intensity \n Star and satellite spots')
plt.subplot(223)
plt.imshow(image_star_horizontal - image_star_horizontal_without_satellite_spots)
plt.colorbar()
plt.title('0 degree polarized light intensity \n Star substracted')
plt.subplot(224)
plt.imshow(image_star_vertical - image_star_vertical_without_satellite_spot)
plt.colorbar()
plt.title('90 degree polarized light intensity \n Star substracted')

fig.suptitle('Polarimetry mode')
fig.subplots_adjust(top=0.88)

Spectroscopy mode

In order to keep the execution time of this notebook short, the demonstration of satellite spots in spectroscopy mode is at the end of spec_slit_prism_demo.