00001
00002
00003
00004
00005
00006
00007 #ifndef MOON_HXX
00008 #define MOON_HXX
00009
00010 extern "C"
00011 {
00012 #include "slalib.h"
00013 }
00014
00015 #include <cmath>
00016
00017 #include "DataMethod.hxx"
00018 #include "StringOps.h"
00019 #include "Options.h"
00020
00023 template<typename DataType>
00024 class Moon : public DataMethod<DataType>
00025 {
00026 public:
00028 Moon() { DataMethod<DataType>::methodName = "Moon"; }
00029
00031 virtual ~Moon() { }
00032
00034 void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00035 { mMJDObsCol = cols[0]; }
00036
00038 void setTarget(IntMap& cols, StringMap units)
00039 {
00040 mRACol = cols[0];
00041 mDecCol = cols[1];
00042 mIllumCol = cols[2];
00043 mTopoRaCol = cols[3];
00044 mTopoDecCol = cols[4];
00045 }
00046
00047 const double& getMoonRA() const { return mRAMoon; }
00048 const double& getMoonDec() const { return mDecMoon; }
00049 const double& getMoonPhase() const { return mIllum; }
00050 const double& getMoonRATopo() const { return mRAMoonTopo; }
00051 const double& getMoonDecTopo() const { return mDecMoonTopo; }
00052
00053 void calcMoonPosAndPhase(const double& aMjdObs)
00054 {
00055 Options options;
00056 const double dpi = M_PI;
00057 const double dd2r = dpi / 180.0;
00058 const double ds2r = dd2r / 3600;
00059 const double dr2d = 180.0 / dpi;
00060 const double ds2d = 1.0 / 86400.0;
00061 const double c_aups = 173.14 * ds2d;
00062 const double dr2h = 12.0 / dpi;
00063 const double dm2au = 1.0 / 149597870691.0;
00064 const double dau2m = 149597870691.0;
00065 double d0 = aMjdObs + 2400000.5 - 2451545.0;
00066 double deqx = 2000.0e0;
00067 double ea = 6378137.0;
00068 double latobs, lonobs, hobs;
00069
00070 if (StringOps::ToLower(options.getArchive()) == "vsa")
00071 {
00072 latobs = -24.61583333 * dd2r;
00073 lonobs = -70.3975 *dd2r;
00074 hobs = 2530.0;
00075 }
00076 else
00077 {
00078 latobs = 19.82561111 * dd2r;
00079 lonobs = -155.47322222 *dd2r;
00080 hobs = 4194.0;
00081 }
00082
00084 int leap = (aMjdObs < 53736) ? 32 : 33;
00085 double gRad = (357.53 + 0.9856003 * d0)*dd2r;
00086 double tt = aMjdObs + (leap + 32.184) * ds2d;
00087 double tdbobs = tt + (0.001658*sin(gRad) + 0.000014*sin(2.*gRad)) * ds2d;
00088
00090 double pvmoon[6];
00091 cslaDmoon(tdbobs, pvmoon);
00092 double xm = pvmoon[0] * dau2m;
00093 double ym = pvmoon[1] * dau2m;
00094 double zm = pvmoon[2] * dau2m;
00095 double raMoonRad = atan2(ym, xm);
00096 double decMoonRad = atan2(zm, sqrt(xm * xm + ym * ym));
00097 double distMoon = sqrt(xm * xm + ym * ym + zm * zm);
00098 mRAMoon = raMoonRad * dr2h;
00099 mDecMoon = decMoonRad * dr2d;
00100 if (mRAMoon < 0.0)
00101 {
00102 mRAMoon += 24.0;
00103 }
00104 else if (mRAMoon > 24.0)
00105 {
00106 mRAMoon -= 24.0;
00107 }
00108
00110 double dvb[3], dpb[3], dvh[3], dph[3];
00111 cslaEvp(tdbobs, deqx, dvb, dpb, dvh, dph);
00112 double pvsun[6];
00113 pvsun[0] = -dph[0];
00114 pvsun[1] = -dph[1];
00115 pvsun[2] = -dph[2];
00116 pvsun[3] = dvh[0];
00117 pvsun[4] = dvh[1];
00118 pvsun[5] = dvh[2];
00119 double xe = pvsun[0];
00120 double ye = pvsun[1];
00121 double ze = pvsun[2];
00122 double vxe = pvsun[3];
00123 double vye = pvsun[4];
00124 double vze = pvsun[5];
00125 double raSunRad = atan2(ye, xe);
00126 double decSunRad = atan2(ze, sqrt(xe * xe + ye * ye));
00127 double distSun = sqrt(xe * xe + ye * ye + ze * ze);
00128
00129
00130 double fd, djm;
00131 int iy, im, id, status;
00132 cslaDjcl(aMjdObs, &iy, &im, &id, &fd, &status);
00133 cslaCldj(iy-1, 12, 31, &djm, &status);
00134
00135 double t0 = (djm - 2415020.0) / 36525.0;
00136 double t = (aMjdObs - djm) / 36525.0;
00137 double t2 = t * t;
00138 double t3 = t2 * t;
00139
00140 double pa = ((2304.25 + 1.396 * t0) * t + 0.302 * t2 + 0.018 * t3) * ds2r;
00141 double pb = (pa + 0.791 * t2) * ds2r;
00142 double pc = ((2004.682 - 0.853 * t0) * t - 0.426 * t2 - 0.042 * t3) * ds2r;
00143
00144 double q = sin(pc) * (tan(decMoonRad) + tan(0.5*pc) * cos(raMoonRad + pa));
00145 double d_ra_p = pa + pb + atan2(q * sin(raMoonRad + pa),
00146 1.0 - q * cos(raMoonRad + pa));
00147 double d_dec_p = 2.0 * atan(
00148 tan(0.5 * pc) * (cos(raMoonRad + pa)
00149 - sin(raMoonRad + pa) * tan(0.5*(d_ra_p-pa-pb))));
00150
00151
00152 t = (aMjdObs - djm) / 36525.0;
00153 double omega = (259.183275 - 1934.142 * t) * dd2r;
00154 double d_long = -8.3597e-05 * sin(omega);
00155 double d_obliq = 4.4678e-05 * cos(omega);
00156 double d_ra_n = d_long * (0.917450511
00157 + 0.39784993 * sin(raMoonRad) * tan(decMoonRad))
00158 - cos(raMoonRad) * tan(decMoonRad) *d_obliq;
00159 double d_dec_n = 0.39784993 * cos(raMoonRad) * d_long
00160 + sin(raMoonRad) * d_obliq;
00161
00162
00163 double d_ra_a = (-vxe * sin(raMoonRad) + vye * cos(raMoonRad))
00164 / (c_aups * cos(decMoonRad));
00165 double d_dec_a = (-vxe * cos(raMoonRad) * sin(decMoonRad)
00166 - vye * sin(raMoonRad) * sin(decMoonRad)
00167 + vze * cos(decMoonRad)) / c_aups;
00168
00169 double raMoonApp = raMoonRad + d_ra_p + d_ra_n + d_ra_a;
00170 double decMoonApp = decMoonRad + d_dec_p + d_dec_n + d_dec_a;
00171
00172
00173 double sinphi = sin(latobs);
00174 double cosphi = cos(latobs);
00175 double gast = cslaGmst(aMjdObs) + cslaEqeqx(tdbobs) + lonobs;
00176 double robs = ea + hobs;
00177 double xobs = robs * cosphi * cos(gast);
00178 double yobs = robs * cosphi * sin(gast);
00179 double zobs = robs * sinphi;
00180
00181
00182 double xmtopo = distMoon * cos(decMoonApp) * cos(raMoonApp) - xobs;
00183 double ymtopo = distMoon * cos(decMoonApp) * sin(raMoonApp) - yobs;
00184 double zmtopo = distMoon * sin(decMoonApp) - zobs;
00185 double rmtopo = sqrt(xmtopo * xmtopo + ymtopo * ymtopo + zmtopo * zmtopo);
00186 double raMoonTopoRad = atan2(ymtopo, xmtopo);
00187 double decMoonTopoRad = asin(zmtopo / rmtopo);
00188 mRAMoonTopo = raMoonTopoRad * dr2h;
00189 mDecMoonTopo = decMoonTopoRad * dr2d;
00190 if (mRAMoonTopo < 0.0)
00191 {
00192 mRAMoonTopo += 24.0;
00193 }
00194 if (mRAMoonTopo > 24.0)
00195 {
00196 mRAMoonTopo -= 24.0;
00197 }
00198
00200 double psi = acos(sin(decSunRad)*sin(decMoonRad)
00201 + cos(decSunRad)*cos(decMoonRad)*cos(raSunRad-raMoonRad));
00202 double phase = atan2(distSun * sin(psi),
00203 distMoon*dm2au - distSun * cos(psi));
00204 mIllum = 0.5 * (1.0 + cos(phase));
00205 }
00206
00207 void doit(TableData<DataType>& data, int row1, int row2)
00208 {
00209 for (int rowNo = row1; rowNo < row2; ++rowNo)
00210 {
00212 double mjdobs;
00213 data.value(mMJDObsCol, rowNo, mjdobs);
00214
00215 calcMoonPosAndPhase(mjdobs);
00216
00218 data.assign(mRACol, rowNo, mRAMoon);
00219 data.assign(mDecCol, rowNo, mDecMoon);
00220 data.assign(mTopoRaCol, rowNo, mRAMoonTopo);
00221 data.assign(mTopoDecCol, rowNo, mDecMoonTopo);
00222 data.assign(mIllumCol, rowNo, mIllum);
00223 }
00224 }
00225
00226
00227 void doit2(TableData<DataType>& rawData, TableData<DataType>& data, int row1, int row2) { }
00228
00229 private:
00230 int mMJDObsCol;
00231 int mIllumCol;
00232 int mRACol;
00233 int mDecCol;
00234 int mTopoRaCol;
00235 int mTopoDecCol;
00236 double mIllum;
00237 double mRAMoon;
00238 double mDecMoon;
00239 double mRAMoonTopo;
00240 double mDecMoonTopo;
00241 };
00242
00243 #endif
00244
00245
00246
00247
00248
00249
00250
00251
00252
