| 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 | |
|---|
| 12 | #include <casa/Exceptions/Error.h> |
|---|
| 13 | |
|---|
| 14 | using namespace casa; |
|---|
| 15 | |
|---|
| 16 | int main() { |
|---|
| 17 | try { |
|---|
| 18 | // observation frequencies at the peak of the maser for three epochs |
|---|
| 19 | Vector<Double> freqs(3); |
|---|
| 20 | freqs(0) = 4.76503417969; |
|---|
| 21 | freqs(1) = 4.76469824219; |
|---|
| 22 | freqs(2) = 4.76493066406; |
|---|
| 23 | // the corresponding epochs |
|---|
| 24 | Vector<Double> epochs(3); |
|---|
| 25 | epochs(0) = 54468.493923611109; |
|---|
| 26 | epochs(1) = 54542.278935185182; |
|---|
| 27 | epochs(2) = 54614.149143518516; |
|---|
| 28 | |
|---|
| 29 | Vector<Double> outfreqs(3); |
|---|
| 30 | |
|---|
| 31 | // Hobart using VLBI position |
|---|
| 32 | MPosition obs(MVPosition(-3950236.7341, |
|---|
| 33 | 2522347.553, |
|---|
| 34 | -4311562.5434),MPosition::ITRF); |
|---|
| 35 | |
|---|
| 36 | // Monr2 maser at 4.7 GHz checked against literature |
|---|
| 37 | MDirection coord(Quantity(1.6048229750035694,"rad"), |
|---|
| 38 | Quantity(-0.11139370025381365,"rad"), |
|---|
| 39 | MDirection::J2000); |
|---|
| 40 | for (uInt i=0; i < freqs.nelements(); ++i) { |
|---|
| 41 | //UTC epoch varies |
|---|
| 42 | MEpoch epo0(Quantity(epochs(i),"d")); |
|---|
| 43 | MeasFrame frame(obs,coord,epo0); |
|---|
| 44 | MFrequency mf0(Quantity(freqs(i), "GHz"), MFrequency::TOPO); |
|---|
| 45 | MFrequency mfout0 = MFrequency::Convert(mf0, |
|---|
| 46 | MFrequency::Ref(MFrequency::LSRK, |
|---|
| 47 | frame))(); |
|---|
| 48 | outfreqs[i] = mfout0.get(Unit("Hz")).getValue(); |
|---|
| 49 | } |
|---|
| 50 | cout << setprecision(10) << outfreqs << endl; |
|---|
| 51 | // The frequency axis increment |
|---|
| 52 | Double deltaf = 976.5625; |
|---|
| 53 | // This should be < 1.0 |
|---|
| 54 | cout << (outfreqs[2] - outfreqs[0])/deltaf << endl; |
|---|
| 55 | } catch (AipsError &e) { |
|---|
| 56 | cout << e.what() << endl; |
|---|
| 57 | } |
|---|
| 58 | return 0; |
|---|
| 59 | } |
|---|
