00001
00002
00003
00004
00005
00006
00007 #ifndef HLRADII_H
00008 #define HLRADII_H
00009
00010 #include <float.h>
00011
00012 #include "Constants.h"
00013 #include "DataMethod.hxx"
00014 #include "FitsFile.h"
00015 #include "StringOps.h"
00016 #include "Logger.h"
00017 #include "DataOps.h"
00018 #include "Converter.h"
00019
00020
00025 template<typename DataType>
00026 class HlRadii : public DataMethod<DataType>
00027 {
00028 public:
00030 HlRadii(const FitsFile& aCatFile, const FitsFile& aFile)
00031 : mCatFile(aCatFile),mFile(aFile)
00032 {
00033 DataMethod<DataType>::methodName = "HlRadii";
00034 }
00035
00037 virtual ~HlRadii() { }
00038
00039
00041 void setSource(IntMap& tbls, IntMap& cols, StringMap units)
00042 {
00043 hlcColNo = cols[0];
00044 ellColNo = cols[1];
00045 extColNo = cols[2];
00046
00047
00048
00049
00050
00051
00052 }
00053
00055 void setTarget(IntMap& cols, StringMap units)
00056 {
00057 hlGeoColNo = cols[0];
00058 hlSMnColNo = cols[1];
00059 hlSMjColNo = cols[2];
00060 hlCorSMnColNo = cols[3];
00061 hlCorSMjColNo = cols[4];
00062 }
00063
00065 void doit(TableData<DataType>& data, int row1, int row2)
00066 {
00067
00068 Logger log;
00069
00070
00071
00072 double seeing[mCatFile.getNumHdus()];
00073 double crval1;
00074 double crval2;
00075 int npixx,npixy;
00076
00077 double centre[2];
00078 double sizeFact[2];
00079
00080
00081
00082 astBegin;
00083 AstFitsChan *fitschan = astFitsChan(NULL, NULL, "%.*s", 0, "NULL");
00084
00085
00086 Converter conv;
00087 sizeFact[0] = conv.getFactor("RADIANS", "Arcseconds");
00088 sizeFact[1] = conv.getFactor("RADIANS", "Arcseconds");
00089
00090
00091 for (int hduNo = 2; hduNo <= mCatFile.getNumHdus(); ++hduNo)
00092 {
00093 int extNo = hduNo - 2;
00094
00095 AstFrameSet *wcsinfo = 0;
00096
00097
00098
00099 mCatFile.movabsHdu(hduNo);
00100 mFile.movabsHdu(hduNo);
00101 mCatFile.readKey("SEEING", seeing[extNo]);
00102
00103
00104
00105 int ncards = mFile.getNumCards();
00106 for (int n = 0; n < ncards; ++n) {
00107 char card[80];
00108 mFile.readCard(n, card);
00109 astPutFits(fitschan, card, 0);
00110 }
00111
00112
00113 astClear( fitschan, "Card" );
00114
00115 mFile.readKey(FitsFile::CRVAL1, crval1);
00116 mFile.readKey(FitsFile::CRVAL2, crval2);
00117 if (crval1 == 0 && crval2 == 0) {
00118 log.addError("Could not read WCS info for " + mCatFile.getFileName());
00119 } else {
00120 wcsinfo = (AstFrameSet*)astRead(fitschan);
00121 }
00122 mCatFile.readKey("NXOUT", npixx);
00123 mCatFile.readKey("NYOUT", npixy);
00124 centre[0] = (double)npixx / 2;
00125 centre[1] = (double)npixy / 2;
00126 Numeric xPixSize,yPixSize;
00127
00128 DataOps::calcPixSize(wcsinfo,centre,sizeFact,xPixSize,yPixSize);
00129 seeing[extNo]*=sqrt(xPixSize*yPixSize);
00130
00131 }
00132
00133
00134 for (int rowNo = row1; rowNo < row2; ++rowNo)
00135 {
00136 double extNum;
00137 data.value(extColNo, rowNo, extNum);
00138 int extNo=(int)(extNum-2);
00139
00140 mCatFile.movabsHdu((int)extNum);
00141
00142 double hlCirc;
00143 double ell;
00144
00145 data.value(hlcColNo, rowNo, hlCirc);
00146
00147
00148
00149
00150
00151 data.value(ellColNo, rowNo, ell);
00152
00153 double hlSMj=DefaultValues::realdef;
00154 double hlSMn=DefaultValues::realdef;
00155 double hlGeo=DefaultValues::realdef;
00156 double hlCorSMj=DefaultValues::realdef;
00157 double hlCorSMn=DefaultValues::realdef;
00158 if (hlCirc>0. && ell>0.)
00159 {
00160 hlSMj=hlCirc*circToEll(ell);
00161 hlSMn=(1-ell)*hlSMj;
00162 hlGeo=sqrt(hlSMj*hlSMn);
00163
00164 double smjSq=hlSMj*hlSMj-seeFact*seeing[extNo]*seeing[extNo];
00165 if (smjSq>0) {
00166 hlCorSMj=sqrt(smjSq);
00167 } else {
00168 hlCorSMj=0.;
00169 }
00170 double smnSq=hlSMn*hlSMn-seeFact*seeing[extNo]*seeing[extNo];
00171 if (smnSq>0) {
00172 hlCorSMn=sqrt(smnSq);
00173 } else {
00174 hlCorSMn=0.;
00175 }
00176
00177 }
00178
00179
00180 data.assign(hlGeoColNo, rowNo, hlGeo);
00181 data.assign(hlSMnColNo, rowNo, hlSMn);
00182 data.assign(hlSMjColNo, rowNo, hlSMj);
00183 data.assign(hlCorSMnColNo, rowNo, hlCorSMn);
00184 data.assign(hlCorSMjColNo, rowNo, hlCorSMj);
00185
00186 }
00187 }
00188
00189
00190 void doit2(TableData<DataType>& rawData, TableData<DataType>& data, int row1, int row2) { }
00191
00192
00193
00194 double circToEll(double ell) {
00195
00196
00197 const double mofPars[3]={1.8243302,0.30914569,0.24304360};
00198 double axisRatio=1.-ell;
00199 axisRatio /= mofPars[1];
00200 return mofPars[0]/pow((1+axisRatio*axisRatio),mofPars[2]);
00201
00202 }
00203 private:
00204 const static double seeFact=0.50;
00205 FitsFile mCatFile;
00206 FitsFile mFile;
00207 int hlGeoColNo;
00208 int hlSMnColNo;
00209 int hlSMjColNo;
00210 int hlCorSMnColNo;
00211 int hlCorSMjColNo;
00212 int hlcColNo;
00213 int extColNo;
00214 int ellColNo;
00215
00216
00217
00218 };
00219
00220 #endif
00221
00222
00223
00224