source: Sophya/trunk/SophyaExt/FitsIOServer/fitsspherethetaphi.cc@ 2050

Last change on this file since 2050 was 1752, checked in by lemeur, 24 years ago

interface avec spheres theta-phi
  • Property svn:executable set to *
File size: 6.6 KB
Line 
1#include "pexceptions.h"
2#include "fitsspherethetaphi.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// SphereThetaPhi
11////////////////////////////////////////////////////////////////////
12
13/*!
14 \class SOPHYA::FITS_SphereThetaPhi
15 \ingroup FitsIOServer
16 FITS format I/O handler for SOPHYA::SphereThetaPhi objects.
17*/
18
19template <class T>
20FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi()
21{
22 dobj_= new SphereThetaPhi<T>;
23 ownobj_= true;
24}
25
26template <class T>
27FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi(char inputfile[],int hdunum)
28{
29 dobj_= new SphereThetaPhi<T>;
30 ownobj_= true;
31 Read(inputfile,hdunum);
32}
33
34template <class T>
35FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi(const SphereThetaPhi<T>& obj)
36{
37 dobj_= new SphereThetaPhi<T>(obj);
38 ownobj_= true;
39}
40
41template <class T>
42FITS_SphereThetaPhi<T>::FITS_SphereThetaPhi(SphereThetaPhi<T>* obj)
43{
44 dobj_= obj;
45 ownobj_= false;
46}
47
48template <class T>
49FITS_SphereThetaPhi<T>::~FITS_SphereThetaPhi()
50{
51 if (ownobj_ && dobj_) delete dobj_;
52}
53
54template <class T>
55AnyDataObj* FITS_SphereThetaPhi<T>::DataObj()
56{
57 return(dobj_);
58}
59
60template <class T>
61void FITS_SphereThetaPhi<T>::SetDataObj(AnyDataObj & o)
62{
63 SphereThetaPhi<T> * po = dynamic_cast< SphereThetaPhi<T> * >(&o);
64 if (po == NULL) return;
65 if (ownobj_ && dobj_) delete dobj_;
66 dobj_ = po;
67 ownobj_ = false;
68}
69
70template <class T>
71void FITS_SphereThetaPhi<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"] = "TETAFI";
86 dvl.SetComment("PIXTYPE", "Theta-Phi-LAL Pixelization");
87 dvl["ORDERING"]= dobj_->TypeOfMap();
88 dvl.SetComment("ORDERING", " (unrelevant) ");
89
90 int IndTheta= dobj_->SizeIndex();
91 dvl["INDTHMAX"]= (int_4) IndTheta;
92 dvl.SetComment("INDTHMAX","Resolution parameter for Theta-Phi pixelization sch." );
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"]= "SphereThetaPhi";
100 dvl.SetComment("Content", "name of SOPHYA object");
101
102 // On ecrit les dataBlocks
103 vector<string> Noms(4);
104 Noms[0] = dvl.GetS("Content");
105 Noms[1] = "NphiParBande";
106 Noms[2] = "CumulPixParBande";
107 Noms[3] = "ThetaBande";
108 string extname("SIMULATION");
109
110
111 string Type;
112 if (typeid(T) == typeid(r_8) ) Type+='D';
113 else
114 if (typeid(T) == typeid(r_4) ) Type+='E';
115 else
116 {
117 cout << " type de la sphere= " << typeid(T).name() << endl;
118 throw IOExc("FITS_SphereThetaPhi:: unknown type");
119 }
120
121 Type += 'J';
122 Type += 'J';
123 Type += 'D';
124
125 vector<int> dummy;
126 os.makeHeaderBntblOnFits(Type, Noms, nPix, 4, &dvl, extname, dummy);
127 os.PutColToFits(0, nPix, dobj_->pixels_.Data());
128 os.PutColToFits(1, IndTheta+1, dobj_->NPhi_.Data());
129 os.PutColToFits(2, IndTheta+1, dobj_->TNphi_.Data());
130 os.PutColToFits(3, IndTheta+1, dobj_->Theta_.Data());
131}
132
133template <class T>
134void FITS_SphereThetaPhi<T>::ReadFromFits(FitsInFile& is)
135{
136
137 if(dobj_ == NULL)
138 {
139 dobj_= new SphereThetaPhi<T>;
140 ownobj_= true;
141 }
142 int nbcols, nbentries;
143
144 nbcols = is.NbColsFromFits();
145 if (nbcols != 4)
146 {
147 throw IOExc("le fichier fits n'est pas une sphere ThetaPhi");
148 }
149 DVList dvl=is.DVListFromFits();
150 nbentries = is.NentriesFromFits(0);
151 int lastpix=dvl.GetI("LASTPIX");
152 if (lastpix>0)
153 {
154 if (nbentries!=lastpix+1)
155 {
156 nbentries=lastpix+1;
157 }
158 }
159
160 // Let's Read the SphereCoordSys object
161 int id= SphereCoordSys_NEUTRAL;
162 string description= "NEUTRAL SphereCoordSystem";
163 string coordsys= dvl.GetS("COORDSYS");
164 if(coordsys.substr(0,8) == "ROTATION")
165 {
166 id= SphereCoordSys_ROTATION;
167 description= "ROTATION SphereCoordSystem";
168 }
169 else if(coordsys.substr(0,5) == "OTHER" )
170 {
171 id= SphereCoordSys_OTHER;
172 description= "OTHER SphereCoordSystem";
173 }
174 dobj_->SetCoordSys(new SphereCoordSys(id,description));
175
176 int IndThetaMax= dvl.GetI("INDTHMAX");
177 int nPix = 0;
178 if(IndThetaMax <= 0)
179 {
180 throw IOExc("FITS_SphereTheta:: No resol parameter INDTHMAX");
181 }
182 nPix = nbentries;
183 double Omega= 4.0*Pi/nPix;
184 dobj_->setParameters(IndThetaMax,nPix,Omega);
185
186 // On lit les DataBlocks;
187 //
188 dobj_->pixels_.ReSize(nPix);
189 is.GetBinTabFCol(dobj_->pixels_.Data(),nPix,0);
190 dobj_->NPhi_.ReSize(IndThetaMax+1);
191 is.GetBinTabFCol(dobj_->NPhi_.Data(),IndThetaMax+1,1);
192 dobj_->TNphi_.ReSize(IndThetaMax+1);
193 is.GetBinTabFCol(dobj_->TNphi_.Data(),IndThetaMax+1,2);
194 dobj_->Theta_.ReSize(IndThetaMax+1);
195 is.GetBinTabFCol(dobj_->Theta_.Data(),IndThetaMax+1,3);
196}
197
198template <class T>
199void FITS_SphereThetaPhi<T>::Mollweide_picture_projection(char filename[])
200{
201 int ni = 300;
202 int nj = 600;
203 TMatrix<float> map(ni, nj);
204
205 int i,j;
206 for(i=0; i < ni; i++)
207 {
208 double yd = (i+0.5)/ni-0.5;
209 double facteur=2.*Pi/sin(acos(yd*2));
210 double theta = (0.5-yd)*Pi;
211 for(j=0; j < nj; j++)
212 {
213 double xa = (j+0.5)/nj-0.5;
214 double phi = xa*facteur+Pi;
215 float th=float(theta);
216 float ff=float(phi);
217 if (phi<2*Pi && phi>= 0)
218 {
219 map(i,j) = (float)dobj_->PixValSph(th, ff);
220 }
221 }
222 }
223 FitsOutFile fout(filename);
224 fout.firstImageOnPrimaryHeader(true);
225 FITS_TArray<r_4> fta(map);
226 fta.Write(fout);
227
228}
229template <class T>
230void FITS_SphereThetaPhi<T>::sinus_picture_projection(char filename[])
231{
232 int ni = 300;
233 int nj = 600;
234 TMatrix<float> map(ni, nj);
235
236 int i,j;
237 for(i=0; i < ni; i++)
238 {
239 double yd = (i+0.5)/ni-0.5;
240 double theta = (0.5-yd)*Pi;
241 double facteur=1./sin(theta);
242 for(j=0; j < nj; j++)
243 {
244 double xa = (j+0.5)/nj-0.5;
245 double phi = 2.*Pi*xa*facteur+Pi;
246 float th=float(theta);
247 float ff=float(phi);
248 if (phi<2*Pi && phi>= 0)
249 {
250 map(i,j) = (float)dobj_->PixValSph(th, ff);
251 }
252 }
253 }
254 FitsOutFile fout(filename);
255 fout.firstImageOnPrimaryHeader(true);
256 FITS_TArray<r_4> fta(map);
257 fta.Write(fout);
258
259
260}
261
262
263#ifdef __CXX_PRAGMA_TEMPLATES__
264#pragma define_template FITS_SphereThetaPhi<r_8>
265#pragma define_template FITS_SphereThetaPhi<r_4>
266//#pragma define_template FITS_SphereThetaPhi<int_4>
267#endif
268#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
269template class FITS_SphereThetaPhi<r_8>;
270template class FITS_SphereThetaPhi<r_4>;
271//template class FITS_SphereThetaPhi<int_4>;
272#endif
Note: See TracBrowser for help on using the repository browser.