#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 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

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;

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

 }
 
 POutPersist pos(ppfrn);
 pos<<PPFNameTag("cube")<<Cube;

 return 0;
}

/*
openppf cmvgsm2cube.ppf

print cube

objaoper cube sliceyz 0
objaoper cube sliceyz 128
objaoper cube sliceyz 255
objaoper cube slicexy 90

 */