In [1]:

from nbtemplate import display_header, get_path
display_header('SelectRowlandParameters.ipynb', status='reviewed for SV')

Revision status: reviewed for SV¶

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 ')

Goal¶

The goal of this study is to optimize the parameters of the Rowland torus. There are several pieces of information that constrain our choice of parameters for the Rowland torus, most importantly engineering constraints that limit the area available to place the SPO channels in the focal plane. Moving the channels too close together will lead to overlapping mounting structures, moving them too far apart requires a larger spacecraft boom.

Here, I present ray-traces for Arcus that explore a reasonable parameter space in the following parameters:

channel spacing (measured center to center and given in mm)
blaze angle (measured in degrees)
Radius of Rowland Torus (measured in mm)

The remaining parameter is the distance between an on-axis grating and the focal point. Here, I have chosen a fixed value of 11800 mm. This will place the actual Arcus gratings (which are not located on-axis) a little more than 200 mm from the nodal plane of the mirrors. Adjusting this by a few cm has very little influence on the instrument performance and thus is can be fine-tuned after the remaining parameters are chosen.

Caveats¶

The analysis presented here is based on full ray-trace calculations but there are a few simplifications for speed and ease of set-up compared to the full Arcus instrument. I list the most important ones here:

No misalignments: The simulations below assume an ideal placement of SPOs and CAT gratings with zero mis-alignments. Typically, misalignments impact the spectral resolving power more than the effective area so it seems wise to chose a configuration that performs better than the mission baseline requirements to leave room for an error budget.
No chip gaps: The calculations of the effective area assume two continuous strips of CCDs. In reality, there will be ~ 2 mm gaps between the CCDs, leading to drops in the effective area. This can be somewhat mitigated by shifting one of the two CCDs strips by 3-5 mm compared to the position simulated here. This will have a very small impact on the effective area, but it avoids the problem that an important wavelength might fall into a chip gap in all channels at the same time. The other strategy that can be employed is dithering, which smears out this gap over a larger wavelength range. Still, the final effective area will a little lower then the curves shown here (probably ~10 % for narrow regions of the spectrum or a few % for wider regions with dithering).
CCDs are curved: In these simulations I employ a cylindrical detector that follows the Rowland circle exactly. The real Arcus detector will consist of flat CCDs that are tiled and deviate slightly from the Rowland circle, thus reducing the spectral resolving power a little. This is a very small effect.

Channel spacing¶

The grid of models that I have run includes a range of channel spacings, but from an engineering point of view a small distance between channels is strongly preferred, because it allows the use of a smaller boom. I did some initial exploration of this, but it seems that we can match the mission baseline requirements with a channel spacing that allows us to use a relatively small boom, so I will concentrate on that in the following and almost all plots are done for a channel spacing of 600 mm (measured center-to-center).

In [2]:

import os

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from astropy.table import Table

from marxs.analysis import gratings as anagrat
from arcus import analyze_design as ad

%matplotlib inline

In [3]:

filepath = get_path('grid2designtorus')
sumtab = Table.read(os.path.join(filepath, 'summary.fits'), hdu=1)
wave = Table.read(os.path.join(filepath, 'summary.fits'), hdu=2)['wave']
orders = Table.read(os.path.join(filepath, 'summary.fits'), hdu=3)['col0']
sumtab['aeff_per_order'] = sumtab['aeff'][..., None] * (sumtab['prob_per_order'] / 
                                                        sumtab['prob_per_order'].sum(axis=2)[..., None])

In [4]:

def res_aeff_in_region(tab, wave, wavelim, min_res):
    '''Summarize Aeff and resolving power in a given spectral region because this is the way the requirements
    in the Science Traceability Matrix are defined.
    
    This function identifies all orders that have a spectral resolving power > min_res and then adds up the 
    effective areas over all orders and averages them over the wavelength range.
    The spectral resolving power over all orders for a given wavelength is averaged weighted by the appropriate
    effective area and then averaged over all wavelengths in the selected region.
    '''
    indwave = (wave >= wavelim[0]) & (wave <= wavelim[1])
    indres = np.isfinite(tab['res_per_order']) & (tab['res_per_order'] > min_res)
    res_per_order = np.ma.array(sumtab['res_per_order'], mask= ~indwave[None, :, None] | ~indres)
    aeff_per_order = np.ma.array(sumtab['aeff_per_order'], mask= ~indwave[None, :, None] | ~indres)
    aeff = aeff_per_order.sum(axis=2).mean(axis=1)
    res = ((res_per_order * aeff_per_order).sum(axis=2) / aeff_per_order.sum(axis=2)).mean(axis=1)
    return res, aeff

In [5]:

# Note that the traceability matrix is not very ordered.
# Sometimes the Aeff required given and then the resolving power and sometimes the other way around.
res, aeff = res_aeff_in_region(sumtab, wave, [21.1, 28.], 100.)
sumtab['G1-1 (a)'] = res
sumtab['G1-1 (b)'] = aeff

res, aeff = res_aeff_in_region(sumtab, wave, [33.7, 40.], 100.)
sumtab['G1-1 (d)'] = res
sumtab['G1-1 (c)'] = aeff

res, aeff = res_aeff_in_region(sumtab, wave, [16., 21.6], 2000.)
sumtab['G1-2 (b)'] = res
sumtab['G1-2 (a)'] = aeff

res, aeff = res_aeff_in_region(sumtab, wave, [12., 50.], 1500.)
sumtab['G1-2 (b)'] = res
sumtab['G1-2 (c)'] = aeff

/nfs/melkor/d1/guenther/soft/anaconda/envs/mayavi46/lib/python3.8/site-packages/astropy/table/column.py:1020: RuntimeWarning: invalid value encountered in greater
  result = getattr(super(), op)(other)

In [6]:

def plot_phi_hist(ax, p):
    p['order_key'] = '     '
    p['order_key'][p['order'] > 0] = '> 0'
    for o in np.arange(-9, 1):
         p['order_key'][p['order'] == o] = str(o)
    p['order_key'][p['order'] < -9] = '< -9'
    p_order = p['order_key', 'phi_folded', 'probability'].group_by('order_key')
    bins = np.linspace(0, .19, 41)
    sl = slice(-1, 2, -1)
    phis = [g['phi_folded'] for g in p_order.groups]
    probs = [g['probability'] for g in p_order.groups]
    ax.hist(phis[sl], weights=probs[sl], bins=bins, stacked=True, 
            label=[f'{k[0]}' for k in p_order.groups.keys[sl]], rwidth=1, ec='none')
    ax.legend(title='Order')
    ax.set_xlabel('dispersion angle [rad]')
    ax.set_ylabel('# photons')

In [7]:

def plot_aeff(ax, row, wave, order0, **kwargs):
    label = kwargs.pop('label', '__no_legend__')
    ax.plot(wave, row['aeff_per_order'][:, ~order0].sum(axis=1), label=label, **kwargs)
    ax.plot(wave, row['aeff_per_order'][:, order0], ':', label='__no_legend__', **kwargs)

In [8]:

from matplotlib import transforms

def plot_everything(p, rows, wave, orders, color='kbgrymcy'):
    
    order0 = (orders == 0)
    
    fig = plt.figure(figsize=(18, 4))
    ax1 = fig.add_subplot(141)
    trans1 = transforms.blended_transform_factory(ax1.transData, ax1.transAxes)
    
    ax2 = fig.add_subplot(142)
    ax3 = fig.add_subplot(143)
    ax4 = fig.add_subplot(144)
    
    plot_phi_hist(ax1, p)
               
    cax4 = ax4.hist2d(p['circ_phi'], p['circ_y'], weights=p['probability'], 
                      bins=[np.arange(-.15, .2,.01), np.arange(-8, -4, .1)],
                     cmap=plt.get_cmap('gist_heat_r'))
    colbar = fig.colorbar(cax4[3])
    colbar.set_label('counts / bin')
    trans4 = transforms.blended_transform_factory(ax4.transData, ax4.transAxes)
    
    
    for i, row in enumerate(rows):
        # This plot is in terms of phi_folded
        ax1.plot([row['phi_start'], row['phi_stop']],
                     [.98 - i * 0.03, .98 - i * 0.03],  
                     transform=trans1, color=color[i], lw=4)
        # This plot is in term of phi, not phi_folded
        phim = p.meta['phi_m']    
        ax4.plot([phim + row['phi_start'], phim + row['phi_stop']],
                     [.98 - i * 0.03, .98 - i * 0.03],
                     transform=trans4, color=color[i], lw=4, label=row['scenario_name'])
        ax4.plot([phim - row['phi_start'], phim - row['phi_stop']],
                     [.98 - i * 0.03, .98 - i * 0.03],
                     transform=trans4, color=color[i], lw=4, label='__no_legend__')
        
        pdet = p[ad.between(p['phi_folded'], [row['phi_start'], row['phi_stop']])]
        plot_aeff(ax2, row, wave, order0, color=color[i])
        
        if i == 0:
            ax3.plot(wave, row['res_per_order'])
        ax3.plot(wave, row['res'], color=color[i], lw=4)

    ax1.set_xlim([0, None])
    ax2.set_ylabel('Aeff [cm$^2$]')
    ax2.set_xlabel('wavelength [$\AA$]')
    ax3.set_ylabel('Resolving power')
    ax3.set_xlabel('wavelength [$\AA$]')
    ax4.set_ylabel('cross-dispersion [mm]')
    ax4.set_xlabel('dispersion angle [rad]')
    ax4.legend(loc=(.4, .5))
    #ax4.set_ylim([-5, 1])
    #ax4.set_xlim([-.17, .17])
    fig.subplots_adjust(wspace=0.3)
    print('Blaze: {:3.1f} deg - d_channels: {:4.0f} mm - max_f: {:6.0f} mm'.format(p.meta['BLAZE'],
                                                                                  p.meta['D_CHAN'], 
                                                                                  p.meta['MAX_F']))
    print('circle r: {:6.0f} mm - torus R: {:6.0f} mm'.format(p.meta['CIRCLE_R'], p.meta['TORUS_R']))
    return fig

The result of each ray-trace is a file with detected rays, where each ray has some wavelength, probability of survival, and position on the detector. The ray-files simulate only one out of four channels. The location of the photons from the remaining channels can be calculated from symmetry. This requires that two strips of 8 detectors (for a total of 16 CCDs) are placed symmetrically around the symmetry plane.

Different parameters can be optimized when placing the CCDs. I require that the zeroth order is detected and then optimize one of the following: Total number of detected photons, number of photons around 22 $\unicode{x212B}$ (see G1-1 (b) in the science traceability matrix), or number of photons around 35 $\unicode{x212B}$ (G1-1(c) in the science traceability matrix). I will add different strategies at a later point for verification, but my experience is that G1-1 (b) and G1-1 (a) set the strictest requirements. The different CCD placement strategies that I tried give very similar answers, so I do not expect this to change the conclusions significantly.

For each simulation and each CCD placement strategy I obtain one set of detected photons and from these photons I derive effective area curves and resolving power curves, which are shown below. The numbers I derive are all "per spectral order" as in "resolving power per order" and "effective area per order". In practice though, only two or three order contribute significantly at any wavelengths, so looking an a summed effective area and resolving power that is averaged weighted by the effective area in each order is good enough for now.

The following figure displays results for one simulation. Several similar plots are shown below for different parameters.

leftmost panel: This is a histrogram of the dispersion angle measured from the plane of symmetry. This diagram is useful to understand the detector placement. In this particular case, the zero-order photons are at $\theta = 0.05$ radian. (Really, they are at $\theta = -0.05$ but since we have two symmetric detector strips, it is easier to show on the absolute value of $\theta$). The low orders are close to the plane of symmetry and the higher orders further away. The black, blue, and green bars in the top indicate the range of $\theta$ values that are covered by the detectors. The green line is a detector large enough to catch all photons and the other colors are scenarios with two strips of 8 CCDs with slightly different positions.

middle left panel: Effective area for the black, blue, red, and green scenario. Solid line is for dispersed photons (all orders added up) and the dotted line is for the zeroth order.

middle right panel: Resolving power. The thick lines show the averaged resolving power for all dispersed orders, weighted by the number of photons in them. The thin, colored lines are resolving powers for individual orders for the first case (black scenario). At any one wavelength, most photons concentrate in only two or three orders, so these lines are mostly useful to diagnose problems with the ray-trace or analysis.

rightmost panel: Detector image. Note that this is not to scale, the y-axis is enlarged. Again, the colored bars on the top indicate the range of dispersion angle covered by the CCDs. The simulations are for a single channel. The distribution of photons becomes wider further away from the zeroths order. This is due to the Rowland geometry, which optimizes the focus in dispersion direction at the cost of artifacts in the cross-dispersion direction. This also explains the banana shape.

In [9]:

p = ad.load_prepare_table(os.path.join(filepath, os.path.basename(sumtab[30]['filename'])))
fig = plot_everything(p, sumtab[30: 34], wave, orders)

Blaze: 1.8 deg - d_channels:  500 mm - max_f:  11880 mm
circle r:   5945 mm - torus R:   6000 mm

In [10]:

fig = plt.figure(figsize=(12, 4))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)



for ax, j in zip([ax1, ax2], [14, 30]):
    for m, s in zip(['o','s', 'p', '3'], ['all_photons', '8 CCDs', 'G1-1 (a/b)', 'G1-1 (c/d)']):
        t = sumtab[(sumtab['D_CHAN'] == 600) & (sumtab['scenario_name'] == s)]
        cax = ax.scatter(t['aeff_per_order'][:, j, orders != 0].sum(axis=1), t['res'][:, j], c=t['BLAZE'],
                         marker=m, cmap=plt.get_cmap('jet'), label=s)

    col = fig.colorbar(cax, ax=ax)
    ax.set_title('$\lambda = {}$ Ang'.format(wave[j]))
    col.set_label('Blaze angle [deg]')
    ax.set_xlabel('effective area [cm$^2$]')
    ax.set_ylabel('Resolving power')
    ax.legend(loc='lower left', scatterpoints=1)

In [11]:

from matplotlib import patches

fig = plt.figure(figsize=(12, 4))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

t = sumtab[((sumtab['D_CHAN'] == 600) | (sumtab['D_CHAN'] == 600) | (sumtab['D_CHAN'] == 600)) &
          (sumtab['scenario_name'] != 'all_photons') & (sumtab['TORUS_R'] > 5700)]

cax1 = ax1.scatter(t['G1-1 (b)'], t['G1-1 (a)'], s=(t['D_CHAN'] - 400) / 5, c=t['BLAZE'], cmap=plt.get_cmap('jet'))
col = fig.colorbar(cax1, ax=ax1)
ax1.set_title('G1-1 (a/b)')

cax2 = ax2.scatter(t['G1-1 (c)'], t['G1-1 (d)'], s=(t['D_CHAN'] - 400) / 5, c=t['BLAZE'], cmap=plt.get_cmap('jet'))
col = fig.colorbar(cax2, ax=ax2)
ax2.set_title('G1-1 (c/d)')

col.set_label('Blaze angle [deg]')
for ax in [ax1, ax2]:
    ax.set_xlabel('effective area [cm$^2$]')
    ax.set_ylabel('Resolving power')
 
def mark_requirements(ax, xmin, ymin, color):
    # add a rectangle
    ylim = ax.get_ylim()
    xlim = ax.get_xlim()
    rect = patches.Rectangle([xmin, ymin], 1000, 10000, ec="none", fc=color, alpha=0.3, zorder=-5)
    ax.set_ylim(ylim)
    ax.set_xlim(xlim)
    ax.add_patch(rect)
    
mark_requirements(ax1, 256, 2500, 'y')
mark_requirements(ax2, 188, 2250, 'r')
mark_requirements(ax1, 180, 2000, 'y')
mark_requirements(ax2, 125, 2000, 'r')

In [12]:

indblaze = (t['BLAZE'] > 1.55) & (t['BLAZE'] < 1.85)
indchan = (t['D_CHAN'] == 600.)
indtor = (t['TORUS_R'] == 5900.)
ind11ab = (t['scenario_name'] == 'G1-1 (a/b)')
ind11cd = (t['scenario_name'] == 'G1-1 (c/d)')
t18 = t[indblaze & indchan & indtor & (ind11ab | ind11cd)]

In [13]:

t18

Out[13]:

Table length=4

BLAZE	CIRCLE_R	D_CHAN	MAX_F	TORUS_R	aeff [43]	filename	phi_start	phi_stop	prob_per_order [43,19]	res [43]	res_per_order [43,19]	scenario_name	aeff_per_order [43,19]	G1-1 (a)	G1-1 (b)	G1-1 (d)	G1-1 (c)	G1-2 (b)	G1-2 (a)	G1-2 (c)
float64	float64	float64	float64	float64	float64	bytes20	float64	float64	float64	float64	float64	bytes11	float64	float64	float64	float64	float64	float64	float64	float64
1.6	5947.570932742206	600.0	11880.0	5900.0	817.1362813363127 .. 105.39685946546189	0600_01.6_05900.fits	0.02050036237725282	0.08708216602081507	1.0091076703654998 .. 0.0	247.5474623002404 .. 2089.2826981277303	nan .. nan	G1-1 (a/b)	1.2579158037864486 .. 0.0	2621.909537243063	312.5207159630367	2838.006695409179	164.78907516447725	2884.4295713913	252.77089617936815	197.1056142690579
1.6	5947.570932742206	600.0	11880.0	5900.0	817.505721506991 .. 105.39685946546189	0600_01.6_05900.fits	0.0030226389208177287	0.06960444256437998	1.0091076703654998 .. 0.0	247.43559277527297 .. 2089.2826981277303	nan .. nan	G1-1 (c/d)	1.2579158037864488 .. 0.0	2468.5744239003866	270.02227419025115	2390.9007786609186	201.70978293469867	2677.2879309004015	221.45139407007855	169.69702165843157
1.8	5947.570932742206	600.0	11880.0	5900.0	799.4918848218655 .. 68.57232975239238	0600_01.8_05900.fits	0.026326270196064518	0.09290807383962676	0.7183910164144742 .. 0.0	183.98502151411242 .. 2252.03920312321	nan .. nan	G1-1 (a/b)	0.8957669978303312 .. 0.0	3176.058031778364	265.3003015447283	3082.7167072086077	146.04146245946507	3182.165480421398	245.03084777678518	178.7003550375701
1.8	5947.570932742206	600.0	11880.0	5900.0	472.8185068231434 .. 144.51613644472596	0600_01.8_05900.fits	0.04297172110695508	0.10955352475051733	0.7183910164144742 .. 0.0	106.96947724439154 .. 3106.053586642954	nan .. nan	G1-1 (c/d)	0.8959382014151932 .. 0.0	3331.4485850600654	272.56555568970094	3342.4128632945876	187.31935816782934	3365.4064651678736	252.91175480495988	195.93535937548018

In [14]:

p = ad.load_prepare_table(os.path.join(filepath, os.path.basename(t18[0]['filename'])))
fig = plot_everything(p, t18, wave, orders)

Blaze: 1.6 deg - d_channels:  600 mm - max_f:  11880 mm
circle r:   5948 mm - torus R:   5900 mm

In [15]:

p = ad.load_prepare_table(os.path.join(filepath, os.path.basename(t18[2]['filename'])))
fig = plot_everything(p, t18, wave, orders)

Blaze: 1.8 deg - d_channels:  600 mm - max_f:  11880 mm
circle r:   5948 mm - torus R:   5900 mm

In [16]:

p = ad.load_prepare_table(os.path.join(filepath, os.path.basename(t18[1]['filename'])))
fig = plot_everything(p[(p['wave'] > 21.6) & (p['wave'] <= 28.)], t18, wave, orders)

Blaze: 1.6 deg - d_channels:  600 mm - max_f:  11880 mm
circle r:   5948 mm - torus R:   5900 mm

In [17]:

fig = plot_everything(p, t18[3::4], wave, orders)

Blaze: 1.6 deg - d_channels:  600 mm - max_f:  11880 mm
circle r:   5948 mm - torus R:   5900 mm

This figure shows the effective area and resolving power in the wavelength range of science requirements G1-1 (a/b) (left) and G1-1 (c/d) (right). (Strictly speaking the requirements call for averages over a range of wavelengths, while these graphs are for a specific wavelength, but the functions are smooth enough that that does not make a major difference.)

The color of the dots shows the blaze angle of this particular simulation. For each blaze angle, there are simulations with different rotation axes for the torus (and therefor different torus R values) and different CCD placements. According to the proposal the baseline requirements for the left figure are $A_\mathrm{eff} > 450$ cm$^2$ and $\lambda/\Delta\lambda =2500$ and $A_\mathrm{eff} > 250$ cm$^2$ and $\lambda/\Delta\lambda =2000$ for the right figure. Blaze angles in the range of 1.7-1.9 deg reach these requirements with enough margin to account for the caveats listed above.

The following plots show several of the scenarios in this range in more detail. In these cases, the black, blue, and green scenarios in each figure differ by the radius of the Rowland torus (or the tilt angle); they all use the same number of CCDs. (The histogram in the leftmost figure only shows the distribution for the black scenario, this can be slightly different for the other Rowland torus radii.)

In [18]:

tab600 = sumtab[(sumtab['D_CHAN'] == 600) & (sumtab['scenario_name'] == 'G1-1 (a/b)')]

In [19]:

p = ad.load_prepare_table(os.path.join(filepath, os.path.basename(tab600[4]['filename'])))
fig = plot_everything(p, tab600[4:], wave, orders)

Blaze: 1.6 deg - d_channels:  600 mm - max_f:  11880 mm
circle r:   5948 mm - torus R:   5900 mm

In [20]:

tab600 = sumtab[(sumtab['D_CHAN'] == 600) & #(sumtab['scenario_name'] == 'G1-1 (c/d)') & 
                (sumtab['BLAZE'] > 1.75)  & (sumtab['BLAZE'] < 1.85)]

Plots for publications¶

Below are a few plots from the notebook above that are optimized for print in publications with e.g. manual modification of plot limits or font sizes.

In [21]:

p = ad.load_prepare_table(os.path.join(filepath, os.path.basename(t18[2]['filename'])))
fig = plot_everything(p, t18[2:], wave, orders)
ax4 = fig.axes[3]
leg = ax4.get_legend()
leg.get_texts()[0].set_text('O VII')
leg.get_texts()[1].set_text('C VI')
#leg.get_texts()[2].set_text('O VII (1.8$^\circ$)')
#leg.get_texts()[3].set_text('C VI (1.8$^\circ$)')

fig.savefig(os.path.join(get_path('figures'), 'selectrowland.pdf'), bbox_inches='tight') 

Blaze: 1.8 deg - d_channels:  600 mm - max_f:  11880 mm
circle r:   5948 mm - torus R:   5900 mm

In [22]:

fig = plt.figure(figsize=(10, 3))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

t = sumtab[(sumtab['D_CHAN'] == 600) & (sumtab['scenario_name'] != 'all_photons') & (sumtab['TORUS_R'] == 5900)]

cax1 = ax1.scatter(t['G1-1 (b)'], t['G1-1 (a)'], s=(t['D_CHAN'] - 400) / 5, c=t['BLAZE'], cmap=plt.get_cmap('jet'))
col = fig.colorbar(cax1, ax=ax1)
ax1.set_title('O VII')

cax2 = ax2.scatter(t['G1-1 (c)'], t['G1-1 (d)'], s=(t['D_CHAN'] - 400) / 5, c=t['BLAZE'], cmap=plt.get_cmap('jet'))
col = fig.colorbar(cax2, ax=ax2)
ax2.set_title('C VI')

col.set_label('Blaze angle [deg]')
for ax in [ax1, ax2]:
    ax.set_xlabel('effective area [cm$^2$]')
    ax.set_ylabel('Resolving power')
    
fig.savefig(os.path.join(get_path('figures'), 'selectrowland2.pdf'), bbox_inches='tight') 