PSF to Cosmology projects

Here I hope to collect relevant results, and papers regarding PSF's and how they affect shape/shear measurements of galaxies for weak lensing.

• https://arxiv.org/abs/0711.4886
• http://faculty.washington.edu/ivezic/Teaching/Astr511/LSST_SNRdoc.pdf
• https://github.com/lsst/sims_seeingModel/files/4037941/simsee_short_summary.pdf 

relevant code:

• Txpipe: https://github.com/LSSTDESC/TXPipe
  
  • https://github.com/LSSTDESC/TXPipe/blob/master/txpipe/input_cats.py#L1036
• seeing simulations https://github.com/lsst/rubin_sim/blob/main/rubin_sim/site_models/seeingModel.py  

Rubin 10 year survey simulations (courtesy of Lynne Jones at UW)

• http://astro-lsst-01.astro.washington.edu:8080/allMetricResults?runId=1#Seeingeff 

Feb 12 2022 version of paper, PDF: RubinSeeingImpact.pdf

Feb 12 2022, Stubbs. 

The FWHM in a given image arises from four factors:

1) the atmospheric 'seeing', defined at looking to the zenith at a wavelength of 500 nm, path-integrated from the telescope slit upwards to the top of the atmosphere,

2) contributions from the Rubin system image quality budget, that includes mirror figure errors, misalignments, sensor effects and (importantly for this work) dome seeing, and

3) The zenith angle (also denoted airmass) at which an observation is obtained. In general we expect worse seeing at larger zenith angles, since we're looking through more atmosphere. There is a simplifying assumption that this effect is axisymmetric. 

4) The optical bandpass through which the image is taken. Shorter wavelengths have poorer FWHM.  

We don't need to use all bands for what we're doing. I suggest we focus on the r band. Seeing in y and z isn't much better than in r, and sky is in general brighter. So shape measurements tend to best in r band and this will tell us what we need to know. 

Here is a figure taken from the MAF output page (http://astro-lsst-01.astro.washington.edu:8080/allMetricResults?runId=1#Seeingeff ) for the histogram of FWHM across all surveyed parts of the sky for the r band. This accounts for all the contributions listed above. The expected site seeing at zenith at 500 nm (from site survey data) is converted to what we'd expect at a given zenith angle for the simulated survey, and is also converted into what we expect for the r band. The Rubin system contribution is estimated at 0.39 arcsec, and is added in quadrature (FWHM^2=(r-band-at-zenith-angle-FWHM)^2 + 0.39^2). 

The tool I used to use for digitizing things doesn't run on my 64 bit Mac so I tried https://automeris.io/WebPlotDigitizer/
That worked well! Here's the digitized data (FWHM, Nimages) . Need to click on 'View Data' to get this. 
0.295, 77.5
0.425, 2.60e+3
0.460, 5.04e+3
0.495, 7.44e+3
0.538, 1.00e+4
0.564, 1.26e+4
0.616, 1.49e+4
0.634, 1.58e+4
0.668, 1.70e+4
0.686, 1.74e+4
0.729, 1.76e+4
0.781, 1.74e+4
0.799, 1.69e+4
0.825, 1.66e+4
0.851, 1.60e+4
0.894, 1.47e+4
0.946, 1.32e+4
0.998, 1.16e+4
1.05, 1.02e+4
1.10, 8.91e+3
1.17, 7.05e+3
1.28, 5.08e+3
1.34, 3.91e+3
1.48, 2.64e+3
1.73, 1.24e+3
2.03, 504
2.44, 233
2.84, 155

Let's compute delivered FWHM histogram as a function of Rubin system contribution to seeing budget. All units are in effective FWHM. 

Vertical axis has been normalized to be fraction of the total. 

Cumulative fractions:

So the next step is to attempt to estimate the effective FWHM that would be generated from an optimal stacking of a set of images with these different seeing profiles. 
We're going to assume that all images in the survey have essentially this distribution. That won't quite be true (ones at declination -30 should be observed at lower airmass than others) but it'll have to do for now. 

As FWHM degrades, three things happen: shapes are circularized, SNR for galaxies goes down, and star-galaxy discrimination and deblending both degrade. For sky-dominated imaging of barely resolved galaxies, SNR goes as 1/FWHM (more sky). Generically let's assume a frame-by-frame weighting for image stacking that goes as 1/FWHM^gamma, where we'll let gamma vary in this study with 0<=gamma <=3. That should bracket the possibilities I think. Gamma = 0 is equal weighting of all images. 

So for each value of system FWHM and value of gamma, we can compute an effective FWHM for the image stack by taking an appropriately weighted average, 

FWHM_stacked(RubinFWHM, gamma) = sum_i  FWHM(RubinFWHM)_i*(image fraction)_i]*FWHM^(-gamma)_i / sum_i  (imagefraction_i)*FWHM^(-gamma)_i   

If I've done this correctly (and it surely needs to be checked!!) then I get the following Table for effective stacked r-band FWHM in arcsec for stacked images from the entire survey:

This does not take into account the change in effective image depth for these different weighting schemes. For comparison, the r band value assumed in the Ivezic et al document linked above is 0.70 arcsec, so we aren't very far off....

To obtain equivalent tables in other bands, multiply the r-band numbers above by the factors given in the Ivezic et al paper above, i.e. by FWHM_ugrizy= r band FWHM * (1.10    1.04    1.00   0.96    0.93    0.90)

We can also compute the derivative d(FWHM_stacked)/d(Rubin_FWHM) from the table above, and obtain (for different weighting schemes):

Unsurprisingly, improving Rubin performance matters more when the best images are weighted most strongly, at gamma=3. If we improve Rubin system FWHM by X, stacked image improves by (0.53 to 0.86) * X depending on weighting scheme, evaluated at baseline expectation of 0.399 arcsec Rubin FWHM. 

MATLAB code for all this: 

------

%% RubinSeeing.m adjusts OpSim seeing histogram to other values of system contribution

% C. Stubbs Feb 12, 2022. 

%% input data from MAF metrics page. This is delivered r band FWHM histogram for a 10 year survey, baseline cadence

% digitized from histogram at http://astro-lsst-01.astro.washington.edu:8080/allMetricResults?runId=1#Seeingeff

% columns are effective FWHM in arcsec and number of images

% Approach here is to compute effective FWHM for a stacked image as a function

% of two parameters: i) the Rubin system contribution to seeing, with

% values of 0.3, 0.4, and 0.5 arcsec, and ii) a weighting factor gamma

% where weight = 1/FWHM^gamma, with gamma ranging from 0 to 3. Gamma=0 is

% equal weighting, Gamma=1 is SNR for a point source, and higher values

% take into account factors like shape transfer function, blending, and

% star-galaxy discrimination. Down the road someone needs to determine best

% value for gamma, or some other weighting scheme. 

OpSimSystem=0.399;  % value added in quadrature for OpSim results

OpSim=[

0.295, 77.5

0.425, 2.60e+3

0.460, 5.04e+3

0.495, 7.44e+3

0.538, 1.00e+4

0.564, 1.26e+4

0.616, 1.49e+4

0.634, 1.58e+4

0.668, 1.70e+4

0.686, 1.74e+4

0.729, 1.76e+4

0.781, 1.74e+4

0.799, 1.69e+4

0.825, 1.66e+4

0.851, 1.60e+4

0.894, 1.47e+4

0.946, 1.32e+4

0.998, 1.16e+4

1.05, 1.02e+4

1.10, 8.91e+3

1.17, 7.05e+3

1.28, 5.08e+3

1.34, 3.91e+3

1.48, 2.64e+3

1.73, 1.24e+3

2.03, 504

2.44, 233

2.84, 155];

FWHM04=OpSim(:,1);

% clip at 0.41 for sanity check, we can't ever have seeing better than 0.4

% given the Rubin system contribution. This avoids artefacts from

% digitization

FWHM04=max(0.41,FWHM04);

N04=OpSim(:,2)/sum(OpSim(:,2)); % normalize to get fraction of images at this seeing value. This is a bit of a fudge since bin widths aren't equal. 

AtmosFWHM=sqrt(FWHM04.^2-OpSimSystem.^2); % back out the system contribution to arrive at bare-atmosphere contribution in r band at these zenith angles. 

% now compute deliverd FWHM if we had different values of Rubin system

% seeing. 

FWHM03=sqrt(0.3.^2 + AtmosFWHM.^2);

FWHM05=sqrt(0.5.^2 + AtmosFWHM.^2);

%% Generate effective FWHM for stacked images, for different values of weighting factor Gamma. 

gamma=[0 1 2 3];

WeightVector03_0=N04.*FWHM03.^(0);

WeightVector03_1=N04.*FWHM03.^(-1);

WeightVector03_2=N04.*FWHM03.^(-2);

WeightVector03_3=N04.*FWHM03.^(-3);

StackedFWHM03_0=sum(FWHM03.*WeightVector03_0)/sum(WeightVector03_0)

StackedFWHM03_1=sum(FWHM03.*WeightVector03_1)/sum(WeightVector03_1)

StackedFWHM03_2=sum(FWHM03.*WeightVector03_2)/sum(WeightVector03_2)

StackedFWHM03_3=sum(FWHM03.*WeightVector03_3)/sum(WeightVector03_3)

WeightVector04_0=N04.*FWHM04.^(0);

WeightVector04_1=N04.*FWHM04.^(-1);

WeightVector04_2=N04.*FWHM04.^(-2);

WeightVector04_3=N04.*FWHM04.^(-3);

StackedFWHM04_0=sum(FWHM04.*WeightVector04_0)/sum(WeightVector04_0)

StackedFWHM04_1=sum(FWHM04.*WeightVector04_1)/sum(WeightVector04_1)

StackedFWHM04_2=sum(FWHM04.*WeightVector04_2)/sum(WeightVector04_2)

StackedFWHM04_3=sum(FWHM04.*WeightVector04_3)/sum(WeightVector04_3)

WeightVector05_0=N04.*FWHM05.^(0);

WeightVector05_1=N04.*FWHM05.^(-1);

WeightVector05_2=N04.*FWHM05.^(-2);

WeightVector05_3=N04.*FWHM05.^(-3);

StackedFWHM05_0=sum(FWHM05.*WeightVector05_0)/sum(WeightVector05_0)

StackedFWHM05_1=sum(FWHM05.*WeightVector05_1)/sum(WeightVector05_1)

StackedFWHM05_2=sum(FWHM05.*WeightVector05_2)/sum(WeightVector05_2)

StackedFWHM05_3=sum(FWHM05.*WeightVector05_3)/sum(WeightVector05_3)

%% plots

figure(10)

plot(AtmosFWHM,N04,'ko-')

hold on

plot(FWHM03,N04,'bo-')

plot(FWHM04,N04,'ro-')

plot(FWHM05,N04,'mo-')

hold off

grid on

legend('Atmos only','0.3" system','0.4" system','0.5" system','FontSize',14)

xlabel('FWHM, arcsec','FontSize',14)

ylabel('Fraction of r band images','FontSize',14)

shg

figure(20)

plot(AtmosFWHM,cumsum(N04),'ko-')

hold on

plot(FWHM03,cumsum(N04),'bo-')

plot(FWHM04,cumsum(N04),'ro-')

plot(FWHM05,cumsum(N04),'mo-')

hold off

legend('Atmos only','0.3" system','0.4" system','0.5" system','FontSize',14)

xlabel('FWHM, arcsec','FontSize',14)

ylabel('Cumulative fraction of r band images','FontSize',14)

grid on

shg

------