Changeset 3983 in Sophya for trunk/Cosmo/SimLSS/cmvgsm2cube.cc

Timestamp:

May 4, 2011, 7:13:16 PM (14 years ago)

Author:

cmv

Message:

amelioration des reconstructions et fits du GSM, cmv 4/5/2011

File:

• trunk/Cosmo/SimLSS/cmvgsm2cube.cc (modified) (4 diffs)

trunk/Cosmo/SimLSS/cmvgsm2cube.cc

@@ -11 +11 @@
 #include "skymap.h"
 
-// set fname = ( `cat gsm820948.list | awk '{print $3}' | sed 's/data/ppf/'` )
-// cmvgsm2cube -o gsm820948_cube.ppf -p 360,0.1666667,160. -t 360,0.1666667,30. $fname
+// set rname = gsm820948
+// set fname = ( `cat ${rname}.list | awk '{print $3}' | sed 's/data/ppf/'` )
+// cmvgsm2cube -o ${rname}_cube.ppf -p 360,0.1666667,160. -t 360,0.1666667,30. $fname
 
 void usage(void)
@@ -19 +20 @@
 <<"cmvgsm2cube [options] liste_des_spheres_healpix"<<endl
 <<" -o outname.ppf : ppf cube output file name (def=cmvgsm2cube.ppf)"<<endl
-<<" -p nphi,dphi,phi0 : phi0 -> phi0+nphi*dphi (degres)"<<endl
-<<" -t ntheta,dtheta,theta0 : theta0 -> theta0+ntheta*dtheta (degres)"<<endl
+<<" -p nphi,dphi,phi0 : phi0 -> phi0 + nphi*dphi (degres)"<<endl
+<<" -t ntheta,dtheta,theta0 : theta0 -> theta0 + ntheta*dtheta (degres)"<<endl
 <<endl;
 }
@@ -66 +67 @@
  cout<<"Number of spheres "<<sphname.size()<<endl;
 
- TArray<r_4> Cube(sphname.size(),nphi,ntheta);
- Cube = 0.;
+ TVector<r_4> Theta(ntheta);
+ int nthgood = 0;
+ float thmin = 1.e30, thmax = -1.e30;
+ for(int it=0;it<ntheta;it++) {
+   Theta(it) = -9999.;
+   double t = theta0 + it*dtheta; t *= M_PI/180.;
+   if(t<0. || t>M_PI) continue;  // [0,Pi]
+   Theta(it) = t;
+   nthgood += 1;
+   if(t<thmin) thmin = t;
+   if(t>thmax) thmax = t;
+ }
+ cout<<"Theta: ngood="<<nthgood<<" ["<<thmin*180/M_PI<<","<<thmax*180/M_PI
+     <<"] middle="<<0.5*(thmin+thmax)*180/M_PI<<endl;
+ 
+ TVector<r_4> Phi(nphi);
+ float phmin = 1.e30, phmax = -1.e30;
+ for(int ip=0;ip<nphi;ip++) {
+   double p = phi0 + ip*dphi; p *= M_PI/180.;
+   while(p<0.) p += 2.*M_PI;
+   while(p>=2.*M_PI) p -= 2.*M_PI;
+   Phi(ip) = p;
+   if(p<phmin) phmin = p;
+   if(p>phmax) phmax = p;
+ }
+ cout<<"Phi: ["<<phmin*180/M_PI<<","<<phmax*180/M_PI
+     <<"] middle="<<0.5*(phmin+phmax)*180/M_PI<<endl;
+
+
+ TVector<r_4> Freq(sphname.size()); Freq = -1.;
+ TArray<r_4> Cube(sphname.size(),nphi,ntheta); Cube = 0.;
  for(unsigned int is=0;is<sphname.size();is++) {
-
    cout<<"... "<<sphname[is]<<endl;
    PInPersist pis(sphname[is]);
    SphereHEALPix<r_4> sph;
    pis>>PPFNameTag("sph")>>sph;
-
-   for(int it=0;it<ntheta;it++) {
-     double t = theta0 + it*dtheta; t *= M_PI/180.;
-     if(t<0. || t>M_PI) continue;
-   for(int ip=0;ip<nphi;ip++) {
-     double p = phi0 + ip*dphi; p *= M_PI/180.;
-     while(p<0.) p += 2.*M_PI;
-     while(p>=2.*M_PI) p -= 2.*M_PI;
-     int_4 ksph = sph.PixIndexSph(t,p);
-     Cube(is,ip,it) = sph(ksph);
+   Freq(is) = (r_4)sph.Info()["FMHz"];
+   for(int it=0;it<Theta.Size();it++) {
+     double t = Theta(it);
+     if(t<0.) continue;
+     for(int ip=0;ip<Phi.Size();ip++) {
+       double p = Phi(ip);
+       int_4 ksph = sph.PixIndexSph(t,p);
+       Cube(is,ip,it) = sph(ksph);
+     }
    }
-   }
-
  }
  
  POutPersist pos(ppfrn);
  pos<<PPFNameTag("cube")<<Cube;
+ pos<<PPFNameTag("theta")<<Theta;
+ pos<<PPFNameTag("phi")<<Phi;
+ pos<<PPFNameTag("freq")<<Freq;
 
  return 0;
@@ -100 +129 @@
 print cube
 
+n/plot phi.val*180./M_PI%n ! ! "cpts"
+n/plot theta.val*180./M_PI%n ! ! "cpts"
+disp freq
+
 objaoper cube sliceyz 0
 objaoper cube sliceyz 128