Create a model for the rest wavelength for a given date by wave type

[mgalloy]

Compute the rest wavelength for CoMP the same way we did for UCoMP, i.e., compute the rest wavelength for each median file, plot and fit, and then create a model that we apply for rest wavelengths.

Compute the rest wavelength three ways:

1. full frame annulus
2. E and W
3. device E and W

Do this for 1074 and 1079, using the synoptic median average file in both cases.

Tasks

• use code from Giuliana
• use method "22" to calculate the RSTWVL and to apply to velocity for both the dynamics file and the quick invert file
• remove "device" method calculations and FITS keywords
• fix issue of undefined RSTWVL

Questions

• Do a single fit? or for each "good" pixel?

[detoma]

I think we want to do it both ways. I started writing the basic code and will pass it to you for checking and integration.

[detoma]

Use code from Giuliana to analyze rest wavelength from median synoptic files

[mgalloy]

The results from the analysis are in comp.restwvl.synoptic.txt. The columns of the file are:

date,
doppler_shift_0, velocity_0, line_width_0, peak_intensity_0,
doppler_shift_1, velocity_1, line_width_1, peak_intensity_1,
doppler_shift_2, velocity_2, line_width_2, peak_intensity_2,
doppler_shift_22, velocity_22, line_width_22, peak_intensity_22,
doppler_shift_3, velocity_3, line_width_3, peak_intensity_3

[mgalloy]

Here are the rest wavelengths from the various methods:

[mgalloy]

Here are the plots as above, but using means instead of medians. Data file is comp.restwvl.mean.synoptic.txt.

[mgalloy]

Here are the old median plots, but on the same fixed y-range for easier comparison.

[detoma]

The code indicate that the east and west method that use these condition works best:

 east = where(mask gt 0 and velocity ne 0 and velocity gt -30 and velocity lt 20 and x lt 0.0  $
                                    and intensity[*,*,0] gt 0.5 and intensity[*,*,1] gt 2 and intensity[*,*,2] gt 0.5 $
                                    and intensity[*,*,0] lt 60   and intensity[*,*,1] lt 60 and intensity[*,*,2] lt 60 $
                                    and line_width gt 15.0 and line_width lt 50.0 , n_east)
                       
 west = where(mask gt 0 and velocity ne 0 and velocity gt -30 and velocity lt 20  and x gt 0.0  $
                                    and intensity[*,*,0] gt 0.5 and intensity[*,*,1] gt 2 and intensity[*,*,2] gt 0.5 $
                                    and intensity[*,*,0] lt 60   and intensity[*,*,1] lt 60 and intensity[*,*,2] lt 60 $
                                    and line_width gt 15.0 and line_width lt 50.0 , n_west)

where mask is a geometric mask that overmask the occulter by 4 pixel and under mask the field stop by 10 pixels. line-width is sigma defined as W/sqrt(2) and W is the output of the gaussian fit. If W is used instead of the sigma then the min and max would be: 22 and 71.

The rest wavelength will be defined as the average of the median velocity(east) and median velocity(west).

Please note that albeit the time series of the mean velocity and median velocity show similar noise level, the mean and median differ in value, especially early in the mission. I analyzed the velocity distribution of 25 median images between 2013 and 2018. The distributions are non gaussian and the median is closer to the peak of the distribution, therefore I decided to us the median velocity to determine the rest wavelength.

[detoma]

If there are not enough points in the east(west) limb to compute a median wavelength, the rest wavelength should not be set equal to the west(east) value. These unlikely cases must be handled differently.

The condition temp_velo ne 0 below is not correct:

 temp_velo = (temp_velo - rest) * c / nominal

 good_vel_indices = where(mask gt 0 $
                            and temp_velo ne 0 $
                            and abs(temp_velo) lt 100 $

At this point temp_velo eq 0 is actually temp_velo eq rest. It needs to be written differently or done before subtracting rest.

[detoma]

The wrong thresholds are used in comp_compute_rest_wavelength.pro because line width has already been converted to FWHM.

I will re-write comp_compute_rest_wavelength.pro and comp_l2_analytical.pro and will commit to GitHub when they cleaned up. Since I cannot test the pipeline code, it will require some testing after I commit it.

Finally, the code that creates the quick invert will need to be made consistent with the new comp_l2_analytical.pro

[detoma]

Save the east median, west median, and the average of the two which is the rest wavelength RSTWVL in the header and database. Possible names are ERESTWVL and WRESTWVL which Ike has already defined.

[mgalloy]

OK, the above change allows the comp_compute_rest_wavelength to find good pixels to find a rest wavelength. The change is to not add the rest wavelength, i.e., replacing

temp_velo = temp_data[*, *, 1] + rest

with

temp_velo = temp_data[*, *, 1]