Index: trunk/FACT++/src/calcsource.cc =================================================================== --- trunk/FACT++/src/calcsource.cc (revision 19053) +++ trunk/FACT++/src/calcsource.cc (revision 19053) @@ -0,0 +1,882 @@ +#include + +#include "Database.h" + +#include "pal.h" +#include "nova.h" +#include "Time.h" +#include "Configuration.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +// ------------------------------------------------------------------------ + +struct Map : pair +{ + Map() { } +}; + +std::istream &operator>>(std::istream &in, Map &m) +{ + string txt; + in >> txt; + + const auto p = txt.find_first_of('/'); + if (p==string::npos) + throw runtime_error("Missing / in map argument."); + if (p!=txt.find_last_of('/')) + throw runtime_error("Only one / allowed in map argument."); + + m.first = txt.substr(0, p); + m.second = txt.substr(p+1); + + return in; +} + + +void SetupConfiguration(Configuration &conf) +{ + po::options_description control("Root to Database"); + control.add_options() + ("uri,u", var()->required(), "Database link as in\n\tuser:password@server[:port]/database.") + ("ra", var()/*5.575539)*/, "") + ("dec", var()/*22.014500)*/, "") + ("focal-dist", var(4889.), "") + //("source-key", var(5), "") + //("source-name", var(""), "") + ("file", var(171011115)->required(), "") + ("drop", po_switch(), "Drop the table (implied create)") + ("list-files", po_switch(), "") + ("tree,t", var("Events"), "Name of the root tree to convert") + ("table", var("Events"), "") + ("update", var(""), "") + ("create", po_switch(), "") + ("verbose,v", var(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)") + ("no-insert", po_switch(), "Does not insert any data to the table") + //("print-insert", po_switch(), "Print the INSERT query (note that it contains all data)") + //("print-create", po_switch(), "Print the CREATE query") + ; + + po::positional_options_description p; + p.add("file", 1); // The 1st positional options (n=1) + + conf.AddOptions(control); + conf.SetArgumentPositions(p); +} + +void PrintUsage() +{ + cout << + "calcsource - Fill a table with source positions\n" +// "\n" +// "Usage: rootifysql [rootify.sql [rootify.root]] [-u URI] [-q query|-f file] [-i] [-o out] [-f] [-cN] [-t tree] [-vN]\n" + "\n" + ; + cout << endl; +} + +enum BasicType_t +{ + kNone = 0, + kFloat, + kDouble, + kInt16, + kUInt16, + kInt32, + kUInt32, + kInt64, + kUInt64, +}; + +static const map> ConvRoot = +{ + { "Float_t", { kFloat, "FLOAT" } }, + { "Double_t", { kDouble, "DOUBLE" } }, + { "ULong64_t", { kUInt64, "BIGINT UNSIGNED" } }, + { "Long64_t", { kInt64, "BIGINT" } }, + { "UInt_t", { kUInt32, "INT UNSIGNED" } }, + { "Int_t", { kInt32, "INT" } }, + { "UShort_t", { kUInt16, "SMALLINT UNSIGNED" } }, + { "Short_t", { kInt16, "SMALLINT" } }, +}; + +struct Container +{ + string branch; // branch name + string column; // column name + BasicType_t type; + void *ptr; + + Container(const string &b, const string &c, const BasicType_t &t) : branch(b), column(c), type(t), ptr(0) + { + switch (t) + { + case kFloat: ptr = new Float_t; break; + case kDouble: ptr = new Double_t; break; + case kInt16: ptr = new Short_t; break; + case kUInt16: ptr = new UShort_t; break; + case kInt32: ptr = new Int_t; break; + case kUInt32: ptr = new UInt_t; break; + case kInt64: ptr = new Long64_t; break; + case kUInt64: ptr = new ULong64_t; break; + case kNone: + break; + } + } + ~Container() + { + //::operator delete(ptr); // It seems root is deleting it already + } + + string fmt() const + { + ostringstream str; + + switch (type) + { + case kFloat: str << setprecision(8) << *reinterpret_cast(ptr); break; + case kDouble: str << setprecision(16) << *reinterpret_cast(ptr); break; + case kInt16: str << *reinterpret_cast(ptr); break; + case kUInt16: str << *reinterpret_cast(ptr); break; + case kInt32: str << *reinterpret_cast(ptr); break; + case kUInt32: str << *reinterpret_cast(ptr); break; + case kInt64: str << *reinterpret_cast(ptr); break; + case kUInt64: str << *reinterpret_cast(ptr); break; + case kNone: + break; + } + + return str.str(); + } +}; + +class MRotation : public TRotation +{ +public: + MRotation() : TRotation(1, 0, 0, 0, -1, 0, 0, 0, 1) + { + } +}; + + +void ErrorHandlerAll(Int_t level, Bool_t abort, const char *location, const char *msg) +{ + if (string(msg).substr(0,24)=="no dictionary for class ") + return; + + DefaultErrorHandler(level, abort, location, msg); +} + +int main(int argc, const char* argv[]) +{ + Time start; + + gROOT->SetBatch(); + SetErrorHandler(ErrorHandlerAll); + + Configuration conf(argv[0]); + conf.SetPrintUsage(PrintUsage); + SetupConfiguration(conf); + + if (!conf.DoParse(argc, argv)) + return 127; + + // ----------------------------- Evaluate options -------------------------- + const bool has_radec = conf.Has("ra") && conf.Has("dec"); + + const string uri = conf.Get("uri"); + const uint32_t file = conf.Get("file"); + const string tree = conf.Get("tree"); + const string table = conf.Get("table"); + //string source_name = conf.Get("source-name"); + //uint32_t source_key = conf.Has("source-key") ? conf.Get("source-key") : 0; + double source_ra = conf.Has("ra") ? conf.Get("ra") : 0; + double source_dec = conf.Has("dec") ? conf.Get("dec") : 0; + double focal_dist = conf.Get("focal-dist"); + //const bool print_create = conf.Get("print-create"); + //const bool print_insert = conf.Get("print-insert"); + const bool drop = conf.Get("drop"); + const bool create = conf.Get("create") || drop; + const string update = conf.Get("update"); + const bool list_files = conf.Get("list-files"); + const uint16_t verbose = conf.Get("verbose"); + + // ------------------------------------------------------------------------- + + if (list_files) + { + const string query = + "SELECT FileId FROM Events GROUP BY FileId"; + + cout << query << endl; + + const mysqlpp::StoreQueryResult res = + Database(uri).query(query).store(); + + for (size_t i=0; i0) + cout << "\n------------------------- Evaluating source ------------------------" << endl; + + cout << "Requesting coordinates from RunInfo/Source table for 20" << file/1000 << "/" << file%1000 << endl; + + if (!has_radec) + { + const string query = + "SELECT Source.fRightAscension, Source.fDeclination, Source.fSourceName" + " FROM RunInfo" + " LEFT JOIN Source" + " USING (fSourceKey)" + " WHERE fNight=20"+to_string(file/1000)+ + " AND fRunID="+to_string(file%1000); + + cout << query << endl; + + const mysqlpp::StoreQueryResult res = + Database(uri).query(query).store(); + + if (res.num_rows()!=1) + { + cerr << "No coordinates from RunInfo for " << file << endl; + return 1; + } + + source_ra = res[0][0]; + source_dec = res[0][1]; + + cout << "Using coordinates " << source_ra << "h / " << source_dec << " deg for " << res[0][2] << endl; + } + else + cout << "Using coordinates " << source_ra << "h / " << source_dec << " deg from resources." << endl; + +/* + if (!source_name.empty()) + { + cout << "Requesting coordinates for '" << source_name << "'" << endl; + + const mysqlpp::StoreQueryResult res = + Database(uri).query("SELECT `Ra`, `Dec` WHERE fSourceName='"+source_name+"'").store(); + + if (res.num_rows()!=1) + { + cerr << "No " << (res.num_rows()>1?"unique ":"") << "coordinates found for '" << source_name << "'" << endl; + return 1; + } + + source_ra = res[0][0]; + source_dec = res[0][1]; + } +*/ + + if (verbose>0) + cout << "\n-------------------------- Evaluating file ------------------------" << endl; + + + Database connection(uri); // Keep alive while fetching rows + + const string query = + "SELECT `Ra`, `Dec`, MJD, MilliSec, NanoSec, Zd, Az, EvtNumber" + " FROM `"+table+"`" + " WHERE FileId="+to_string(file); + + cout << query << endl; + + const mysqlpp::UseQueryResult res1 = + connection.query(query).use(); + + Nova::RaDecPosn source(source_ra, source_dec); + + source_ra *= M_PI/12; + source_dec *= M_PI/180; + + auto obs = Nova::kORM; + + obs.lng *= M_PI/180; + obs.lat *= M_PI/180; + + //const double mm2deg = 1.17193246260285378e-02; + + ostringstream ins; + ins << setprecision(16); + + ostringstream upd; + upd << setprecision(16); + + size_t count = 0; + while (auto row=res1.fetch_row()) + { + count++; + + double point_ra = row[0]; + double point_dec = row[1]; + uint32_t mjd = row[2]; + int64_t millisec = row[3]; + //uint32_t nanosec = row[4]; + //double zd = row[5]; + //double az = row[6]; + uint32_t event = row[7]; + + Nova::RaDecPosn point(point_ra, point_dec); + + point_ra *= M_PI/12; + point_dec *= M_PI/180; + + /* + // ============================ Mars ================================ + + TVector3 pos; // pos: source position + TVector3 pos0; // pos: source position + + pos.SetMagThetaPhi(1, M_PI/2-source_dec, source_ra); + pos0.SetMagThetaPhi(1, M_PI/2-point_dec, point_ra); + + const double ut = (nanosec/1e6+millisec)/(24*3600000); + + // Julian centuries since J2000. + const double t = (ut -(51544.5-mjd)) / 36525.0; + + // GMST at this UT1 + const double r1 = 24110.54841+(8640184.812866+(0.093104-6.2e-6*t)*t)*t; + const double r2 = 86400.0*ut; + + const double sum = (r1+r2)/(3600*24); + + double gmst = fmod(sum, 1) * 2*M_PI; + + MRotation conv; + conv.RotateZ(gmst + obs.lng); + conv.RotateY(obs.lat-M_PI/2); + conv.RotateZ(M_PI); + + pos *= conv; + pos0 *= conv; + + pos.RotateZ(-pos0.Phi()); + pos.RotateY(-pos0.Theta()); + pos.RotateZ(-M_PI/2); // exchange x and y + pos *= -focal_dist/pos.Z(); + + TVector2 v = pos.XYvector(); + + //if (fDeviation) + // v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg(); + + //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n'; + */ + + // ================================= Nova =============================== + + Nova::ZdAzPosn ppos = Nova::GetHrzFromEqu(source, 2400000.5+mjd+millisec/1000./3600/24); + Nova::ZdAzPosn ppos0 = Nova::GetHrzFromEqu(point, 2400000.5+mjd+millisec/1000./3600/24); + + TVector3 pos; + TVector3 pos0; + pos.SetMagThetaPhi( 1, ppos.zd *M_PI/180, ppos.az *M_PI/180); + pos0.SetMagThetaPhi(1, ppos0.zd*M_PI/180, ppos0.az*M_PI/180); + + pos.RotateZ(-pos0.Phi()); + pos.RotateY(-pos0.Theta()); + pos.RotateZ(-M_PI/2); // exchange x and y + pos *= -focal_dist/pos.Z(); + + TVector2 v = pos.XYvector(); + + //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n'; + + /* + // =============================== Slalib ========================== + const double height = 2200; + const double temp = 10; + const double hum = 0.25; + const double press = 780; + + const double _mjd = mjd+millisec/1000./3600/24; + + const double dtt = palDtt(_mjd); // 32.184 + 35 + + const double tdb = _mjd + dtt/3600/24; + const double dut = 0; + + // prepare calculation: Mean Place to geocentric apperent + // (UTC would also do, except for the moon?) + double fAmprms[21]; + palMappa(2000.0, tdb, fAmprms); // Epoche, TDB + + // prepare: Apperent to observed place + double fAoprms[14]; + palAoppa(_mjd, dut, // mjd, Delta UT=UT1-UTC + obs.lng, obs.lat, height, // long, lat, height + 0, 0, // polar motion x, y-coordinate (radians) + 273.155+temp, press, hum, // temp, pressure, humidity + 0.40, 0.0065, // wavelength, tropo lapse rate + fAoprms); + + // ---- Mean to apparent ---- + double r=0, d=0; + palMapqkz(point_ra, point_dec, fAmprms, &r, &d); + + double _zd, _az, ha, ra, dec; + // -- apparent to observed -- + palAopqk(r, d, fAoprms, + &_az, // observed azimuth (radians: N=0,E=90) [-pi, pi] + &_zd, // observed zenith distance (radians) [-pi/2, pi/2] + &ha, // observed hour angle (radians) + &dec, // observed declination (radians) + &ra); // observed right ascension (radians) + + //cout << _zd*180/M_PI << " " << _az*180/M_PI << endl; + + pos0.SetMagThetaPhi(1, _zd, _az); + + r=0, d=0; + palMapqkz(source_ra, source_dec, fAmprms, &r, &d); + + // -- apparent to observed -- + palAopqk(r, d, fAoprms, + &_az, // observed azimuth (radians: N=0,E=90) [-pi, pi] + &_zd, // observed zenith distance (radians) [-pi/2, pi/2] + &ha, // observed hour angle (radians) + &dec, // observed declination (radians) + &ra); // observed right ascension (radians) + + pos.SetMagThetaPhi(1, _zd, _az); + + //cout << _zd*180/M_PI << " " << _az*180/M_PI << endl; + + pos.RotateZ(-pos0.Phi()); + pos.RotateY(-pos0.Theta()); + pos.RotateZ(-M_PI/2); // exchange x and y + pos *= -focal_dist/pos.Z(); + + v = pos.XYvector(); + + //cout << v.X() << " " << v.Y() << " " << v.Mod()*mm2deg << '\n'; + */ + + if (1/*insert*/) + ins << "( " << file << ", " << event << ", " << v.X() << ", " << v.Y() << " ),\n"; + + if (!update.empty()) + upd << "UPDATE " << update << " SET X=" << v.X() << ", Y=" << v.Y() << " WHERE FileId=" << file << " AND EvtNumber=" << event <<";\n"; + }; + + if (connection.errnum()) + { + cerr << "SQL error fetching row: " << connection.error() << endl; + return 7; + } + + + if (drop) + { + cout << "Drop table Position." << endl; + const mysqlpp::SimpleResult res2 = + Database(uri).query("DROP TABLE Position").execute(); + + //cout << res.rows() << " rows affected." << res.info() << endl; + } + + + if (create) + { + cout << "Create table Position." << endl; + const mysqlpp::SimpleResult res3 = + Database(uri).query("CREATE TABLE IF NOT EXISTS Position\n" + "(\n" + " FileId INT UNSIGNED NOT NULL,\n" + " EvtNumber INT UNSIGNED NOT NULL,\n" + " X FLOAT NOT NULL,\n" + " Y FLOAT NOT NULL,\n" + " PRIMARY KEY (FileId, EvtNumber)\n" + ")\n" + "DEFAULT CHARSET=latin1 COLLATE=latin1_general_ci\n" + "ENGINE=MyISAM\n" + "COMMENT='created by root2db'\n").execute(); + + //cout << res.rows() << " rows affected." << res.info() << endl; + + } + else + { + // FIXME: Only if entries exist? + cout << "Delete old entries from Position." << endl; + const mysqlpp::SimpleResult res3 = + Database(uri).query("DELETE FROM Position WHERE FileId="+to_string(file)).execute(); + + cout << res3.rows() << " rows affected." << endl; + } + + if (1/*insert*/) + { + cout << "Insert data into table Position." << endl; + const string query3 = + "INSERT `Position` (FileId, EvtNumber, X, Y) VALUES\n"+ + ins.str().substr(0, ins.str().size()-2)+ + "\n"; + + if (0/*print_insert*/) + cout << query3 << endl; + + const mysqlpp::SimpleResult res3 = + Database(uri).query(query3).execute(); + + cout << res3.info() << endl; + } + + if (!update.empty()) + { + cout << upd.str() << endl; + const mysqlpp::StoreQueryResult res3 = + connection.query(upd.str()).store(); + + /* + mysqlpp::Connection con(db, server, user, pass); + // Set option to allow multiple queries to be issued. + mysqlpp::MultiStatementsOption* opt = new mysqlpp::MultiStatementsOption(true); + con.set_option(opt); + // Build either single queries or multiple queries strung together. + mysqlpp::Query query = con.query(); + .... + // Issue query + query.store(); + // Do NOT delete 'opt'; it will be destroyed by the 'con' object when it + // falls out of scope. + */ + } + + //cout << ins.str().substr(0, ins.str().size()-2) << endl; + //cout << count << endl; + + return 0; +} + +/* + TRotation & TRotation::RotateX(Double_t a) { + //rotate around x + Double_t c = TMath::Cos(a); + Double_t s = TMath::Sin(a); + Double_t x = fyx, y = fyy, z = fyz; + fyx = c*x - s*fzx; + fyy = c*y - s*fzy; + fyz = c*z - s*fzz; + fzx = s*x + c*fzx; + fzy = s*y + c*fzy; + fzz = s*z + c*fzz; + return *this; +} + +TRotation & TRotation::RotateY(Double_t a){ + //rotate around y + Double_t c = TMath::Cos(a); + Double_t s = TMath::Sin(a); + Double_t x = fzx, y = fzy, z = fzz; + fzx = c*x - s*fxx; + fzy = c*y - s*fxy; + fzz = c*z - s*fxz; + fxx = s*x + c*fxx; + fxy = s*y + c*fxy; + fxz = s*z + c*fxz; + return *this; +} + +TRotation & TRotation::RotateZ(Double_t a) { + //rotate around z + Double_t c = TMath::Cos(a); + Double_t s = TMath::Sin(a); + Double_t x = fxx, y = fxy, z = fxz; + fxx = c*x - s*fyx; + fxy = c*y - s*fyy; + fxz = c*z - s*fyz; + fyx = s*x + c*fyx; + fyy = s*y + c*fyy; + fyz = s*z + c*fyz; + return *this; +} + + */ + + +// Reuired: +// pointing_ra[8] +// pointing_dec[8] +// fNanoSec[4], fTime[fMiliSec,8], fMjd[4] +// source_ra -> rc +// source_dec -> rc +// fLong, fLat -> rc +// fCameraDist -> rc +// 32 byte per row, 1 Monat = 4mio rows => 120 MB + +/* + +void TVector3::SetThetaPhi(Double_t theta, Double_t phi) +{ + //setter with mag, theta, phi + fX = TMath::Sin(theta) * TMath::Cos(phi); + fY = TMath::Sin(theta) * TMath::Sin(phi); + fZ = TMath::Cos(theta); +} + + + // Set Sky coordinates of source, taken from container "MSourcePos" + // of type MPointingPos. The sky coordinates must be J2000, as the + // sky coordinates of the camera center that we get from the container + // "MPointingPos" filled by the Drive. + //MVector3 pos; // pos: source position + //pos.SetRaDec(fSourcePos->GetRaRad(), fSourcePos->GetDecRad()); + + TVector3 pos; // pos: source position + pos.SetMagThetaPhi(1, TMath::Pi()/2-source_dec, source_ra); + + // Set sky coordinates of camera center in pos0 (for details see below) + //MVector3 pos0; // pos0: camera center + //pos0.SetRaDec(fPointPos->GetRaRad(), fPointPos->GetDecRad()); + + TVector3 pos0; // pos: source position + pos0.SetMagThetaPhi(1, TMath::Pi()/2-pointing_dec, pointing_ra); + + // Convert sky coordinates of source to local coordinates. Warning! These are not the "true" local + // coordinates, since this transformation ignores precession and nutation effects. + + //const MAstroSky2Local conv(*fTime, *fObservatory); + //pos *= conv; + + const Double_t ut = (Double_t)(fNanoSec/1e6+(Long_t)fTime)/kDay; + + // Julian centuries since J2000. + const Double_t t = (ut -(51544.5-fMjd)) / 36525.0; + + // GMST at this UT1 + const Double_t r1 = 24110.54841+(8640184.812866+(0.093104-6.2e-6*t)*t)*t; + const Double_t r2 = 86400.0*ut; + + const Double_t sum = (r1+r2)/kDaySec; + + double gmst = fmod(sum, 1)*TMath::TwoPi(); + + // TRotation::TRotation(Double_t mxx, Double_t mxy, Double_t mxz, + // Double_t myx, Double_t myy, Double_t myz, + // Double_t mzx, Double_t mzy, Double_t mzz) + // : fxx(mxx), fxy(mxy), fxz(mxz), + // fyx(myx), fyy(myy), fyz(myz), + // fzx(mzx), fzy(mzy), fzz(mzz) {} + + TRotation conv(1, 0, 0, 0, -1, 0, 0, 0, 1); + conv.RotateZ(gmst + obs.GetElong()); + conv.RotateY(obs.GetPhi()-TMath::Pi()/2); + conv.RotateZ(TMath::Pi()); + + // TVector3 & TVector3::operator *= (const TRotation & m){ + // return *this = m * (*this); + // } + + // inline TVector3 TRotation::operator * (const TVector3 & p) const { + // return TVector3(fxx*p.X() + fxy*p.Y() + fxz*p.Z(), + // fyx*p.X() + fyy*p.Y() + fyz*p.Z(), + // fzx*p.X() + fzy*p.Y() + fzz*p.Z()); + // } + + // Convert sky coordinates of camera center convert to local. + // Same comment as above. These coordinates differ from the true + // local coordinates of the camera center that one could get from + // "MPointingPos", calculated by the Drive: the reason is that the Drive + // takes into account precession and nutation corrections, while + // MAstroSky2Local (as of Jan 27 2005 at least) does not. Since we just + // want to get the source position on the camera from the local + // coordinates of the center and the source, it does not matter that + // the coordinates contained in pos and pos0 ignore precession and + // nutation... since the shift would be the same in both cases. What + // would be wrong is to set in pos0 directly the local coordinates + // found in MPointingPos! + pos0 *= conv; + + if (fDeviation) + { + // Position at which the starguider camera is pointing in real: + // pointing position = nominal position - dev + // + //vx = CalcXYinCamera(pos0, pos)*fGeom->GetCameraDist()*1000; + pos0.SetZdAz(pos0.Theta()-fDeviation->GetDevZdRad(), + pos0.Phi() -fDeviation->GetDevAzRad()); + } + + // Calculate source position in camera, and convert to mm: + TVector2 v = MAstro::GetDistOnPlain(pos0, pos, -fGeom->GetCameraDist()*1000); + + pos.RotateZ(-pos0.Phi()); + pos.RotateY(-pos0.Theta()); + pos.RotateZ(-TMath::Pi()/2); // exchange x and y + pos *= -fGeom->GetCameraDist()*1000/pos.Z(); + + TVector2 v = pos.XYvector(); + + if (fDeviation) + v -= fDeviation->GetDevXY()/fGeom->GetConvMm2Deg(); + + SetSrcPos(v); + + // ================================================ + + if (fMode==kWobble) + { + // The trick here is that the anti-source position in case + // of the off-source regions is always the on-source positon + // thus a proper anti-source cut is possible. + fSrcPosAnti->SetXY(v); + if (fCallback) + v = Rotate(v, fCallback->GetNumPass(), fCallback->GetNumPasses()); + else + v *= -1; + fSrcPosCam->SetXY(v); + } + else + { + // Because we don't process this three times like in the + // wobble case we have to find another way to determine which + // off-source region is the right anti-source position + // We do it randomly. + fSrcPosCam->SetXY(v); + if (fNumRandomOffPositions>1) + v = Rotate(v, gRandom->Integer(fNumRandomOffPositions), fNumRandomOffPositions); + else + v *= -1; + fSrcPosAnti->SetXY(v); + } + + + + */ + + + + +/* + PointingData CalcPointingPos(const PointingSetup &setup, double _mjd, const Weather &weather, uint16_t timeout, bool tpoint=false) + { + PointingData out; + out.mjd = _mjd; + + const double elong = Nova::kORM.lng * M_PI/180; + const double lat = Nova::kORM.lat * M_PI/180; + const double height = 2200; + + const bool valid = weather.time+boost::posix_time::seconds(timeout) > Time(); + + const double temp = valid ? weather.temp : 10; + const double hum = valid ? weather.hum : 0.25; + const double press = valid ? weather.press : 780; + + const double dtt = palDtt(_mjd); // 32.184 + 35 + + const double tdb = _mjd + dtt/3600/24; + const double dut = 0; + + // prepare calculation: Mean Place to geocentric apperent + // (UTC would also do, except for the moon?) + double fAmprms[21]; + palMappa(2000.0, tdb, fAmprms); // Epoche, TDB + + // prepare: Apperent to observed place + double fAoprms[14]; + palAoppa(_mjd, dut, // mjd, Delta UT=UT1-UTC + elong, lat, height, // long, lat, height + 0, 0, // polar motion x, y-coordinate (radians) + 273.155+temp, press, hum, // temp, pressure, humidity + 0.40, 0.0065, // wavelength, tropo lapse rate + fAoprms); + + out.source.ra = setup.source.ra * M_PI/ 12; + out.source.dec = setup.source.dec * M_PI/180; + + if (setup.planet!=kENone) + { + // coordinates of planet: topocentric, equatorial, J2000 + // One can use TT instead of TDB for all planets (except the moon?) + double ra, dec, diam; + palRdplan(tdb, setup.planet, elong, lat, &ra, &dec, &diam); + + // ---- apparent to mean ---- + palAmpqk(ra, dec, fAmprms, &out.source.ra, &out.source.dec); + } + + if (setup.wobble_offset<=0 || tpoint) + { + out.pointing.dec = out.source.dec; + out.pointing.ra = out.source.ra; + } + else + { + const double dphi = + setup.orbit_period==0 ? 0 : 2*M_PI*(_mjd-setup.start)/setup.orbit_period; + + const double phi = setup.wobble_angle + dphi; + + const double cosdir = cos(phi); + const double sindir = sin(phi); + const double cosoff = cos(setup.wobble_offset); + const double sinoff = sin(setup.wobble_offset); + const double cosdec = cos(out.source.dec); + const double sindec = sin(out.source.dec); + + const double sintheta = sindec*cosoff + cosdec*sinoff*cosdir; + + const double costheta = sintheta>1 ? 0 : sqrt(1 - sintheta*sintheta); + + const double cosdeltara = (cosoff - sindec*sintheta)/(cosdec*costheta); + const double sindeltara = sindir*sinoff/costheta; + + out.pointing.dec = asin(sintheta); + out.pointing.ra = atan2(sindeltara, cosdeltara) + out.source.ra; + } + + // ---- Mean to apparent ---- + double r=0, d=0; + palMapqkz(out.pointing.ra, out.pointing.dec, fAmprms, &r, &d); + + // + // Doesn't work - don't know why + // + // slaMapqk (radec.Ra(), radec.Dec(), rdpm.Ra(), rdpm.Dec(), + // 0, 0, (double*)fAmprms, &r, &d); + // + + // -- apparent to observed -- + palAopqk(r, d, fAoprms, + &out.sky.az, // observed azimuth (radians: N=0,E=90) [-pi, pi] + &out.sky.zd, // observed zenith distance (radians) [-pi/2, pi/2] + &out.apparent.ha, // observed hour angle (radians) + &out.apparent.dec, // observed declination (radians) + &out.apparent.ra); // observed right ascension (radians) + + // ----- fix ambiguity ----- + if (out.sky.zd<0) + { + out.sky.zd = -out.sky.zd; + out.sky.az += out.sky.az<0 ? M_PI : -M_PI; + } + + // Star culminating behind zenith and Az between ~90 and ~180deg + if (out.source.dec0) + out.sky.az -= 2*M_PI; + + out.mount = SkyToMount(out.sky); + + return out; + } +}; + +*/ +