| 1 | #include "sopnamsp.h"
|
|---|
| 2 | #include "pexceptions.h"
|
|---|
| 3 | #include "fitsspherehealpix.h"
|
|---|
| 4 | #include "tarray.h"
|
|---|
| 5 | #include "fitstarray.h"
|
|---|
| 6 |
|
|---|
| 7 |
|
|---|
| 8 |
|
|---|
| 9 | ///////////////////////////////////////////////////////////////////////
|
|---|
| 10 | // ----objets delegues pour la gestion de persistance I/O format fits--
|
|---|
| 11 | // SphereHealpix
|
|---|
| 12 | ////////////////////////////////////////////////////////////////////
|
|---|
| 13 |
|
|---|
| 14 | /*!
|
|---|
| 15 | \class SOPHYA::FITS_SphereHEALPix
|
|---|
| 16 | \ingroup FitsIOServer
|
|---|
| 17 | FITS format I/O handler for SOPHYA::SphereHEALPix objects.
|
|---|
| 18 | */
|
|---|
| 19 |
|
|---|
| 20 |
|
|---|
| 21 | template <class T>
|
|---|
| 22 | FITS_SphereHEALPix<T>::FITS_SphereHEALPix()
|
|---|
| 23 | {
|
|---|
| 24 | dobj_= new SphereHEALPix<T>;
|
|---|
| 25 | ownobj_= true;
|
|---|
| 26 | }
|
|---|
| 27 |
|
|---|
| 28 | template <class T>
|
|---|
| 29 | FITS_SphereHEALPix<T>::FITS_SphereHEALPix(char inputfile[],int hdunum)
|
|---|
| 30 | {
|
|---|
| 31 | dobj_= new SphereHEALPix<T>;
|
|---|
| 32 | ownobj_= true;
|
|---|
| 33 | Read(inputfile,hdunum);
|
|---|
| 34 | }
|
|---|
| 35 |
|
|---|
| 36 |
|
|---|
| 37 | template <class T>
|
|---|
| 38 | FITS_SphereHEALPix<T>::FITS_SphereHEALPix(const SphereHEALPix<T>& obj)
|
|---|
| 39 | {
|
|---|
| 40 | dobj_= new SphereHEALPix<T>(obj);
|
|---|
| 41 | ownobj_= true;
|
|---|
| 42 | }
|
|---|
| 43 |
|
|---|
| 44 | template <class T>
|
|---|
| 45 | FITS_SphereHEALPix<T>::FITS_SphereHEALPix(SphereHEALPix<T>* obj)
|
|---|
| 46 | {
|
|---|
| 47 | dobj_= obj;
|
|---|
| 48 | ownobj_= false;
|
|---|
| 49 | }
|
|---|
| 50 |
|
|---|
| 51 | template <class T>
|
|---|
| 52 | FITS_SphereHEALPix<T>::~FITS_SphereHEALPix()
|
|---|
| 53 | {
|
|---|
| 54 | if (ownobj_ && dobj_) delete dobj_;
|
|---|
| 55 | }
|
|---|
| 56 |
|
|---|
| 57 | template <class T>
|
|---|
| 58 | AnyDataObj* FITS_SphereHEALPix<T>::DataObj()
|
|---|
| 59 | {
|
|---|
| 60 | return(dobj_);
|
|---|
| 61 | }
|
|---|
| 62 |
|
|---|
| 63 | template <class T>
|
|---|
| 64 | void FITS_SphereHEALPix<T>::SetDataObj(AnyDataObj & o)
|
|---|
| 65 | {
|
|---|
| 66 | SphereHEALPix<T> * po = dynamic_cast< SphereHEALPix<T> * >(&o);
|
|---|
| 67 | if (po == NULL) return;
|
|---|
| 68 | if (ownobj_ && dobj_) delete dobj_;
|
|---|
| 69 | dobj_ = po;
|
|---|
| 70 | ownobj_ = false;
|
|---|
| 71 | }
|
|---|
| 72 |
|
|---|
| 73 | template <class T>
|
|---|
| 74 | int FITS_SphereHEALPix<T>::CheckHandling(AnyDataObj & o)
|
|---|
| 75 | {
|
|---|
| 76 | if ( typeid(o) == typeid( SphereHEALPix< T > ) ) return 2;
|
|---|
| 77 | SphereHEALPix<T> * po = dynamic_cast< SphereHEALPix<T> * >(&o);
|
|---|
| 78 | if (po == NULL) return 0;
|
|---|
| 79 | else return 1;
|
|---|
| 80 | }
|
|---|
| 81 |
|
|---|
| 82 | template <class T>
|
|---|
| 83 | int FITS_SphereHEALPix<T>::CheckReadability(FitsInOutFile& is)
|
|---|
| 84 | {
|
|---|
| 85 | if (is.CurrentHDUType() != BINARY_TBL ) return 0;
|
|---|
| 86 | string key;
|
|---|
| 87 | key = "PIXTYPE";
|
|---|
| 88 | if ( is.KeyValue(key) != "HEALPIX") return 0;
|
|---|
| 89 |
|
|---|
| 90 | bool nosk = false;
|
|---|
| 91 | key = "ORDERING";
|
|---|
| 92 | is.KeyValue(key, nosk);
|
|---|
| 93 | bool nosk2 = false;
|
|---|
| 94 | key = "NSIDE";
|
|---|
| 95 | is.KeyValue(key, nosk2);
|
|---|
| 96 | if (nosk || nosk2) return 0;
|
|---|
| 97 |
|
|---|
| 98 | vector<string> colnames;
|
|---|
| 99 | vector<int> coltypes;
|
|---|
| 100 | vector<long> repcnt;
|
|---|
| 101 | vector<long> width;
|
|---|
| 102 | long ncols = is.GetColInfo(colnames, coltypes, repcnt, width);
|
|---|
| 103 | if (ncols < 1) return 0;
|
|---|
| 104 | T x = 0;
|
|---|
| 105 | if (coltypes[0] == FitsTypes::DataType(x)) return 2 ;
|
|---|
| 106 | else return 1;
|
|---|
| 107 | }
|
|---|
| 108 |
|
|---|
| 109 | template <class T>
|
|---|
| 110 | FitsHandlerInterface* FITS_SphereHEALPix<T>::Clone()
|
|---|
| 111 | {
|
|---|
| 112 | return new FITS_SphereHEALPix<T>() ;
|
|---|
| 113 | }
|
|---|
| 114 |
|
|---|
| 115 |
|
|---|
| 116 | template <class T>
|
|---|
| 117 | void FITS_SphereHEALPix<T>::WriteToFits(FitsOutFile& os)
|
|---|
| 118 | {
|
|---|
| 119 | if(dobj_ == NULL)
|
|---|
| 120 | {
|
|---|
| 121 | cout << " WriteToFits:: dobj_= null " << endl;
|
|---|
| 122 | return;
|
|---|
| 123 | }
|
|---|
| 124 |
|
|---|
| 125 | DVList dvl( dobj_->Info() );
|
|---|
| 126 |
|
|---|
| 127 | SphereCoordSys* cs= dobj_->GetCoordSys();
|
|---|
| 128 | string description= cs->description();
|
|---|
| 129 | dvl["COORDSYS"]= description;
|
|---|
| 130 |
|
|---|
| 131 | dvl["PIXTYPE"] = "HEALPIX";
|
|---|
| 132 | dvl.SetComment("PIXTYPE", "HEALPIX Pixelization");
|
|---|
| 133 | dvl["ORDERING"]= dobj_->TypeOfMap();
|
|---|
| 134 | dvl.SetComment("ORDERING", "Pixel ordering scheme, either RING or NESTED");
|
|---|
| 135 |
|
|---|
| 136 | int nSide= dobj_->SizeIndex();
|
|---|
| 137 | dvl["NSIDE"]= (int_4) nSide;
|
|---|
| 138 | dvl.SetComment("NSIDE","Resolution parameter for HEALPIX" );
|
|---|
| 139 |
|
|---|
| 140 | int nPix= dobj_->NbPixels();
|
|---|
| 141 | dvl["FIRSTPIX"]= (int_4) 0;
|
|---|
| 142 | dvl.SetComment("FIRSTPIX", "First pixel # (0 based)");
|
|---|
| 143 | dvl["LASTPIX"]= (int_4) nPix-1;
|
|---|
| 144 | dvl.SetComment("LASTPIX", "Last pixel # (0 based)");
|
|---|
| 145 | dvl["Content"]= "SphereHEALPix";
|
|---|
| 146 | dvl.SetComment("Content", "name of SOPHYA object");
|
|---|
| 147 | dvl["SOPCLSNM"]= "SOPHYA::SphereHEALPix<T>";
|
|---|
| 148 | dvl.SetComment("SOPCLSNM", "SOPHYA class name");
|
|---|
| 149 |
|
|---|
| 150 | // On ecrit les dataBlocks
|
|---|
| 151 | vector<string> Noms;
|
|---|
| 152 | Noms.push_back(dvl.GetS("Content"));
|
|---|
| 153 | // string extname("SIMULATION");
|
|---|
| 154 | string extname = os.NextExtensionName();
|
|---|
| 155 |
|
|---|
| 156 | string Type;
|
|---|
| 157 | if (typeid(T) == typeid(r_8) ) Type+='D';
|
|---|
| 158 | else if (typeid(T) == typeid(r_4) ) Type+='E';
|
|---|
| 159 | else if (typeid(T) == typeid(int_4) ) Type+='J';
|
|---|
| 160 | else
|
|---|
| 161 | {
|
|---|
| 162 | cout << "FITS_SphereHEALPix<T>::WriteToFits()- Unsupported data type ("
|
|---|
| 163 | << typeid(T).name() << ") throwing IOExc " << endl;
|
|---|
| 164 | throw IOExc("FITS_SphereHEALPix<T>::WriteToFits() unsupported data type");
|
|---|
| 165 | }
|
|---|
| 166 | vector<int> dummy;
|
|---|
| 167 | os.makeHeaderBntblOnFits(Type, Noms, nPix, 1, &dvl, extname, dummy);
|
|---|
| 168 | os.PutColToFits(0, nPix, dobj_->pixels_.Data());
|
|---|
| 169 | }
|
|---|
| 170 |
|
|---|
| 171 | template <class T>
|
|---|
| 172 | void FITS_SphereHEALPix<T>::ReadFromFits(FitsInFile& is)
|
|---|
| 173 | {
|
|---|
| 174 | if(dobj_ == NULL)
|
|---|
| 175 | {
|
|---|
| 176 | dobj_= new SphereHEALPix<T>;
|
|---|
| 177 | ownobj_= true;
|
|---|
| 178 | }
|
|---|
| 179 | int nbcols, nbentries;
|
|---|
| 180 | //
|
|---|
| 181 |
|
|---|
| 182 | nbcols = is.NbColsFromFits();
|
|---|
| 183 | if (nbcols != 1)
|
|---|
| 184 | cout << "FITS_SphereHEALPix<T>::ReadFromFits()/Warning Fits::NbCols="
|
|---|
| 185 | << nbcols << " not equal 1" << endl;
|
|---|
| 186 |
|
|---|
| 187 | DVList dvl=is.DVListFromFits();
|
|---|
| 188 | nbentries = is.NentriesFromFits(0);
|
|---|
| 189 | int lastpix=dvl.GetI("LASTPIX");
|
|---|
| 190 | if (lastpix>0)
|
|---|
| 191 | {
|
|---|
| 192 | if (nbentries!=lastpix+1)
|
|---|
| 193 | {
|
|---|
| 194 | nbentries=lastpix+1;
|
|---|
| 195 | }
|
|---|
| 196 | }
|
|---|
| 197 | // Let's Read the SphereCoordSys object
|
|---|
| 198 | int id= SphereCoordSys::NEUTRAL;
|
|---|
| 199 | string description= "NEUTRAL SphereCoordSystem";
|
|---|
| 200 | string coordsys= dvl.GetS("COORDSYS");
|
|---|
| 201 | // $CHECK$ - Reza 2/05/2000 - compare(size_t, size_t,...) passe pas sur g++
|
|---|
| 202 | // if(coordsys.compare(0,7,"ROTATION",0,7) == 0) $CHECK$
|
|---|
| 203 | if(coordsys.substr(0,8) == "ROTATION")
|
|---|
| 204 | {
|
|---|
| 205 | id= SphereCoordSys::ROTATION;
|
|---|
| 206 | description= "ROTATION SphereCoordSystem";
|
|---|
| 207 | }
|
|---|
| 208 | // else if(coordsys.compare(0,4,"OTHER",0,4) == 0) $CHECK$ Reza 2/05/2000
|
|---|
| 209 | else if(coordsys.substr(0,5) == "OTHER" )
|
|---|
| 210 | {
|
|---|
| 211 | id= SphereCoordSys::OTHER;
|
|---|
| 212 | description= "OTHER SphereCoordSystem";
|
|---|
| 213 | }
|
|---|
| 214 | dobj_->SetCoordSys(new SphereCoordSys(id,description));
|
|---|
| 215 |
|
|---|
| 216 | string ordering= dvl.GetS("ORDERING");
|
|---|
| 217 | if ( (ordering.substr(0,4) != "RING" ) && (ordering.substr(0,6) != "NESTED" ) ) {
|
|---|
| 218 | cerr << "FITS_SphereHEALPix/Error Not supported ORDERING= " << ordering
|
|---|
| 219 | << " != RING/NESTED " << endl;
|
|---|
| 220 | throw FitsIOException("FITS_SphereHEALPix/Error ordering != RING/NESTED ");
|
|---|
| 221 | }
|
|---|
| 222 |
|
|---|
| 223 | bool fgring = true;
|
|---|
| 224 | if (ordering.substr(0,4) != "RING" ) fgring = false;
|
|---|
| 225 |
|
|---|
| 226 | int nside= dvl.GetI("NSIDE");
|
|---|
| 227 | int nPix = 0;
|
|---|
| 228 | if(nside <= 0)
|
|---|
| 229 | {
|
|---|
| 230 | // throw IOExc("FITS_SphereHEALPix:: No resol parameter nside");
|
|---|
| 231 | if (lastpix>0)
|
|---|
| 232 | {
|
|---|
| 233 | nside = sqrt((double)(nbentries/12));
|
|---|
| 234 | }
|
|---|
| 235 | else
|
|---|
| 236 | {
|
|---|
| 237 | // dvl.Print();
|
|---|
| 238 | // int naxis2 = dvl.GetI("NAXIS2");
|
|---|
| 239 | // if (naxis2 > 0 )
|
|---|
| 240 | // {
|
|---|
| 241 | nside = sqrt((double)(nbentries/12));
|
|---|
| 242 | // }
|
|---|
| 243 | // else
|
|---|
| 244 | // throw IOExc("FITS_SphereHEALPix:: No NSIDE nor LASTPIX nor NAXIS2 on file");
|
|---|
| 245 | }
|
|---|
| 246 | }
|
|---|
| 247 | nPix = 12*nside*nside;
|
|---|
| 248 | if (nPix != nbentries)
|
|---|
| 249 | {
|
|---|
| 250 | cout << "WARNING: le nombre de pixels relu " << nbentries << " est incoherent avec la pixelisation Healpix qui donne nPix= " << nPix << endl;
|
|---|
| 251 | }
|
|---|
| 252 | double Omega= 4.0*Pi/nPix;
|
|---|
| 253 | dobj_->setParameters(nside,nPix,Omega,fgring);
|
|---|
| 254 | // On lit les DataBlocks;
|
|---|
| 255 | //
|
|---|
| 256 | dobj_->pixels_.ReSize(nPix);
|
|---|
| 257 | is.GetSingleColumn(dobj_->pixels_.Data(),nPix);
|
|---|
| 258 | //
|
|---|
| 259 |
|
|---|
| 260 | // on effectue le decoupage en tranches
|
|---|
| 261 | dobj_->SetThetaSlices();
|
|---|
| 262 | //
|
|---|
| 263 | dobj_->Info() = is.DVListFromFits();
|
|---|
| 264 | }
|
|---|
| 265 |
|
|---|
| 266 | template <class T>
|
|---|
| 267 | void FITS_SphereHEALPix<T>::Mollweide_picture_projection(char filename[])
|
|---|
| 268 | {
|
|---|
| 269 | int ni = 300;
|
|---|
| 270 | int nj = 600;
|
|---|
| 271 | TMatrix<float> map(ni, nj);
|
|---|
| 272 |
|
|---|
| 273 | int i,j;
|
|---|
| 274 | for(i=0; i < ni; i++)
|
|---|
| 275 | {
|
|---|
| 276 | double yd = (i+0.5)/ni-0.5;
|
|---|
| 277 | double facteur=2.*Pi/sin(acos(yd*2));
|
|---|
| 278 | double theta = (0.5-yd)*Pi;
|
|---|
| 279 | for(j=0; j < nj; j++)
|
|---|
| 280 | {
|
|---|
| 281 | double xa = (j+0.5)/nj-0.5;
|
|---|
| 282 | double phi = xa*facteur+Pi;
|
|---|
| 283 | float th=float(theta);
|
|---|
| 284 | float ff=float(phi);
|
|---|
| 285 | if (phi<2*Pi && phi>= 0)
|
|---|
| 286 | {
|
|---|
| 287 | map(i,j) = (float)dobj_->PixValSph(th, ff);
|
|---|
| 288 | }
|
|---|
| 289 | }
|
|---|
| 290 | }
|
|---|
| 291 | FitsOutFile fout(filename);
|
|---|
| 292 | fout.firstImageOnPrimaryHeader(true);
|
|---|
| 293 | FITS_TArray<r_4> fta(map);
|
|---|
| 294 | fta.Write(fout);
|
|---|
| 295 |
|
|---|
| 296 | }
|
|---|
| 297 | template <class T>
|
|---|
| 298 | void FITS_SphereHEALPix<T>::sinus_picture_projection(char filename[])
|
|---|
| 299 | {
|
|---|
| 300 | int ni = 300;
|
|---|
| 301 | int nj = 600;
|
|---|
| 302 | TMatrix<float> map(ni, nj);
|
|---|
| 303 |
|
|---|
| 304 | int i,j;
|
|---|
| 305 | for(i=0; i < ni; i++)
|
|---|
| 306 | {
|
|---|
| 307 | double yd = (i+0.5)/ni-0.5;
|
|---|
| 308 | double theta = (0.5-yd)*Pi;
|
|---|
| 309 | double facteur=1./sin(theta);
|
|---|
| 310 | for(j=0; j < nj; j++)
|
|---|
| 311 | {
|
|---|
| 312 | double xa = (j+0.5)/nj-0.5;
|
|---|
| 313 | double phi = 2.*Pi*xa*facteur+Pi;
|
|---|
| 314 | float th=float(theta);
|
|---|
| 315 | float ff=float(phi);
|
|---|
| 316 | if (phi<2*Pi && phi>= 0)
|
|---|
| 317 | {
|
|---|
| 318 | map(i,j) = (float)dobj_->PixValSph(th, ff);
|
|---|
| 319 | }
|
|---|
| 320 | }
|
|---|
| 321 | }
|
|---|
| 322 | FitsOutFile fout(filename);
|
|---|
| 323 | fout.firstImageOnPrimaryHeader(true);
|
|---|
| 324 | FITS_TArray<r_4> fta(map);
|
|---|
| 325 | fta.Write(fout);
|
|---|
| 326 |
|
|---|
| 327 |
|
|---|
| 328 | }
|
|---|
| 329 |
|
|---|
| 330 |
|
|---|
| 331 | #ifdef __CXX_PRAGMA_TEMPLATES__
|
|---|
| 332 | #pragma define_template FITS_SphereHEALPix<r_8>
|
|---|
| 333 | #pragma define_template FITS_SphereHEALPix<r_4>
|
|---|
| 334 | #pragma define_template FITS_SphereHEALPix<int_4>
|
|---|
| 335 | #endif
|
|---|
| 336 | #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
|
|---|
| 337 | namespace SOPHYA {
|
|---|
| 338 | template class FITS_SphereHEALPix<r_8>;
|
|---|
| 339 | template class FITS_SphereHEALPix<r_4>;
|
|---|
| 340 | template class FITS_SphereHEALPix<int_4>;
|
|---|
| 341 | }
|
|---|
| 342 | #endif
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
