old_datasources/healpix/healpix.cpp

0001 /***************************************************************************
0002  *                                                                         *
0003  *   copyright : (C) 2007 The University of Toronto                        *
0004  *                   netterfield@astro.utoronto.ca                         *
0005  *                                                                         *
0006  *   This program is free software; you can redistribute it and/or modify  *
0007  *   it under the terms of the GNU General Public License as published by  *
0008  *   the Free Software Foundation; either version 2 of the License, or     *
0009  *   (at your option) any later version.                                   *
0010  *                                                                         *
0011  ***************************************************************************/
0012
0013 #include "healpix.h"
0014
0015 #include <assert.h>
0016 #include <QXmlStreamWriter>
0017 #include <math.h>
0018
0019 #include "ui_healpixconfig.h"
0020
0021 #define DEFAULT_XDIM 800
0022 #define DEFAULT_YDIM 600
0023 #define HPUNIT_RAD 0
0024 #define HPUNIT_DEG 1
0025 #define HPUNIT_RADEC 2
0026 #define HPUNIT_LATLON 3
0027
0028 #include "kst_i18n.h"
0029
0030 static const QString healpixTypeString = I18N_NOOP("HEALPIX image");
0031
0032 class HealpixSource::Config {
0033   public:
0034     Config() {
0035       _nX = 800;
0036       _nY = 600;
0037       _autoTheta = true;
0038       _autoPhi = true;
0039       //FIXME switch to radians default later
0040       _thetaUnits = HPUNIT_RADEC;
0041       _phiUnits = HPUNIT_RADEC;
0042       _vecUnits = HPUNIT_RADEC;
0043       _autoMag = true;
0044       _vecDegrade = 0;
0045       _vecTheta = 0;
0046       _vecPhi = 0;
0047       _vecQU = false;
0048     }
0049
0050     void read(QSettings *cfg, const QString& fileName = QString::null) {
0051       Q_UNUSED(fileName);
0052       cfg->beginGroup(healpixTypeString);
0053       double confThetaMin;
0054       double confThetaMax;
0055       double confPhiMin;
0056       double confPhiMax;
0057       int tempdegrade;
0058
0059       if (!fileName.isEmpty()) {
0060         cfg->endGroup();
0061         cfg->beginGroup(fileName);
0062       }
0063       _nX = cfg->value("Matrix X Dimension", DEFAULT_XDIM).toInt();
0064       _nY = cfg->value("Matrix Y Dimension", DEFAULT_YDIM).toInt();
0065       _autoTheta = cfg->value("Theta Autoscale", true).toBool();
0066       _thetaUnits = cfg->value("Theta Units", HPUNIT_RADEC).toInt();
0067       confThetaMin = cfg->value("Theta Min", 0).toDouble();
0068       confThetaMax = cfg->value("Theta Max", 0).toDouble();
0069       _autoPhi = cfg->value("Phi Autoscale", true).toBool();
0070       _phiUnits = cfg->value("Phi Units", HPUNIT_RADEC).toInt();
0071       confPhiMin = cfg->value("Phi Min", 0).toDouble();
0072       confPhiMax = cfg->value("Phi Max", 0).toDouble();
0073       _vecTheta = cfg->value("Vector Theta", 0).toInt();
0074       _vecPhi = cfg->value("Vector Phi", 0).toInt();
0075       tempdegrade = cfg->value("Vector Degrade Factor", 1).toInt();
0076       _autoMag = cfg->value("Vector Magnitude Autoscale", true).toBool();
0077       _maxMag = cfg->value("Vector Max Magnitude", 0).toDouble();
0078       _vecQU = cfg->value("Vector is QU", false).toBool();
0079
0080       // check degrade factor
0081       checkDegrade(tempdegrade);
0082       _vecDegrade = tempdegrade;
0083
0084       // convert the entered range values into radians and
0085       // force them to the correct range.
0086       theta2Internal(_thetaUnits, confThetaMin);
0087       theta2Internal(_thetaUnits, confThetaMax);
0088       phi2Internal(_phiUnits, confPhiMin);
0089       phi2Internal(_phiUnits, confPhiMax);
0090
0091       // swap theta min/max if coordinate system requires it
0092       if (confThetaMax < confThetaMin) {
0093         double temp = confThetaMax;
0094         confThetaMax = confThetaMin;
0095         confThetaMin = temp;
0096       }
0097
0098       _thetaMin = confThetaMin;
0099       _thetaMax = confThetaMax;
0100       _phiMin = confPhiMin;
0101       _phiMax = confPhiMax;
0102
0103       cfg->endGroup();
0104     }
0105
0106     void save(QXmlStreamWriter& s) {
0107       double confThetaMin = _thetaMin;
0108       double confThetaMax = _thetaMax;
0109       double confPhiMin = _phiMin;
0110       double confPhiMax = _phiMax;
0111
0112       // export the internal range (in radians) to the
0113       // selected coordinate system
0114       theta2External(_thetaUnits, confThetaMin);
0115       theta2External(_thetaUnits, confThetaMax);
0116       phi2External(_phiUnits, confPhiMin);
0117       phi2External(_phiUnits, confPhiMax);
0118
0119       // swap theta min/max if coordinate system requires it
0120       if (confThetaMax < confThetaMin) {
0121         double temp = confThetaMax;
0122         confThetaMax = confThetaMin;
0123         confThetaMin = temp;
0124       }
0125       s.writeStartElement("properties");
0126       s.writeAttribute("dim-x", QString::number(_nX));
0127       s.writeAttribute("dim-y", QString::number(_nY));
0128       s.writeAttribute("theta-auto", QVariant(_autoTheta).toString());
0129       s.writeAttribute("theta-units", QString::number(_thetaUnits));
0130       s.writeAttribute("theta-min", QString::number(confThetaMin));
0131       s.writeAttribute("theta-max", QString::number(confThetaMax));
0132       s.writeAttribute("phi-auto", QVariant(_autoPhi).toString());
0133       s.writeAttribute("phi-units", QString::number(_phiUnits));
0134       s.writeAttribute("phi-min", QString::number(confPhiMin));
0135       s.writeAttribute("phi-max", QString::number(confPhiMax));
0136       s.writeAttribute("vector-theta", QString::number(_vecTheta));
0137       s.writeAttribute("vector-phi", QString::number(_vecPhi));
0138       s.writeAttribute("vector-degrade", QVariant(_vecDegrade).toString());
0139       s.writeAttribute("vector-auto", QVariant(_autoMag).toString());
0140       s.writeAttribute("vector-max", QString::number(_maxMag));
0141       s.writeAttribute("vector-QU", QVariant(_vecQU).toString());
0142       s.writeEndElement();
0143     }
0144
0145     void parseProperties(QXmlStreamAttributes &properties) {
0146       double confThetaMin;
0147       double confThetaMax;
0148       double confPhiMin;
0149       double confPhiMax;
0150       int tempdegrade;
0151
0152       _nX = properties.value("dim-x").toString().toInt();
0153       _nY = properties.value("dim-y").toString().toInt();
0154       _autoTheta = QVariant(properties.value("theta-auto").toString()).toBool();
0155       _thetaUnits = properties.value("theta-units").toString().toInt();
0156       confThetaMin = properties.value("theta-min").toString().toDouble();
0157       confThetaMax = properties.value("theta-max").toString().toDouble();
0158       _autoPhi = QVariant(properties.value("phi-auto").toString()).toBool();
0159       _phiUnits = properties.value("phi-units").toString().toInt();
0160       confPhiMin = properties.value("phi-min").toString().toDouble();
0161       confPhiMax = properties.value("phi-max").toString().toDouble();
0162       _vecTheta = properties.value("vector-theta").toString().toInt();
0163       _vecPhi = properties.value("vector-phi").toString().toInt();
0164       tempdegrade = QVariant(properties.value("vector-degrade").toString()).toBool();
0165       _autoMag = QVariant(properties.value("vector-auto").toString()).toBool();
0166       _vecQU = QVariant(properties.value("vector-QU").toString()).toBool();
0167       _maxMag = properties.value("vector-max").toString().toDouble();
0168
0169       checkDegrade(tempdegrade);
0170       _vecDegrade = tempdegrade;
0171
0172       // convert the entered range values into radians and
0173       // force them to the correct range.
0174       theta2Internal(_thetaUnits, confThetaMin);
0175       theta2Internal(_thetaUnits, confThetaMax);
0176       phi2Internal(_phiUnits, confPhiMin);
0177       phi2Internal(_phiUnits, confPhiMax);
0178
0179       // swap theta min/max if coordinate system requires it
0180       if (confThetaMax < confThetaMin) {
0181         double temp = confThetaMax;
0182         confThetaMax = confThetaMin;
0183         confThetaMin = temp;
0184       }
0185
0186       _thetaMin = confThetaMin;
0187       _thetaMax = confThetaMax;
0188       _phiMin = confPhiMin;
0189       _phiMax = confPhiMax;
0190     }
0191
0192
0193     void load(const QDomElement& e) {
0194       double confThetaMin;
0195       double confThetaMax;
0196       double confPhiMin;
0197       double confPhiMax;
0198       int tempdegrade;
0199       QDomNode n = e.firstChild();
0200       while (!n.isNull()) {
0201         QDomElement e = n.toElement();
0202         if (!e.isNull()) {
0203           if (e.tagName() == "dim") {
0204             if (e.hasAttribute("x")) {
0205               _nX = e.attribute("x").toInt();
0206             }
0207             if (e.hasAttribute("y")) {
0208               _nX = e.attribute("y").toInt();
0209             }
0210           } else if (e.tagName() == "theta") {
0211             if (e.hasAttribute("auto")) {
0212               _autoTheta = e.attribute("auto").toInt();
0213             }
0214             if (e.hasAttribute("units")) {
0215               _thetaUnits = e.attribute("units").toInt();
0216             }
0217             if (e.hasAttribute("min")) {
0218               confThetaMin = e.attribute("min").toDouble();
0219             }
0220             if (e.hasAttribute("max")) {
0221               confThetaMax = e.attribute("max").toDouble();
0222             }
0223           } else if (e.tagName() == "phi") {
0224             if (e.hasAttribute("auto")) {
0225               _autoPhi = e.attribute("auto").toInt();
0226             }
0227             if (e.hasAttribute("units")) {
0228               _phiUnits = e.attribute("units").toInt();
0229             }
0230             if (e.hasAttribute("min")) {
0231               confPhiMin = e.attribute("min").toDouble();
0232             }
0233             if (e.hasAttribute("max")) {
0234               confPhiMax = e.attribute("max").toDouble();
0235             }
0236           } else if (e.tagName() == "vector") {
0237             if (e.hasAttribute("auto")) {
0238               _autoMag = e.attribute("auto").toInt();
0239             }
0240             if (e.hasAttribute("degrade")) {
0241               tempdegrade = e.attribute("degrade").toInt();
0242               checkDegrade(tempdegrade);
0243               _vecDegrade = tempdegrade;
0244             }
0245             if (e.hasAttribute("theta")) {
0246               _vecTheta = e.attribute("theta").toInt();
0247             }
0248             if (e.hasAttribute("phi")) {
0249               _vecPhi = e.attribute("phi").toInt();
0250             }
0251             if (e.hasAttribute("QU")) {
0252               _vecQU = e.attribute("QU").toInt();
0253             }
0254             if (e.hasAttribute("max")) {
0255               _maxMag = e.attribute("max").toDouble();
0256             }
0257           }
0258         }
0259         n = n.nextSibling();
0260       }
0261
0262       // convert the entered range values into radians and
0263       // force them to the correct range.
0264       theta2Internal(_thetaUnits, confThetaMin);
0265       theta2Internal(_thetaUnits, confThetaMax);
0266       phi2Internal(_phiUnits, confPhiMin);
0267       phi2Internal(_phiUnits, confPhiMax);
0268
0269       // swap theta min/max if coordinate system requires it
0270       if (confThetaMax < confThetaMin) {
0271         double temp = confThetaMax;
0272         confThetaMax = confThetaMin;
0273         confThetaMin = temp;
0274       }
0275
0276       _thetaMin = confThetaMin;
0277       _thetaMax = confThetaMax;
0278       _phiMin = confPhiMin;
0279       _phiMax = confPhiMax;
0280     }
0281
0282     // data range
0283     int _nX;
0284     int _nY;
0285     double _thetaMin;
0286     double _phiMin;
0287     double _thetaMax;
0288     double _phiMax;
0289     bool _autoTheta;
0290     bool _autoPhi;
0291
0292     //values are
0293     // 0 : Radians
0294     // 1 : Degrees
0295     // 2 : Degrees RA/DEC
0296     // 3 : Degrees Lat/Long
0297     int _thetaUnits;
0298     int _phiUnits;
0299     int _vecUnits;
0300
0301     // vector field
0302     int _vecDegrade;
0303     int _vecTheta;
0304     int _vecPhi;
0305     bool _autoMag;
0306     double _maxMag;
0307     bool _vecQU;
0308
0309     size_t _mapNside;
0310
0311     void checkDegrade(int& degrade) {
0312       int nside = _mapNside;
0313       if (!degrade) {
0314         return;
0315       }
0316       if ((nside == 1) && (degrade != 0)) {
0317         degrade = 0;
0318         return;
0319       }
0320       for (int i = 0; i < degrade; i++) {
0321         nside = (int)(nside/2);
0322         if (nside == 1) {
0323           degrade = i+1;
0324           return;
0325         }
0326       }
0327       return;
0328     }
0329
0330     void theta2Internal(int units, double& theta) {
0331       switch (units) {
0332         case HPUNIT_RAD:
0333           break;
0334         case HPUNIT_DEG:
0335           theta *= HEALPIX_PI/180.0;
0336           break;
0337         case HPUNIT_RADEC: case HPUNIT_LATLON:
0338           theta = (90.0 - theta) * HEALPIX_PI/180.0;
0339           break;
0340         default:
0341           break;
0342       }
0343       while (theta < 0.0) {
0344         theta += HEALPIX_PI;
0345       }
0346       while (theta > HEALPIX_PI) {
0347         theta -= HEALPIX_PI;
0348       }
0349       return;
0350     }
0351
0352     void theta2External(int units, double& theta) {
0353       switch (units) {
0354         case HPUNIT_RAD:
0355           break;
0356         case HPUNIT_DEG:
0357           theta *= 180.0/HEALPIX_PI;
0358           break;
0359         case HPUNIT_RADEC: case HPUNIT_LATLON:
0360           theta = 90.0 - theta*180.0/HEALPIX_PI;
0361           break;
0362         default:
0363           break;
0364       }
0365       return;
0366     }
0367
0368     void phi2Internal(int units, double& phi) {
0369       switch (units) {
0370         case HPUNIT_RAD:
0371           break;
0372         case HPUNIT_DEG: case HPUNIT_RADEC: case HPUNIT_LATLON:
0373           phi *= HEALPIX_PI/180.0;
0374           break;
0375         default:
0376           break;
0377       }
0378       while (phi < 0.0) {
0379         phi += 2.0*HEALPIX_PI;
0380       }
0381       while (phi > 2.0*HEALPIX_PI) {
0382         phi -= 2.0*HEALPIX_PI;
0383       }
0384       return;
0385     }
0386
0387     void phi2External(int units, double& phi) {
0388       switch (units) {
0389         case HPUNIT_RAD:
0390           break;
0391         case HPUNIT_DEG: case HPUNIT_RADEC:
0392           phi *= 180.0/HEALPIX_PI;
0393           break;
0394         case HPUNIT_LATLON:
0395           phi *= 180.0/HEALPIX_PI;
0396           if (phi >= 180.0) {
0397             phi -= 360.0;
0398           }
0399           break;
0400         default:
0401           break;
0402       }
0403       return;
0404     }
0405 };
0406
0407
0408 HealpixSource::HealpixSource(Kst::ObjectStore *store, QSettings *cfg, const QString& filename, const QString& type, const QDomElement& e)
0409 : Kst::DataSource(store, cfg, filename, type), _config(0L) {
0410   _valid = false;
0411
0412   if (!type.isEmpty() && type != "HEALPIX Source") {
0413     return;
0414   }
0415
0416   _config = new HealpixSource::Config;
0417   _config->read(cfg, filename);
0418   if (!e.isNull()) {
0419     _config->load(e);
0420   }
0421
0422   if (init()) {
0423     _valid = true;
0424   }
0425
0426   update();
0427 }
0428
0429
0430 HealpixSource::~HealpixSource() {
0431   if (_keys) {
0432     healpix_keys_free(_keys);
0433   }
0434   if (_names) {
0435     healpix_strarr_free(_names, HEALPIX_FITS_MAXCOL);
0436   }
0437   if (_units) {
0438     healpix_strarr_free(_units, HEALPIX_FITS_MAXCOL);
0439   }
0440 }
0441
0442
0443 void HealpixSource::reset() {
0444   return;
0445 }
0446
0447
0448 bool HealpixSource::init() {
0449   size_t poff;
0450
0451   strncpy(_healpixfile, _filename.toLatin1().data(), HEALPIX_STRNL);
0452
0453   _names = healpix_strarr_alloc(HEALPIX_FITS_MAXCOL);
0454   _units = healpix_strarr_alloc(HEALPIX_FITS_MAXCOL);
0455   _keys = healpix_keys_alloc();
0456
0457   if (!healpix_fits_map_info(_healpixfile, &_mapNside, &_mapOrder, &_mapCoord, &_mapType, &_nMaps, _creator, _extname, _names, _units, _keys)) {
0458     healpix_keys_free(_keys);
0459     healpix_strarr_free(_names, HEALPIX_FITS_MAXCOL);
0460     healpix_strarr_free(_units, HEALPIX_FITS_MAXCOL);
0461
0462     return false;
0463   }
0464
0465   _config->_mapNside = _mapNside;
0466   if (_mapType == HEALPIX_FITS_CUT) {
0467     poff = 1;
0468   } else {
0469     poff = 0;
0470   }
0471
0472   _mapNpix = 12 * _mapNside * _mapNside;
0473
0474   // populate the metadata
0475   QString metaVal;
0476   QString metaName;
0477
0478   metaName = "NSIDE";
0479   metaVal.sprintf("%lu", (long unsigned int)_mapNside);
0480   _metaData.insert(metaName, metaVal);
0481
0482   metaName = "NPIX";
0483   metaVal.sprintf("%lu", (long unsigned int)_mapNpix);
0484   _metaData.insert(metaName, metaVal);
0485
0486   metaName = "ORDER";
0487   metaVal.sprintf("%d", _mapOrder);
0488   _metaData.insert(metaName, metaVal);
0489
0490   metaName = "COORD";
0491   metaVal.sprintf("%d", _mapCoord);
0492   _metaData.insert(metaName, metaVal);
0493
0494   metaName = "TYPE";
0495   metaVal.sprintf("%d", _mapType);
0496   _metaData.insert(metaName, metaVal);
0497
0498   metaName = "NMAPS";
0499   metaVal.sprintf("%lu", (long unsigned int)_nMaps);
0500   _metaData.insert(metaName, metaVal);
0501
0502   metaName = "CREATOR";
0503   metaVal.sprintf("%s", _creator);
0504   _metaData.insert(metaName, metaVal);
0505
0506   metaName = "EXTNAME";
0507   metaVal.sprintf("%s", _extname);
0508   _metaData.insert(metaName, metaVal);
0509
0510   for (size_t j = 0; j < _keys->nskeys; j++) {
0511     metaName.sprintf("%s", _keys->skeynames[j]);
0512     metaVal.sprintf("%s", _keys->skeyvals[j]);
0513     _metaData.insert(metaName, metaVal);
0514   }
0515
0516   for (size_t j = 0; j < _keys->nikeys; j++) {
0517     metaName.sprintf("%s", _keys->ikeynames[j]);
0518     metaVal.sprintf("%d", _keys->ikeyvals[j]);
0519     _metaData.insert(metaName, metaVal);
0520   }
0521
0522   for (size_t j = 0; j < _keys->nfkeys; j++) {
0523     metaName.sprintf("%s", _keys->fkeynames[j]);
0524     metaVal.sprintf("%e", _keys->fkeyvals[j]);
0525     _metaData.insert(metaName, metaVal);
0526   }
0527
0528   // populate the field list
0529   QString mapName;
0530   for (size_t i = 0; i < _nMaps; i++) {
0531     if (strlen(_names[i+poff]) == 0) {
0532       mapName.sprintf("%d - %s",(int)(i+1),"MAP");
0533     } else {
0534       mapName.sprintf("%d - %s",(int)(i+1),_names[i+poff]);
0535     }
0536     if (strlen(_units[i+poff]) == 0) {
0537       mapName.sprintf("%s (%s)",mapName.toLatin1().data(),"Unknown Units");
0538     } else {
0539       mapName.sprintf("%s (%s)",mapName.toLatin1().data(),_units[i+poff]);
0540     }
0541     _matrixList.append(mapName);
0542   }
0543
0544   if (_mapType == HEALPIX_FITS_CUT) {
0545     if (strlen(_names[_nMaps+1]) == 0) {
0546       mapName.sprintf("%s","HITS");
0547     } else {
0548       mapName.sprintf("%s",_names[_nMaps+1]);
0549     }
0550     _matrixList.append(mapName);
0551     if (strlen(_names[_nMaps+2]) == 0) {
0552       mapName.sprintf("%s","ERRORS");
0553     } else {
0554       mapName.sprintf("%s",_names[_nMaps+2]);
0555     }
0556     if (strlen(_units[_nMaps+2]) == 0) {
0557       mapName.sprintf("%s (%s)",mapName.toLatin1().data(),"Unknown Units");
0558     } else {
0559       mapName.sprintf("%s (%s)",mapName.toLatin1().data(),_units[_nMaps+2]);
0560     }
0561     _matrixList.append(mapName);
0562   }
0563
0564   _fieldList.append("1 - Vector Field Head Theta");
0565   _fieldList.append("2 - Vector Field Head Phi");
0566   _fieldList.append("3 - Vector Field Tail Theta");
0567   _fieldList.append("4 - Vector Field Tail Phi");
0568
0569   return true;
0570 }
0571
0572
0573 Kst::Object::UpdateType HealpixSource::update() {
0574   // updates not supported yet
0575   // we should check to see if the FITS file has changed on disk
0576   return Kst::Object::NO_CHANGE;
0577 }
0578
0579
0580 bool HealpixSource::matrixDimensions( const QString& matrix, int* xDim, int* yDim) {
0581   Q_UNUSED(matrix)
0582   if (_valid) {
0583     (*xDim) = _config->_nX;
0584     (*yDim) = _config->_nY;
0585     return true;
0586   }
0587   return false;
0588 }
0589
0590
0591 int HealpixSource::readMatrix(Kst::MatrixData* data, const QString& field, int xStart,
0592                                      int yStart, int xNumSteps,
0593                                      int yNumSteps) {
0594   // If the file is _valid, then we already
0595   // have all the header information- no need to read it again.
0596   // We also know that the matrix index is not out-of-range.
0597   if (_valid && isValidMatrix(field)) {
0598     fitsfile *fp;
0599     int colnum;
0600     int fieldnum;
0601     int ret = 0;
0602     int ncol;
0603     int ttype;
0604     long nrows;
0605     long pcount;
0606     int tfields;
0607     char extname[HEALPIX_STRNL];
0608     char comment[HEALPIX_STRNL];
0609     float *datavec;
0610     int *pixvec;
0611     float *mapdata;
0612     float nullval = 0.0;
0613     int nnull = 0;
0614     int keynpix;
0615     int keyfirst;
0616     int ischunk;
0617     long nelem;
0618
0619     if (_matrixList.contains(field)) {
0620       fieldnum = _matrixList.findIndex(field);
0621     } else {
0622       fieldnum = field.toInt();
0623       fieldnum--;
0624     }
0625
0626     // check range
0627     if ((xStart >= _config->_nX) || (yStart >= _config->_nY)) {
0628       return -1;
0629     }
0630     int nxread = xNumSteps;
0631     if (nxread < 0) {
0632       nxread = 1;
0633     }
0634     if (xStart + nxread > _config->_nX) {
0635       nxread = _config->_nX - xStart;
0636     }
0637     int nyread = yNumSteps;
0638     if (nyread < 0) {
0639       nyread = 1;
0640     }
0641     if (yStart + nyread > _config->_nY) {
0642       nyread = _config->_nY - yStart;
0643     }
0644
0645     if (_mapType == HEALPIX_FITS_CUT) {
0646       // cut-sphere files have extra pixel/hits/error columns
0647       colnum = fieldnum + 2;
0648       ncol = (int)_nMaps + 3;
0649     } else {
0650       colnum = fieldnum + 1;
0651       ncol = (int)_nMaps;
0652     }
0653
0654     // open file and move to second header unit
0655     if (fits_open_file(&fp, _healpixfile, READONLY, &ret)) {
0656       return -1;
0657     }
0658
0659     if (fits_movabs_hdu(fp, 2, &ttype, &ret)) {
0660       return -1;
0661     }
0662
0663     // read the number of rows
0664     if (fits_read_btblhdr(fp, ncol, &nrows, &tfields, NULL, NULL, NULL, extname, &pcount, &ret)) {
0665       ret = 0;
0666       fits_close_file(fp, &ret);
0667       return -1;
0668     }
0669
0670     mapdata = (float*)calloc(_mapNpix, sizeof(float));
0671     //initialize data to HEALPIX_NULL
0672     for (size_t j = 0; j < _mapNpix; j++) {
0673       mapdata[j] = HEALPIX_NULL;
0674     }
0675
0676     if (_mapType == HEALPIX_FITS_CUT) {
0677       // For a cut-sphere file, we must read the entire
0678       // file and then re-map the data onto a full-sphere
0679       // vector.
0680
0681       datavec = (float*)calloc((size_t) nrows, sizeof(float));
0682       pixvec = (int*)calloc((size_t) nrows, sizeof(int));
0683
0684       if (fits_read_col(fp, TINT, 1, 1, 1, nrows, &nullval, pixvec, &nnull, &ret)) {
0685         free(pixvec);
0686         free(datavec);
0687         free(mapdata);
0688         ret = 0;
0689         fits_close_file(fp, &ret);
0690         return -1;
0691       }
0692
0693       if (fits_read_col(fp, TFLOAT, colnum, 1, 1, nrows, &nullval, datavec, &nnull, &ret)) {
0694         free(pixvec);
0695         free(datavec);
0696         free(mapdata);
0697         ret = 0;
0698         fits_close_file(fp, &ret);
0699         return -1;
0700       }
0701
0702       for (long j = 0; j < nrows; j++) {
0703         if ((pixvec[j] >= 0) && (pixvec[j] < (int)_mapNpix)) {
0704           mapdata[pixvec[j]] = datavec[j];
0705         }
0706       }
0707       free(pixvec);
0708       free(datavec);
0709     } else {
0710       /* is this a chunk? */
0711       if ((nrows != (long)(_mapNpix))&&(1024*nrows != (long)(_mapNpix))) {
0712         /*this must be a chunk file*/
0713         char charFirstPix[] = "FIRSTPIX";
0714         if (fits_read_key(fp, TLONG, charFirstPix, &keyfirst, comment, &ret)) {
0715           /*must at least have FIRSTPIX key*/
0716           fits_close_file(fp, &ret);
0717           return 0;
0718         } else {
0719           char charNPix[] = "NPIX";
0720           if (fits_read_key(fp, TLONG, charNPix, &keynpix, comment, &ret)) {
0721             ret = 0;
0722             /*might be using LASTPIX instead*/
0723             char charLastPix[] = "LASTPIX";
0724             if (fits_read_key(fp, TLONG, charLastPix, &keynpix, comment, &ret)) {
0725               fits_close_file(fp, &ret);
0726               return 0;
0727             } else {
0728               keynpix = keynpix - keyfirst + 1;
0729               ischunk = 1;
0730             }
0731           } else {
0732             ischunk = 1;
0733           }
0734         }
0735       } else {
0736         ischunk = 0;
0737       }
0738       if (ischunk) {
0739         datavec = (float*)calloc((size_t)keynpix, sizeof(float));
0740         nelem = (long)keynpix;
0741       } else {
0742         datavec = (float*)calloc(_mapNpix, sizeof(float));
0743         nelem = (long)(_mapNpix);
0744       }
0745       if (fits_read_col(fp, TFLOAT, colnum, 1, 1, nelem, &nullval, datavec, &nnull, &ret)) {
0746         free(datavec);
0747         free(mapdata);
0748         ret = 0;
0749         fits_close_file(fp, &ret);
0750         return -1;
0751       }
0752       if (ischunk) {
0753         for (long j = 0; j < nelem; j++) {
0754           mapdata[j+keyfirst] = datavec[j];
0755         }
0756       } else {
0757         for (long j = 0; j < nelem; j++) {
0758           mapdata[j] = datavec[j];
0759         }
0760       }
0761       free(datavec);
0762     }
0763
0764     fits_close_file(fp, &ret);
0765
0766     // set the matrix to NULL
0767     for (int i = 0; i < nxread*nyread; i++) {
0768       data->z[i] = NAN;
0769     }
0770
0771     double theta, phi;
0772     size_t ppix;
0773
0774     // compute autorange parameters if necessary
0775     if (_config->_autoTheta || _config->_autoPhi) {
0776       double mapMinTheta = HEALPIX_PI;
0777       double mapMaxTheta = 0.0;
0778       double mapMinPhi = 2.0*HEALPIX_PI;
0779       double mapMaxPhi = 0.0;
0780
0781       for (size_t i = 0; i < _mapNpix; i++) {
0782         if (!healpix_is_fnull(mapdata[i])) {
0783           if (_mapOrder == HEALPIX_RING) {
0784             healpix_pix2ang_ring(_mapNside, i, &theta, &phi);
0785           } else {
0786             healpix_pix2ang_nest(_mapNside, i, &theta, &phi);
0787           }
0788           if (theta < mapMinTheta) {
0789             mapMinTheta = theta;
0790           }
0791           if (theta > mapMaxTheta) {
0792             mapMaxTheta = theta;
0793           }
0794           if (phi < mapMinPhi) {
0795             mapMinPhi = phi;
0796           }
0797           if (phi > mapMaxPhi) {
0798             mapMaxPhi = phi;
0799           }
0800         }
0801       }
0802       if (mapMaxTheta < mapMinTheta) { // no valid data in map
0803         mapMaxTheta = HEALPIX_PI;
0804         mapMinTheta = 0.0;
0805         mapMaxPhi = 2.0*HEALPIX_PI;
0806         mapMinPhi = 0.0;
0807       }
0808       if (_config->_autoTheta) {
0809         _config->_thetaMin = mapMinTheta;
0810         _config->_thetaMax = mapMaxTheta;
0811       }
0812       if (_config->_autoPhi) {
0813         _config->_phiMin = mapMinPhi;
0814         _config->_phiMax = mapMaxPhi;
0815       }
0816       //qDebug() << "HEALPIX autorange is Theta=[" << mapMinTheta << "..." << mapMaxTheta << "] Phi=[" << mapMinPhi << "..." << mapMaxPhi << "]";
0817     }
0818     //qDebug() << "HEALPIX using range Theta=[" << _thetaMin << "..." << _thetaMax << "] Phi=[" << _phiMin << "..." << _phiMax << "]";
0819
0820     // copy sphere data to matrix.
0821     for (int i = xStart; i < nxread; i++) {
0822       for (int j = yStart; j < nyread; j++) {
0823         healpix_proj_rev_car(_config->_thetaMin, _config->_thetaMax, _config->_phiMin, _config->_phiMax, (double)_config->_nX, (double)_config->_nY, (double)i, (double)j, &theta, &phi);
0824         if ((!healpix_is_dnull(theta)) && (!healpix_is_dnull(phi))) {
0825           if (_mapOrder == HEALPIX_RING) {
0826             healpix_ang2pix_ring(_mapNside, theta, phi, &ppix);
0827           } else {
0828             healpix_ang2pix_nest(_mapNside, theta, phi, &ppix);
0829           }
0830           if (!healpix_is_fnull(mapdata[ppix])) {
0831             data->z[i*nyread+j] = (double)mapdata[ppix];
0832           }
0833         }
0834       }
0835     }
0836     free(mapdata);
0837
0838     // FIXME
0839     // Eventually, we can just always use radians for the
0840     // matrix range, since actual display units will be
0841     // handled by the Kst2DRenderer.  For now, set the matrix
0842     // range to the same units used in the config dialog.
0843     // Note that this behaviour is broken, since 2 maps
0844     // from different files could be put on the same plot
0845     // in different units.  It's the best we can do for now.
0846     //
0847     //data->yMin = _thetaMax;
0848     //data->yStepSize = -(_thetaMax - _thetaMin)/(double)_nY;
0849     //data->xMin = _phiMin;
0850     //if (_phiMin > _phiMax) {
0851     //  data->xStepSize = (_phiMax + (2.0*HEALPIX_PI - _phiMin)) / (double)_nX;
0852     //} else {
0853     //  data->xStepSize = ((_phiMax - _phiMin) / (double)_nX;
0854     //}
0855
0856     double tMin = _config->_thetaMin;
0857     double tMax = _config->_thetaMax;
0858     double pMin = _config->_phiMin;
0859     double pMax = _config->_phiMax;
0860     _config->theta2External(_config->_thetaUnits, tMin);
0861     _config->theta2External(_config->_thetaUnits, tMax);
0862     _config->phi2External(_config->_phiUnits, pMin);
0863     _config->phi2External(_config->_phiUnits, pMax);
0864     //qDebug() << "HEALPIX exported range Theta=[" << tMin << "..." << tMax << "] Phi=[" << pMin << "..." << pMax << "]";
0865
0866     switch (_config->_thetaUnits) {
0867       case HPUNIT_RAD: case HPUNIT_DEG:
0868         data->yMin = tMax;
0869         data->yStepSize = (tMin - tMax)/(double)_config->_nY;
0870         break;
0871       case HPUNIT_RADEC: case HPUNIT_LATLON:
0872         data->yMin = tMax;
0873         data->yStepSize = (tMin - tMax)/(double)_config->_nY;
0874         break;
0875       default:
0876         break;
0877     }
0878
0879     data->xMin = pMin;
0880     switch (_config->_phiUnits) {
0881       case HPUNIT_RAD:
0882         if (pMin > pMax) {
0883           data->xStepSize = (pMax + (2.0*HEALPIX_PI-pMin)) / (double)_config->_nX;
0884         } else {
0885           data->xStepSize = (pMax - pMin) / (double)_config->_nX;
0886         }
0887         break;
0888       case HPUNIT_DEG: case HPUNIT_RADEC:
0889         if (pMin > pMax) {
0890           data->xStepSize = (pMax + (360.0-pMin)) / (double)_config->_nX;
0891         } else {
0892           data->xStepSize = (pMax - pMin) / (double)_config->_nX;
0893         }
0894         break;
0895       case HPUNIT_LATLON:
0896         if (pMin > pMax) {
0897           data->xStepSize = ((180.0+pMax) + (180.0-pMin)) / (double)_config->_nX;
0898         } else {
0899           data->xStepSize = (pMax - pMin) / (double)_config->_nX;
0900         }
0901         break;
0902       default:
0903         break;
0904     }
0905     return nxread*nyread;
0906   } else {
0907     return -1;
0908   }
0909 }
0910
0911
0912 int HealpixSource::readField(double *v, const QString& field, int s, int n) {
0913   Q_UNUSED(s)
0914   if (_valid && isValidField(field)) {
0915     fitsfile *fp;
0916     int coltheta;
0917     int colphi;
0918     int ret = 0;
0919     int ncol;
0920     int ttype;
0921     long nrows;
0922     long pcount;
0923     int tfields;
0924     char extname[HEALPIX_STRNL];
0925     char comment[HEALPIX_STRNL];
0926     float *datavec;
0927     int *pixvec;
0928     int *hits;
0929     float *mapdata;
0930     float *comptheta;
0931     float *compphi;
0932     float nullval;
0933     int nnull;
0934     size_t vecNside;
0935     size_t temppix;
0936     int fieldnum;
0937     double theta, phi;
0938     long nearest[8];
0939     int keynpix;
0940     int keyfirst;
0941     int ischunk;
0942     long nelem;
0943
0944     if (_fieldList.contains(field)) {
0945       fieldnum = _fieldList.findIndex(field);
0946     } else {
0947       fieldnum = field.toInt();
0948       fieldnum--;
0949     }
0950
0951     // check range
0952     if (n <= 0) {
0953       return -1;
0954     }
0955
0956     vecNside = _mapNside;
0957     for (int i = 0; i < _config->_vecDegrade; i++) {
0958       vecNside = (size_t)(vecNside/2);
0959     }
0960
0961     if (_mapType == HEALPIX_FITS_CUT) {
0962       // cut-sphere files have extra pixel/hits/error columns
0963       coltheta = _config->_vecTheta + 2;
0964       colphi = _config->_vecPhi + 2;
0965       ncol = (int)_nMaps + 3;
0966     } else {
0967       coltheta = _config->_vecTheta + 1;
0968       colphi = _config->_vecPhi + 1;
0969       ncol = (int)_nMaps;
0970     }
0971
0972     // open file and move to second header unit
0973     if (fits_open_file(&fp, _healpixfile, READONLY, &ret)) {
0974       return -1;
0975     }
0976
0977     if (fits_movabs_hdu(fp, 2, &ttype, &ret)) {
0978       return -1;
0979     }
0980
0981     // read the number of rows
0982     if (fits_read_btblhdr(fp, ncol, &nrows, &tfields, NULL, NULL, NULL, extname, &pcount, &ret)) {
0983       ret = 0;
0984       fits_close_file(fp, &ret);
0985       return -1;
0986     }
0987
0988     mapdata = (float*)calloc(_mapNpix, sizeof(float));
0989     hits = (int*)calloc(12*vecNside*vecNside, sizeof(int));
0990
0991     //-----------------------------------------
0992     // Read in the first component and degrade
0993     //-----------------------------------------
0994     for (size_t j = 0; j < _mapNpix; j++) {
0995       mapdata[j] = HEALPIX_NULL;
0996     }
0997     if (_mapType == HEALPIX_FITS_CUT) {
0998       datavec = (float*)calloc((size_t) nrows, sizeof(float));
0999       pixvec = (int*)calloc((size_t) nrows, sizeof(int));
1000       if (fits_read_col(fp, TINT, 1, 1, 1, nrows, &nullval, pixvec, &nnull, &ret)) {
1001         free(pixvec);
1002         free(datavec);
1003         free(mapdata);
1004         ret = 0;
1005         fits_close_file(fp, &ret);
1006         return -1;
1007       }
1008       if (fits_read_col(fp, TFLOAT, coltheta, 1, 1, nrows, &nullval, datavec, &nnull, &ret)) {
1009         free(pixvec);
1010         free(datavec);
1011         free(mapdata);
1012         ret = 0;
1013         fits_close_file(fp, &ret);
1014         return -1;
1015       }
1016       for (long j = 0; j < nrows; j++) {
1017         if ((pixvec[j] >= 0) && (pixvec[j] < (int)_mapNpix)) {
1018           mapdata[pixvec[j]] = datavec[j];
1019         }
1020       }
1021       free(pixvec);
1022       free(datavec);
1023     } else {
1024       /* is this a chunk? */
1025       if ((nrows != (long)(_mapNpix))&&(1024*nrows != (long)(_mapNpix))) {
1026         /*this must be a chunk file*/
1027         char charFirstPix[] = "FIRSTPIX";
1028         if (fits_read_key(fp, TLONG, charFirstPix, &keyfirst, comment, &ret)) {
1029           /*must at least have FIRSTPIX key*/
1030           fits_close_file(fp, &ret);
1031           return 0;
1032         } else {
1033           char charNPix[] = "NPIX";
1034           if (fits_read_key(fp, TLONG, charNPix, &keynpix, comment, &ret)) {
1035             ret = 0;
1036             /*might be using LASTPIX instead*/
1037             char charLastPix[] = "LASTPIX";
1038             if (fits_read_key(fp, TLONG, charLastPix, &keynpix, comment, &ret)) {
1039               fits_close_file(fp, &ret);
1040               return 0;
1041             } else {
1042               keynpix = keynpix - keyfirst + 1;
1043               ischunk = 1;
1044             }
1045           } else {
1046             ischunk = 1;
1047           }
1048         }
1049       } else {
1050         ischunk = 0;
1051       }
1052       if (ischunk) {
1053         datavec = (float*)calloc((size_t)keynpix, sizeof(float));
1054         nelem = (long)keynpix;
1055       } else {
1056         datavec = (float*)calloc(_mapNpix, sizeof(float));
1057         nelem = (long)(_mapNpix);
1058       }
1059       if (fits_read_col(fp, TFLOAT, coltheta, 1, 1, nelem, &nullval, datavec, &nnull, &ret)) {
1060         free(datavec);
1061         free(mapdata);
1062         ret = 0;
1063         fits_close_file(fp, &ret);
1064         return -1;
1065       }
1066       if (ischunk) {
1067         for (long j = 0; j < nelem; j++) {
1068           mapdata[j+keyfirst] = datavec[j];
1069         }
1070       } else {
1071         for (long j = 0; j < nelem; j++) {
1072           mapdata[j] = datavec[j];
1073         }
1074       }
1075       free(datavec);
1076     }
1077     comptheta = (float*)calloc(12*vecNside*vecNside, sizeof(float));
1078     for (size_t i = 0; i < 12*vecNside*vecNside; i++) {
1079       comptheta[i] = HEALPIX_NULL;
1080       hits[i] = 0;
1081     }
1082     for (size_t i = 0; i < _mapNpix; i++) {
1083       if (!healpix_is_fnull(mapdata[i])) {
1084         if (_mapOrder == HEALPIX_NEST) {
1085           healpix_degrade_nest(_mapNside, i, vecNside, &temppix);
1086         } else {
1087           healpix_degrade_ring(_mapNside, i, vecNside, &temppix);
1088         }
1089         if (!healpix_is_fnull(comptheta[temppix])) {
1090           comptheta[temppix] += mapdata[i];
1091         } else {
1092           comptheta[temppix] = mapdata[i];
1093         }
1094         hits[temppix] += 1;
1095       }
1096     }
1097     for (size_t i = 0; i < 12*vecNside*vecNside; i++) {
1098       if (!healpix_is_fnull(comptheta[i])) {
1099         comptheta[i] = comptheta[i] / (double)hits[i];
1100       }
1101     }
1102
1103     //-----------------------------------------
1104     // Read in the second component and degrade
1105     //-----------------------------------------
1106     for (size_t j = 0; j < _mapNpix; j++) {
1107       mapdata[j] = HEALPIX_NULL;
1108     }
1109     if (_mapType == HEALPIX_FITS_CUT) {
1110       datavec = (float*)calloc((size_t) nrows, sizeof(float));
1111       pixvec = (int*)calloc((size_t) nrows, sizeof(int));
1112       if (fits_read_col(fp, TINT, 1, 1, 1, nrows, &nullval, pixvec, &nnull, &ret)) {
1113         free(pixvec);
1114         free(datavec);
1115         free(mapdata);
1116         ret = 0;
1117         fits_close_file(fp, &ret);
1118         return -1;
1119       }
1120       if (fits_read_col(fp, TFLOAT, colphi, 1, 1, nrows, &nullval, datavec, &nnull, &ret)) {
1121         free(pixvec);
1122         free(datavec);
1123         free(mapdata);
1124         ret = 0;
1125         fits_close_file(fp, &ret);
1126         return -1;
1127       }
1128       for (long j = 0; j < nrows; j++) {
1129         if ((pixvec[j] >= 0) && (pixvec[j] < (int)_mapNpix)) {
1130           mapdata[pixvec[j]] = datavec[j];
1131         }
1132       }
1133       free(pixvec);
1134       free(datavec);
1135     } else {
1136       /* is this a chunk? */
1137       if ((nrows != (long)(_mapNpix))&&(1024*nrows != (long)(_mapNpix))) {
1138         /*this must be a chunk file*/
1139         char charFirstPix[] = "FIRSTPIX";
1140         if (fits_read_key(fp, TLONG, charFirstPix, &keyfirst, comment, &ret)) {
1141           /*must at least have FIRSTPIX key*/
1142           fits_close_file(fp, &ret);
1143           return 0;
1144         } else {
1145           char charNPix[] = "NPIX";
1146           if (fits_read_key(fp, TLONG, charNPix, &keynpix, comment, &ret)) {
1147             ret = 0;
1148             /*might be using LASTPIX instead*/
1149             char charLastPix[] = "LASTPIX";
1150             if (fits_read_key(fp, TLONG, charLastPix, &keynpix, comment, &ret)) {
1151               fits_close_file(fp, &ret);
1152               return 0;
1153             } else {
1154               keynpix = keynpix - keyfirst + 1;
1155               ischunk = 1;
1156             }
1157           } else {
1158             ischunk = 1;
1159           }
1160         }
1161       } else {
1162         ischunk = 0;
1163       }
1164       if (ischunk) {
1165         datavec = (float*)calloc((size_t)keynpix, sizeof(float));
1166         nelem = (long)keynpix;
1167       } else {
1168         datavec = (float*)calloc(_mapNpix, sizeof(float));
1169         nelem = (long)(_mapNpix);
1170       }
1171       if (fits_read_col(fp, TFLOAT, colphi, 1, 1, nelem, &nullval, datavec, &nnull, &ret)) {
1172         free(datavec);
1173         free(mapdata);
1174         ret = 0;
1175         fits_close_file(fp, &ret);
1176         return -1;
1177       }
1178       if (ischunk) {
1179         for (long j = 0; j < nelem; j++) {
1180           mapdata[j+keyfirst] = datavec[j];
1181         }
1182       } else {
1183         for (long j = 0; j < nelem; j++) {
1184           mapdata[j] = datavec[j];
1185         }
1186       }
1187       free(datavec);
1188     }
1189
1190     compphi = (float*)calloc(12*vecNside*vecNside, sizeof(float));
1191     for (size_t i = 0; i < 12*vecNside*vecNside; i++) {
1192       compphi[i] = HEALPIX_NULL;
1193       hits[i] = 0;
1194     }
1195     for (size_t i = 0; i < _mapNpix; i++) {
1196       if (!healpix_is_fnull(mapdata[i])) {
1197         if (_mapOrder == HEALPIX_NEST) {
1198           healpix_degrade_nest(_mapNside, i, vecNside, &temppix);
1199         } else {
1200           healpix_degrade_ring(_mapNside, i, vecNside, &temppix);
1201         }
1202         if (!healpix_is_fnull(compphi[temppix])) {
1203           compphi[temppix] += mapdata[i];
1204         } else {
1205           compphi[temppix] = mapdata[i];
1206         }
1207         hits[temppix] += 1;
1208       }
1209     }
1210     for (size_t i = 0; i < 12*vecNside*vecNside; i++) {
1211       if (!healpix_is_fnull(compphi[i])) {
1212         compphi[i] = compphi[i] / (double)hits[i];
1213       }
1214     }
1215
1216     free(hits);
1217     free(mapdata);
1218     fits_close_file(fp, &ret);
1219
1220     //----------------------------------------------
1221     // compute head and tail coordinates in radians
1222     //----------------------------------------------
1223
1224     // autoscale magnitude if necessary
1225     double vecMag = 0.0;
1226     if (_config->_autoMag) {
1227       for (size_t i = 0; i < 12*vecNside*vecNside; i++) {
1228         if (sqrt((double)(comptheta[i]*comptheta[i] + compphi[i]*compphi[i])) > vecMag) {
1229           vecMag = sqrt((double)(comptheta[i]*comptheta[i] + compphi[i]*compphi[i]));
1230         }
1231       }
1232     } else {
1233       vecMag = _config->_maxMag;
1234     }
1235
1236     // find distance to 4 nearest neighbor pixels and
1237     // use the average to compute how "long" the max
1238     // valued vector should be for a given pixel.
1239     // This is *expensive*, but we only load vector
1240     // fields infrequently, so should be ok.
1241
1242     // only actually compute the component of the
1243     // vectorfield that is requested.
1244
1245     double maxang;
1246     double thetapart;
1247     double phipart;
1248     double headtheta;
1249     double tailtheta;
1250     double headphi;
1251     double tailphi;
1252     double pmag;
1253     double alpha;
1254     // only compute as many elements as requested
1255     for (size_t i = 0; i < (size_t)n; i++) {
1256       if (!healpix_is_fnull(comptheta[i]) && !healpix_is_fnull(compphi[i])) {
1257         if (_config->_vecQU) { //theta is really Q, phi really U
1258           pmag = sqrt((double)(comptheta[i]*comptheta[i] + compphi[i]*compphi[i]));
1259           alpha = 0.5 * atan((double)compphi[i] / (double)comptheta[i]);
1260           thetapart = pmag * cos(alpha);
1261           phipart = pmag * sin(alpha);
1262         } else {
1263           thetapart = (double)comptheta[i];
1264           phipart = (double)compphi[i];
1265         }
1266         if (_mapOrder == HEALPIX_RING) {
1267           healpix_pix2ang_ring(vecNside, i, &theta, &phi);
1268         } else {
1269           healpix_pix2ang_nest(vecNside, i, &theta, &phi);
1270         }
1271         //find max angle for this pixel
1272         healpix_neighbors(vecNside, _mapOrder, i, nearest);
1273         maxang = healpix_loc_dist(vecNside, _mapOrder, i, nearest[0]);
1274         maxang += healpix_loc_dist(vecNside, _mapOrder, i, nearest[1]);
1275         maxang += healpix_loc_dist(vecNside, _mapOrder, i, nearest[2]);
1276         maxang += healpix_loc_dist(vecNside, _mapOrder, i, nearest[3]);
1277         maxang /= 4.0;
1278         thetapart *= maxang / vecMag;
1279         phipart *= maxang / vecMag;
1280         //find coordinates of head and tail
1281         headtheta = theta + 0.5 * thetapart;
1282         tailtheta = theta - 0.5 * thetapart;
1283         headphi = phi + 0.5 * phipart;
1284         tailphi = phi - 0.5 * phipart;
1285         while (headtheta > HEALPIX_PI) {
1286           headtheta -= HEALPIX_PI;
1287         }
1288         while (headtheta < 0.0) {
1289           headtheta += HEALPIX_PI;
1290         }
1291         while (tailtheta > HEALPIX_PI) {
1292           tailtheta -= HEALPIX_PI;
1293         }
1294         while (tailtheta < 0.0) {
1295           tailtheta += HEALPIX_PI;
1296         }
1297         while (headphi > 2.0 * HEALPIX_PI) {
1298           headphi -= 2.0 * HEALPIX_PI;
1299         }
1300         while (headphi < 0.0) {
1301           headphi += 2.0 * HEALPIX_PI;
1302         }
1303         while (tailphi > 2.0 * HEALPIX_PI) {
1304           tailphi -= 2.0 * HEALPIX_PI;
1305         }
1306         while (tailphi < 0.0) {
1307           tailphi += 2.0 * HEALPIX_PI;
1308         }
1309         switch (fieldnum) {
1310           case 0: //head theta
1311             v[i] = headtheta;
1312             break;
1313           case 1: //head phi
1314             v[i] = headphi;
1315             break;
1316           case 2: //tail theta
1317             v[i] = tailtheta;
1318             break;
1319           case 3: //tail phi
1320             v[i] = tailphi;
1321             break;
1322           default:
1323             break;
1324         }
1325       }
1326     }
1327     return n;
1328   } else {
1329     return -1;
1330   }
1331 }
1332
1333
1334 bool HealpixSource::isValidField(const QString& field) const {
1335   // need to allow for referencing fields by number, so
1336   // that command line options work.
1337   if (_fieldList.contains(field)) {
1338     return true;
1339   } else {
1340     bool ok = false;
1341     int num = field.toInt(&ok);
1342     if (ok && num <= (int)_fieldList.count() && num != 0) {
1343       return true;
1344     } else {
1345       return false;
1346     }
1347   }
1348
1349   return _fieldList.contains(field);
1350 }
1351
1352
1353 bool HealpixSource::isValidMatrix(const QString& field) const {
1354   // need to allow for referencing fields by number, so
1355   // that command line options work.
1356   if (_matrixList.contains(field)) {
1357     return true;
1358   } else {
1359     bool ok = false;
1360     int num = field.toInt(&ok);
1361     if (ok && num <= (int)_matrixList.count() && num != 0) {
1362       return true;
1363     } else {
1364       return false;
1365     }
1366   }
1367 }
1368
1369
1370 int HealpixSource::samplesPerFrame(const QString &field) {
1371   Q_UNUSED(field)
1372   return 1;
1373 }
1374
1375
1376 int HealpixSource::frameCount(const QString& field) const {
1377   Q_UNUSED(field)
1378   if (_valid) {
1379     // return the number of pixels in the sky for
1380     // the degraded vectorfield
1381     size_t vecNside = _mapNside;
1382     for (int i = 0; i < _config->_vecDegrade; i++) {
1383       vecNside = (size_t)(vecNside/2);
1384     }
1385     return (int)(12*vecNside*vecNside);
1386   } else {
1387     return 0;
1388   }
1389 }
1390
1391
1392 bool HealpixSource::isEmpty() const {
1393   return !_valid;
1394 }
1395
1396
1397 QString HealpixSource::fileType() const {
1398   return healpixTypeString;
1399 }
1400
1401
1402 void HealpixSource::save(QXmlStreamWriter &streamWriter) {
1403   Kst::DataSource::save(streamWriter);
1404 }
1405
1406
1407 void HealpixSource::parseProperties(QXmlStreamAttributes &properties) {
1408   _config->parseProperties(properties);
1409 }
1410
1411
1412 int HealpixSource::readScalar(double &S, const QString& scalar) {
1413   if (scalar == "FRAMES") {
1414     S = 1;
1415     return 1;
1416   }
1417   return 0;
1418 }
1419
1420
1421 int HealpixSource::readString(QString &S, const QString& string) {
1422   if (string == "FILE") {
1423     S = _filename;
1424     return 1;
1425   } else if (_metaData.contains(string)) {
1426     S = _metaData[string];
1427     return 1;
1428   }
1429   return 0;
1430 }
1431
1432
1433 class ConfigWidgetHealpixInternal : public QWidget, public Ui_HealpixConfig {
1434   public:
1435     ConfigWidgetHealpixInternal(QWidget *parent) : QWidget(parent), Ui_HealpixConfig() { setupUi(this); }
1436 };
1437
1438
1439 class ConfigWidgetHealpix : public Kst::DataSourceConfigWidget {
1440   public:
1441     ConfigWidgetHealpix() : Kst::DataSourceConfigWidget() {
1442       QGridLayout *layout = new QGridLayout(this);
1443       _hc = new ConfigWidgetHealpixInternal(this);
1444       layout->addWidget(_hc, 0, 0);
1445       layout->activate();
1446     }
1447
1448     ~ConfigWidgetHealpix() {}
1449
1450     void setConfig(QSettings *cfg) {
1451       Kst::DataSourceConfigWidget::setConfig(cfg);
1452     }
1453
1454     void load() {
1455        QStringList unitList;
1456       //FIXME add back in when we have a renderer
1457       //unitList.append("(Radians)");
1458       //unitList.append("(Degrees)");
1459       unitList.append("(RA/DEC)");
1460       //unitList.append("(Lat/Long)");
1461       _cfg->beginGroup("Healpix General");
1462       _hc->matThetaUnits->clear();
1463       _hc->matPhiUnits->clear();
1464       _hc->vecTheta->clear();
1465       _hc->vecPhi->clear();
1466       _hc->matThetaUnits->insertStringList(unitList);
1467       _hc->matPhiUnits->insertStringList(unitList);
1468       //FIXME change to HPUNIT_RAD
1469       _hc->matThetaUnits->setCurrentItem(0);
1470       _hc->matPhiUnits->setCurrentItem(0);
1471       _hc->matDimX->setValue(_cfg->value("Matrix X Dimension", DEFAULT_XDIM).toInt());
1472       _hc->matDimY->setValue(_cfg->value("Matrix Y Dimension", DEFAULT_YDIM).toInt());
1473       _hc->matThetaAuto->setChecked(_cfg->value("Theta Autoscale", true).toBool());
1474       //FIXME Only one combobox item for now...
1475       //_hc->matThetaUnits->setCurrentItem(_cfg->readEntry("Theta Units", HPUNIT_RAD));
1476       _hc->matThetaUnits->setCurrentItem(0);
1477       _hc->matThetaMin->setText(_cfg->value("Theta Min").toString());
1478       _hc->matThetaMax->setText(_cfg->value("Theta Max").toString());
1479       _hc->matPhiAuto->setChecked(_cfg->value("Phi Autoscale", true).toBool());
1480       //FIXME Only one combobox item for now...
1481       //_hc->matPhiUnits->setCurrentItem(_cfg->readEntry("Phi Units", HPUNIT_RAD));
1482       _hc->matPhiUnits->setCurrentItem(0);
1483       _hc->matPhiMin->setText(_cfg->value("Phi Min").toString());
1484       _hc->matPhiMax->setText(_cfg->value("Phi Max").toString());
1485       _hc->vecTheta->setCurrentItem(_cfg->value("Vector Theta", 0).toInt());
1486       _hc->vecPhi->setCurrentItem(_cfg->value("Vector Phi", 0).toInt());
1487       _hc->vecDegrade->setValue(_cfg->value("Vector Degrade Factor", 1).toInt());
1488       _hc->vecMagAuto->setChecked(_cfg->value("Vector Magnitude Autoscale", true).toBool());
1489       _hc->vecMag->setText(_cfg->value("Vector Max Magnitude").toString());
1490       _hc->vecQU->setChecked(_cfg->value("Vector is QU", false).toBool());
1491
1492       bool hasInstance = (_instance != 0L);
1493
1494       if (hasInstance) {
1495         _hc->vecTheta->insertStringList(_instance->matrixList());
1496         _hc->vecPhi->insertStringList(_instance->matrixList());
1497         Kst::SharedPtr<HealpixSource> src = Kst::kst_cast<HealpixSource>(_instance);
1498         assert(src);
1499         _cfg->endGroup();
1500         _cfg->beginGroup(src->fileName());
1501         _hc->matDimX->setValue(_cfg->value("Matrix X Dimension", DEFAULT_XDIM).toInt());
1502         _hc->matDimY->setValue(_cfg->value("Matrix Y Dimension", DEFAULT_YDIM).toInt());
1503         _hc->matThetaAuto->setChecked(_cfg->value("Theta Autoscale", true).toBool());
1504         //FIXME Only one combobox item for now...
1505         //_hc->matThetaUnits->setCurrentItem(_cfg->readEntry("Theta Units", HPUNIT_RAD));
1506         _hc->matThetaUnits->setCurrentItem(0);
1507         _hc->matThetaMin->setText(_cfg->value("Theta Min").toString());
1508         _hc->matThetaMax->setText(_cfg->value("Theta Max").toString());
1509         _hc->matPhiAuto->setChecked(_cfg->value("Phi Autoscale", true).toBool());
1510         //FIXME Only one combobox item for now...
1511         //_hc->matPhiUnits->setCurrentItem(_cfg->readEntry("Phi Units", HPUNIT_RAD));
1512         _hc->matPhiUnits->setCurrentItem(0);
1513         _hc->matPhiMin->setText(_cfg->value("Phi Min").toString());
1514         _hc->matPhiMax->setText(_cfg->value("Phi Max").toString());
1515         _hc->vecTheta->setCurrentItem(_cfg->value("Vector Theta", 0).toInt());
1516         _hc->vecPhi->setCurrentItem(_cfg->value("Vector Phi", 0).toInt());
1517         _hc->vecDegrade->setValue(_cfg->value("Vector Degrade Factor", 1).toInt());
1518         _hc->vecMagAuto->setChecked(_cfg->value("Vector Magnitude Autoscale", true).toBool());
1519         _hc->vecMag->setText(_cfg->value("Vector Max Magnitude").toString());
1520         _hc->vecQU->setChecked(_cfg->value("Vector is QU", false).toBool());
1521         _cfg->endGroup();
1522       }
1523     }
1524
1525     void save() {
1526       assert(_cfg);
1527       _cfg->beginGroup("Healpix General");
1528       // If we have an instance, save settings for that instance only
1529       Kst::SharedPtr<HealpixSource> src = Kst::kst_cast<HealpixSource>(_instance);
1530       if (src) {
1531         _cfg->endGroup();
1532         _cfg->beginGroup(src->fileName());
1533       }
1534       _cfg->writeEntry("Matrix X Dimension", _hc->matDimX->value());
1535       _cfg->writeEntry("Matrix Y Dimension", _hc->matDimY->value());
1536       _cfg->writeEntry("Theta Autoscale", _hc->matThetaAuto->isChecked());
1537       //FIXME override for now
1538       //_cfg->writeEntry("Theta Units", _hc->matThetaUnits->currentItem());
1539       _cfg->writeEntry("Theta Units", HPUNIT_RADEC);
1540       _cfg->writeEntry("Theta Min", _hc->matThetaMin->text());
1541       _cfg->writeEntry("Theta Max", _hc->matThetaMax->text());
1542       _cfg->writeEntry("Phi Autoscale", _hc->matPhiAuto->isChecked());
1543       //FIXME override for now
1544       //_cfg->writeEntry("Phi Units", _hc->matPhiUnits->currentItem());
1545       _cfg->writeEntry("Phi Units", HPUNIT_RADEC);
1546       _cfg->writeEntry("Phi Min", _hc->matPhiMin->text());
1547       _cfg->writeEntry("Phi Max", _hc->matPhiMax->text());
1548       _cfg->writeEntry("Vector Theta", _hc->vecTheta->currentItem());
1549       _cfg->writeEntry("Vector Phi", _hc->vecPhi->currentItem());
1550       _cfg->writeEntry("Vector Degrade Factor", _hc->vecDegrade->value());
1551       _cfg->writeEntry("Vector Magnitude Autoscale", _hc->vecMagAuto->isChecked());
1552       _cfg->writeEntry("Vector Max Magnitude", _hc->vecMag->text());
1553       _cfg->writeEntry("Vector is QU", _hc->vecQU->isChecked());
1554       _cfg->endGroup();
1555
1556       // Update the instance from our new settings
1557       if (src && src->reusable()) {
1558         src->_config->read(_cfg, src->fileName());
1559         src->reset();
1560       }
1561     }
1562
1563     ConfigWidgetHealpixInternal *_hc;
1564 };
1565
1566
1567 QString HealpixPlugin::pluginName() const { return "HEALPIX Source Reader"; }
1568 QString HealpixPlugin::pluginDescription() const { return "HEALPIX Source Reader"; }
1569
1570
1571 Kst::DataSource *HealpixPlugin::create(Kst::ObjectStore *store,
1572                                             QSettings *cfg,
1573                                             const QString &filename,
1574                                             const QString &type,
1575                                             const QDomElement &element) const {
1576
1577   return new HealpixSource(store, cfg, filename, type, element);
1578 }
1579
1580
1581
1582 QStringList HealpixPlugin::matrixList(QSettings *cfg,
1583                                              const QString& filename,
1584                                              const QString& type,
1585                                              QString *typeSuggestion,
1586                                              bool *complete) const {
1587   Q_UNUSED(type)
1588   QStringList matrixList;
1589
1590   if (complete) {
1591     *complete = true;
1592   }
1593
1594   if (typeSuggestion) {
1595     *typeSuggestion = healpixTypeString;
1596   }
1597   if ( understands(cfg, filename) ) {
1598     size_t poff;
1599
1600     int mapType;
1601     int mapOrder;
1602     int mapCoord;
1603     size_t mapNside;
1604     size_t nMaps;
1605     char healpixfile[HEALPIX_STRNL];
1606     healpix_keys *keys;
1607     char creator[HEALPIX_STRNL];
1608     char extname[HEALPIX_STRNL];
1609     char **names;
1610     char **units;
1611
1612     strncpy(healpixfile, filename.toLatin1().data(), HEALPIX_STRNL);
1613
1614     names = healpix_strarr_alloc(HEALPIX_FITS_MAXCOL);
1615     units = healpix_strarr_alloc(HEALPIX_FITS_MAXCOL);
1616     keys = healpix_keys_alloc();
1617
1618     if (healpix_fits_map_info(healpixfile, &mapNside, &mapOrder, &mapCoord, &mapType, &nMaps, creator, extname, names, units, keys)) {
1619       if (mapType == HEALPIX_FITS_CUT) {
1620         poff = 1;
1621       } else {
1622         poff = 0;
1623       }
1624
1625       // populate the field list
1626       QString mapName;
1627       for (size_t i = 0; i < nMaps; i++) {
1628         if (strlen(names[i+poff]) == 0) {
1629           mapName.sprintf("%d - %s",(int)(i+1),"MAP");
1630         } else {
1631           mapName.sprintf("%d - %s",(int)(i+1),names[i+poff]);
1632         }
1633         if (strlen(units[i+poff]) == 0) {
1634           mapName.sprintf("%s (%s)",mapName.toLatin1().data(),"Unknown Units");
1635         } else {
1636           mapName.sprintf("%s (%s)",mapName.toLatin1().data(),units[i+poff]);
1637         }
1638         matrixList.append(mapName);
1639       }
1640
1641       if (mapType == HEALPIX_FITS_CUT) {
1642         if (strlen(names[nMaps+1]) == 0) {
1643           mapName.sprintf("%s","HITS");
1644         } else {
1645           mapName.sprintf("%s",names[nMaps+1]);
1646         }
1647         matrixList.append(mapName);
1648         if (strlen(names[nMaps+2]) == 0) {
1649           mapName.sprintf("%s","ERRORS");
1650         } else {
1651           mapName.sprintf("%s",names[nMaps+2]);
1652         }
1653         if (strlen(units[nMaps+2]) == 0) {
1654           mapName.sprintf("%s (%s)",mapName.toLatin1().data(),"Unknown Units");
1655         } else {
1656           mapName.sprintf("%s (%s)",mapName.toLatin1().data(),units[nMaps+2]);
1657         }
1658         matrixList.append(mapName);
1659       }
1660     }
1661
1662     healpix_keys_free(keys);
1663     healpix_strarr_free(names, HEALPIX_FITS_MAXCOL);
1664     healpix_strarr_free(units, HEALPIX_FITS_MAXCOL);
1665   }
1666
1667   return matrixList;
1668
1669 }
1670
1671
1672 QStringList HealpixPlugin::scalarList(QSettings *cfg,
1673                                             const QString& filename,
1674                                             const QString& type,
1675                                             QString *typeSuggestion,
1676                                             bool *complete) const {
1677
1678   QStringList scalarList;
1679
1680   if ((!type.isEmpty() && !provides().contains(type)) || 0 == understands(cfg, filename)) {
1681     if (complete) {
1682       *complete = false;
1683     }
1684     return QStringList();
1685   }
1686
1687   if (typeSuggestion) {
1688     *typeSuggestion = healpixTypeString;
1689   }
1690
1691   scalarList.append("FRAMES");
1692   return scalarList;
1693
1694 }
1695
1696
1697 QStringList HealpixPlugin::stringList(QSettings *cfg,
1698                                       const QString& filename,
1699                                       const QString& type,
1700                                       QString *typeSuggestion,
1701                                       bool *complete) const {
1702
1703   QStringList stringList;
1704
1705   if ((!type.isEmpty() && !provides().contains(type)) || 0 == understands(cfg, filename)) {
1706     if (complete) {
1707       *complete = false;
1708     }
1709     return QStringList();
1710   }
1711
1712   if (typeSuggestion) {
1713     *typeSuggestion = healpixTypeString;
1714   }
1715
1716   stringList.append("FILENAME");
1717
1718   int mapType;
1719   int mapOrder;
1720   int mapCoord;
1721   size_t mapNside;
1722   size_t nMaps;
1723   char healpixfile[HEALPIX_STRNL];
1724   healpix_keys *keys;
1725   char creator[HEALPIX_STRNL];
1726   char extname[HEALPIX_STRNL];
1727   char **names;
1728   char **units;
1729
1730   strncpy(healpixfile, filename.toLatin1().data(), HEALPIX_STRNL);
1731
1732   names = healpix_strarr_alloc(HEALPIX_FITS_MAXCOL);
1733   units = healpix_strarr_alloc(HEALPIX_FITS_MAXCOL);
1734   keys = healpix_keys_alloc();
1735
1736   if (healpix_fits_map_info(healpixfile, &mapNside, &mapOrder, &mapCoord, &mapType, &nMaps, creator, extname, names, units, keys)) {
1737     // populate the metadata
1738     QString metaVal;
1739     QString metaName;
1740
1741     metaName = "NSIDE";
1742     metaVal.sprintf("%lu", (long unsigned int)mapNside);
1743     stringList.append(metaName);
1744
1745     metaName = "NPIX";
1746     size_t mapNpix = 12 * mapNside * mapNside;
1747     metaVal.sprintf("%lu", (long unsigned int)mapNpix);
1748     stringList.append(metaName);
1749
1750     metaName = "ORDER";
1751     metaVal.sprintf("%d", mapOrder);
1752     stringList.append(metaName);
1753
1754     metaName = "COORD";
1755     metaVal.sprintf("%d", mapCoord);
1756     stringList.append(metaName);
1757
1758     metaName = "TYPE";
1759     metaVal.sprintf("%d", mapType);
1760     stringList.append(metaName);
1761
1762     metaName = "NMAPS";
1763     metaVal.sprintf("%lu", (long unsigned int)nMaps);
1764     stringList.append(metaName);
1765
1766     metaName = "CREATOR";
1767     metaVal.sprintf("%s", creator);
1768     stringList.append(metaName);
1769
1770     metaName = "EXTNAME";
1771     metaVal.sprintf("%s", extname);
1772     stringList.append(metaName);
1773
1774     for (size_t j = 0; j < keys->nskeys; j++) {
1775       metaName.sprintf("%s", keys->skeynames[j]);
1776       metaVal.sprintf("%s", keys->skeyvals[j]);
1777       stringList.append(metaName);
1778     }
1779
1780     for (size_t j = 0; j < keys->nikeys; j++) {
1781       metaName.sprintf("%s", keys->ikeynames[j]);
1782       metaVal.sprintf("%d", keys->ikeyvals[j]);
1783       stringList.append(metaName);
1784     }
1785
1786     for (size_t j = 0; j < keys->nfkeys; j++) {
1787       metaName.sprintf("%s", keys->fkeynames[j]);
1788       metaVal.sprintf("%e", keys->fkeyvals[j]);
1789       stringList.append(metaName);
1790     }
1791   }
1792
1793   healpix_keys_free(keys);
1794   healpix_strarr_free(names, HEALPIX_FITS_MAXCOL);
1795   healpix_strarr_free(units, HEALPIX_FITS_MAXCOL);
1796
1797   return stringList;
1798
1799 }
1800
1801 QStringList HealpixPlugin::fieldList(QSettings *cfg,
1802                                             const QString& filename,
1803                                             const QString& type,
1804                                             QString *typeSuggestion,
1805                                             bool *complete) const {
1806   Q_UNUSED(type)
1807   QStringList fieldList;
1808
1809   if (complete) {
1810     *complete = true;
1811   }
1812
1813   if (typeSuggestion) {
1814     *typeSuggestion = healpixTypeString;
1815   }
1816   if ( understands(cfg, filename) ) {
1817     fieldList.append("1 - Vector Field Head Theta");
1818     fieldList.append("2 - Vector Field Head Phi");
1819     fieldList.append("3 - Vector Field Tail Theta");
1820     fieldList.append("4 - Vector Field Tail Phi");
1821   }
1822   return fieldList;
1823 }
1824
1825
1826 int HealpixPlugin::understands(QSettings *cfg, const QString& filename) const {
1827   Q_UNUSED(cfg)
1828
1829   int ret;
1830   char thealpixfile[HEALPIX_STRNL];
1831   int tOrder;
1832   int tCoord;
1833   int tType;
1834   size_t tNside;
1835   size_t tMaps;
1836
1837   strncpy(thealpixfile, filename.toLatin1().data(), HEALPIX_STRNL);
1838   ret = healpix_fits_map_test(thealpixfile, &tNside, &tOrder, &tCoord, &tType, &tMaps);
1839
1840   if (ret) {
1841     // MUST return 100, since LFIIO datasource returns 99
1842     // for *all* valid FITS files
1843     return 100;
1844   } else {
1845     return 0;
1846   }
1847 }
1848
1849
1850 bool HealpixPlugin::supportsTime(QSettings *cfg, const QString& filename) const {
1851   //FIXME
1852   Q_UNUSED(cfg)
1853   Q_UNUSED(filename)
1854   return true;
1855 }
1856
1857
1858 QStringList HealpixPlugin::provides() const {
1859   QStringList rc;
1860   rc += "HEALPIX Source";
1861   return rc;
1862 }
1863
1864
1865 Kst::DataSourceConfigWidget *HealpixPlugin::configWidget(QSettings *cfg, const QString& filename) const {
1866   Q_UNUSED(filename)
1867   ConfigWidgetHealpix *config = new ConfigWidgetHealpix;
1868   config->setConfig(cfg);
1869   config->load();
1870   return config;
1871 }
1872
1873 Q_EXPORT_PLUGIN2(kstdata_qimagesource, HealpixPlugin)
1874
1875
1876 // vim: ts=2 sw=2 et