source: Sophya/trunk/ArchTOIPipe/ProcWSophya/toi2map.cc@ 1984

Last change on this file since 1984 was 1809, checked in by cmv, 24 years ago
  • map2toi et toi2map avec nouvelle interface sys coor et unites.
  • les modifs map2toi (LocalMap) de EA sont introduites + certains bugs corriges.

cmv 3/12/01
File size: 7.4 KB
Line 
1// ArchTOIPipe (C) CEA/DAPNIA/SPP IN2P3/LAL
2// Eric Aubourg
3// Christophe Magneville
4// Reza Ansari
5// $Id: toi2map.cc,v 1.24 2001-12-03 20:38:02 cmv Exp $
6
7#include "machdefs.h"
8#include "toimanager.h"
9#include "pexceptions.h"
10#include "ctimer.h"
11#include "toi2map.h"
12// La valeur "Pi" doit etre celle de smathconst.h a cause du test sur theta
13#include "smathconst.h"
14
15////////////////////////////////////////////////////////////////////////
16TOI2Map::TOI2Map(PixelMap<r_8>* map,PixelMap<r_8>* wmap)
17 : mMap(map), mWMap(wmap), mWMapInternal(false)
18{
19 SetEquinox();
20 SetCoorIn();
21 SetCoorMap();
22 SetTestFlag();
23 SetTestMin();
24 SetTestMax();
25
26 if(mMap->NbPixels()<1) {
27 cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mMap "
28 <<mMap->NbPixels()<<endl;
29 throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
30 }
31 mMap->SetPixels(0.);
32 int nlat = mMap->SizeIndex();
33 string typmap = mMap->TypeOfMap();
34
35 if(mWMap==NULL) {
36 // We would need a cloning function in maps. $$TBD$$ $CHECK$
37 if (typmap == "LOCAL") {
38 mWMap = new LocalMap<r_8>(*(LocalMap<r_8>*)mMap);
39 } else if (typmap == "RING") {
40 mWMap = new SphereHEALPix<r_8>(nlat);
41 } else {
42 cout << "TOI2Map::TOI2Map() cannot handle map of type " << typmap << endl;
43 throw ParmError("TOI2Map::TOI2Map() - bad type of map");
44 }
45 mWMapInternal = true;
46 } else {
47 mWMapInternal = false;
48 if(nlat != mWMap->SizeIndex()) {
49 cout<<"TOI2Map::TOI2Map() Bad size for sphere mWMap, does not "
50 <<"correspond to mMap : "<<mMap->SizeIndex()<<", "<<mWMap->SizeIndex()<<endl;
51 if (typmap == "LOCAL") {
52 throw ParmError("TOI2Map::TOI2Map() - Different sizes for map and wmap");
53 } else if(typmap == "RING") {
54 ((SphereHEALPix<r_8> *)mWMap)->Resize(nlat);
55 cout<<"Resize have been done..."<<endl;
56 } else {
57 cout << "TOI2Map::TOI2Map() cannot handle map of type " << typmap << endl;
58 throw ParmError("TOI2Map::TOI2Map() - bad type of map");
59 }
60 }
61 }
62 if(mWMap->NbPixels()<1) {
63 cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mWMap "
64 <<mWMap->NbPixels()<<endl;
65 throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
66 }
67 mWMap->SetPixels(0);
68}
69
70TOI2Map::~TOI2Map()
71{
72 if(mWMap && mWMapInternal) delete mWMap;
73}
74
75////////////////////////////////////////////////////////////////////////
76void TOI2Map::Print(::ostream & os)
77{
78 os<<"TOI2Map::Print -- Map type " << mMap->TypeOfMap() << " SizeIndex = "<<mMap->SizeIndex()<<endl
79 <<" - Equinoxe="<<mActualYear<<endl
80 <<" - TypCoorIn: "<<mTypCoorIn<<" = "<<DecodeTypAstro(mTypCoorIn)<<endl
81 <<" - TypCoorMap: "<<mTypCoorMap<<" = "<<DecodeTypAstro(mTypCoorMap)<<endl
82 <<" - Tests: Flag("<<mTFlag<<") bad="<<mBadFlag
83 <<" / Value Min("<<mTMin<<")="<<mValMin
84 <<" , Max("<<mTMax<<")="<<mValMax<<endl;
85}
86
87////////////////////////////////////////////////////////////////////////
88void TOI2Map::init() {
89 cout << "TOI2Map::init" << endl;
90 declareInput("Coord1In"); // input index 0
91 declareInput("Coord2In"); // input index 1
92 declareInput("BoloIn"); // input index 2
93}
94
95////////////////////////////////////////////////////////////////////////
96// define SANS_BUFFER
97void TOI2Map::run()
98{
99long snb = getMinIn();
100long sne = getMaxIn();
101
102if(snb>sne) {
103 cout<<"TOI2Map::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
104 throw ParmError("TOI2Map::run() - Bad sample interval");
105}
106if(!checkInputTOIIndex(0) || !checkInputTOIIndex(1) || !checkInputTOIIndex(2)) {
107 cout<<"TOI2Map::run() - Input TOI (Coord1In or Coord2In or BoloIn) not connected! "<<endl;
108 throw ParmError("TOI2Map::run() Output TOI (Coord1In or Coord2In or BoloIn) not connected!");
109}
110if( !(mTypCoorIn&TypCoordEq || mTypCoorIn&TypCoordGal) ) {
111 cout<<"TOI2Map::run() - Input Coordinates not Eq or Gal! "<<endl;
112 throw ParmError("TOI2Map::run() - Input Coordinates not Eq or Gal!");
113}
114if( !(mTypCoorMap&TypCoordEq || mTypCoorMap&TypCoordGal) ) {
115 cout<<"TOI2Map::run() - Output Coordinates not Eq or Gal! "<<endl;
116 throw ParmError("TOI2Map::run() - Output Coordinates not Eq or Gal!");
117}
118
119//---------------------------------------------------------
120#define NFILL 25
121try {
122
123int ii;
124uint_4 mNSnFill=0, mNpixFill=0, NFill[NFILL], BadCoorRange=0;
125for(ii=0;ii<NFILL;ii++) NFill[ii]=0;
126double mjd = MJDfrYear(mActualYear);
127
128cout<<"TOI2Map::run() from "<<snb<<" to "<<sne;
129#ifndef SANS_BUFFER
130int bufsz = 100;
131uint_8* bfgbolo = new uint_8[bufsz];
132double* bbolo = new double[bufsz];
133double* bc1 = new double[bufsz];
134double* bc2 = new double[bufsz];
135int i0 = -1;
136cout<<" (getData() bufferized)"<<endl;
137#else
138cout<<" (getData() not bufferized)"<<endl;
139#endif
140
141// Remplissage des spheres
142for(int s=snb;s<=sne;s++) {
143 uint_8 fgbolo = 0;
144 double bolo,coord1,coord2;
145 // Equatoriales / Galactiques
146 // coord1,2 = alpha,delta / gLon,gLat
147
148#ifndef SANS_BUFFER
149 if(i0<0 || s<i0 || s>=i0+bufsz) {
150 i0 = s;
151 int nget = (sne-s+1<bufsz)? nget=sne-s+1: bufsz;
152 getData(0,i0,nget,bc1);
153 getData(1,i0,nget,bc2);
154 getData(2,i0,nget,bbolo,bfgbolo);
155 }
156 bolo = bbolo[s-i0]; fgbolo = bfgbolo[s-i0];
157 coord1 = bc1[s-i0]; coord2 = bc2[s-i0];
158#else
159 getData(2,s,bolo,fgbolo);
160 coord1 = getData(0,s);
161 coord2 = getData(1,s);
162#endif
163
164 // Comme il n'y a pas de toi en sortie, il faut dire
165 // aux processeur/toi que l'on a plus besoin des donnees.
166 if (s%100 == 0) wontNeedBefore(s-1);
167 if(mTFlag && fgbolo&mBadFlag) continue;
168 if(mTMin && bolo<mValMin) continue;
169 if(mTMax && bolo>mValMax) continue;
170
171 // sphere phi entre [0,2*Pi] en radians
172 // sphere theta entre [0,Pi] en radians
173 double phi=-1.;
174 CoordConvertToStd(mTypCoorIn,&coord1,&coord2);
175
176 if(mTypCoorIn&TypCoordEq && mTypCoorMap&TypCoordGal) { // Eq -> Gal
177 EqtoGal(mjd,coord1,coord2,&coord1,&coord2);
178 phi = coord1 * Pi/180.;
179 } else if(mTypCoorIn&TypCoordGal && mTypCoorMap&TypCoordEq) { // Gal -> Eq
180 GaltoEq(mjd,coord1,coord2,&coord1,&coord2);
181 phi = coord1 * Pi/12.;
182 } else if(mTypCoorMap&TypCoordGal) { // Gal -> Gal
183 phi = coord1 * Pi/180.;
184 } else if(mTypCoorMap&TypCoordEq) { // Eq -> Eq
185 phi = coord1 * Pi/12.;
186 }
187 ToCoLat(&coord2,TypUniteD);
188 double theta = coord2 * Pi/180.;
189
190 if(phi<0. | phi>=2*Pi || theta<0. || theta>Pi)
191 {BadCoorRange++; continue;}
192
193 int_4 ipix = mMap->PixIndexSph(theta,phi);
194 (*mMap)(ipix) += bolo;
195 ((*mWMap)(ipix)) += 1;
196 mNSnFill++;
197}
198
199 cout<<"TOI2Map::run(): Fin de boucle sur les sampleNum"<<endl;
200
201// Remplissage des spheres
202 for(int_4 i=0;i<mMap->NbPixels();i++) {
203 r_8 wf = (*mWMap)(i);
204 if(wf>0.) {mNpixFill++; (*mMap)(i) /= wf;}
205 int_4 nf = int_4(wf);
206 if(nf>=NFILL) nf=NFILL-1; NFill[nf]++;
207 }
208
209 cout<<"TOI2Map::run(): mNpixTot="<<mMap->NbPixels()
210 <<" mNpixFill="<<mNpixFill
211 <<" mNSnFill="<<mNSnFill<<endl
212 <<" --> FracSky="<<mNpixFill*100./(double)mMap->NbPixels()<<"%"
213 <<" NFill["<<NFILL<<"] ="<<endl;
214 for(ii=0;ii<NFILL;ii++) {cout<<NFill[ii]<<" "; if(ii%10==9) cout<<endl;}
215 cout<<endl;
216 cout<<" BadCoorRange="<<BadCoorRange<<endl;
217
218#ifndef SANS_BUFFER
219delete [] bbolo; delete [] bfgbolo;
220delete [] bc1; delete [] bc2;
221#endif
222
223//---------------------------------------------------------
224} catch (PException & exc) {
225 cout<<"TOI2Map: Catched Exception "<<(string)typeid(exc).name()
226 <<"\n .... Msg= "<<exc.Msg()<<endl;
227}
228
229return;
230}
Note: See TracBrowser for help on using the repository browser.