| 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 | }
|
|---|
