Ticket #150: dTOPO2LSRK2.cc

File dTOPO2LSRK2.cc, 4.1 KB (added by Max Voronkov, 17 years ago)

Test of the casacore conversions between TOPO and LSRK

Line 
1#include <casa/aips.h>
2#include <casa/iostream.h>
3#include <casa/iomanip.h>
4
5#include <measures/Measures.h>
6#include <measures/Measures/MeasFrame.h>
7#include <measures/Measures/MCEpoch.h>
8#include <measures/Measures/MCDirection.h>
9#include <measures/Measures/MCPosition.h>
10#include <measures/Measures/MCFrequency.h>
11#include <measures/Measures/MCRadialVelocity.h>
12#include <measures/Measures/MCDoppler.h>
13
14
15#include <casa/Exceptions/Error.h>
16
17using namespace casa;
18
19int main() {
20 try {
21 // Observed peak topocentric frequencies of the maser at 3 given epochs
22 Vector<Double> freqs(3);
23 freqs(0) = 4.76503417969;
24 freqs(1) = 4.76469824219;
25 freqs(2) = 4.76493066406;
26 // epochs of observations (mjd)
27 Vector<Double> epochs(3);
28 epochs(0) = 54468.493923611109;
29 epochs(1) = 54542.278935185182;
30 epochs(2) = 54614.149143518516;
31
32 // rest frequency in Hz (needed for velocity-based conversion)
33 const Double restfreq = 4.765562e9;
34
35 // buffer for LSRK frequencies and velocity
36 Vector<Double> outfreqs(3);
37 Vector<Double> outvels(3);
38 Vector<Double> outfreqs_viavel(3);
39
40
41 // The frequency axis increment in Hz, given here just to show the error in terms of the
42 // spectral channel width. The main reason for this is to give a right feeling how big the
43 // error is. All conversions are done for single frequencies.
44 Double deltaf = 976.5625;
45
46 // Hobart using VLBI position
47 MPosition obs(MVPosition(-3950236.7341,
48 2522347.553,
49 -4311562.5434),MPosition::ITRF);
50
51 // Monr2 maser at 4.7 GHz checked against literature
52 MDirection coord(Quantity(1.6048229750035694,"rad"),
53 Quantity(-0.11139370025381365,"rad"),
54 MDirection::J2000);
55
56 for (uInt i=0; i < epochs.nelements(); ++i) {
57 //UTC epoch varies
58 MEpoch epo0(Quantity(epochs(i),"d"));
59 MeasFrame frame(obs,coord,epo0);
60 MFrequency mf0(Quantity(freqs(i), "GHz"), MFrequency::TOPO);
61
62 // convert TOPO frequency to LSRK
63 MFrequency mfout0 = MFrequency::Convert(mf0,
64 MFrequency::Ref(MFrequency::LSRK,
65 frame))();
66 outfreqs[i] = mfout0.get(Unit("Hz")).getValue();
67
68 // convert TOPO frequency to TOPO velocity
69 MRadialVelocity mv0 = MRadialVelocity::fromDoppler(mf0.toDoppler(restfreq),MRadialVelocity::TOPO);
70
71 // convert TOPO velocity to LSRK velocity
72 MRadialVelocity mvout0 = MRadialVelocity::Convert(mv0,
73 MRadialVelocity::Ref(MRadialVelocity::LSRK, frame))();
74 outvels[i] = mvout0.get("km/s").getValue();
75
76 // convert LSRK velocity to LSRK frequency
77 mfout0 = MFrequency::fromDoppler(mvout0.toDoppler(), restfreq, MFrequency::LSRK);
78 outfreqs_viavel[i] = mfout0.get("Hz").getValue();
79
80 }
81 cout << "In freqs: " << setprecision(10) << freqs << endl;
82 cout << "Out freqs: " << setprecision(10) << outfreqs << endl;
83 cout << "... via velocity: " << setprecision(10) << outfreqs_viavel << endl;
84 cout << "Out vel: " << setprecision(10) << outvels << endl;
85
86
87 // This should be < 1.0
88 cout << "Difference 2-0" << endl;
89 cout << "LSR via freq (chan): " << (outfreqs[2] - outfreqs[0])/deltaf << endl;
90 cout << "LSR via vel (chan): " << (outfreqs_viavel[2] - outfreqs_viavel[0])/deltaf << endl;
91 cout << "LSR velocities (km/s): " << (outvels[2] - outvels[0]) << endl;
92 cout << "TOPO (in): " << (freqs[2] - freqs[0])*1e9/deltaf << endl;
93 cout << "Difference 1-0" << endl;
94 cout << "LSR via freq (chan): " << (outfreqs[1] - outfreqs[0])/deltaf << endl;
95 cout << "LSR via vel (chan): " << (outfreqs_viavel[1] - outfreqs_viavel[0])/deltaf << endl;
96 cout << "LSR velocities (km/s): " << (outvels[1] - outvels[0]) << endl;
97 cout << "TOPO (in): " << (freqs[1] - freqs[0])*1e9/deltaf << endl;
98 } catch (AipsError &e) {
99 cout << e.what() << endl;
100 }
101 return 0;
102}