source: tags/release-1.3.2/src/Outputs/columns.cc @ 1327

Last change on this file since 1327 was 1192, checked in by MatthewWhiting, 11 years ago

Ticket #187 - Fixing the problem with Z_PEAK - now have better way of dealing with the default columns.

File size: 35.5 KB
Line 
1// -----------------------------------------------------------------------
2// columns.cc: Define a set of columns for Duchamp output.
3// -----------------------------------------------------------------------
4// Copyright (C) 2006, Matthew Whiting, ATNF
5//
6// This program is free software; you can redistribute it and/or modify it
7// under the terms of the GNU General Public License as published by the
8// Free Software Foundation; either version 2 of the License, or (at your
9// option) any later version.
10//
11// Duchamp is distributed in the hope that it will be useful, but WITHOUT
12// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13// FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14// for more details.
15//
16// You should have received a copy of the GNU General Public License
17// along with Duchamp; if not, write to the Free Software Foundation,
18// Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
19//
20// Correspondence concerning Duchamp may be directed to:
21//    Internet email: Matthew.Whiting [at] atnf.csiro.au
22//    Postal address: Dr. Matthew Whiting
23//                    Australia Telescope National Facility, CSIRO
24//                    PO Box 76
25//                    Epping NSW 1710
26//                    AUSTRALIA
27// -----------------------------------------------------------------------
28#include <iostream>
29#include <sstream>
30
31#include <math.h>
32#include <duchamp/duchamp.hh>
33#include <duchamp/fitsHeader.hh>
34#include <vector>
35#include <string>
36#include <duchamp/Detection/detection.hh>
37#include <duchamp/Outputs/columns.hh>
38
39namespace duchamp
40{
41
42  namespace Catalogues
43  {
44
45
46    Column::Column()
47    {
48      /// @brief Set the default values for the parameters.
49      this->itsWidth=1;
50      this->itsPrecision=0;
51      this->itsName=" ";
52      this->itsUnits=" ";
53      this->itsType="";
54      this->itsUCD = "";
55      this->itsDatatype = "";
56      this->itsExtraInfo="";
57      this->itsColID="";
58    }
59
60    Column::Column(std::string type)
61    {
62      operator=(defaultColumn(type));
63    }
64
65    Column::Column(const Column& c)
66    {
67      operator=(c);
68    }
69
70    Column& Column::operator=(const Column& c)
71    {
72      if(this==&c) return *this;
73      this->itsWidth = c.itsWidth;
74      this->itsPrecision = c.itsPrecision;
75      this->itsName = c.itsName;
76      this->itsUnits = c.itsUnits;
77      this->itsUCD = c.itsUCD;
78      this->itsDatatype = c.itsDatatype;
79      this->itsExtraInfo = c.itsExtraInfo;
80      this->itsColID = c.itsColID;
81      this->itsType = c.itsType;
82      return *this;
83    }
84
85    Column::Column(std::string type, std::string name, std::string units, int width, int prec, std::string ucd, std::string datatype, std::string colID, std::string extraInfo):
86      itsType(type), itsName(name), itsUnits(units), itsWidth(width), itsPrecision(prec), itsUCD(ucd), itsDatatype(datatype), itsColID(colID), itsExtraInfo(extraInfo)
87    {
88      /// A generic constructor designed to make a Column other than
89      /// the predefined ones listed in the Column namespace
90      /// \param name The title of the column
91      /// \param units What units the column takes
92      /// \param width The starting width
93      /// \param prec The starting precision
94      /// \param ucd The UCD relating to the quantity (VOTable use only)
95      /// \param datatype The datatype (float/char/int/etc) of the quantity (VOTable)
96      /// \param colID The column ID string used in the FIELD entry of a VOTable
97      /// \param extraInfo Other information used in the VOTable (eg. a reference to a COOSYS)
98   }
99
100    //------------------------------------------------------------
101
102    void Column::changePrec(int p)
103    {
104      /// A direct way to alter the precision without having to use
105      /// Column::upPrec(). If the precision increases, the width will
106      /// increase by the same amount. If it decreases, no change is
107      /// made to the width.
108      /// \param p The new precision value.
109
110      if(p > this->itsPrecision) this->itsWidth += (p-this->itsPrecision);
111      this->itsPrecision = p;
112
113    }
114   
115    //------------------------------------------------------------
116    void Column::printTitle(std::ostream &stream)
117    {
118      stream << std::setw(this->itsWidth)
119             << std::setfill(' ')
120             << this->itsName;
121    }
122
123    void Column::printUnits(std::ostream &stream)
124    {
125      stream << std::setw(this->itsWidth)
126             << std::setfill(' ')
127             << this->itsUnits;
128    }
129 
130    void Column::printDash (std::ostream &stream)
131    {
132      stream << std::setw(this->itsWidth)
133             << std::setfill('-')
134             << ""
135             << std::setfill(' ');
136    }
137
138    void Column::printBlank(std::ostream &stream)
139    {
140      stream << std::setw(this->itsWidth)
141             << std::setfill(' ')
142             << "";
143    }
144 
145    template <class T> void Column::printEntry(std::ostream &stream, T value)
146    {
147      ///  \param stream Where the printing is done.
148      ///  \param value  The value to be printed.
149
150      stream << std::setprecision(this->itsPrecision)
151             << std::setw(this->itsWidth)
152             << std::setfill(' ')
153             << value;
154    }
155    template void Column::printEntry<int>(std::ostream &stream, int value);
156    template void Column::printEntry<long>(std::ostream &stream, long value);
157    template void Column::printEntry<unsigned int>(std::ostream &stream, unsigned int value);
158    template void Column::printEntry<unsigned long>(std::ostream &stream, unsigned long value);
159    template void Column::printEntry<float>(std::ostream &stream, float value);
160    template void Column::printEntry<double>(std::ostream &stream, double value);
161    template void Column::printEntry<std::string>(std::ostream &stream, std::string value);
162
163 
164    bool Column::doCol(DESTINATION tableType, bool flagFint)
165    {
166      ///  @details Determines whether a given column is used for a
167      ///  given table type, based on the Column::type string.
168      /// \param tableType The type of table: one of FILE, SCREEN, LOG, VOTABLE.
169      /// \param flagFint Whether to use FINT (true) or FTOT
170      /// (false). This defaults to true, so need not be given. It only
171      /// applies to the screen and votable cases -- both are written for the
172      /// results file case.
173      /// \return True if column is used for given table type. False
174      /// otherwise. False if tableType not one of four listed.
175     
176 
177      const size_t sizeFileList=39;
178      std::string FileList[sizeFileList]={"NUM","NAME","X","Y","Z","RA","DEC","VEL","MAJ","MIN","PA","WRA","WDEC",
179                                          "W50","W20","WVEL","FINT", "FINTERR", "FTOT", "FTOTERR", "FPEAK","SNRPEAK",
180                                          "X1","X2","Y1","Y2","Z1","Z2","NPIX","FLAG","XAV","YAV",
181                                          "ZAV","XCENT","YCENT","ZCENT","XPEAK","YPEAK","ZPEAK"};
182      const size_t sizeScreenList=22;
183      std::string ScreenList[sizeScreenList]={"NUM","NAME","X","Y","Z","RA","DEC","VEL","MAJ","MIN","PA",
184                                              "W50","FPEAK","SNRPEAK","X1","X2","Y1",
185                                              "Y2","Z1","Z2","NPIX","FLAG"};
186      const size_t sizeLogList=16;
187      std::string LogList[sizeLogList]={"NUM","X","Y","Z","FTOT","FPEAK","SNRPEAK",
188                                        "X1","X2","Y1","Y2","Z1","Z2","NPIX","NUMCH","SPATSIZE"};
189      const size_t sizeVOList=23;
190      std::string VOList[sizeVOList]={"NUM","NAME","RAJD","DECJD","VEL","MAJ","MIN","PA",
191                                      "W50","W20","WVEL","FPEAK","SNRPEAK","FLAG","XAV","YAV",
192                                      "ZAV","XCENT","YCENT","ZCENT","XPEAK","YPEAK","ZPEAK"};
193
194      bool doIt=false;
195      if(tableType == FILE){
196        for(size_t i=0;i<sizeFileList && !doIt;i++) doIt = doIt || (this->itsType==FileList[i]);
197      }
198      else if(tableType == SCREEN){
199        if(this->itsType=="FTOT") doIt = !flagFint;
200        else if(this->itsType == "FINT") doIt = flagFint;
201        else for(size_t i=0;i<sizeScreenList && !doIt;i++) doIt = doIt || (this->itsType==ScreenList[i]);
202      }
203      else if(tableType == LOG){
204        for(size_t i=0;i<sizeLogList && !doIt;i++) doIt = doIt || (this->itsType==LogList[i]);
205      }
206      else if(tableType == VOTABLE){
207        if(this->itsType=="FTOT" || this->itsType=="FTOTERR") doIt = !flagFint;
208        else if(this->itsType == "FINT" || this->itsType=="FINTERR") doIt = flagFint;
209        else for(size_t i=0;i<sizeVOList && !doIt;i++) doIt = doIt || (this->itsType==VOList[i]);
210      }
211      return doIt;
212       
213    }
214
215    Column defaultColumn(std::string type)
216    {
217      // type|name|units|width|prec|ucd|datatype|extraInfo
218      if (type == "NUM") return Column(type,"ObjID","",6,0,"meta.id","int","col_objnum","");
219      else if(type=="NAME") return Column(type,"Name","",8,0,"meta.id;meta.main","char","col_name","");
220      else if(type=="X") return Column(type,"X","",6,prXYZ,"pos.cartesian.x","float","col_x","");
221      else if(type=="Y") return Column(type,"Y","",6,prXYZ,"pos.cartesian.y","float","col_y","");
222      else if(type=="Z") return Column(type,"Z","",6,prXYZ,"pos.cartesian.z","float","col_z","");
223      else if(type=="RA") return Column(type,"RA","",11,0,"","char","col_ra",coordRef);
224      else if(type=="DEC") return Column(type,"DEC","",10,0,"","char","col_dec",coordRef);
225      else if(type=="RAJD") return Column(type,"RAJD","[deg]",11,prPOS,"","float","col_rajd",coordRef);
226      else if(type=="DECJD") return Column(type,"DECJD","[deg]",11,prPOS,"","float","col_decjd",coordRef);
227      else if(type=="VEL") return Column(type,"VEL","",7,prVEL,"","float","col_vel","");
228      else if(type=="MAJ") return Column(type,"MAJ","[deg]",6,prWPOS,"phys.angSize.smajAxis","float","col_maj",coordRef);
229      else if(type=="MIN") return Column(type,"MIN","[deg]",6,prWPOS,"phys.angSize.sminAxis","float","col_min",coordRef);
230      else if(type=="PA") return Column(type,"PA","[deg]",7,prWPOS,"pos.posAng;phys.angSize","float","col_pa",coordRef);
231      else if(type=="WRA") return Column(type,"w_RA","[arcmin]",9,prWPOS,"","float","col_wra",coordRef);
232      else if(type=="WDEC") return Column(type,"w_DEC","[arcmin]",9,prWPOS,"","float","col_wdec",coordRef);
233      else if(type=="W50") return Column(type,"w_50","",7,prVEL,"","float","col_w50","");
234      else if(type=="W20") return Column(type,"w_20","",7,prVEL,"","float","col_w20","");
235      else if(type=="WVEL") return Column(type,"w_VEL","",7,prVEL,"","float","col_wvel","");
236      else if(type=="FINT") return Column(type,"F_int","",10,prFLUX,"phot.flux.density.integrated","float","col_fint","");
237      else if(type=="FINTERR") return Column(type,"eF_int","",10,prFLUX,"phot.flux.density.integrated;stat.err","float","col_efint","");
238      else if(type=="FTOT") return Column(type,"F_tot","",10,prFLUX,"phot.flux.density","float","col_ftot","");
239      else if(type=="FTOTERR") return Column(type,"eF_tot","",10,prFLUX,"phot.flux.density;stat.err","float","col_eftot","");
240      else if(type=="FPEAK") return Column(type,"F_peak","",9,prFLUX,"phot.flux;stat.max","float","col_fpeak","");
241      else if(type=="SNRPEAK") return Column(type,"S/Nmax","",7,prSNR,"stat.snr;phot.flux","float","col_snrmax","");
242      else if(type=="X1") return Column(type,"X1","",4,0,"pos.cartesian.x;stat.min","int","col_x1","");
243      else if(type=="Y1") return Column(type,"Y1","",4,0,"pos.cartesian.y;stat.min","int","col_y1","");
244      else if(type=="Z1") return Column(type,"Z1","",4,0,"pos.cartesian.z;stat.min","int","col_z1","");
245      else if(type=="X2") return Column(type,"X2","",4,0,"pos.cartesian.x;stat.max","int","col_x2","");
246      else if(type=="Y2") return Column(type,"Y2","",4,0,"pos.cartesian.y;stat.max","int","col_y2","");
247      else if(type=="Z2") return Column(type,"Z2","",4,0,"pos.cartesian.z;stat.max","int","col_z2","");
248      else if(type=="NPIX") return Column(type,"Npix","[pix]",6,0,"","int","col_npix","");
249      else if(type=="FLAG") return Column(type,"Flag","",5,0,"meta.code.qual","char","col_flag","");
250      else if(type=="XAV") return Column(type,"X_av","",6,prXYZ,"pos.cartesian.x;stat.mean","float","col_xav","");
251      else if(type=="YAV") return Column(type,"Y_av","",6,prXYZ,"pos.cartesian.y;stat.mean","float","col_yav","");
252      else if(type=="ZAV") return Column(type,"Z_av","",6,prXYZ,"pos.cartesian.z;stat.mean","float","col_zav","");
253      else if(type=="XCENT") return Column(type,"X_cent","",7,prXYZ,"pos.cartesian.x;stat.centroid","float","col_xcent","");
254      else if(type=="YCENT") return Column(type,"Y_cent","",7,prXYZ,"pos.cartesian.y;stat.centroid","float","col_ycent","");
255      else if(type=="ZCENT") return Column(type,"Z_cent","",7,prXYZ,"pos.cartesian.z;stat.centroid","float","col_zcent","");
256      else if(type=="XPEAK") return Column(type,"X_peak","",7,prXYZ,"pos.cartesian.x;phot.flux;stat.max","int","col_xpeak","");
257      else if(type=="YPEAK") return Column(type,"Y_peak","",7,prXYZ,"pos.cartesian.y;phot.flux;stat.max","int","col_ypeak","");
258      else if(type=="ZPEAK") return Column(type,"Z_peak","",7,prXYZ,"pos.cartesian.z;phot.flux;stat.max","int","col_zpeak","");
259      else if(type=="NUMCH") return Column(type,"Nch","",6,0,"em.bin;stat.sum","int","col_nch","");
260      else if(type=="SPATSIZE") return Column(type,"Sp_size","",8,0,"instr.pixel;stat.sum","int","col_spsize","");
261      else {
262        Column col;
263        col.setType(type);
264        return col;
265      }
266
267    }
268
269    CatalogueSpecification getFullColSet(std::vector<Detection> &objectList, FitsHeader &head)
270    {
271      ///  @details A function that returns a catalogue specification
272      ///  containing information on the columns necessary for output
273      ///  to the results file:
274      ///  Obj#,NAME,X,Y,Z,RA,DEC,VEL,w_RA,w_DEC,w_VEL,F_tot,F_int,F_peak,
275      ///  X1,X2,Y1,Y2,Z1,Z2,Npix,Flag,
276      ///  XAV,YAV,ZAV,XCENT,YCENT,ZCENT,XPEAK,YPEAK,ZPEAK
277      ///
278      ///   Each object in the provided objectList is checked to see if it
279      ///   requires any column to be widened, or for that column to have
280      ///   its precision increased.
281      ///
282      ///   Both Ftot and Fint are provided -- it is up to the calling
283      ///   function to determine which to use.
284      ///
285      /// \param objectList A std::vector list of Detection objects that the
286      /// columns need to fit.
287      /// \param head The FitsHeader object defining the World Coordinate
288      /// System.
289      /// \return A CatalogueSpecification
290
291      CatalogueSpecification newset;
292
293      // desired precisions for fluxes, velocities and SNR value
294      int precVel,precFlux,precSNR;
295      if(objectList.size()>0){
296        precVel=objectList[0].getVelPrec();
297        precFlux=objectList[0].getFpeakPrec(); // same as FintPrec at this point.
298        precSNR=objectList[0].getSNRPrec();
299      }
300      else {
301        precVel = prVEL;
302        precFlux = prFLUX;
303        precSNR = prSNR;
304      }
305
306      // set up the default columns
307      newset.addColumn( Column("NUM") );
308      newset.addColumn( Column("NAME") );
309      newset.addColumn( Column("X") );
310      newset.addColumn( Column("Y") );
311      newset.addColumn( Column("Z") );
312      newset.addColumn( Column("RA") );
313      newset.addColumn( Column("DEC") );
314      newset.addColumn( Column("RAJD") );
315      newset.addColumn( Column("DECJD") );
316      newset.addColumn( Column("VEL") );
317      newset.addColumn( Column("MAJ") );
318      newset.addColumn( Column("MIN") );
319      newset.addColumn( Column("PA") );
320      newset.addColumn( Column("WRA") );
321      newset.addColumn( Column("WDEC") );
322      newset.addColumn( Column("W50") );
323      newset.addColumn( Column("W20") );
324      newset.addColumn( Column("WVEL") );
325      newset.addColumn( Column("FINT") );
326      newset.addColumn( Column("FINTERR") );
327      newset.addColumn( Column("FTOT") );
328      newset.addColumn( Column("FTOTERR") );
329      newset.addColumn( Column("FPEAK") );
330      newset.addColumn( Column("SNRPEAK") );
331      newset.addColumn( Column("X1") );
332      newset.addColumn( Column("X2") );
333      newset.addColumn( Column("Y1") );
334      newset.addColumn( Column("Y2") );
335      newset.addColumn( Column("Z1") );
336      newset.addColumn( Column("Z2") );
337      newset.addColumn( Column("NPIX") );
338      newset.addColumn( Column("FLAG") );
339      newset.addColumn( Column("XAV") );
340      newset.addColumn( Column("YAV") );
341      newset.addColumn( Column("ZAV") );
342      newset.addColumn( Column("XCENT") );
343      newset.addColumn( Column("YCENT") );
344      newset.addColumn( Column("ZCENT") );
345      newset.addColumn( Column("XPEAK") );
346      newset.addColumn( Column("YPEAK") );
347      newset.addColumn( Column("ZPEAK") );
348      newset.addColumn( Column("NUMCH") );
349      newset.addColumn( Column("SPATSIZE") );
350
351
352      // Now test each object against each new column
353      std::vector<Detection>::iterator obj;
354      for(obj = objectList.begin(); obj < objectList.end(); obj++){
355        std::string tempstr;
356        int tempwidth;
357        float val,minval;
358        double valD,minvalD;
359
360        // Obj#
361        tempwidth = int( log10(obj->getID()) + 1) + 1;
362        for(int i=newset.column("NUM").getWidth();i<tempwidth;i++) newset.column("NUM").widen();
363
364        // Name
365        tempwidth = obj->getName().size() + 1;
366        for(int i=newset.column("NAME").getWidth();i<tempwidth;i++) newset.column("NAME").widen();
367
368        // X position
369        val = obj->getXcentre() + obj->getXOffset();
370        if((val<1.)&&(val>0.)){
371          minval = pow(10, -1. * (newset.column("X").getPrecision()+1));
372          if(val < minval) newset.column("X").upPrec();
373        }
374        tempwidth = int( log10(val) + 1) + newset.column("X").getPrecision() + 2;
375        for(int i=newset.column("X").getWidth();i<tempwidth;i++) newset.column("X").widen();
376        // Y position
377        val = obj->getYcentre() + obj->getYOffset();
378        if((val<1.)&&(val>0.)){
379          minval = pow(10, -1. * (newset.column("Y").getPrecision()+1));
380          if(val < minval) newset.column("Y").upPrec();
381        }
382        tempwidth = int( log10(val) + 1) + newset.column("Y").getPrecision() + 2;
383        for(int i=newset.column("Y").getWidth();i<tempwidth;i++) newset.column("Y").widen();
384        // Z position
385        val = obj->getZcentre() + obj->getZOffset();
386        if((val<1.)&&(val>0.)){
387          minval = pow(10, -1. * (newset.column("Z").getPrecision()+1));
388          if((val>0.)&&(val < minval)) newset.column("Z").upPrec();
389        }
390        tempwidth = int( log10(val) + 1) + newset.column("Z").getPrecision() + 2;
391        for(int i=newset.column("Z").getWidth();i<tempwidth;i++) newset.column("Z").widen();
392
393        if(head.isWCS()){ 
394          // RA -- assign correct title. Check width but should be ok by definition
395          tempstr = head.lngtype();
396          newset.column("RA").setName(tempstr);
397          tempwidth = obj->getRAs().size() + 1;
398          for(int i=newset.column("RA").getWidth();i<tempwidth;i++) newset.column("RA").widen();
399     
400          // Dec -- assign correct title. Check width, should be ok by definition
401          tempstr = head.lattype();
402          newset.column("DEC").setName(tempstr);
403          tempwidth = obj->getDecs().size() + 1;
404          for(int i=newset.column("DEC").getWidth();i<tempwidth;i++) newset.column("DEC").widen();
405
406          // RA decimal degrees -- assign correct title. Check width but should be OK
407          tempstr = head.lngtype();
408          newset.column("RAJD").setName(tempstr);
409          valD = obj->getRA();
410          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("RAJD").getPrecision() + 2;
411          for(int i=newset.column("RAJD").getWidth();i<tempwidth;i++) newset.column("RAJD").widen();
412     
413          // Dec decimal degrees -- assign correct title. Check width but should be OK
414          tempstr = head.lattype();
415          newset.column("DECJD").setName(tempstr);
416          valD = obj->getDec();
417          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("DECJD").getPrecision() + 2;
418          for(int i=newset.column("DECJD").getWidth();i<tempwidth;i++) newset.column("DECJD").widen();
419
420          // Vel -- check width, title and units.
421          if(head.canUseThirdAxis()){
422            if(head.WCS().spec < 0)  // if it's not a spectral axis
423              newset.column("VEL").setName( head.WCS().ctype[2] );
424            else
425              newset.column("VEL").setName(head.getSpectralType());
426            tempwidth = newset.column("VEL").getName().size() + 1;
427            for(int i=newset.column("VEL").getWidth();i<tempwidth;i++) newset.column("VEL").widen();
428            if(head.getSpectralUnits().size()>0)
429              newset.column("VEL").setUnits("[" + head.getSpectralUnits() + "]");
430            tempwidth = newset.column("VEL").getUnits().size() + 1;
431            for(int i=newset.column("VEL").getWidth();i<tempwidth;i++) newset.column("VEL").widen();
432       
433            valD = obj->getVel();
434            if((fabs(valD) < 1.)&&(valD>0.)){
435              minvalD = pow(10, -1. * (newset.column("VEL").getPrecision()+1));
436              if(valD < minvalD) newset.column("VEL").upPrec();
437            }
438            tempwidth = int(log10(fabs(valD)) + 1) + newset.column("VEL").getPrecision() + 2;
439            if(valD<0) tempwidth++;
440            for(int i=newset.column("VEL").getWidth();i<tempwidth;i++) newset.column("VEL").widen();
441          }
442
443          // MAJ -- check width & title. leave units for the moment.
444          newset.column("MAJ").setUnits("["+head.getShapeUnits()+"]");
445          tempwidth = newset.column("MAJ").getUnits().size() + 1;
446          for(int i=newset.column("MAJ").getWidth();i<tempwidth;i++) newset.column("MAJ").widen();
447          valD = obj->getMajorAxis();
448          if((fabs(valD) < 1.)&&(valD>0.)){
449            minvalD = pow(10, -1. * (newset.column("MAJ").getPrecision()));
450            if(valD < minvalD) newset.column("MAJ").upPrec();
451          }
452          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("MAJ").getPrecision() + 2;
453          if(valD<0) tempwidth++;
454          for(int i=newset.column("MAJ").getWidth();i<tempwidth;i++) newset.column("MAJ").widen();
455
456          // MIN -- check width & title. leave units for the moment.
457          newset.column("MIN").setUnits("["+head.getShapeUnits()+"]");
458          tempwidth = newset.column("MIN").getUnits().size() + 1;
459          for(int i=newset.column("MIN").getWidth();i<tempwidth;i++) newset.column("MIN").widen();
460          valD = obj->getMinorAxis();
461          if((fabs(valD) < 1.)&&(valD>0.)){
462            minvalD = pow(10, -1. * (newset.column("MIN").getPrecision()));
463            if(valD < minvalD) newset.column("MIN").upPrec();
464          }
465          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("MIN").getPrecision() + 2;
466          if(valD<0) tempwidth++;
467          for(int i=newset.column("MIN").getWidth();i<tempwidth;i++) newset.column("MIN").widen();
468
469          // PA -- check width & title. leave units for the moment.
470          tempwidth = newset.column("PA").getUnits().size() + 1;
471          for(int i=newset.column("PA").getWidth();i<tempwidth;i++) newset.column("PA").widen();
472          valD = obj->getPositionAngle();
473          if((fabs(valD) < 1.)&&(valD>0.)){
474            minvalD = pow(10, -1. * (newset.column("PA").getPrecision()+1));
475            if(valD < minvalD) newset.column("PA").upPrec();
476          }
477          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("PA").getPrecision() + 2;
478          if(valD<0) tempwidth++;
479          for(int i=newset.column("PA").getWidth();i<tempwidth;i++) newset.column("PA").widen();
480
481          // w_RA -- check width & title. leave units for the moment.
482          tempstr = "w_" + head.lngtype();
483          newset.column("WRA").setName(tempstr);
484          tempwidth = newset.column("RA").getUnits().size() + 1;
485          for(int i=newset.column("WRA").getWidth();i<tempwidth;i++) newset.column("WRA").widen();
486          valD = obj->getRAWidth();
487          if((fabs(valD) < 1.)&&(valD>0.)){
488            minvalD = pow(10, -1. * (newset.column("WRA").getPrecision()+1));
489            if(valD < minvalD) newset.column("WRA").upPrec();
490          }
491          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("WRA").getPrecision() + 2;
492          if(valD<0) tempwidth++;
493          for(int i=newset.column("WRA").getWidth();i<tempwidth;i++) newset.column("WRA").widen();
494
495          // w_DEC -- check width & title. leave units for the moment.
496          tempstr = "w_" + head.lattype();
497          newset.column("WDEC").setName(tempstr);
498          tempwidth = newset.column("DEC").getUnits().size() + 1;
499          for(int i=newset.column("WDEC").getWidth();i<tempwidth;i++) newset.column("WDEC").widen();
500          valD = obj->getDecWidth();
501          if((fabs(valD) < 1.)&&(valD>0.)){
502            minvalD = pow(10, -1. * (newset.column("WDEC").getPrecision()+1));
503            if(valD < minvalD) newset.column("WDEC").upPrec();
504          }
505          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("WDEC").getPrecision() + 2;
506          if(valD<0) tempwidth++;
507          for(int i=newset.column("WDEC").getWidth();i<tempwidth;i++) newset.column("WDEC").widen();
508
509          // w_50 -- check width, title and units.
510          if(head.canUseThirdAxis()){
511            if(head.getSpectralUnits().size()>0)
512              newset.column("W50").setUnits("[" + head.getSpectralUnits() + "]");
513            tempwidth = newset.column("W50").getUnits().size() + 1;
514            for(int i=newset.column("W50").getWidth();i<tempwidth;i++) newset.column("W50").widen();
515            valD = obj->getW50();
516            if((fabs(valD) < 1.)&&(valD>0.)){
517              minvalD = pow(10, -1. * (newset.column("W50").getPrecision()+1));
518              if(valD < minvalD) newset.column("W50").upPrec();
519            }
520            tempwidth = int( log10(fabs(valD)) + 1) + newset.column("W50").getPrecision() + 2;
521            if(valD<0) tempwidth++;
522            for(int i=newset.column("W50").getWidth();i<tempwidth;i++) newset.column("W50").widen();
523          }
524
525          // w_20 -- check width, title and units.
526          if(head.canUseThirdAxis()){
527            if(head.getSpectralUnits().size()>0)
528              newset.column("W20").setUnits("[" + head.getSpectralUnits() + "]");
529            tempwidth = newset.column("W20").getUnits().size() + 1;
530            for(int i=newset.column("W20").getWidth();i<tempwidth;i++)newset.column("W20").widen();
531            valD = obj->getW20();
532            if((fabs(valD) < 1.)&&(valD>0.)){
533              minvalD = pow(10, -1. * (newset.column("W20").getPrecision()+1));
534              if(valD < minvalD) newset.column("W20").upPrec();
535            }
536            tempwidth = int( log10(fabs(valD)) + 1) + newset.column("W20").getPrecision() + 2;
537            if(valD<0) tempwidth++;
538            for(int i=newset.column("W20").getWidth();i<tempwidth;i++) newset.column("W20").widen();
539          }
540
541          // w_Vel -- check width, title and units.
542          if(head.canUseThirdAxis()){
543            if(head.WCS().spec < 0) // if it's not a spectral axis
544              newset.column("WVEL").setName( std::string("w_") + head.WCS().ctype[2] );
545            else
546              newset.column("WVEL").setName(std::string("w_")+head.getSpectralType());
547            if(head.getSpectralUnits().size()>0)
548              newset.column("WVEL").setUnits("[" + head.getSpectralUnits() + "]");
549            tempwidth = newset.column("WVEL").getUnits().size() + 1;
550            for(int i=newset.column("WVEL").getWidth();i<tempwidth;i++)newset.column("WVEL").widen();
551            tempwidth = newset.column("WVEL").getName().size() + 1;
552            for(int i=newset.column("WVEL").getWidth();i<tempwidth;i++) newset.column("WVEL").widen();
553            valD = obj->getVelWidth();
554            if((fabs(valD) < 1.)&&(valD>0.)){
555              minvalD = pow(10, -1. * (newset.column("WVEL").getPrecision()+1));
556              if(valD < minvalD) newset.column("WVEL").upPrec();
557            }
558            tempwidth = int( log10(fabs(valD)) + 1) + newset.column("WVEL").getPrecision() + 2;
559            if(valD<0) tempwidth++;
560            for(int i=newset.column("WVEL").getWidth();i<tempwidth;i++) newset.column("WVEL").widen();
561          }
562
563          // F_int -- check width & units
564          if(head.getIntFluxUnits().size()>0)
565            newset.column("FINT").setUnits("[" + head.getIntFluxUnits() + "]");
566          tempwidth = newset.column("FINT").getUnits().size() + 1;
567          for(int i=newset.column("FINT").getWidth();i<tempwidth;i++) newset.column("FINT").widen();
568          valD = obj->getIntegFlux();
569          if((fabs(valD) < 1.)// &&(valD>0.)
570             ){
571            int minprec = int(fabs(log10(fabs(valD))))+2;
572            for(int i=newset.column("FINT").getPrecision();i<minprec;i++) newset.column("FINT").upPrec();
573          }
574          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("FINT").getPrecision() + 2;
575          if(valD<0) tempwidth++;
576          for(int i=newset.column("FINT").getWidth();i<tempwidth;i++) newset.column("FINT").widen();
577     
578          // eF_int -- check width & units
579          if(obj->getIntegFluxError()>0.){
580              if(head.getIntFluxUnits().size()>0)
581                  newset.column("FINTERR").setUnits("[" + head.getIntFluxUnits() + "]");
582              tempwidth = newset.column("FINTERR").getUnits().size() + 1;
583              for(int i=newset.column("FINTERR").getWidth();i<tempwidth;i++) newset.column("FINTERR").widen();
584              valD = obj->getIntegFluxError();
585              if((fabs(valD) < 1.)// &&(valD>0.)
586                  ){
587                  int minprec = int(fabs(log10(fabs(valD))))+2;
588                  for(int i=newset.column("FINTERR").getPrecision();i<minprec;i++) newset.column("FINTERR").upPrec();
589              }
590              tempwidth = int( log10(fabs(valD)) + 1) + newset.column("FINTERR").getPrecision() + 2;
591              if(valD<0) tempwidth++;
592              for(int i=newset.column("FINTERR").getWidth();i<tempwidth;i++) newset.column("FINTERR").widen();
593          }
594
595        }
596
597        // F_tot
598        newset.column("FTOT").setUnits("[" + head.getFluxUnits() + "]");
599        tempwidth = newset.column("FTOT").getUnits().size() + 1;
600        for(int i=newset.column("FTOT").getWidth();i<tempwidth;i++) newset.column("FTOT").widen();
601        val = obj->getTotalFlux();
602        //     std::cerr << val << "\n";
603        if((fabs(val) < 1.)// &&(val>0.)
604           ){
605          int minprec = int(fabs(log10(fabs(val))))+2;
606          for(int i=newset.column("FTOT").getPrecision();i<minprec;i++) newset.column("FTOT").upPrec();
607        }
608        tempwidth = int( log10(fabs(val)) + 1) + newset.column("FTOT").getPrecision() + 2;
609        if(val<0) tempwidth++;
610        for(int i=newset.column("FTOT").getWidth();i<tempwidth;i++) newset.column("FTOT").widen();
611
612        // eF_tot -- check width & units
613        if(obj->getIntegFluxError()>0.){
614            newset.column("FTOTERR").setUnits("[" + head.getFluxUnits() + "]");
615            tempwidth = newset.column("FTOTERR").getUnits().size() + 1;
616            for(int i=newset.column("FTOTERR").getWidth();i<tempwidth;i++) newset.column("FTOTERR").widen();
617            valD = obj->getTotalFluxError();
618            if((fabs(valD) < 1.)// &&(valD>0.)
619                ){
620                int minprec = int(fabs(log10(fabs(valD))))+2;
621                for(int i=newset.column("FTOTERR").getPrecision();i<minprec;i++) newset.column("FTOTERR").upPrec();
622              }
623            tempwidth = int( log10(fabs(valD)) + 1) + newset.column("FTOTERR").getPrecision() + 2;
624            if(valD<0) tempwidth++;
625            for(int i=newset.column("FTOTERR").getWidth();i<tempwidth;i++) newset.column("FTOTERR").widen();
626        }
627
628        // F_peak
629        newset.column("FPEAK").setUnits("[" + head.getFluxUnits() + "]");
630        tempwidth = newset.column("FPEAK").getUnits().size() + 1;
631        for(int i=newset.column("FPEAK").getWidth();i<tempwidth;i++) newset.column("FPEAK").widen();
632        val = obj->getPeakFlux();
633        if((fabs(val) < 1.)// &&(val>0.)
634           ){
635          int minprec = int(fabs(log10(fabs(val))))+2;
636          for(int i=newset.column("FPEAK").getPrecision();i<minprec;i++) newset.column("FPEAK").upPrec();
637        }
638        tempwidth = int( log10(fabs(val)) + 1) + newset.column("FPEAK").getPrecision() + 2;
639        if(val<0) tempwidth++;
640        for(int i=newset.column("FPEAK").getWidth();i<tempwidth;i++) newset.column("FPEAK").widen();
641
642        // S/N_peak
643        if(obj->getPeakSNR()>0.){
644            val = obj->getPeakSNR();
645            if((fabs(val) < 1.)&&(val>0.)){
646                minval = pow(10, -1. * (newset.column("SNRPEAK").getPrecision()+1));
647                if(val < minval) newset.column("SNRPEAK").upPrec();
648            }
649            tempwidth = int( log10(fabs(val)) + 1) + newset.column("SNRPEAK").getPrecision() +2;
650            if(val<0) tempwidth++;
651            for(int i=newset.column("SNRPEAK").getWidth();i<tempwidth;i++) newset.column("SNRPEAK").widen();
652        }
653
654        // X1 position
655        val = obj->getXmin() + obj->getXOffset();
656        tempwidth = int( log10(val) + 1) + newset.column("X1").getPrecision() + 1;
657        for(int i=newset.column("X1").getWidth();i<tempwidth;i++) newset.column("X1").widen();
658        // X2 position
659        val = obj->getXmax() + obj->getXOffset();
660        tempwidth = int( log10(val) + 1) + newset.column("X2").getPrecision() + 1;
661        for(int i=newset.column("X2").getWidth();i<tempwidth;i++) newset.column("X2").widen();
662        // Y1 position
663        val = obj->getYmin() + obj->getYOffset();
664        tempwidth = int( log10(val) + 1) + newset.column("Y1").getPrecision() + 1;
665        for(int i=newset.column("Y1").getWidth();i<tempwidth;i++) newset.column("Y1").widen();
666        // Y2 position
667        val = obj->getYmax() + obj->getYOffset();
668        tempwidth = int( log10(val) + 1) + newset.column("Y2").getPrecision() + 1;
669        for(int i=newset.column("Y2").getWidth();i<tempwidth;i++) newset.column("Y2").widen();
670        // Z1 position
671        val = obj->getZmin() + obj->getZOffset();
672        tempwidth = int( log10(val) + 1) + newset.column("Z1").getPrecision() + 1;
673        for(int i=newset.column("Z1").getWidth();i<tempwidth;i++) newset.column("Z1").widen();
674        // Z2 position
675        val = obj->getZmax() + obj->getZOffset();
676        tempwidth = int( log10(val) + 1) + newset.column("Z2").getPrecision() + 1;
677        for(int i=newset.column("Z2").getWidth();i<tempwidth;i++) newset.column("Z2").widen();
678
679        // Npix
680        tempwidth = int( log10(obj->getSize()) + 1) +
681          newset.column("NPIX").getPrecision() + 1;
682        for(int i=newset.column("NPIX").getWidth();i<tempwidth;i++) newset.column("NPIX").widen();
683   
684        // average X position
685        val = obj->getXaverage() + obj->getXOffset();
686        if((val<1.)&&(val>0.)){
687          minval = pow(10, -1. * (newset.column("XAV").getPrecision()+1));
688          if(val < minval) newset.column("XAV").upPrec();
689        }
690        tempwidth = int( log10(val) + 1) + newset.column("XAV").getPrecision() + 2;
691        for(int i=newset.column("XAV").getWidth();i<tempwidth;i++) newset.column("XAV").widen();
692        // average Y position
693        val = obj->getYaverage() + obj->getYOffset();
694        if((val<1.)&&(val>0.)){
695          minval = pow(10, -1. * (newset.column("YAV").getPrecision()+1));
696          if(val < minval) newset.column("YAV").upPrec();
697        }
698        tempwidth = int( log10(val) + 1) + newset.column("YAV").getPrecision() + 2;
699        for(int i=newset.column("YAV").getWidth();i<tempwidth;i++) newset.column("YAV").widen();
700        // average Z position
701        val = obj->getZaverage() + obj->getZOffset();
702        if((val<1.)&&(val>0.)){
703          minval = pow(10, -1. * (newset.column("ZAV").getPrecision()+1));
704          if((val>0.)&&(val < minval)) newset.column("ZAV").upPrec();
705        }
706        tempwidth = int( log10(val) + 1) + newset.column("ZAV").getPrecision() + 2;
707        for(int i=newset.column("ZAV").getWidth();i<tempwidth;i++) newset.column("ZAV").widen();
708   
709        // X position of centroid
710        val = obj->getXCentroid() + obj->getXOffset();
711        if((val<1.)&&(val>0.)){
712          minval = pow(10, -1. * (newset.column("XCENT").getPrecision()+1));
713          if(val < minval) newset.column("XCENT").upPrec();
714        }
715        tempwidth = int( log10(val) + 1) + newset.column("XCENT").getPrecision() + 2;
716        for(int i=newset.column("XCENT").getWidth();i<tempwidth;i++) newset.column("XCENT").widen();
717        // Y position of centroid
718        val = obj->getYCentroid() + obj->getYOffset();
719        if((val<1.)&&(val>0.)){
720          minval = pow(10, -1. * (newset.column("YCENT").getPrecision()+1));
721          if(val < minval) newset.column("YCENT").upPrec();
722        }
723        tempwidth = int( log10(val) + 1) + newset.column("YCENT").getPrecision() + 2;
724        for(int i=newset.column("YCENT").getWidth();i<tempwidth;i++) newset.column("YCENT").widen();
725        // Z position of centroid
726        val = obj->getZCentroid() + obj->getZOffset();
727        if((val<1.)&&(val>0.)){
728          minval = pow(10, -1. * (newset.column("ZCENT").getPrecision()+1));
729          if((val>0.)&&(val < minval)) newset.column("ZCENT").upPrec();
730        }
731        tempwidth = int( log10(val) + 1) + newset.column("ZCENT").getPrecision() + 2;
732        for(int i=newset.column("ZCENT").getWidth();i<tempwidth;i++) newset.column("ZCENT").widen();
733   
734        // X position of peak flux
735        val = obj->getXPeak() + obj->getXOffset();
736        if((val<1.)&&(val>0.)){
737          minval = pow(10, -1. * (newset.column("XPEAK").getPrecision()+1));
738          if(val < minval) newset.column("XPEAK").upPrec();
739        }
740        tempwidth = int( log10(val) + 1) + newset.column("XPEAK").getPrecision() + 2;
741        for(int i=newset.column("XPEAK").getWidth();i<tempwidth;i++) newset.column("XPEAK").widen();
742        // Y position of peak flux
743        val = obj->getYPeak() + obj->getYOffset();
744        if((val<1.)&&(val>0.)){
745          minval = pow(10, -1. * (newset.column("YPEAK").getPrecision()+1));
746          if(val < minval) newset.column("YPEAK").upPrec();
747        }
748        tempwidth = int( log10(val) + 1) + newset.column("YPEAK").getPrecision() + 2;
749        for(int i=newset.column("YPEAK").getWidth();i<tempwidth;i++) newset.column("YPEAK").widen();
750        // Z position of peak flux
751        val = obj->getZPeak() + obj->getZOffset();
752        if((val<1.)&&(val>0.)){
753          minval = pow(10, -1. * (newset.column("ZPEAK").getPrecision()+1));
754          if((val>0.)&&(val < minval)) newset.column("ZPEAK").upPrec();
755        }
756        tempwidth = int( log10(val) + 1) + newset.column("ZPEAK").getPrecision() + 2;
757        for(int i=newset.column("ZPEAK").getWidth();i<tempwidth;i++) newset.column("ZPEAK").widen();
758
759        //Number of contiguous channels
760        tempwidth = int( log10(obj->getMaxAdjacentChannels()) + 1) +
761          newset.column("NUMCH").getPrecision() + 1;
762        for(int i=newset.column("NUMCH").getWidth();i<tempwidth;i++) newset.column("NUMCH").widen();
763        //Spatial size
764        tempwidth = int( log10(obj->getSpatialSize()) + 1) +
765          newset.column("SPATSIZE").getPrecision() + 1;
766        for(int i=newset.column("SPATSIZE").getWidth();i<tempwidth;i++) newset.column("SPATSIZE").widen();
767       
768      }
769     
770      return newset;
771
772    }
773
774  }
775
776}
Note: See TracBrowser for help on using the repository browser.