source: trunk/src/Outputs/columns.cc @ 1101

Last change on this file since 1101 was 1101, checked in by MatthewWhiting, 12 years ago

Ticket #168: Updating the UCDs to make them unique across the catalogue, and to make the flux UCDs conditional on there being a spectral dimension (in which case we add spect.line.intensity). VOField interface has been updated as well.

File size: 30.8 KB
RevLine 
[300]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// -----------------------------------------------------------------------
[144]28#include <iostream>
29#include <sstream>
30
[155]31#include <math.h>
[393]32#include <duchamp/duchamp.hh>
33#include <duchamp/fitsHeader.hh>
[93]34#include <vector>
35#include <string>
[393]36#include <duchamp/Detection/detection.hh>
[1061]37#include <duchamp/Outputs/columns.hh>
[93]38
[378]39namespace duchamp
[232]40{
[221]41
[1061]42  namespace Catalogues
[272]43  {
[221]44
[439]45
[1061]46    Column::Column()
[378]47    {
[528]48      /// @brief Set the default values for the parameters.
[1061]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="";
[419]58    }
[232]59
[1061]60    Column::Column(std::string type)
61    {
62      operator=(defaultColumn(type));
63    }
[365]64
[1061]65    Column::Column(const Column& c)
[378]66    {
67      operator=(c);
68    }
[365]69
[1061]70    Column& Column::operator=(const Column& c)
[378]71    {
72      if(this==&c) return *this;
[1061]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;
[378]82      return *this;
[232]83    }
[378]84
[1064]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):
[1067]86      itsType(type), itsName(name), itsUnits(units), itsWidth(width), itsPrecision(prec), itsUCD(ucd), itsDatatype(datatype), itsColID(colID), itsExtraInfo(extraInfo)
[540]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
[1066]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)
[1061]98   }
[540]99
[378]100    //------------------------------------------------------------
[540]101
[1061]102    void Column::changePrec(int p)
[540]103    {
104      /// A direct way to alter the precision without having to use
[1061]105      /// Column::upPrec(). If the precision increases, the width will
[540]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
[1061]110      if(p > this->itsPrecision) this->itsWidth += (p-this->itsPrecision);
111      this->itsPrecision = p;
[540]112
113    }
114   
115    //------------------------------------------------------------
[1061]116    void Column::printTitle(std::ostream &stream)
[378]117    {
[1061]118      stream << std::setw(this->itsWidth)
[378]119             << std::setfill(' ')
[1061]120             << this->itsName;
[378]121    }
[93]122
[1061]123    void Column::printUnits(std::ostream &stream)
[378]124    {
[1061]125      stream << std::setw(this->itsWidth)
[378]126             << std::setfill(' ')
[1061]127             << this->itsUnits;
[378]128    }
[272]129 
[1061]130    void Column::printDash (std::ostream &stream)
[378]131    {
[1061]132      stream << std::setw(this->itsWidth)
[378]133             << std::setfill('-')
134             << ""
135             << std::setfill(' ');
136    }
[232]137
[1061]138    void Column::printBlank(std::ostream &stream)
[378]139    {
[1061]140      stream << std::setw(this->itsWidth)
[378]141             << std::setfill(' ')
142             << "";
143    }
[272]144 
[1061]145    template <class T> void Column::printEntry(std::ostream &stream, T value)
[378]146    {
[528]147      ///  \param stream Where the printing is done.
148      ///  \param value  The value to be printed.
149
[1061]150      stream << std::setprecision(this->itsPrecision)
151             << std::setw(this->itsWidth)
[378]152             << std::setfill(' ')
153             << value;
154    }
[1061]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);
[272]162
[440]163 
[1061]164    bool Column::doCol(DESTINATION tableType, bool flagFint)
[440]165    {
[1066]166      ///  @details Determines whether a given column is used for a
167      ///  given table type, based on the Column::type string.
[1047]168      /// \param tableType The type of table: one of FILE, SCREEN, LOG, VOTABLE.
[528]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.
[440]175     
[1066]176 
[1064]177      std::string FileList[34]={"NUM","NAME","X","Y","Z","RA","DEC","VEL","WRA","WDEC",
[1061]178                                "W50","W20","WVEL","FINT","FTOT","FPEAK","SNRPEAK",
179                                "X1","X2","Y1","Y2","Z1","Z2","NPIX","FLAG","XAV","YAV",
180                                "ZAV","XCENT","YCENT","ZCENT","XPEAK","YPEAK","ZPEAK"};
[1064]181      std::string ScreenList[21]={"NUM","NAME","X","Y","Z","RA","DEC","VEL","WRA","WDEC",
[1061]182                                  "W50","FPEAK","SNRPEAK","X1","X2","Y1",
183                                  "Y2","Z1","Z2","NPIX","FLAG"};
184      std::string LogList[16]={"NUM","X","Y","Z","FTOT","FPEAK","SNRPEAK",
185                               "X1","X2","Y1","Y2","Z1","Z2","NPIX","NUMCH","SPATSIZE"};
[1064]186      std::string VOList[22]={"NUM","NAME","RAJD","DECJD","VEL","WRA","WDEC",
[1061]187                            "W50","W20","WVEL","FPEAK","SNRPEAK","FLAG","XAV","YAV",
188                            "ZAV","XCENT","YCENT","ZCENT","XPEAK","YPEAK","ZPEAK"};
189
190      bool doIt=false;
191      if(tableType == FILE){
192        for(int i=0;i<34 && !doIt;i++) doIt = doIt || (this->itsType==FileList[i]);
193      }
[1047]194      else if(tableType == SCREEN){
[1061]195        if(this->itsType=="FTOT") doIt = !flagFint;
196        else if(this->itsType == "FINT") doIt = flagFint;
197        else for(int i=0;i<21 && !doIt;i++) doIt = doIt || (this->itsType==ScreenList[i]);
[440]198      }
[1061]199      else if(tableType == LOG){
200        for(int i=0;i<16 && !doIt;i++) doIt = doIt || (this->itsType==LogList[i]);
201      }
[1047]202      else if(tableType == VOTABLE){
[1061]203        if(this->itsType=="FTOT") doIt = !flagFint;
204        else if(this->itsType == "FINT") doIt = flagFint;
[1064]205        else for(int i=0;i<22 && !doIt;i++) doIt = doIt || (this->itsType==VOList[i]);
[440]206      }
[1061]207      return doIt;
208       
[440]209    }
210
[1061]211    Column defaultColumn(std::string type)
212    {
213      // type|name|units|width|prec|ucd|datatype|extraInfo
[1064]214      if (type == "NUM") return Column(type,"Obj#","",5,0,"meta.id","int","col_objnum","");
215      else if(type=="NAME") return Column(type,"Name","",8,0,"meta.id;meta.main","char","col_name","");
216      else if(type=="X") return Column(type,"X","",6,prXYZ,"pos.cartesian.x","float","col_x","");
217      else if(type=="Y") return Column(type,"Y","",6,prXYZ,"pos.cartesian.y","float","col_y","");
218      else if(type=="Z") return Column(type,"Z","",6,prXYZ,"pos.cartesian.z","float","col_z","");
219      else if(type=="RA") return Column(type,"RA","",11,0,"","char","col_ra",coordRef);
220      else if(type=="DEC") return Column(type,"DEC","",10,0,"","char","col_dec",coordRef);
221      else if(type=="RAJD") return Column(type,"RAJD","[deg]",11,prPOS,"","float","col_rajd",coordRef);
222      else if(type=="DECJD") return Column(type,"DECJD","[deg]",11,prPOS,"","float","col_decjd",coordRef);
223      else if(type=="VEL") return Column(type,"VEL","",7,prVEL,"","float","col_vel","");
224      else if(type=="WRA") return Column(type,"w_RA","[arcmin]",9,prWPOS,"","float","col_wra",coordRef);
225      else if(type=="WDEC") return Column(type,"w_DEC","[arcmin]",9,prWPOS,"","float","col_wdec",coordRef);
226      else if(type=="W50") return Column(type,"w_50","",7,prVEL,"","float","col_w50","");
227      else if(type=="W20") return Column(type,"w_20","",7,prVEL,"","float","col_w20","");
228      else if(type=="WVEL") return Column(type,"w_VEL","",7,prVEL,"","float","col_wvel","");
[1101]229      else if(type=="FINT") return Column(type,"F_int","",10,prFLUX,"phot.flux.density.integrated","float","col_fint","");
230      else if(type=="FTOT") return Column(type,"F_tot","",10,prFLUX,"phot.flux.density","float","col_ftot","");
231      else if(type=="FPEAK") return Column(type,"F_peak","",9,prFLUX,"phot.flux;stat.max","float","col_fpeak","");
232      else if(type=="SNRPEAK") return Column(type,"S/Nmax","",7,prSNR,"stat.snr;phot.flux","float","col_snrmax","");
[1064]233      else if(type=="X1") return Column(type,"X1","",4,0,"pos.cartesian.x;stat.min","int","col_x1","");
234      else if(type=="Y1") return Column(type,"Y1","",4,0,"pos.cartesian.y;stat.min","int","col_y1","");
235      else if(type=="Z1") return Column(type,"Z1","",4,0,"pos.cartesian.z;stat.min","int","col_z1","");
236      else if(type=="X2") return Column(type,"X2","",4,0,"pos.cartesian.x;stat.max","int","col_x2","");
237      else if(type=="Y2") return Column(type,"Y2","",4,0,"pos.cartesian.y;stat.max","int","col_y2","");
238      else if(type=="Z2") return Column(type,"Z2","",4,0,"pos.cartesian.z;stat.max","int","col_z2","");
239      else if(type=="NPIX") return Column(type,"Npix","[pix]",6,0,"","int","col_npix","");
240      else if(type=="FLAG") return Column(type,"Flag","",5,0,"meta.code.qual","char","col_flag","");
[1101]241      else if(type=="XAV") return Column(type,"X_av","",6,prXYZ,"pos.cartesian.x;stat.mean","float","col_xav","");
242      else if(type=="YAV") return Column(type,"Y_av","",6,prXYZ,"pos.cartesian.y;stat.mean","float","col_yav","");
243      else if(type=="ZAV") return Column(type,"Z_av","",6,prXYZ,"pos.cartesian.z;stat.mean","float","col_zav","");
244      else if(type=="XCENT") return Column(type,"X_cent","",7,prXYZ,"pos.cartesian.x;stat.centroid","float","col_xcent","");
245      else if(type=="YCENT") return Column(type,"Y_cent","",7,prXYZ,"pos.cartesian.y;stat.centroid","float","col_ycent","");
246      else if(type=="ZCENT") return Column(type,"Z_cent","",7,prXYZ,"pos.cartesian.z;stat.centroid","float","col_zcent","");
247      else if(type=="XPEAK") return Column(type,"X_peak","",7,prXYZ,"pos.cartesian.x;phot.flux;stat.max","int","col_xpeak","");
248      else if(type=="YPEAK") return Column(type,"Y_peak","",7,prXYZ,"pos.cartesian.y;phot.flux;stat.max","int","col_ypeak","");
249      else if(type=="ZPEAK") return Column(type,"Z_peak","",7,prXYZ,"pos.cartesian.z;phot.flux;stat.max","int","col_zpeak","");
250      else if(type=="NUMCH") return Column(type,"Nch","em.bin;stat.sum",6,0,"","int","col_nch","");
251      else if(type=="SPATSIZE") return Column(type,"Sp_size","instr.pixel;stat.sum",8,0,"","int","col_spsize","");
[1061]252      else {
253        Column col;
254        col.setType(type);
255        return col;
256      }
[221]257
[1061]258    }
259
[1064]260    CatalogueSpecification getFullColSet(std::vector<Detection> &objectList, FitsHeader &head)
[378]261    {
[1064]262      ///  @details A function that returns a catalogue specification
263      ///  containing information on the columns necessary for output
264      ///  to the results file:
265      ///  Obj#,NAME,X,Y,Z,RA,DEC,VEL,w_RA,w_DEC,w_VEL,F_tot,F_int,F_peak,
266      ///  X1,X2,Y1,Y2,Z1,Z2,Npix,Flag,
267      ///  XAV,YAV,ZAV,XCENT,YCENT,ZCENT,XPEAK,YPEAK,ZPEAK
[528]268      ///
269      ///   Each object in the provided objectList is checked to see if it
270      ///   requires any column to be widened, or for that column to have
271      ///   its precision increased.
272      ///
273      ///   Both Ftot and Fint are provided -- it is up to the calling
274      ///   function to determine which to use.
275      ///
276      /// \param objectList A std::vector list of Detection objects that the
277      /// columns need to fit.
278      /// \param head The FitsHeader object defining the World Coordinate
279      /// System.
[1064]280      /// \return A CatalogueSpecification
[93]281
[1064]282      CatalogueSpecification newset;
[93]283
[438]284      // desired precisions for fluxes, velocities and SNR value
[461]285      int precVel,precFlux,precSNR;
286      if(objectList.size()>0){
287        precVel=objectList[0].getVelPrec();
288        precFlux=objectList[0].getFpeakPrec(); // same as FintPrec at this point.
[438]289        precSNR=objectList[0].getSNRPrec();
[461]290      }
291      else {
292        precVel = prVEL;
293        precFlux = prFLUX;
294        precSNR = prSNR;
295      }
[438]296
[378]297      // set up the default columns
[1064]298      newset.addColumn( Column("NUM") );
299      newset.addColumn( Column("NAME") );
300      newset.addColumn( Column("X") );
301      newset.addColumn( Column("Y") );
302      newset.addColumn( Column("Z") );
303      newset.addColumn( Column("RA") );
304      newset.addColumn( Column("DEC") );
305      newset.addColumn( Column("RAJD") );
306      newset.addColumn( Column("DECJD") );
307      newset.addColumn( Column("VEL") );
308      newset.addColumn( Column("WRA") );
309      newset.addColumn( Column("WDEC") );
310      newset.addColumn( Column("W50") );
311      newset.addColumn( Column("W20") );
312      newset.addColumn( Column("WVEL") );
313      newset.addColumn( Column("FINT") );
314      newset.addColumn( Column("FTOT") );
315      newset.addColumn( Column("FPEAK") );
316      newset.addColumn( Column("SNRPEAK") );
317      newset.addColumn( Column("X1") );
318      newset.addColumn( Column("X2") );
319      newset.addColumn( Column("Y1") );
320      newset.addColumn( Column("Y2") );
321      newset.addColumn( Column("Z1") );
322      newset.addColumn( Column("Z2") );
323      newset.addColumn( Column("NPIX") );
324      newset.addColumn( Column("FLAG") );
325      newset.addColumn( Column("XAV") );
326      newset.addColumn( Column("YAV") );
327      newset.addColumn( Column("ZAV") );
328      newset.addColumn( Column("XCENT") );
329      newset.addColumn( Column("YCENT") );
330      newset.addColumn( Column("ZCENT") );
331      newset.addColumn( Column("XPEAK") );
332      newset.addColumn( Column("YPEAK") );
333      newset.addColumn( Column("ZPEAK") );
334      newset.addColumn( Column("NUMCH") );
335      newset.addColumn( Column("SPATSIZE") );
[93]336
[438]337
[378]338      // Now test each object against each new column
[623]339      std::vector<Detection>::iterator obj;
340      for(obj = objectList.begin(); obj < objectList.end(); obj++){
[378]341        std::string tempstr;
342        int tempwidth;
343        float val,minval;
[953]344        double valD,minvalD;
[93]345
[378]346        // Obj#
[623]347        tempwidth = int( log10(obj->getID()) + 1) + 1;
[1064]348        for(int i=newset.column("NUM").getWidth();i<tempwidth;i++) newset.column("NUM").widen();
[93]349
[378]350        // Name
[623]351        tempwidth = obj->getName().size() + 1;
[1064]352        for(int i=newset.column("NAME").getWidth();i<tempwidth;i++) newset.column("NAME").widen();
[93]353
[378]354        // X position
[623]355        val = obj->getXcentre() + obj->getXOffset();
[378]356        if((val<1.)&&(val>0.)){
[1064]357          minval = pow(10, -1. * (newset.column("X").getPrecision()+1));
358          if(val < minval) newset.column("X").upPrec();
[378]359        }
[1064]360        tempwidth = int( log10(val) + 1) + newset.column("X").getPrecision() + 2;
361        for(int i=newset.column("X").getWidth();i<tempwidth;i++) newset.column("X").widen();
[378]362        // Y position
[623]363        val = obj->getYcentre() + obj->getYOffset();
[378]364        if((val<1.)&&(val>0.)){
[1064]365          minval = pow(10, -1. * (newset.column("Y").getPrecision()+1));
366          if(val < minval) newset.column("Y").upPrec();
[378]367        }
[1064]368        tempwidth = int( log10(val) + 1) + newset.column("Y").getPrecision() + 2;
369        for(int i=newset.column("Y").getWidth();i<tempwidth;i++) newset.column("Y").widen();
[378]370        // Z position
[623]371        val = obj->getZcentre() + obj->getZOffset();
[378]372        if((val<1.)&&(val>0.)){
[1064]373          minval = pow(10, -1. * (newset.column("Z").getPrecision()+1));
374          if((val>0.)&&(val < minval)) newset.column("Z").upPrec();
[378]375        }
[1064]376        tempwidth = int( log10(val) + 1) + newset.column("Z").getPrecision() + 2;
377        for(int i=newset.column("Z").getWidth();i<tempwidth;i++) newset.column("Z").widen();
[93]378
[378]379        if(head.isWCS()){ 
380          // RA -- assign correct title. Check width but should be ok by definition
[963]381          tempstr = head.lngtype();
[1064]382          newset.column("RA").setName(tempstr);
[623]383          tempwidth = obj->getRAs().size() + 1;
[1064]384          for(int i=newset.column("RA").getWidth();i<tempwidth;i++) newset.column("RA").widen();
[144]385     
[378]386          // Dec -- assign correct title. Check width, should be ok by definition
[963]387          tempstr = head.lattype();
[1064]388          newset.column("DEC").setName(tempstr);
[623]389          tempwidth = obj->getDecs().size() + 1;
[1064]390          for(int i=newset.column("DEC").getWidth();i<tempwidth;i++) newset.column("DEC").widen();
[93]391
[439]392          // RA decimal degrees -- assign correct title. Check width but should be OK
[963]393          tempstr = head.lngtype();
[1064]394          newset.column("RAJD").setName(tempstr);
[953]395          valD = obj->getRA();
[1064]396          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("RAJD").getPrecision() + 2;
397          for(int i=newset.column("RAJD").getWidth();i<tempwidth;i++) newset.column("RAJD").widen();
[439]398     
399          // Dec decimal degrees -- assign correct title. Check width but should be OK
[963]400          tempstr = head.lattype();
[1064]401          newset.column("DECJD").setName(tempstr);
[953]402          valD = obj->getDec();
[1064]403          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("DECJD").getPrecision() + 2;
404          for(int i=newset.column("DECJD").getWidth();i<tempwidth;i++) newset.column("DECJD").widen();
[439]405
[378]406          // Vel -- check width, title and units.
407          if(head.canUseThirdAxis()){
[479]408            if(head.WCS().spec < 0)  // if it's not a spectral axis
[1064]409              newset.column("VEL").setName( head.WCS().ctype[2] );
[947]410            else
[1064]411              newset.column("VEL").setName(head.getSpectralType());
412            tempwidth = newset.column("VEL").getName().size() + 1;
413            for(int i=newset.column("VEL").getWidth();i<tempwidth;i++) newset.column("VEL").widen();
[378]414            if(head.getSpectralUnits().size()>0)
[1064]415              newset.column("VEL").setUnits("[" + head.getSpectralUnits() + "]");
416            tempwidth = newset.column("VEL").getUnits().size() + 1;
417            for(int i=newset.column("VEL").getWidth();i<tempwidth;i++) newset.column("VEL").widen();
[271]418       
[953]419            valD = obj->getVel();
420            if((fabs(valD) < 1.)&&(valD>0.)){
[1064]421              minvalD = pow(10, -1. * (newset.column("VEL").getPrecision()+1));
422              if(valD < minvalD) newset.column("VEL").upPrec();
[378]423            }
[1064]424            tempwidth = int(log10(fabs(valD)) + 1) + newset.column("VEL").getPrecision() + 2;
[953]425            if(valD<0) tempwidth++;
[1064]426            for(int i=newset.column("VEL").getWidth();i<tempwidth;i++) newset.column("VEL").widen();
[378]427          }
[93]428
[378]429          // w_RA -- check width & title. leave units for the moment.
[963]430          tempstr = "w_" + head.lngtype();
[1064]431          newset.column("WRA").setName(tempstr);
432          tempwidth = newset.column("RA").getUnits().size() + 1;
433          for(int i=newset.column("WRA").getWidth();i<tempwidth;i++) newset.column("WRA").widen();
[953]434          valD = obj->getRAWidth();
435          if((fabs(valD) < 1.)&&(valD>0.)){
[1064]436            minvalD = pow(10, -1. * (newset.column("WRA").getPrecision()+1));
437            if(valD < minvalD) newset.column("WRA").upPrec();
[378]438          }
[1064]439          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("WRA").getPrecision() + 2;
[953]440          if(valD<0) tempwidth++;
[1064]441          for(int i=newset.column("WRA").getWidth();i<tempwidth;i++) newset.column("WRA").widen();
[136]442
[378]443          // w_DEC -- check width & title. leave units for the moment.
[963]444          tempstr = "w_" + head.lattype();
[1064]445          newset.column("WDEC").setName(tempstr);
446          tempwidth = newset.column("DEC").getUnits().size() + 1;
447          for(int i=newset.column("WDEC").getWidth();i<tempwidth;i++) newset.column("WDEC").widen();
[953]448          valD = obj->getDecWidth();
449          if((fabs(valD) < 1.)&&(valD>0.)){
[1064]450            minvalD = pow(10, -1. * (newset.column("WDEC").getPrecision()+1));
451            if(valD < minvalD) newset.column("WDEC").upPrec();
[378]452          }
[1064]453          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("WDEC").getPrecision() + 2;
[953]454          if(valD<0) tempwidth++;
[1064]455          for(int i=newset.column("WDEC").getWidth();i<tempwidth;i++) newset.column("WDEC").widen();
[136]456
[463]457          // w_50 -- check width, title and units.
458          if(head.canUseThirdAxis()){
459            if(head.getSpectralUnits().size()>0)
[1064]460              newset.column("W50").setUnits("[" + head.getSpectralUnits() + "]");
461            tempwidth = newset.column("W50").getUnits().size() + 1;
462            for(int i=newset.column("W50").getWidth();i<tempwidth;i++) newset.column("W50").widen();
[953]463            valD = obj->getW50();
464            if((fabs(valD) < 1.)&&(valD>0.)){
[1064]465              minvalD = pow(10, -1. * (newset.column("W50").getPrecision()+1));
466              if(valD < minvalD) newset.column("W50").upPrec();
[463]467            }
[1064]468            tempwidth = int( log10(fabs(valD)) + 1) + newset.column("W50").getPrecision() + 2;
[953]469            if(valD<0) tempwidth++;
[1064]470            for(int i=newset.column("W50").getWidth();i<tempwidth;i++) newset.column("W50").widen();
[463]471          }
472
473          // w_20 -- check width, title and units.
474          if(head.canUseThirdAxis()){
475            if(head.getSpectralUnits().size()>0)
[1064]476              newset.column("W20").setUnits("[" + head.getSpectralUnits() + "]");
477            tempwidth = newset.column("W20").getUnits().size() + 1;
478            for(int i=newset.column("W20").getWidth();i<tempwidth;i++)newset.column("W20").widen();
[953]479            valD = obj->getW20();
480            if((fabs(valD) < 1.)&&(valD>0.)){
[1064]481              minvalD = pow(10, -1. * (newset.column("W20").getPrecision()+1));
482              if(valD < minvalD) newset.column("W20").upPrec();
[463]483            }
[1064]484            tempwidth = int( log10(fabs(valD)) + 1) + newset.column("W20").getPrecision() + 2;
[953]485            if(valD<0) tempwidth++;
[1064]486            for(int i=newset.column("W20").getWidth();i<tempwidth;i++) newset.column("W20").widen();
[463]487          }
488
[378]489          // w_Vel -- check width, title and units.
490          if(head.canUseThirdAxis()){
[479]491            if(head.WCS().spec < 0) // if it's not a spectral axis
[1064]492              newset.column("WVEL").setName( std::string("w_") + head.WCS().ctype[2] );
[947]493            else
[1064]494              newset.column("WVEL").setName(std::string("w_")+head.getSpectralType());
[378]495            if(head.getSpectralUnits().size()>0)
[1064]496              newset.column("WVEL").setUnits("[" + head.getSpectralUnits() + "]");
497            tempwidth = newset.column("WVEL").getUnits().size() + 1;
498            for(int i=newset.column("WVEL").getWidth();i<tempwidth;i++)newset.column("WVEL").widen();
499            tempwidth = newset.column("WVEL").getName().size() + 1;
500            for(int i=newset.column("WVEL").getWidth();i<tempwidth;i++) newset.column("WVEL").widen();
[953]501            valD = obj->getVelWidth();
502            if((fabs(valD) < 1.)&&(valD>0.)){
[1064]503              minvalD = pow(10, -1. * (newset.column("WVEL").getPrecision()+1));
504              if(valD < minvalD) newset.column("WVEL").upPrec();
[378]505            }
[1064]506            tempwidth = int( log10(fabs(valD)) + 1) + newset.column("WVEL").getPrecision() + 2;
[953]507            if(valD<0) tempwidth++;
[1064]508            for(int i=newset.column("WVEL").getWidth();i<tempwidth;i++) newset.column("WVEL").widen();
[378]509          }
[136]510
[378]511          // F_int -- check width & units
512          if(head.getIntFluxUnits().size()>0)
[1064]513            newset.column("FINT").setUnits("[" + head.getIntFluxUnits() + "]");
514          tempwidth = newset.column("FINT").getUnits().size() + 1;
515          for(int i=newset.column("FINT").getWidth();i<tempwidth;i++) newset.column("FINT").widen();
[953]516          valD = obj->getIntegFlux();
517          if((fabs(valD) < 1.)// &&(valD>0.)
[378]518             ){
[953]519            int minprec = int(fabs(log10(fabs(valD))))+2;
[1064]520            for(int i=newset.column("FINT").getPrecision();i<minprec;i++) newset.column("FINT").upPrec();
[378]521          }
[1064]522          tempwidth = int( log10(fabs(valD)) + 1) + newset.column("FINT").getPrecision() + 2;
[953]523          if(valD<0) tempwidth++;
[1064]524          for(int i=newset.column("FINT").getWidth();i<tempwidth;i++) newset.column("FINT").widen();
[365]525     
[378]526        }
[136]527
[378]528        // F_tot
[1064]529        newset.column("FTOT").setUnits("[" + head.getFluxUnits() + "]");
530        tempwidth = newset.column("FTOT").getUnits().size() + 1;
531        for(int i=newset.column("FTOT").getWidth();i<tempwidth;i++) newset.column("FTOT").widen();
[623]532        val = obj->getTotalFlux();
[378]533        //     std::cerr << val << "\n";
534        if((fabs(val) < 1.)// &&(val>0.)
535           ){
536          int minprec = int(fabs(log10(fabs(val))))+2;
[1064]537          for(int i=newset.column("FTOT").getPrecision();i<minprec;i++) newset.column("FTOT").upPrec();
[378]538        }
[1064]539        tempwidth = int( log10(fabs(val)) + 1) + newset.column("FTOT").getPrecision() + 2;
[378]540        if(val<0) tempwidth++;
[1064]541        for(int i=newset.column("FTOT").getWidth();i<tempwidth;i++) newset.column("FTOT").widen();
[136]542
[378]543        // F_peak
[1064]544        newset.column("FPEAK").setUnits("[" + head.getFluxUnits() + "]");
545        tempwidth = newset.column("FPEAK").getUnits().size() + 1;
546        for(int i=newset.column("FPEAK").getWidth();i<tempwidth;i++) newset.column("FPEAK").widen();
[623]547        val = obj->getPeakFlux();
[378]548        if((fabs(val) < 1.)// &&(val>0.)
549           ){
550          int minprec = int(fabs(log10(fabs(val))))+2;
[1064]551          for(int i=newset.column("FPEAK").getPrecision();i<minprec;i++) newset.column("FPEAK").upPrec();
[378]552        }
[1064]553        tempwidth = int( log10(fabs(val)) + 1) + newset.column("FPEAK").getPrecision() + 2;
[378]554        if(val<0) tempwidth++;
[1064]555        for(int i=newset.column("FPEAK").getWidth();i<tempwidth;i++) newset.column("FPEAK").widen();
[136]556
[378]557        // S/N_peak
[623]558        val = obj->getPeakSNR();
[378]559        if((fabs(val) < 1.)&&(val>0.)){
[1064]560          minval = pow(10, -1. * (newset.column("SNRPEAK").getPrecision()+1));
561          if(val < minval) newset.column("SNRPEAK").upPrec();
[378]562        }
[1064]563        tempwidth = int( log10(fabs(val)) + 1) + newset.column("SNRPEAK").getPrecision() +2;
[378]564        if(val<0) tempwidth++;
[1064]565        for(int i=newset.column("SNRPEAK").getWidth();i<tempwidth;i++) newset.column("SNRPEAK").widen();
[191]566
[378]567        // X1 position
[623]568        val = obj->getXmin() + obj->getXOffset();
[1064]569        tempwidth = int( log10(val) + 1) + newset.column("X1").getPrecision() + 1;
570        for(int i=newset.column("X1").getWidth();i<tempwidth;i++) newset.column("X1").widen();
[378]571        // X2 position
[623]572        val = obj->getXmax() + obj->getXOffset();
[1064]573        tempwidth = int( log10(val) + 1) + newset.column("X2").getPrecision() + 1;
574        for(int i=newset.column("X2").getWidth();i<tempwidth;i++) newset.column("X2").widen();
[378]575        // Y1 position
[623]576        val = obj->getYmin() + obj->getYOffset();
[1064]577        tempwidth = int( log10(val) + 1) + newset.column("Y1").getPrecision() + 1;
578        for(int i=newset.column("Y1").getWidth();i<tempwidth;i++) newset.column("Y1").widen();
[378]579        // Y2 position
[623]580        val = obj->getYmax() + obj->getYOffset();
[1064]581        tempwidth = int( log10(val) + 1) + newset.column("Y2").getPrecision() + 1;
582        for(int i=newset.column("Y2").getWidth();i<tempwidth;i++) newset.column("Y2").widen();
[378]583        // Z1 position
[623]584        val = obj->getZmin() + obj->getZOffset();
[1064]585        tempwidth = int( log10(val) + 1) + newset.column("Z1").getPrecision() + 1;
586        for(int i=newset.column("Z1").getWidth();i<tempwidth;i++) newset.column("Z1").widen();
[378]587        // Z2 position
[623]588        val = obj->getZmax() + obj->getZOffset();
[1064]589        tempwidth = int( log10(val) + 1) + newset.column("Z2").getPrecision() + 1;
590        for(int i=newset.column("Z2").getWidth();i<tempwidth;i++) newset.column("Z2").widen();
[136]591
[378]592        // Npix
[623]593        tempwidth = int( log10(obj->getSize()) + 1) +
[1064]594          newset.column("NPIX").getPrecision() + 1;
595        for(int i=newset.column("NPIX").getWidth();i<tempwidth;i++) newset.column("NPIX").widen();
[144]596   
[378]597        // average X position
[623]598        val = obj->getXaverage() + obj->getXOffset();
[378]599        if((val<1.)&&(val>0.)){
[1064]600          minval = pow(10, -1. * (newset.column("XAV").getPrecision()+1));
601          if(val < minval) newset.column("XAV").upPrec();
[378]602        }
[1064]603        tempwidth = int( log10(val) + 1) + newset.column("XAV").getPrecision() + 2;
604        for(int i=newset.column("XAV").getWidth();i<tempwidth;i++) newset.column("XAV").widen();
[378]605        // average Y position
[623]606        val = obj->getYaverage() + obj->getYOffset();
[378]607        if((val<1.)&&(val>0.)){
[1064]608          minval = pow(10, -1. * (newset.column("YAV").getPrecision()+1));
609          if(val < minval) newset.column("YAV").upPrec();
[378]610        }
[1064]611        tempwidth = int( log10(val) + 1) + newset.column("YAV").getPrecision() + 2;
612        for(int i=newset.column("YAV").getWidth();i<tempwidth;i++) newset.column("YAV").widen();
[378]613        // average Z position
[623]614        val = obj->getZaverage() + obj->getZOffset();
[378]615        if((val<1.)&&(val>0.)){
[1064]616          minval = pow(10, -1. * (newset.column("ZAV").getPrecision()+1));
617          if((val>0.)&&(val < minval)) newset.column("ZAV").upPrec();
[378]618        }
[1064]619        tempwidth = int( log10(val) + 1) + newset.column("ZAV").getPrecision() + 2;
620        for(int i=newset.column("ZAV").getWidth();i<tempwidth;i++) newset.column("ZAV").widen();
[265]621   
[378]622        // X position of centroid
[623]623        val = obj->getXCentroid() + obj->getXOffset();
[378]624        if((val<1.)&&(val>0.)){
[1064]625          minval = pow(10, -1. * (newset.column("XCENT").getPrecision()+1));
626          if(val < minval) newset.column("XCENT").upPrec();
[378]627        }
[1064]628        tempwidth = int( log10(val) + 1) + newset.column("XCENT").getPrecision() + 2;
629        for(int i=newset.column("XCENT").getWidth();i<tempwidth;i++) newset.column("XCENT").widen();
[378]630        // Y position of centroid
[623]631        val = obj->getYCentroid() + obj->getYOffset();
[378]632        if((val<1.)&&(val>0.)){
[1064]633          minval = pow(10, -1. * (newset.column("YCENT").getPrecision()+1));
634          if(val < minval) newset.column("YCENT").upPrec();
[378]635        }
[1064]636        tempwidth = int( log10(val) + 1) + newset.column("YCENT").getPrecision() + 2;
637        for(int i=newset.column("YCENT").getWidth();i<tempwidth;i++) newset.column("YCENT").widen();
[378]638        // Z position of centroid
[623]639        val = obj->getZCentroid() + obj->getZOffset();
[378]640        if((val<1.)&&(val>0.)){
[1064]641          minval = pow(10, -1. * (newset.column("ZCENT").getPrecision()+1));
642          if((val>0.)&&(val < minval)) newset.column("ZCENT").upPrec();
[378]643        }
[1064]644        tempwidth = int( log10(val) + 1) + newset.column("ZCENT").getPrecision() + 2;
645        for(int i=newset.column("ZCENT").getWidth();i<tempwidth;i++) newset.column("ZCENT").widen();
[265]646   
[378]647        // X position of peak flux
[623]648        val = obj->getXPeak() + obj->getXOffset();
[378]649        if((val<1.)&&(val>0.)){
[1064]650          minval = pow(10, -1. * (newset.column("XPEAK").getPrecision()+1));
651          if(val < minval) newset.column("XPEAK").upPrec();
[378]652        }
[1064]653        tempwidth = int( log10(val) + 1) + newset.column("XPEAK").getPrecision() + 2;
654        for(int i=newset.column("XPEAK").getWidth();i<tempwidth;i++) newset.column("XPEAK").widen();
[378]655        // Y position of peak flux
[623]656        val = obj->getYPeak() + obj->getYOffset();
[378]657        if((val<1.)&&(val>0.)){
[1064]658          minval = pow(10, -1. * (newset.column("YPEAK").getPrecision()+1));
659          if(val < minval) newset.column("YPEAK").upPrec();
[378]660        }
[1064]661        tempwidth = int( log10(val) + 1) + newset.column("YPEAK").getPrecision() + 2;
662        for(int i=newset.column("YPEAK").getWidth();i<tempwidth;i++) newset.column("YPEAK").widen();
[378]663        // Z position of peak flux
[623]664        val = obj->getZPeak() + obj->getZOffset();
[378]665        if((val<1.)&&(val>0.)){
[1064]666          minval = pow(10, -1. * (newset.column("ZPEAK").getPrecision()+1));
667          if((val>0.)&&(val < minval)) newset.column("ZPEAK").upPrec();
[378]668        }
[1064]669        tempwidth = int( log10(val) + 1) + newset.column("ZPEAK").getPrecision() + 2;
670        for(int i=newset.column("ZPEAK").getWidth();i<tempwidth;i++) newset.column("ZPEAK").widen();
[688]671
672        //Number of contiguous channels
673        tempwidth = int( log10(obj->getMaxAdjacentChannels()) + 1) +
[1064]674          newset.column("NUMCH").getPrecision() + 1;
675        for(int i=newset.column("NUMCH").getWidth();i<tempwidth;i++) newset.column("NUMCH").widen();
[688]676        //Spatial size
677        tempwidth = int( log10(obj->getSpatialSize()) + 1) +
[1064]678          newset.column("SPATSIZE").getPrecision() + 1;
679        for(int i=newset.column("SPATSIZE").getWidth();i<tempwidth;i++) newset.column("SPATSIZE").widen();
[688]680       
[378]681      }
[688]682     
[378]683      return newset;
[136]684
[378]685    }
[136]686
[378]687  }
688
[93]689}
Note: See TracBrowser for help on using the repository browser.