00001
00002
00003
00004
00005
00006
00007 #ifndef QUANTITIES_HXX
00008 #define QUANTITIES_HXX
00009
00010
00011 extern "C"
00012 {
00013 #include "ast.h"
00014 }
00015
00016 #include <cmath>
00017 #include <string>
00018
00019 #include "Constants.h"
00020 #include "Converter.h"
00021 #include "FitsFile.h"
00022 #include "Logger.h"
00023 #include "DataMethod.hxx"
00024 #include "DataOps.h"
00025 #include "StringOps.h"
00026
00034 template<typename DataType>
00035 class QuantitiesByAst : public DataMethod<DataType>
00036 {
00037 public:
00039 QuantitiesByAst(const FitsFile& aFile)
00040 : mFile(aFile)
00041 { DataMethod<DataType>::methodName = "QuantitiesByAst"; }
00042
00044 virtual ~QuantitiesByAst() { }
00045
00048 void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00049 { extnumColno = cols[0]; }
00050
00052 void setTarget(IntMap& cols, StringMap units)
00053 {
00054 raColno = cols[0];
00055 decColno = cols[1];
00056 posangColno = cols[2];
00057 extentColno = cols[3];
00058 xsizeColno = cols[4];
00059 ysizeColno = cols[5];
00060 minRaColno = cols[6];
00061 minDecColno = cols[7];
00062 maxRaColno = cols[8];
00063 maxDecColno = cols[9];
00064
00065
00066
00067 Converter conv;
00068 raFact = conv.getFactor("RADIANS", units[0]);
00069 decFact = conv.getFactor("RADIANS", units[1]);
00070 posangFact = conv.getFactor("RADIANS", units[2]);
00071 extentFact = conv.getFactor("RADIANS", units[3]);
00072 xsizeFact = conv.getFactor("RADIANS", units[4]);
00073 ysizeFact = conv.getFactor("RADIANS", units[5]);
00074 minRaFact = conv.getFactor("RADIANS", units[6]);
00075 minDecFact = conv.getFactor("RADIANS", units[7]);
00076 maxRaFact = conv.getFactor("RADIANS", units[8]);
00077 maxDecFact = conv.getFactor("RADIANS", units[9]);
00078 }
00079
00081 void doit(TableData<DataType>& data, int row1, int row2)
00082 {
00083 Logger log;
00084
00085 std::string filename = mFile.getFileName();
00086 int numHDUs = mFile.getNumHdus();
00087 double crval1[numHDUs],crval2[numHDUs];
00088 int npixx[numHDUs], npixy[numHDUs];
00089
00090
00091 double xpix[5], ypix[5], xsky[5], ysky[5];
00092
00093
00094 astBegin;
00095
00096
00097 AstFitsChan *fitschan = astFitsChan(NULL, NULL, "%.*s", 0, "NULL");
00098 for (int rowNo = row1; rowNo < row2; ++rowNo)
00099 {
00100
00101
00102
00103
00104 Numeric value;
00105 data.value(extnumColno, rowNo, value);
00106 int hduNo = (int)value;
00107
00108 int extNum = hduNo - 2;
00109
00110
00111 AstFrameSet *wcsinfo = 0;
00112
00113
00114 mFile.movabsHdu(hduNo);
00115
00116
00117
00118 int ncards = mFile.getNumCards();
00119 for (int n = 0; n < ncards; ++n) {
00120 char card[80];
00121 mFile.readCard(n, card);
00122 astPutFits(fitschan, card, 0);
00123 }
00124
00125 astClear( fitschan, "Card" );
00126
00127 mFile.readKey(FitsFile::CRVAL1, crval1[extNum]);
00128 mFile.readKey(FitsFile::CRVAL2, crval2[extNum]);
00129
00130 if (crval1[extNum] == 0 && crval2[extNum] == 0) {
00131 if (filename.find("dark") == string::npos) {
00132 log.addError("Could not read WCS info in QBA for " + filename);
00133 }
00134 } else {
00135 wcsinfo = (AstFrameSet*)astRead(fitschan);
00136
00137
00138
00139 if (mFile.hasKey("ZIMAGE"))
00140 {
00141 mFile.readKey(FitsFile::ZNAXIS1, npixx[extNum]);
00142 mFile.readKey(FitsFile::ZNAXIS2, npixy[extNum]);
00143 }
00144 else
00145 {
00146 mFile.readKey(FitsFile::NAXIS1, npixx[extNum]);
00147 mFile.readKey(FitsFile::NAXIS2, npixy[extNum]);
00148 }
00149 }
00150
00151 if (crval1[extNum] != 0 || crval2[extNum] != 0)
00152 {
00153
00154
00155 xpix[0] = (double)npixx[extNum] / 2;
00156 ypix[0] = (double)npixy[extNum] / 2;
00157
00158
00159 xpix[1] = 0;
00160 ypix[1] = 0;
00161 xpix[2] = npixx[extNum];
00162 ypix[2] = 0;
00163 xpix[3] = 0;
00164 ypix[3] = npixy[extNum];
00165 xpix[4] = npixx[extNum];
00166 ypix[4] = npixy[extNum];
00167
00168
00169
00170 astTran2(wcsinfo, 5, xpix, ypix, 1, xsky, ysky);
00171
00172
00173 for (int n = 0; n < 5; ++n)
00174 {
00175 if (xsky[n] < 0)
00176 {
00177 xsky[n] += 2*M_PI;
00178 }
00179 }
00180
00181 Numeric raCen = raFact * xsky[0];
00182 Numeric decCen = decFact * ysky[0];
00183
00184
00185
00186 double xskyMin = 9999.999;
00187 double xskyMax = -9999.999;
00188 double yskyMin = 9999.999;
00189 double yskyMax = -9999.999;
00190 for (int n = 1; n < 5; ++n)
00191 {
00192 if (xsky[n] < xskyMin) xskyMin = xsky[n];
00193 if (xsky[n] > xskyMax) xskyMax = xsky[n];
00194 if (ysky[n] < yskyMin) yskyMin = ysky[n];
00195 if (ysky[n] > yskyMax) yskyMax = ysky[n];
00196 }
00197
00198 double h = 0.0;
00199 if (xskyMax - xskyMin > M_PI)
00200 {
00201 h = xskyMax;
00202 xskyMax = xskyMin + 2*M_PI;
00203 xskyMin = h;
00204 }
00205 Numeric raMin = minRaFact * xskyMin;
00206 Numeric decMin = minDecFact * yskyMin;
00207 Numeric raMax = maxRaFact * xskyMax;
00208 Numeric decMax = maxDecFact * yskyMax;
00209
00210
00211
00212 double maxsep = 0;
00213 for (int n = 0; n < 5; ++n)
00214 {
00215 double ang = DataOps::calcAngsep(xsky[0], ysky[0],
00216 xsky[n], ysky[n]);
00217 if (ang > maxsep) maxsep = ang;
00218 }
00219
00220 Numeric imageExtent = maxsep * extentFact;
00221
00222
00223
00224
00225 xsky[1] = xsky[0];
00226 ysky[1] = ysky[0] + 0.01;
00227
00228 astTran2(wcsinfo, 5, xsky, ysky, 0, xpix, ypix);
00229
00230
00231 Numeric posAng = posangFact * atan2(ypix[1] - ypix[0],
00232 xpix[1] - xpix[0]);
00233
00234
00235
00236 double centre[2]={xpix[0],ypix[0]};
00237 double convfact[2]={xsizeFact,ysizeFact};
00238
00239 Numeric xPixSize,yPixSize;
00240
00241 DataOps::calcPixSize(wcsinfo,centre,convfact,xPixSize,yPixSize);
00242
00243 data.assign(raColno, rowNo, raCen);
00244 data.assign(decColno, rowNo, decCen);
00245 data.assign(minRaColno, rowNo, raMin);
00246 data.assign(minDecColno, rowNo, decMin);
00247 data.assign(maxRaColno, rowNo, raMax);
00248 data.assign(maxDecColno, rowNo, decMax);
00249 data.assign(extentColno, rowNo, imageExtent);
00250 data.assign(posangColno, rowNo, posAng);
00251 data.assign(xsizeColno, rowNo, xPixSize);
00252 data.assign(ysizeColno, rowNo, yPixSize);
00253 }
00254 }
00255
00256 astEnd;
00257 }
00258
00259 void doit2(TableData<DataType>& rawData, TableData<DataType>& data, int row1, int row2) { }
00260
00261 private:
00262 FitsFile mFile;
00263 Logger log;
00264
00265 int extnumColno, raColno, decColno, posangColno;
00266 int extentColno, xsizeColno, ysizeColno;
00267 int minRaColno,maxRaColno,minDecColno,maxDecColno;
00268 double raFact, decFact, posangFact, extentFact, xsizeFact, ysizeFact;
00269 double minRaFact, maxRaFact, minDecFact, maxDecFact;
00270 };
00271
00272 #endif
00273
00274
00275
00276
00277
00278
00279