source: Sophya/trunk/SophyaExt/FitsIOServer/fitsspherehealpix.cc@ 1176

Last change on this file since 1176 was 1176, checked in by ansari, 25 years ago

Mollweide et sinus projections
File size: 6.8 KB
Line 
1#include "pexceptions.h"
2#include "fitsspherehealpix.h"
3#include "tarray.h"
4#include "fitstarray.h"
5
6
7
8///////////////////////////////////////////////////////////////////////
9// ----objets delegues pour la gestion de persistance I/O format fits--
10// SphereHealpix
11////////////////////////////////////////////////////////////////////
12
13
14template <class T>
15FITS_SphereHEALPix<T>::FITS_SphereHEALPix()
16{
17 dobj_= new SphereHEALPix<T>;
18 ownobj_= true;
19}
20
21template <class T>
22FITS_SphereHEALPix<T>::FITS_SphereHEALPix(char inputfile[],int hdunum)
23{
24 dobj_= new SphereHEALPix<T>;
25 ownobj_= true;
26 Read(inputfile,hdunum);
27}
28
29
30template <class T>
31FITS_SphereHEALPix<T>::FITS_SphereHEALPix(const SphereHEALPix<T>& obj)
32{
33 cout << " constructeur &obj " << endl;
34 dobj_= new SphereHEALPix<T>(obj);
35 ownobj_= true;
36}
37
38template <class T>
39FITS_SphereHEALPix<T>::FITS_SphereHEALPix(SphereHEALPix<T>* obj)
40{
41 dobj_= obj;
42 ownobj_= false;
43}
44
45template <class T>
46FITS_SphereHEALPix<T>::~FITS_SphereHEALPix()
47{
48 if (ownobj_ && dobj_) delete dobj_;
49}
50
51template <class T>
52AnyDataObj* FITS_SphereHEALPix<T>::DataObj()
53{
54 return(dobj_);
55}
56
57template <class T>
58void FITS_SphereHEALPix<T>::SetDataObj(AnyDataObj & o)
59{
60 SphereHEALPix<T> * po = dynamic_cast< SphereHEALPix<T> * >(&o);
61 if (po == NULL) return;
62 if (ownobj_ && dobj_) delete dobj_;
63 dobj_ = po;
64 ownobj_ = false;
65}
66
67
68
69
70template <class T>
71void FITS_SphereHEALPix<T>::WriteToFits(FitsOutFile& os)
72{
73 if(dobj_ == NULL)
74 {
75 cout << " WriteToFits:: dobj_= null " << endl;
76 return;
77 }
78
79 DVList dvl( dobj_->Info() );
80
81 SphereCoordSys* cs= dobj_->GetCoordSys();
82 string description= cs->description();
83 dvl["COORDSYS"]= description;
84
85 dvl["PIXTYPE"] = "HEALPIX";
86 dvl.SetComment("PIXTYPE", "HEALPIX Pixelization");
87 dvl["ORDERING"]= dobj_->TypeOfMap();
88 dvl.SetComment("ORDERING", "Pixel ordering scheme, either RING or NESTED");
89
90 int nSide= dobj_->SizeIndex();
91 dvl["NSIDE"]= (int_4) nSide;
92 dvl.SetComment("NSIDE","Resolution parameter for HEALPIX" );
93
94 int nPix= dobj_->NbPixels();
95 dvl["FIRSTPIX"]= (int_4) 0;
96 dvl.SetComment("FIRSTPIX", "First pixel # (0 based)");
97 dvl["LASTPIX"]= (int_4) nPix-1;
98 dvl.SetComment("LASTPIX", "Last pixel # (0 based)");
99 dvl["Content"]= "SphereHEALPix";
100
101 // On ecrit les dataBlocks
102 char** Noms = new char*[1];
103 Noms[0]= new char[15];
104 strncpy(Noms[0],dvl.GetS("Content").c_str(),15);
105 char extname[15] = "SIMULATION";
106
107 char Type[2];
108 if (typeid(T) == typeid(r_8) ) Type[0]='D';
109 else
110 if (typeid(T) == typeid(r_4) ) Type[0]='E';
111 else
112 {
113 cout << " type de la sphere= " << typeid(T).name() << endl;
114 throw IOExc("FITS_SphereHEALPix:: unknown type");
115 }
116 Type[1]='\0';
117 vector<int> dummy;
118 os.makeHeaderBntblOnFits(Type, Noms, nPix, 1, dvl, extname, dummy);
119 delete [] Noms[0];
120 delete [] Noms;
121 os.putColToFits(0, nPix, dobj_->pixels_.Data());
122}
123
124template <class T>
125void FITS_SphereHEALPix<T>::ReadFromFits(FitsInFile& is)
126{
127 if(dobj_ == NULL)
128 {
129 dobj_= new SphereHEALPix<T>;
130 ownobj_= true;
131 }
132 int nbcols, nbentries;
133 //
134
135 nbcols = is.NbColsFromFits();
136 if (nbcols != 1)
137 {
138 throw IOExc("le fichier fits n'est pas une sphere Healpix");
139 }
140 DVList dvl=is.DVListFromFits();
141 nbentries = is.NentriesFromFits(0);
142 int lastpix=dvl.GetI("LASTPIX");
143 if (lastpix>0)
144 {
145 if (nbentries!=lastpix+1)
146 {
147 nbentries=lastpix+1;
148 }
149 }
150 // Let's Read the SphereCoordSys object
151 int id= SphereCoordSys_NEUTRAL;
152 string description= "NEUTRAL SphereCoordSystem";
153 string coordsys= dvl.GetS("COORDSYS");
154 // $CHECK$ - Reza 2/05/2000 - compare(size_t, size_t,...) passe pas sur g++
155 // if(coordsys.compare(0,7,"ROTATION",0,7) == 0) $CHECK$
156 if(coordsys.substr(0,8) == "ROTATION")
157 {
158 id= SphereCoordSys_ROTATION;
159 description= "ROTATION SphereCoordSystem";
160 }
161 // else if(coordsys.compare(0,4,"OTHER",0,4) == 0) $CHECK$ Reza 2/05/2000
162 else if(coordsys.substr(0,5) == "OTHER" )
163 {
164 id= SphereCoordSys_OTHER;
165 description= "OTHER SphereCoordSystem";
166 }
167 dobj_->SetCoordSys(new SphereCoordSys(id,description));
168
169 string ordering= dvl.GetS("ORDERING");
170 // if(ordering.compare(0,3,"RING",0,3) != 0) $CHECK$ Reza 2/05/2000
171 if(ordering.substr(0,4) != "RING" )
172 {
173 cerr << "FITS_SphereHEALPix/Error Not supported ORDERING= "
174 << ordering << " substr(0,4)=[" << ordering.substr(0,4) << "]" << endl;
175 throw IOExc(" FITS_SphereHEALPix:: numerotation non RING");
176 }
177
178 int nside= dvl.GetI("NSIDE");
179 if(nside <= 0)
180 throw IOExc("FITS_SphereHEALPix:: No resol parameter nside");
181 int nPix = 12*nside*nside;
182 if (nPix != nbentries)
183 {
184 cout << "WARNING: le nombre de pixels relu " << nbentries << " est incoherent avec la pixelisation Healpix qui donne nPix= " << nPix << endl;
185 }
186 double Omega= 4.0*Pi/nPix;
187 dobj_->setParameters(nside,nPix,Omega);
188 // On lit les DataBlocks;
189 //
190 dobj_->pixels_.ReSize(nPix);
191 is.GetSingleColumn(dobj_->pixels_.Data(),nPix);
192 //
193
194 // on effectue le decoupage en tranches
195 dobj_->SetThetaSlices();
196 //
197 dobj_->Info() = is.DVListFromFits();
198 cout << " end of file reading " << endl;
199}
200
201template <class T>
202void FITS_SphereHEALPix<T>::Mollweide_picture_projection(char filename[])
203{
204 int ni = 300;
205 int nj = 600;
206 TMatrix<float> map(ni, nj);
207
208 int i,j;
209 for(i=0; i < ni; i++)
210 {
211 double yd = (i+0.5)/ni-0.5;
212 double facteur=2.*Pi/sin(acos(yd*2));
213 double theta = (0.5-yd)*Pi;
214 for(j=0; j < nj; j++)
215 {
216 double xa = (j+0.5)/nj-0.5;
217 double phi = xa*facteur+Pi;
218 float th=float(theta);
219 float ff=float(phi);
220 if (phi<2*Pi && phi>= 0)
221 {
222 map(i,j) = (float)dobj_->PixValSph(th, ff);
223 }
224 }
225 }
226 FitsOutFile fout(filename);
227 fout.firstImageOnPrimaryHeader();
228 FITS_TArray<r_4> fta(map);
229 fta.Write(fout);
230
231}
232template <class T>
233void FITS_SphereHEALPix<T>::sinus_picture_projection(char filename[])
234{
235 int ni = 300;
236 int nj = 600;
237 TMatrix<float> map(ni, nj);
238
239 int i,j;
240 for(i=0; i < ni; i++)
241 {
242 double yd = (i+0.5)/ni-0.5;
243 double theta = (0.5-yd)*Pi;
244 double facteur=1./sin(theta);
245 for(j=0; j < nj; j++)
246 {
247 double xa = (j+0.5)/nj-0.5;
248 double phi = 2.*Pi*xa*facteur+Pi;
249 float th=float(theta);
250 float ff=float(phi);
251 if (phi<2*Pi && phi>= 0)
252 {
253 map(i,j) = (float)dobj_->PixValSph(th, ff);
254 }
255 }
256 }
257 FitsOutFile fout(filename);
258 fout.firstImageOnPrimaryHeader();
259 FITS_TArray<r_4> fta(map);
260 fta.Write(fout);
261
262
263}
264
265
266#ifdef __CXX_PRAGMA_TEMPLATES__
267#pragma define_template FITS_SphereHEALPix<r_8>
268#pragma define_template FITS_SphereHEALPix<r_4>
269#endif
270#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
271template class FITS_SphereHEALPix<r_8>;
272template class FITS_SphereHEALPix<r_4>;
273#endif
Note: See TracBrowser for help on using the repository browser.