source: Sophya/trunk/Cosmo/RadioBeam/srcat2cube.cc@ 3964

Last change on this file since 3964 was 3829, checked in by ansari, 15 years ago

Adaptation programme srcat2cube.cc pour utilisation catalogue North20cm avec indice spectral, Reza 03/08/2010
File size: 10.4 KB
Line 
1/* ------------------------ Projet BAORadio --------------------
2 Programme de fabrication d'un cube 3D (angles,fre)
3 a partir du catalogue de source radio (NVSS)
4 R. Ansari , C. Magneville - Juin 2010
5
6 Usage: srccat2cube CatalogFitsName Out3DPPFName [Out2DMapName]
7--------------------------------------------------------------- */
8
9#include "sopnamsp.h"
10#include "machdefs.h"
11#include <math.h>
12#include <iostream>
13#include <typeinfo>
14
15#include "array.h"
16#include "histats.h"
17
18#include "swfitsdtable.h"
19#include "fitshdtable.h"
20
21#include "randr48.h"
22
23#include "xastropack.h" // Pour faire les conversions de coordonnees celestes
24
25#include "radutil.h"
26
27// Pour l'initialisation des modules
28#include "tarrinit.h"
29#include "histinit.h"
30#include "fiosinit.h"
31
32#include "timing.h"
33#include "ctimer.h"
34
35#include "cubedef.h"
36
37//----------------
38int nvssTocube(DataTable& nvss, TArray<r_4>& omap, TArray<r_4>& cube);
39int north20Tocube(DataTable& nor, TArray<r_4>& omap, TArray<r_4>& cube);
40
41//----------------------------------------------------------------------------
42//----------------------------------------------------------------------------
43int main(int narg, char* arg[])
44{
45 // Sophya modules initialization
46 TArrayInitiator _inia;
47 HiStatsInitiator _inih;
48 FitsIOServerInitiator _inif;
49 //------- AU LIEU DE ------> SophyaInit();
50
51 InitTim(); // Initializing the CPU timer
52 Timer tm("srcat2cube");
53
54 if (narg < 4) {
55 cout << "Usage: srccat2cube -nvss/-north20 CatalogFitsName Out3DPPFName [Out2DMapName]\n" << endl;
56 return 1;
57 }
58
59
60 // decodage arguments
61
62 string copt=arg[1];
63 string outname=arg[3];
64 string inname=arg[2];
65 int rc = 91;
66
67 cout << " ====== srccat2cube : Input catalog name= " << inname << " OutName=" << outname;
68 bool fginmap=true;
69 try {
70 DataTable cat;
71 cout << "srccat2cube[1]: reading source catalog from " << inname << endl;
72 {
73 FitsInOutFile fis(inname, FitsInOutFile::Fits_RO);
74 fis >> cat;
75 }
76 cout << cat ;
77 TArray<r_4> omap(NPhi,NTheta);
78 TArray<r_4> ocube(NPhi,NTheta,NFreq);
79 if (copt=="-nvss") nvssTocube(cat, omap, ocube);
80 else north20Tocube(cat, omap, ocube);
81
82 { // On sauve le cube de sortie
83 cout << " srccat2cube[7]: Saving output cube to -> " << outname << endl;
84 POutPersist poc(outname);
85 poc << ocube;
86 }
87
88 if (narg > 4) {
89 string ppfname = arg[4];
90 cout << " srccat2cube[8]: saving 2D source map to PPF file-> " << ppfname << endl;
91 POutPersist po(ppfname);
92 po << omap;
93 }
94 rc = 0;
95 }
96 catch (PThrowable& exc) {
97 cerr << " srccat2cube.cc catched Exception " << exc.Msg() << endl;
98 rc = 77;
99 }
100 catch (std::exception& sex) {
101 cerr << "\n srccat2cube.cc std::exception :"
102 << (string)typeid(sex).name() << "\n msg= "
103 << sex.what() << endl;
104 }
105 catch (...) {
106 cerr << " srccat2cube.cc catched unknown (...) exception " << endl;
107 rc = 78;
108 }
109
110 cout << ">>>> srccat2cube[9] ------- FIN ----------- Rc=" << rc << endl;
111 return rc;
112}
113
114
115
116/* -- Fonction -- */
117int nvssTocube(DataTable& nvss, TArray<r_4>& omap, TArray<r_4>& ocube)
118{
119 sa_size_t idxa = nvss.IndexNom("C_RAJ2000");
120 sa_size_t idxd = nvss.IndexNom("C_DEJ2000");
121 sa_size_t idxf = nvss.IndexNom("S1_4");
122 sa_size_t idxmajax = nvss.IndexNom("MajAxis");
123 sa_size_t idxminax = nvss.IndexNom("MinAxis");
124
125 cout << " NVSS catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf
126 << " MajAxis: " << idxmajax << " MajAxis: " << idxminax << endl;
127
128 double tet0 = Theta0Degre;
129 double phi0 = Phi0Degre;
130 double tetmax = tet0+ThetaSizeDegre;
131 double phimax = phi0+PhiSizeDegre;
132
133 cout << "srccat2cube/NVSS[2]: projecting sources to map ..." << endl;
134
135 sa_size_t srccnt=0;
136 sa_size_t extendedsrccnt=0;
137
138 double meanflx=0.;
139 double flxmin=9.e99;
140 double flxmax=-9.e99;
141
142 double dtet = ThetaSizeDegre/(double)NTheta;
143 double dphi = PhiSizeDegre/(double)NPhi;
144 double mpixsizarcmin = 0.5*(dtet+dphi)*60.;
145
146 for (sa_size_t n=0; n<nvss.NRows(); n++) {
147 r_8* pline=nvss.GetLineD(n);
148 double alpha=pline[idxa]; // alpha en degre
149 double delta=pline[idxd]; // delta en degre
150 double flx=pline[idxf]*1.e-3; // flux en Jy
151 double srcszarcmin=0.5*(pline[idxmajax]+pline[idxminax])/60.; // taille (extension de la source en arcmin
152 if (srcszarcmin<1.) srcszarcmin=1.;
153 double tet = 90.-delta;
154 double phi = alpha;
155 sa_size_t i = (phi-phi0)/dphi;
156 sa_size_t j = (tet-tet0)/dtet;
157 if ((i<0)||(i>=omap.SizeX())) continue;
158 if ((j<0)||(j>=omap.SizeY())) continue;
159 if (srcszarcmin<(0.5*mpixsizarcmin)) { // Toute l'energie dans un seul pixel
160 omap(i,j) += flx;
161 }
162 else { // on repartit l'energie de la source dans plusieurs pixels
163 extendedsrccnt++;
164 for(int bi=-1;bi<=1;bi++) {
165 for(sa_size_t bj=-1; bj<=1; bj++) {
166 sa_size_t ii = (phi-phi0+bi*srcszarcmin/60.)/dphi;
167 sa_size_t jj = (tet-tet0+bj*srcszarcmin/60.)/dtet;
168 if ((ii<0)||(ii>=omap.SizeX())) continue;
169 if ((jj<0)||(jj>=omap.SizeY())) continue;
170 if ((bi==0)&&(bj==0)) omap(ii,jj) += flx*0.3;
171 else omap(ii,jj) += flx*0.7/8.;
172 }
173 }
174 }
175 srccnt++; meanflx+=flx;
176 if (flx<flxmin) flxmin=flx;
177 if (flx>flxmax) flxmax=flx;
178 }
179
180 cout << "srccat2cube/NVSS[3]: Output rectangular map computed " << endl;
181 meanflx /= (double)srccnt;
182 cout << " SrcCount in map: " << srccnt << " extended=" << extendedsrccnt
183 << " -> meanFlx=" << meanflx << " min=" << flxmin
184 << " max=" << flxmax << " Jy" << endl;
185
186 double mean, sigma;
187 r_4 minjy, maxjy;
188 omap.MinMax(minjy, maxjy);
189 MeanSigma(omap, mean, sigma);
190 cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky ; Sizes:" << endl;
191 omap.Show();
192
193 H21Conversions conv;
194 conv.setRedshift(0.);
195 conv.setOmegaPixDeg2(dphi*dtet);
196 cout << "srccat2cube/NVSS[4] H21Conversions, OmegaPix=" << conv.getOmegaPix() << " srad"
197 << " toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
198 omap *= (r_4)conv.toKelvin(1.);
199 MeanSigma(omap, mean, sigma);
200 r_4 minT, maxT;
201 omap.MinMax(minT, maxT);
202 cout << " NVSS/ After conversion : Mean=" << mean << " Sigma=" << sigma << " Min=" << minT
203 << " Max=" << maxT << " Kelvin " << endl;
204
205 double infreq = 1420.; // frequence de reference du flux des sources
206 double freq0 = Freq0MHz; // Freq0 du cube de sortie
207 double dfreq = FreqSizeMHz/(double)NFreq;
208
209 ThSDR48RandGen rg;
210 for (sa_size_t j=0; j<ocube.SizeY(); j++) {
211 for (sa_size_t i=0; i<ocube.SizeX(); i++) {
212 double freqexpo = rg.Gaussian(sigPLidxSrc,PLidxSrc);
213 for (sa_size_t k=0; k<ocube.SizeZ(); k++) {
214 double rapfreq = pow((freq0+k*dfreq)/infreq, freqexpo);
215 ocube(i,j,k) = AmpPL1*omap(i,j)*rapfreq;
216 }
217 }
218 }
219 cout << "srccat2cube/NVSS[5] data cube created from sources " << endl;
220 ocube.MinMax(minT, maxT);
221 MeanSigma(ocube, mean, sigma);
222 cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
223
224 return srccnt;
225}
226
227/* -- Fonction -- */
228int north20Tocube(DataTable& nor, TArray<r_4>& omap, TArray<r_4>& ocube)
229{
230
231 sa_size_t idxa = nor.IndexNom("ra");
232 sa_size_t idxd = nor.IndexNom("dec");
233 sa_size_t idxf = nor.IndexNom("flux");
234 sa_size_t idxslo = nor.IndexNom("SpLO");
235 sa_size_t idxshi = nor.IndexNom("SpHI");
236
237 cout << " North20cm catalog ... Index Alpha: " << idxa << " Delta: " << idxd << " Flux: " << idxf
238 << " SpLO: " << idxslo << " SpHI: " << idxshi << endl;
239
240 double tet0 = Theta0Degre;
241 double phi0 = Phi0Degre;
242 double tetmax = tet0+ThetaSizeDegre;
243 double phimax = phi0+PhiSizeDegre;
244
245 cout << "srccat2cube/North20[2]: projecting sources to map ..." << endl;
246
247 sa_size_t srccnt=0;
248 sa_size_t lowoksrccnt=0;
249
250 double meanflx=0.;
251 double flxmin=9.e99;
252 double flxmax=-9.e99;
253
254 double dtet = ThetaSizeDegre/(double)NTheta;
255 double dphi = PhiSizeDegre/(double)NPhi;
256
257 double infreq = 1420.; // frequence de reference du flux des sources
258 double freq0 = Freq0MHz; // Freq0 du cube de sortie
259 double dfreq = FreqSizeMHz/(double)NFreq;
260
261 for (sa_size_t n=0; n<nor.NRows(); n++) {
262 r_8* pline=nor.GetLineD(n);
263 double alpha=pline[idxa]; // alpha en degre
264 double delta=pline[idxd]; // delta en degre
265 double flx=pline[idxf]*1.e-3; // flux en Jy
266 double tet = 90.-delta;
267 double phi = alpha*360./24.;
268 sa_size_t i = (phi-phi0)/dphi;
269 sa_size_t j = (tet-tet0)/dtet;
270 if ((i<0)||(i>=omap.SizeX())) continue;
271 if ((j<0)||(j>=omap.SizeY())) continue;
272 omap(i,j) += flx;
273 srccnt++; meanflx+=flx;
274 if (flx<flxmin) flxmin=flx;
275 if (flx>flxmax) flxmax=flx;
276 double slo=pline[idxslo];
277 if (slo<9.) { // source detected at 80 cm
278 lowoksrccnt++;
279 }
280 else slo=5.;
281 for (sa_size_t k=0; k<ocube.SizeZ(); k++) {
282 double rapfreq = pow((freq0+k*dfreq)/infreq, slo);
283 ocube(i,j,k) += flx*rapfreq;
284 }
285 }
286
287 cout << "srccat2cube/North20[3]: Output rectangular map computed " << endl;
288 meanflx /= (double)srccnt;
289 cout << " SrcCount in map: " << srccnt << " SpLowOK=" << lowoksrccnt
290 << " -> meanFlx=" << meanflx << " min=" << flxmin
291 << " max=" << flxmax << " Jy" << endl;
292
293 double mean, sigma;
294 r_4 minjy, maxjy;
295 omap.MinMax(minjy, maxjy);
296 MeanSigma(omap, mean, sigma);
297 cout << " Src Map : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;
298
299 ocube.MinMax(minjy, maxjy);
300 MeanSigma(ocube, mean, sigma);
301 cout << " Cube : Mean=" << mean << " Sigma=" << sigma << " Min=" << minjy << " Max=" << maxjy << " Jansky" << endl;
302
303 H21Conversions conv;
304 conv.setRedshift(0.);
305 conv.setOmegaPixDeg2(dphi*dtet);
306 cout << "srccat2cube/North20[4] H21Conversions, OmegaPix=" << conv.getOmegaPix() << " srad"
307 << " @1400MHz toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
308
309 // Jansky to Kelvin conversion
310 for (sa_size_t k=0; k<ocube.SizeZ(); k++) {
311 conv.setFrequency(freq0+k*dfreq);
312 // cout << " DBG* Freq= " << freq0+k*dfreq << " -> toKelvin(1 Jy)= " << conv.toKelvin(1.) << endl;
313 ocube(Range::all(), Range::all(), Range(k)) *= (r_4)conv.toKelvin(1.);
314 }
315 cout << "srccat2cube/North20[5] data cube in Kelvin computed " << endl;
316
317 r_4 minT, maxT;
318 ocube.MinMax(minT, maxT);
319 MeanSigma(ocube, mean, sigma);
320 cout << "... Mean=" << mean << " Sigma=" << sigma << " Min=" << minT << " Max=" << maxT << " Kelvin" << endl;
321
322 return srccnt;
323}
Note: See TracBrowser for help on using the repository browser.