Stubbs Analysis Notes Sept 24 2022

Ali sent some FITS files on Slack. Images are 

-rw-r--r--@     1 cstubbs  staff    58M Sep 24 06:53 CU2A6058.fits

-rw-r--r--@     1 cstubbs  staff    58M Sep 24 06:55 CU2A6042.fits

-rw-r--r--@     1 cstubbs  staff    58M Sep 24 06:55 CU2A6050.fits

-rw-r--r--@     1 cstubbs  staff    58M Sep 24 06:55 CU2A6081.fits

my first comment is that we need to be sure the analysis pipeline appends a date+time stamp when the FITS files are made. I think the reduction scripts I wrote long ago do that, if not we need to make sure we do that. 

Images look like this: 

This image 'stretch' (mapping of intensity values to greyscale) shows more clearly that the main sources totally stand out in the frame: 

Here is what I see: 

• there is horizontal striping that we can get rid of by subtracting a median column from each column in the frame. 
• There are 7 bright sources and each of them has a ghost to the left, presumably from an internal reflection in the lens
• There is a source of stray light in the bottom of the frame that's bright smudge
• There is what looks like a lens flare ghost as an arc emanating out clockwise from the bright-smudge source, that overlaps one of the sources. 
• In all 4 frames, the sources of interest don't shift much. 

Sources of interest are at (roughly) 

Bad smudge ends at y=500 and we don't have a source below y=500. 

Here is what I would do to process these frames and pull out the centroids of interest: 

1. Find and kill the stray light source. But we can't go back in time and do that, so instead just cut off y<500 from the frames, to make things cleaner
2. compute the median of all columns and subtract that from each column, to get rid of the horizontal streaking
3. Either smooth the image with a Gaussian kernel with sigma~3 pixels to account for residual Bayer pixel artifacts and fragmentation of the PSF
4. run a source detection code with a fairly high threshold of significance, and then select out the bright sources that have centroids near where we expect them. 
5. compute mean of all 7 centroids, call that the 'average position' of all 7
6. subtract that average from each of the 7 positions, which is differential motion of each spot. 
7. I'd also compute the mean intensity for all 7 spots, and take the ratio of each spot's intensity to that mean. This will let us look at scintillation as well as image motion. 
8. Put the results into some sensible data structure(s)

For reasons of time I'm going to skip steps 1,2,3 and run Source Extractor on the four frames, with a tweaked configuration file that requires high significance for detection. 

I ran source extractor with this configuration file: sex.config

and this default.sex file: default.sex (note this has an edited value for DEBLEND_NTHRESH, set to 1 here but was 32 for initial results shown below):

sex CU2A60*1.fits 

And got this output file, called test.cat:

So that seems good. It caught all 7 sources plus 5 more. All look fairly compact. So why the extra ones, and where are they? 

Extract x and y values: Use grep to get rid of comment lines and awk to print column 2. These are x values

cat test.cat | grep -v "#" | awk '{print $2}'

4506.706

772.435

3774.000

792.160

5232.965

5225.302

5239.854

5234.985

1511.580

3739.769

752.534

750.389

Same for y:

cat test.cat | grep -v "#" | awk '{print $3}'

888.634

1602.611

1639.883

1591.234

2393.825

2398.723

2395.239

2402.175

3099.798

3863.995

3849.634

3835.964

Plot x vs y to see where they are:

Ah, OK. So it's finding all 7 sources that we want, but the spiky PSF is making it break ('deblend') them into multiple adjacent sources. We can do one of two things: 1) convolve image with a Gaussian to to smooth out the PSF, or 2) tweak the deblending value that Source Extractor uses to fragment sources. Let's try setting number of deblending levels to 1, edited default.sex file to have 

DEBLEND_NTHRESH  1             # Number of deblending sub-thresholds

Re-ran source extraction with this, got this result: 

Produce x-ordered listing to compare to expectations: 

cat test.cat | grep -v "#" | sort -k 2 -n 

This looks good now! 

Run on all 4 frames: 

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ ls *.fits

CU2A6042.fits  CU2A6050.fits  CU2A6058.fits  CU2A6081.fits

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ sex CU2A6042.fits 

----- SExtractor 2.5.0 started on 2022-09-24 at 09:47:41 with 1 thread

Measuring from: "Unnamed"  / 6744 x 4502 / 16 bits INTEGER data

(M+D) Background: 511.912    RMS: 122.625    / Threshold: 183.937    

Objects: detected 7        / sextracted 7                      

> All done (in 3 s)

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ mv test.cat  CU2A6042.cat 

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ sex CU2A6050.fits 

----- SExtractor 2.5.0 started on 2022-09-24 at 09:48:02 with 1 thread

Measuring from: "Unnamed"  / 6744 x 4502 / 16 bits INTEGER data

(M+D) Background: 511.917    RMS: 122.511    / Threshold: 183.767    

Objects: detected 7        / sextracted 7                      

> All done (in 3 s)

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ mv test.cat CU2A6050.cat

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ sex CU2A6058.fits 

----- SExtractor 2.5.0 started on 2022-09-24 at 09:48:21 with 1 thread

Measuring from: "Unnamed"  / 6744 x 4502 / 16 bits INTEGER data

(M+D) Background: 511.912    RMS: 122.556    / Threshold: 183.835    

Objects: detected 7        / sextracted 7                      

> All done (in 3 s)

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ mv test.cat CU2A6058.cat

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ sex CU2A6081.fits 

----- SExtractor 2.5.0 started on 2022-09-24 at 09:48:49 with 1 thread

Measuring from: "Unnamed"  / 6744 x 4502 / 16 bits INTEGER data

(M+D) Background: 511.919    RMS: 122.474    / Threshold: 183.71     

Objects: detected 8        / sextracted 8                      

> All done (in 3 s)

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ mv test.cat CU2A6081.cat

HUIT-FAS-MacBook-Stubbs:Q42022 cstubbs$ wc *.cat

      25     267    2500 CU2A6042.cat

      25     267    2500 CU2A6050.cat

      25     267    2500 CU2A6058.cat

      26     285    2683 CU2A6081.cat

      30     357    3415 test2.cat

     131    1443   13598 total

Ok, so we have one extra source in frame 81. Huh. 

I think it found a companion ghost. extract flux values: 

cat CU2A6081.cat | grep -v "#" | awk  '{print $4}' 

make a diagnostic plot of flux vs. x

Sure enough, there is a discrepantly low flux point. So if we pick the 7 brightest sources we ought to be able to avoid the ghosts. This command line does that: 

cat CU2A6081.cat | grep -v "#" | sort -k 4 -n -r | head -7 > stripped81.cat

That prints the file to the grep command that strips out comment lines, pipe the output to sort tool using column 4 interpreted as numbers, in descending order piped to head command to keep brightest 7 sources. 

do 'man grep', 'man sort', etc to learn about command line tools. 

What should we expect for sensitivity? RMS width of PSF is around 4 pixels, and FWHM=2.35*RMS ~ 9.5 pixels. Fluxes are around 100K. Assume (guess!) around 1e per ADU so that corresponds to 100K electrons. Signal to Noise ratio is square root of that, or ~300
Sigma of centroid should be FWHM/SNR = 9/300 ~ 0.03 pixels. Cool. 

I made trimmed source catalogs that retained the brightest 7 compact sources in each image. 

We then need to do source association in (x,y) space. I cheated a bit and did this in 1 dimension, using ordering in y. This suggests we arrange the sources and their orientation on the detector to not have overlap in either x or y, to ease source association. 

Computed mean position for each spot, then computed magnitude of difference of (postiion-mean) for each spot pair, and plotted that vs. separation. Looks interesting. The different frames have different amounts of turbulence: 

Motion of center-of-spots: 

Magnitude of differential motion vs. separation between pairs, for all 4 images: 

Looks like fitting differential motion vs. sqrt(separation) might work. Did those fits, here are the results: 

Fit to

differential motion = a + b * sqrt(separation)

MATLAB code DIMSUM.m is here: DIMSUM.m