source: Sophya/trunk/Cosmo/SimLSS/cmvgsm2cube.cc@ 4045

#include "sopnamsp.h"
#include "machdefs.h"
#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <string.h>
#include <unistd.h>
#include <vector>

#include "fioarr.h"
#include "skymap.h"

// 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)
{
cout
<<"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
<<endl;
}


int main(int narg,char *arg[])
{
 string ppfrn = "cmvgsm2cube.ppf";
 int nphi=0, ntheta=0;
 double dphi=-1., phi0=0.;
 double dtheta=-1., theta0=0.;

 if(narg<=1) {usage(); return -1;}

// Read arguments
 char c;
 while ((c = getopt(narg, arg, "ho:p:t:")) != -1) {
  switch (c) {
  case 'o' :
    ppfrn = optarg;
    break;
  case 'p' :
    sscanf(optarg,"%d,%lf,%lf",&nphi,&dphi,&phi0);
    break;
  case 't' :
    sscanf(optarg,"%d,%lf,%lf",&ntheta,&dtheta,&theta0);
    break;
  case 'h' :
  default:
    usage();
    return -1;
  }
 }

 if(nphi<=0 || dphi<=0. || ntheta<=0 || dtheta<=0.) {
   cout<<"bad value for transverse dim."<<endl;
   return -2;
 }
 cout<<"nphi="<<nphi<<" dphi="<<dphi<<" phi0="<<phi0<<endl;
 cout<<"ntheta="<<ntheta<<" dtheta="<<dtheta<<" theta0="<<theta0<<" deg"<<endl;

 vector<string> sphname;
 for(int i=optind;i<narg;i++) sphname.push_back(arg[i]);
 if(sphname.size()==0) {usage(); return -3;}
 cout<<"Number of spheres "<<sphname.size()<<endl;

 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;
   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;
}

/*
openppf cmvgsm2cube.ppf

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
objaoper cube sliceyz 255
objaoper cube slicexy 90

 */