ASAP Velocities incorrect

[Simon.Ellingsen@…]

We have been monitoring a particular spectral line source with Hobart (Mon R2) for more than a year now. We observe with fixed frequency and align in velocity when averaging the data during ASAP processing. When comparing spectra taken months apart we noticed differences of around 1 km/s between the same features in different spectra.

We have data collected with the ATCA and Hobart on the same day. Processing the ATCA data in miriad and the Hobart data in ASAP the velocities are different there as well. The source coordinates and observing frequency are recorded correctly in the RPFITS files (as best we can tell) and ASAP identifies the data as coming for Hobart.

[Malte Marquarding]

Hi Simon,

Max just fixed a defect with respect to frequency alignment (Ticket #142). Can you make the data and script available so I can run it through the fixed version. I'll send you back the results.

Cheers, Malte.

[Malte Marquarding]

Hi Malte,

Attached is 3 epochs of the Mon R2 data that Simon was talking about.  
The script that we would usually use to reduce it is fairly lengthy  
and so below is just a basic text that produces the same problem with  
the velocities. I have also attached an eps file which is a spectrum  
of the three (plus an additional one) epochs of data overlaid.

If you need anything else let me know,

Shari

[Malte Marquarding]

ASAP version 2.3 pre-release result

[Malte Marquarding]

Hi

I can't see a 1 km/s shift in the feature for 0.1 which is close to the channel width anyway.

If you point me at the ATCA data and reduction script I get Max to have a look what the difference could be.

Malte.

[Malte Marquarding]

Hi Malte,

That is exactly the same as what I got - the only difference is that I  
had an extra epoch of data in mine. If you look at the plot that I  
sent you there is a large difference in the data set that I didn't  
send you because the file was bigger. Between the two extremes there  
is 3 channels difference. I could probably find some with even bigger  
differences if I go through all of the data rather than just selecting  
files at random if you would like them?

Shari

[Malte Marquarding]

It would be good if you could send me the ones with the biggest difference.

Can you reply through the ticket interface rather than email. This way we can keep a history.

Cheers, Malte

[Malte Marquarding]

For reference that data/script are on delphinus:

/DATA/DELPHINUS_3/mar637/monr2.tar.bz2

[Malte Marquarding]

Hi,

I had a look at the input data and you are using TOPOcentric frequency, but your Sky frequency is the same for all epochs. This should vary for different epochs.

The only workaround I can see you doing is calculating the shift in frequency manually for each scantable and use scantable.shift_refpix. This only allows you integer shift though.

Long term you need to fix how these frequencies are written into the RPFITS file.

[Max Voronkov]

A reply from Simon is below'

Hi Malte,

Sorry I'm not sending this through trac, but I can't remember my username/password at the moment.

You've got me completely confused with a reference frame for frequencies??!!?? Its not like we have a choice as to which reference frame the observations are made in, they are made in the observatory frame (which I thought was the topocentric frame) and this is true for every observatory, not just Hobart? The sky frequency recorded in the RPFITS file was the frequency at which the observations were made.

What should the reference frame be set to for the frequency in the RPFITS file in our case? An integer channel shift wont be sufficiently accurate for our purposes, so if we can find out what the RPFITS should have to do what we want it to then we can edit the RPFITS prior to processing, either by hand or a program.

[Max Voronkov]

my reply to Simon's e-mail

Hi Simon,

Malte and myself discussed it a bit last Friday. It looks like the issue is not as simple as we originally thought. We'll continue investigation and get back to you. At this stage it looks like nothing is done to the frequency axis by ASAP for your dataset and it is the correct behavior according to the header information and our current understanding. I guess the question about RPFITS frame is largely about how the doppler tracking is done if it is done at all. All three datasets I've seen so far have the same frequency of the first channel, which would be the case for no doppler tracking (i.e. like for ATCA observations). However, I would expect larger velocity shifts than a few channels in this case, if the velocity correction is completely ignored. If doppler tracking is done at the observatory, I'd expect the frequency of the first channel to be different for sufficiiently different dates of observations (as it is for Mopra and Parkes with the doppler tracking on). Either case should be handled well by ASAP, but if the header information in RPFITS doesn't reflect the doppler tracking accurately there could be a problem. How do you handle doppler tracking in Hobart?

Cheers,

Max

[Max Voronkov]

A reply from Simon

We observe with fixed frequency (i.e. no Doppler tracking). ASAP is certainly doing something with the frequency information as the velocity is nearly correct (i.e. its always within about 1 km/s). The way it changes as a function of time suggests to me one of the following is most likely :

• the wrong observatory location is being used in the Doppler

correction calculations

• the wrong source position is being used in the Doppler correction

calculations (I've seen this sort of thing in the past when J2000 coordinates were interpretted as B1950)

[Max Voronkov]

I have spent quite some time today on the ticket, but unfortunately couldn't find the fix. Below there is a list of my findings for the record

• the problem is also present in the frequency (rather than velocity) spectra
• it seems that FrequencyAligner is doing the right thing. First, it doesn't have much to align (if we process individual datasets separately) and the alignment is not necessary if the spectrum is just exported into an ascii file. Second, the spectra look very similar if the alignment is turned off completely.
• When the spectrum is exported, the conversion takes place in Scantable::getAbcissa. The code looks alright from the first glance. The conversion is done with the SpectralCoordinate. All frame information is passed via STFrequencies::getSpectralCoordinate.
• The input spectra are in the TOPO frame and the output are in the LSRK frame. If TOPO frame is chosen for output the velocities/frequencies are notably different (so the conversion roughly works).
• I have checked that the source position is set to 06:07:48 -6:22:57 (J2000), which seems correct
• I have also checked that antenna position (in ITRF) agrees with that used for VLBI to within 3cm, which should be more than enough for frequency conversion purposes.
• The direction/position given above are passed with the frame information in both places where the conversion is done (inside getAbcissa and inside FrequencyAligner)

So at the end I wouldn't be surprised if the issue lies at a higher level in casacore (SpectralCoordinate or measures). It is also worth mentioning that I recently had a problem with casacore's SpectralCoordinate in another code also based on casacore. It appeared recently when the FITS reading code was changed to use WCS. Apparently spectral coordinate didn't receive much attention. Although that problem doesn't look like the subject of this ticket, it is possible that the spectral axis is corrupted somehow at the time an rpfits file is read. It may also be worth checking whether the problem exists for Mopra or Parkes data taken without Doppler tracking or it is just Hobart problem.

[Max Voronkov]

Spectra obtained via miriad for 3 selected dates

[Max Voronkov]

Hi Malte,

I managed to run miriad's uvspec on the 3 rpf files you gave me and attached the result. There is a better alignment indeed. The data reduction was very basic. I just extracted a single ON spectrum from each date, and subtracted the flux of the first pixel to get roughly similar baselines for all three spectra. We can probably take the peak LSR velocity from these spectra for the clean test of casacore's measures.

Cheers,

Max

[Malte Marquarding]

from simon

Hi Maxim,
		We do have an observation of Mon R2 at the ATCA on the same day as  
one of the Hobart observations if that would help in trying to track  
down the problem.  The multiple ATCA observations we have all align  
nicely.  I also have a fortran program based on Mark Calabretta's AT  
ephemeris library which I have always found accurate to 0.1 km/s or  
better.  If you can give me the value the CASA routines are  
calculating for the Doppler shift for a particular date/time on Mon R2  
I can make a comparison using that.  It might also be worth doing some  
comparisons with a public utility like RV for a broader comparison?

[Max Voronkov]

Hobart position from Chris (ITRF)

VLBA SCHED

X= -3950236.7341 Y= 2522347.5530 Z= -4311562.5434

DiFX

HOB -3950236.733 2522347.553 -4311562.537 -0.03861 0.00778 0.03842 1997.0 4 8.1913

[Malte Marquarding]

casacore program to demonstrate shift in frequency

[Malte Marquarding]

I have create a small program for casacore which shows this problem and send it to Wim Brouw (the creator of the measures code). The program is attached to the page.

[Malte Marquarding]

from Wim Brouw:

Onr suggestion: (I never heard any problems from wsrt where it is used 
always): could it be that the online frequencies are not corrected for the 
topo effects, and hence they will be applied twice in the end?

As a result I tried the input frequency as GEO and the values are much closer. The maximum distance between peaks is now only 1.7 x delta f instead of 3.2. Also, the middle scan is now the furthest off not the last one.

Could it be that there is something wrong in the frequency setup for Hobart as suggested by Wim?

Malte.

[Malte Marquarding]

Using GEO as teh base frame instead of TOPO

[Max Voronkov]

two replies from Simon

Hi Malte,

I'm not sure that I put a lot of store in the problem not being observed at WRST. To my knowledge they do _no_ maser work there (I've certainly never seen a paper). This error is small enough that you would never notice it in broad thermal HI lines.

I don't really understand Wim's comment "could it be that the online frequencies are not corrected for the topo effects, and hence they will be applied twice in the end". There are no frequency corrections/ changes during, the observing frequency is fixed at some integral MHz value for the entire observations, I think for Mon R2 we have used the same observing frequency for every observing session (in excess of 50 to date).

Hi Malte,

I've just used Mark Calabretta's AT ephemeris library code to calculate the sky frequency of the Mon R2 peak (10.4 +/- 0.1 km/s) for the three days of data that Shari sent. If you can plot the spectra with frequency on the x-axis, prior to any velocity correction (I assume with the TOPO frame) then I expect the maser peak to be at the following frequencies :

Frequency for 10.4 km/s at 06:07:47.9 -06:22:56.6 (J2000) at Hobart 2008 003/11:51:15 4765.034290 MHz (velocity 33.197247) 2008 077/06:41:40 4764.697714 MHz (velocity 54.370620) 2008 149/03:34:46 4764.929807 MHz (velocity 39.770061)

The times I've used are the times given in the UT08XXX...txt files from the Trac page and I have assume that this is the time used for the velocity alignment within ASAP when it is applied.

[Max Voronkov]

Reply from Malte

Hi Simon,

if you look at the cc file attached to the wiki, those are roughly the frequencies I have determined from the spectrum using max()

So at least these seem to be consistent, which would mean the observed frequencies are TOPOcentric.

I'll get some more data next week, maybe by plotting enough points we can observe a trend. Wim is still going to look at it anyway.

I will also do the same "backwards" calculation as you using casacore.

Cheers, Malte.

[Max Voronkov]

My experiments with miriad showed that it gave the aligned LSR velocity/frequency for all epochs. It is clear from rpfits headers that the frequency of the first channel didn't change at all (i.e. no doppler tracking in any form). I'd say this supports the view that the problem is at our end (probably casacore), not in Hobart.

[Malte Marquarding]

Hi all,

I have tracked this down. casacore preforms frequency-to-frequency conversions (TOPO to LSRK) in the frequency domain. Somewhere it is documented that this is only an approximation for v/c << 1. I have tested doing conversions in the velocity domain and attached the resulting image. It is slow for now but works.

The question is now how we interface to this in asap. One idea is to have an rc parameter which controls the accuracy, default being the current way and optional doing conversions in velocity. I have to think about alignment, where that is done.

Cheers, Malte.

[Malte Marquarding]

monr2 - conversion in velocity domain

[Max Voronkov]

I'd say the conversion in ASAP should be as accurate as possible all the time even if it takes longer to process the data. Otherwise, it would bite again at some stage. In Mon R2 case v/c<<1 is true, but the spectral resolution is high enough to be able to see the difference. The v which causes problems here is the velocity of LSRK with respect to TOPO frame (a few tens of km/s, depending on the direction).

I am surprised that the frequency-to-frequency conversion between these two frames is allowed at all in casacore without the rest frequency supplied (which is needed, strictly speaking, to account for the relative velocity correctly). The code probably uses the actual frequency represented by the measure instead of the rest frequency. I would prefer to supply a wrong rest frequency explicitly if this approximation is suitable, rather than assume it all the time.

[Malte Marquarding]

We will need a new function scantable.vel_align otherwise averaging will broaden the line due to the effects discussed above.

My current implementation for scantable.get_abscissa will do on-the-fly conversion to velocities in velocity. This is a little slow but should work. Simon, if you need to average the final data, then this will require a bit of work (implementing vel_align). If you are happy with just the abscissa being in velocity for display/export purposes, I can implement it very quickly. This will require you to specify an rc parameter at the start of the script and make sure you do NOT use freq_align or average(align=True)

[Max Voronkov]

Hi All,

It seems to me I have a better approach in mind. It is also an approximation, but a better approximation than the one used so far. We can perform an accurate conversion for a single spectral channel (e.g. the centre of the band), compute the difference with the old method using frequency-to-frequency conversion and then add this correction to all channels. The good thing of this approach is that it can be done in frequencies. Therefore, the old frequency aligner can be used (although the code has to be changed to calculate the offset per spectrum) and there will be no two different interfaces giving slightly different results (and therefore less confusion for users).

This approach uses the fact that we usually don't have wide bandwidth with high spectral resolution. If the bandwidth is too wide to bring high order terms into play this can't be used. Any comments?

Cheers,

Max

[Max Voronkov]

I did some back of the envelope calculations of the error introduced by the method suggested above (i.e. calculating the correction between frequency-based and velocity-based conversion for a single channel and adding it to all channels). The upper limit for this error in channels is Nv/2c, where N is the number of channels and v is the relative velocity between two frames involved in the conversion. This formula assumes that the correction is calculated for the centre of the band and the line of interest is near the edge. In the case of TOPO to LSRK correction the velocity is less than 60 km/s (worst case scenario of 30 km/s due to orbital motion and 30 km/s due to Sun's motion with respect to LSR, although I think these two vectors are never exactly parallel). For 4096 channels (Hobart and Mopra zoom mode) the error is less than 0.4 channels. It is also fine with broad band mode with its 8192 channels if BARY frame is used. However, if one uses a faster moving frame (e.g. local group or CMB), the correction can exceed one channel. It may not be a problem as usually the lines are wider in this case. Note that results doesn't seem to depend on the bandwidth or observed frequency.

A completely different approach may be to introduce non-linearity into spectral axis. However, it may have an impact on the frequency aligner code and needs some investigation.

[Max Voronkov]

Test of the casacore conversions between TOPO and LSRK

[Max Voronkov]

my e-mail to Wim is below for the trac record

Hi Wim,

First, I agree with your comment that the frequency conversion doesn't need a rest frequency. I guess the example sent by Malte shows largely that the conversion does a reasonably good job (and we have a topocentric frequency), but there is still some residual effect.

I made a different example (attached), which may clarify what we meant by the velocity route. The "IN" frequencies are just peak TOPO frequencies in the spectrum as observed by Hobart at 3 different epochs (no doppler tracking in any form, i.e. the same frequency of the first spectral channel). I think Malte fitted a Gaussian into a line profile for every epoch.

These values are converted to LSRK, where they ideally should give the same value for all epochs. The spectral resolution is present in the code just to report the difference in spectral channels (so we want it to be less than 1). We convert single spectral channels throughout the example (but you're right that if we converted a band, the channel increment would have to be adjusted as well). The example reports differences of 3.2 and 1.2 channels between the last and the first epochs, and between the first two epochs, respectively.

Then, the same TOPO frequencies were converted to LSRK using a different approach. A topocentric velocity is formed first and then it is converted to LSRK velocity using MRadialVelocity::Convert, which is then converted back to the frequency. In this case the differences are 0.96 and 0.56 channels, which is less than in the case of the pure frequency conversion. The error is still close to 1 channel, but we would need a higher spectral resolution data to make sure that it is not a result of the measurement. The conversion between velocity and frequency requires the rest frequency, which is hard coded in the example. I checked that the frequency differences obtained using the "velocity approach" don't seem to depend on its value as expected. I am a bit puzzled why going via velocity gives a bit better result (chance coincidence). I can't remember whether Malte mentioned this, but I ran the same data through miriad and saw the lines nicely aligned. We can compare miriad's formulae with what is used in casacore.

Cheers,

Max

The output I have:

In freqs:  [4.76503418, 4.764698242, 4.764930664]
Out freqs: [4765394923, 4765396095, 4765398091]
... via velocity: [4765396586, 4765397129, 4765397527]
Out vel: [10.4060516, 10.37191675, 10.34683492]
Difference 2-0
LSR via freq (chan): 3.244485066
LSR via vel (chan): 0.963878665
LSR velocities (km/s): -0.05921667541
TOPO (in):  -106.0000051
Difference 1-0
LSR via freq (chan): 1.200670038
LSR via vel (chan): 0.5556178965
LSR velocities (km/s): -0.03413484205
TOPO (in):  -344

[Max Voronkov]

a reply from Wim is below

Maxim, Simon,

After many manual checks I found the error in the frequency conversion. I used at one stage the wrong coordinate system (apparent) rather than J2000 by applying the forward annual aberration. Was not noted earlier, since it was written and tested around the year 2000, when these differences where small. Also the WSRT observations were ok, since they do an online shift incorporating aberrations using the on-line coordinates.

The frequency result is now identical to the 'using velocities' route, as it should be.

The remaining difference could be due by not taking transversal velocities into account. I will try to check all that when I get back from holidays in April.

Sorry it took so long; and thanks for Simon on finding it.

Wim

ps I checked the corrected MCFrequency.cc into casacore.