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