3 vs 4 sided boom on Arcus¶

In [1]:

from nbtemplate import display_header, get_path
display_header('3vs4sidedboom.ipynb', status='Review Arcus 2020')

Revision status: Review Arcus 2020¶

Last revision in version control:

Author: Hans Moritz Guenther hgunther@mit.edu
commit cec12fa383d7ab490123c6ca02cd0a13286196af
Date: Thu Oct 22 15:35:06 2020 -0400

This document is git version controlled. The repository is available at https://github.com/hamogu/arcus. See git commit log for full revision history.

Code was last run with:

MARXS ray-trace code version 1.2.dev747 (commit hash: 40c68619ea from 2020-05-27)
ARCUS python code version 0.2.dev186 (commit hash: 8c6308a4ec from 2020-10-19)
ARCUS CALDB version hash: b'6ee78f7' (commited on b'2020-10-19 ')

In [2]:

import os
import sys
import numpy as np
import astropy.units as u
import marxs

%matplotlib inline
import matplotlib.pyplot as plt

In [3]:

figureout = get_path('figures')

In [4]:

from marxs.simulator import Sequence
from arcus import Arcus, xyz2zxy
import arcus.boom
from arcus.defaults import DefaultPointing, DefaultSource

The Arcus focal length is too large to fit a fully expanded spacecraft into the launcher. Instead, Arcus will have an extensible boom that is coiled up at launch and will be extended in space. If the boom is large enough, all the photons will be inside, but if the boom diameter is smaller, photons from the outer SPOs will have to pass through the structure of the boom. This is true in particular for a three-sided boom. In this notebook, I study how the shape and size of the boom affects the detected signal.

Arcus is designed with a 4-sided boom that is large enough that almost all photons paths fit entirely inside of the boom. It is easy to see analytically that a 4-sided boom should be oriented parallel to edges of the SPO petals. Only a a small fraction of the outermost SPOs is larger than then boom and the photons have to pass from outside the boom to inside and some photons can be lost when they hit the elements the boom is made of. This is a relatively easy case, and is simulated first. As a design study, we then look at a 3-sided boom with different orientation and diameters.

The simulations assume that every photons that intersects any boom element is absorbed. Given that the diameters of the boom elements are large compared to the wavelengths of the photons, we ignore diffraction.

In [5]:

n_photons = 1e5

# define position and spectrum of source
energy = np.arange(0.25, 1.7, .01)
mysource = DefaultSource(energy={'energy': energy, 'flux': np.ones_like(energy)/energy**2})
jitterpointing = DefaultPointing()
fixedpointing = DefaultPointing(jitter=0. * u.rad)

class ArcusNoBoom(Arcus):
    def add_boom(self, conf):
        return []

instrum = ArcusNoBoom()
# Make an object that holds just the Front Assambly (FA)
instrum_FA = Sequence(elements=instrum.elements[:3])
instrum_FA.postprocess_steps = instrum.postprocess_steps
instrum_BA = Sequence(elements=instrum.elements[3:])

In [6]:

photons = mysource.generate_photons(n_photons)
photons = jitterpointing(photons)
photons = instrum_FA(photons)

In this section I use a single simulation of $10^5$ rays with a flat input spectrum ("flat" here means that bins of constant width in wavelength will give you a flat spectrum - see plot below). The details of the spectrum are not terribly important; the point is that the numbers below are not specific to a single energy, but represent the average over a range of energies. The analysis below is limited to the photons that actually make it onto a detector for the current configuration of 16 CCDs in two strips as above. Changing the detector layout would slightly change some of the numbers below.

I save the path that each rays takes in this simulation. Then, I look at different booms (different sizes and rotation angles) and check which of the photons would have interacted with each boom.

In [7]:

photons['wave'] = photons['energy'].to(u.Angstrom, equivalencies=u.spectral())
fig = plt.figure()
ax1 = fig.add_subplot(111)
out = ax1.hist(photons['wave'], bins=100)
ax1.set_xlabel('wavelength [$\AA$]')
ax1.set_ylabel('Flux (with random errors\ndue to Poisson sampling)')
out = ax1.set_title('Spectrum of photons used for the following simulations')

The 4 sided boom¶

Arcus is designed with a 4-sided boom and relatively easy, so we study this first.

Rotation angle of boom¶

For the four sided boom it is a natural choice to orient it such that the battens are parallel to the $x$ and $y$ axis of the instrument coordinate system. For a relatively large boom, it can be shown analytically that this configuration gives the smallest obscuration. As below for the three-sided boom, the larger the diameter, the smaller the loss of photons will be.

In [8]:

photons = mysource.generate_photons(n_photons)
photons = jitterpointing(photons)
instrum_FA.postprocess_steps[0].data = []
photons = instrum_FA(photons)

fourboom = arcus.boom.FourSidedBoom(orientation=xyz2zxy[:3, :3],
                                   position=[0, 0, 596.])
hit_4 = fourboom(photons.copy())['hitrod']

photons = instrum_BA(photons)
good = (photons['CCD'] > 0) & (photons['probability'] > 0)

In [9]:

pos = instrum_FA.postprocess_steps[0].data[0]
ind =  (photons['CCD'] > 0)
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(111, aspect='equal')
#plt.plot(pos[:, 0], pos[:, 1] , '.')
ax.plot(pos[ind, 0], pos[ind, 1], '.', label='detected on CCD')
ax.plot(pos[ind & hit_4, 0], pos[ind & hit_4, 1],
        'r.', label='hitting boom')
ax.legend(loc='center')

Out[9]:

<matplotlib.legend.Legend at 0x7f77847d8850>

In [10]:

pos = instrum.KeepPos.data[0]
ind =  (photons['CCD'] > 0) & (photons['aperture'] == 0)
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(111, aspect='equal')
ax.plot(pos[ind, 0], pos[ind, 1], '.', label='Photon detected on CCD')
ax.plot(pos[ind & hit_4, 0], pos[ind & hit_4, 1], 'ro', label='Photon hitting boom')
ax.legend(loc='center', numpoints=1)
ax.set_xlim([-530, -80])
ax.set_ylim([270, 850])
ax.set_xlabel('Position in aperture plane [mm]')
ax.set_ylabel('Position in aperture plane [mm]')
ax.set_title('How many photons hit the boom?\n(here: 1.85 m boom diameter)')
fig.savefig(os.path.join(get_path('figures'), '4sidedboom.png'), dpi=600, bbox_inches='tight') 

The figures above show that even with a boom diameter of 1.85 m there are still a few photons which are absorbed by the boom. The inner rows of SPOs are on the inside of the boom and do not suffer from obscuration, but the two outermost rows that are located at a radius larger than the boom loose some of the photons that pass though the them.

In [11]:

ind = (photons['CCD'] >=-6)
ind0 = ind & (photons['order'] == 0)
inddisp = ind & (photons['order'] != 0)
loss_0 = (photons['probability'][ind0 & hit_4]).sum() / (photons['probability'][ind0]).sum()
loss_disp = (photons['probability'][inddisp & hit_4]).sum() / (photons['probability'][inddisp]).sum() 
print(f'Fraction of zero-order photons lost to boom: {loss_0:6.3f}')
print(f'Fraction of diffracted photons lost to boom: {loss_disp:6.3f}')

Fraction of zero-order photons lost to boom:  0.008
Fraction of diffracted photons lost to boom:  0.014

In [12]:

fig, ax = plt.subplots()
orderok = np.isfinite(photons['order'])
out = ax.hist(photons['order'][hit_4 & orderok], 
         bins=np.arange(np.min(photons['order'][orderok]) - 0.5,
                          np.max(photons['order'][orderok]) + 0.5),
         weights=photons['probability'][hit_4 & orderok])
ax.set_xlabel('Diffraction order')
out = ax.set_ylabel('Number of photons lost')

This figure shows how many photons are lost in each order. Since each photon in the simulation has a survival probability in the range 0..1, the sum of the survival probabilities is not an integer. This figure obviously depends strongly on the spectrum of the source, which determines how many photons are seen in each order but the main point is that the boom has an impact on essentially all orders.

The 3 sided boom¶

Next, we will look at the design of a three-sided boom.

Geometry¶

In [13]:

angles = np.arange(0, 2. * np.pi/3., 0.1)
l_batten = np.array([1.6, 1.8, 2.0, 2.25, 2.5, 2.75]) * 1e3
colors = ['k', 'b', 'r', 'g', 'm', '0.6']

In [14]:

import matplotlib.patches as patches
from transforms3d.affines import decompose

fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(111, aspect='equal')

from arcus.arcus import Aperture, DetCamera, defaultconf

for e in Aperture(defaultconf).elements:
    trans, rot, zoom, shear = decompose(e.pos4d)
    aper = patches.Rectangle(trans[:2] - zoom[1:], zoom[1] * 2, zoom[2] * 2, 
                             linewidth=1, edgecolor='0.7', facecolor='0.7', label='SPO')
    ax.add_patch(aper)
for e in DetCamera(defaultconf).elements:
    trans, rot, zoom, shear = decompose(e.pos4d)
    det = patches.Rectangle(trans[:2] - zoom[1:], zoom[1] * 2, zoom[2] * 2, 
                             linewidth=1, edgecolor='orange', facecolor='orange', label='detector')
    ax.add_patch(det)
    
for i in range(len(l_batten)):
    for a in angles[::7]:
        triag = patches.RegularPolygon((0, 0), 3, l_batten[i] * 3**(-0.5), a, 
                                       edgecolor=colors[i], linewidth=1, facecolor='none',
                                       label='boom', lw=1.5)
        ax.add_patch(triag)
for i, chanpos in enumerate(defaultconf['pos_opt_ax'].values()):
    optax = plt.plot(chanpos[0], chanpos[1], 'bp', label='opt. axis' if i == 0 else '__no_legend__')
ax.set_ylim([-1.5e3, 1.5e3])
ax.set_xlabel('distance [mm]')
out = ax.set_ylabel('distance [mm]')
fig.subplots_adjust(left=0.17, right=.95, top = .98, bottom=.05)
fig.savefig(figureout + 'boom_scematic.png')
fig.savefig(figureout + 'boom_scematic.pdf')

The plot above shows a view along the optical axis. The gray boxes are the entrance apertures that are filled with SPOs, the small orange rectangles along the y=0 axis are the 16 CCDs (set in two groups of 8). Remember that in the current layout Arcus has four of SPO channels, where each channel has its own optical axis. The channels are paired and the optical axes in each pair hit the same CCD, thus the symbols overlap on this plot. This allows us to place the CCD detectors in such a way that we catch the zeroths order for each pair as well as the dispersed spectrum in the region where the CAT gratings cast most of the signal.

The left pair of SPO modules focusses on the left blue pentagon, so zeroth order rays from the left pair of SPOs will be located there. The dispersed spectrum is imaged on the right detector strip. Conversely, the right pair of SPO channels places its zeroth order on the right blue pentagon and the dispersed photons on the left detector strip.

Overlayed are several triangles marking the footprint of a three-sided boom with different dimensions. The colors are the same as in the figures in the next section, where the exact dimensions are given. There are many ways to overlay a triangular boom on the rectangular geometry of the SPO channels. With changing rotation angle, rays going through different parts of the SPOs will have a chance of hitting the boom. The figure shows three possible configurations, but in the simulations below I use a much denser grid of rotation angles with a step size of just a few degrees.

How does absorption depend on boom size and boom rotation angle?¶

In [15]:

def boom_abs(photons, angles, boom_dims={}):
    hitrod = []

    for i, angle in enumerate(angles):
        rot = np.array([[np.cos(angle), -np.sin(angle), 0],
                        [np.sin(angle), np.cos(angle), 0],
                        [0, 0, 1]])
        myboom = arcus.boom.ThreeSidedBoom(orientation=np.dot(rot, xyz2zxy[:3, :3]), boom_dimensions=boom_dims)
        photons['hitrod'] = False
        photons = myboom(photons)
        hitrod.append(photons['hitrod'].copy())
    return hitrod

In [16]:

boom_hits = []
for l in l_batten:
    print('Working on l={}'.format(l))
    boom_hits.append(boom_abs(photons.copy(), angles, boom_dims={'l_batten': l}))

Working on l=1600.0
Working on l=1800.0
Working on l=2000.0
Working on l=2250.0
Working on l=2500.0
Working on l=2750.0

In [17]:

def calc_prob_matrix(ind):
    prob = np.empty((len(angles), len(l_batten)), dtype=float)
    for i in range(len(l_batten)):
        for j in range(len(angles)):
            prob[j, i] = photons['probability'][ind & boom_hits[i][j]].sum()
    return prob / photons['probability'][ind].sum()

In [18]:

photons = instrum_BA(photons)

ind = (photons['CCD'] >=-6)
ind0 = ind & (photons['order'] == 0)
prob0 = calc_prob_matrix(ind0)

inddisp = ind & (photons['order'] != 0)
probdisp = calc_prob_matrix(inddisp)
import matplotlib as mpl
from cycler import cycler
with mpl.style.context('fivethirtyeight'):
    #plt.rc('axes', prop_cycle=(cycler('linestyle', ['-', '--', ':']) * 
    #                           cycler('color', ['r', 'b', 'y', 'c', 'g', 'm', '#7A68A6', '#188487', '#E24A33'])))

    fig = plt.figure(figsize=(7,5))
    ax = fig.add_subplot(111)

    for i, l in enumerate(l_batten):
        ax.plot(np.rad2deg(angles), prob0[:, i], label='{:4.2f} m'.format(l / 1e3), color=colors[i], lw=2)
        ax.plot(np.rad2deg(angles), probdisp[:, i], label='__no_legend__', color=colors[i], ls=':', lw=3)
    ax.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
    ax.set_xlabel('rotation angle of boom [deg]')
    ax.set_title('Photon loss for dispersed (dotted line)\nand 0th order (solid line) photons')
    out = ax.set_ylabel('relative photon loss\n(averaged over all orders)')
    fig.savefig(figureout + 'boom_loss.png', bbox_inches='tight')
    fig.savefig(figureout + 'boom_loss.pdf', bbox_inches='tight')

The plot above shows which fraction of photons is lost due to hitting elements of the boom. Lines of different color are for different length of the battens (the boom structures running parallel to the focal plane). Solid lines are for the dispersed spectrum, dotted lines are for the zero order photons. For a 3-sided boom, a rotation of 120$^\circ$ maps the boom on itself, thus only the range up to 120$^\circ$ is shown. The dimension given by the color is the length of one side of the triangle, i.e. the length of a "batten".

There are noticeable differences with the boom rotation angle particularly for intermediate boom sizes. For a small boom (black) essentially the entire circumfence of the boom is under the SPOs (see plot of the geometry above), while for very large booms only the outer corners of the SPO channel are touched by the boom. Thus, only a small fraction of photons is absorbed independent of the rotation angle. For intermediate boom sizes with triangle edges around 2 m a clever choice or rotation angle can change the fraction of absorbed photons from 13% to 9%!

Energy dependence¶

With the caveat in mind that the current mirror implementation does not simulate the angle dependence of the reflection probability in the SPOs and thus ignores that SPOs at different radius actually have a slightly different importance for different photon energies, we will now look at the photon losses as a function of photon energy. In the following, we will pick a particular orientation angle for the boom since the shape of the following plots differs only little with the rotation angle. Also, we will concentrate on the dispsered signal.

In [19]:

photons['wave'] = marxs.energy2wave / photons['energy'] * 1e7 # in ang

In [20]:

spec = plt.hist(photons['wave'][inddisp], 
                weights=photons['probability'][inddisp], bins=20,
                label='no boom', color='y', lw=0)
specabs = []
for i, l in enumerate(l_batten):
    indboom = boom_hits[i][0]
    specabs.append(plt.hist(photons['wave'][inddisp & ~indboom], 
                    weights=photons['probability'][inddisp & ~indboom], 
                       bins=spec[1],
                      histtype='step', color=colors[i],
                      label='{:4.2f} m'.format(l / 1e3)))

plt.xlabel('wavelength [$\AA$]')
plt.ylabel('signal')
plt.title('Change of effective area due to the boom absorption')
out = plt.legend(loc='upper right')

The plot above gives the observed grating signal (with energy on the x-axis in this case). The input spectrum is flat, but the effective area has an energy dependence due to e.g. energy dependend reflection probabilities on the mirror, grating efficiency, filter curves of the optical blocking filter, detector gaps in the focal plane, CCD QE etc. This plot shows us the loss of effective area due to the boom. Below, the same data is plotted again, but this time relative to the signal that would be detected in the absence of a boom.

In [21]:

outn = np.histogram(photons['wave'][inddisp], bins=20)
outnabs = np.histogram(photons['wave'][inddisp & ~indboom], bins=outn[1])

spec = np.histogram(photons['wave'][inddisp], 
                weights=photons['probability'][inddisp], bins=outn[1])
specabs = []

with mpl.style.context('fivethirtyeight'):
    fig = plt.figure(figsize=(7, 5))
    ax = fig.add_subplot(111)
    for i, l in enumerate(l_batten):
        indboom = boom_hits[i][0]
        pboom = photons.copy()
        pboom['probability'] = pboom['probability'] * ~indboom
        specabs.append(np.histogram(pboom['wave'][inddisp], 
                       weights=pboom['probability'][inddisp], 
                       bins=outn[1]))
        ax.plot(0.5 * (outn[1][:-1] + outn[1][1:]), 
                specabs[i][0] / spec[0],
                c=colors[i], lw=2,
                label='{:4.2f} m'.format(l / 1e3))

    plt.xlabel('wavelength [$\AA$]')
    plt.ylabel('relative effective area')
    plt.title('Change of $A_\mathrm{eff}$ due to the boom absorption\n(for the dispersed spectrum)')
    out = plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
    fig.savefig(figureout + 'boom_aeff.png', bbox_inches='tight')
    fig.savefig(figureout + 'boom_aeff.pdf', bbox_inches='tight')

The plot above shows the change in effective area for booms of different size (for a fixed orientation angle). Particularly for small booms, the short wavelength end of the spectrum is hit hardest.

Last, we want to know if the absorption is smooth over the full wavelength range or if it produces line-like features that would be confused with an astrophysical signal. In particular, line-like features wold move with any boom movement or misalignment, so they should be avoided. While the entire wavelength range should be checked with simulations using a larger number of rays eventually, I will just select some small regions here and zoom in to bin them to the Arcus spectral resolving power. Also, I will only show the plot below for one particular boom (rotation = 0, fiducial size 1.6 m).

In [22]:

def simulate_one_boom(energyrange):
    mysource = DefaultSource(energy={'energy': energyrange,
                                     'flux': np.ones(2)})
    myboom = arcus.boom.ThreeSidedBoom(orientation=xyz2zxy[:3, :3])

    photonj = mysource.generate_photons(n_photons)
    photonj = jitterpointing(photonj)
    photonj = instrum(photonj)
    photonj = myboom(photonj)
    return photonj

In [23]:

photonreg = [[0.62, 0.65], [0.8,0.802], [1., 1.02]]
photonsinreg = [simulate_one_boom(enrange) for enrange in photonreg]

In [24]:

def plot_change_in_aeff(photonj, ax):
    photonj['wave'] = marxs.energy2wave / photonj['energy'] * 1e7 # in ang
    indj = np.isfinite(photonj['det_x'])

    for i in [-6, -7, -8]:
        indorder = photonj['order'] == i
        spec = np.histogram(photonj['wave'][indj & indorder], bins=20, 
                            weights=photonj['probability'][indj & indorder])
        specabs = np.histogram(photonj['wave'][indj & indorder & ~indboom], 
                           weights=photonj['probability'][indj & indorder & ~indboom], 
                           bins=spec[1])
        ax.plot(0.5 * (spec[1][:-1] + spec[1][1:]), specabs[0] / spec[0], label='order = {}'.format(int(np.abs(i))))


    ax.set_xlabel('wavelength [$\AA$]')
    ax.set_ylabel('relative effective area')
    out = ax.legend(loc='lower left')

In the dispersed orders, the change in effective area looks relatively smooth. The wiggeling up and down in the plot is consistent with Poisson statistics. If we decide to continue to investigate this in more detail, I will run longer simulations with a larger number of photons, but for now, there are no indications that line-like features appear in the dispersed orders.

In [25]:

fig = plt.figure()
ax = fig.add_subplot(111)
plot_change_in_aeff(photonsinreg[1], ax)