/* 


This is an example for MPE installation, replace with appropriate paths.

setenv HEALPIX_PATH /afs/ipp-garching.mpg.de/home/a/aws/Healpix/Healpix_2.20a_sles11_olga2/Healpix_2.20a/src/cxx/generic_gcc

GNU compiler:
g++ healpix_skymap_convert.cc -I$HEALPIX_PATH/include -L$HEALPIX_PATH/lib -lhealpix_cxx -lcxxsupport -lpsht  -lfftpack  -lc_utils  -lcfitsio -fopenmp -o healpix_skymap_convert

Intel compiler:
 icpc healpix_skymap_convert.cc -I$HEALPIX_PATH/include -L$HEALPIX_PATH/lib -lhealpix_cxx -lcxxsupport -lpsht -lfftpack -lc_utils -lcfitsio -openmp  -o healpix_skymap_convert


NB OpenMP is not actually used here, but is required when HealPix libraries are built with OpenMP, as in the default makefiles.

*/

/*
   converts from GALPROP/GARDIAN Skymap class HealPix to standard HealPix with arrays in columns instead of pixels.
   both original and energy-integrated values are output to separate files. Gaussian-smoothed energy bins and 
   energy-integrated values are   also output.
   If two input files are present, the first is divided by the second, with interpolation of the second in energy, according to the energies
   specified in the second extension (HDU 3) which is ENERGIES in the case of GALPROP and may be EBOUNDS for counts maps.

   Requires HealPix 2.20a available from http://healpix.jpl.nasa.gov/

   usage: ./healpix_skymap_convert filename            FWHM prefix
      or  ./healpix_skymap_convert filename1 filename2 FWHM prefix

                  filename input Skymap-format file (may be .gz, in which case the output file is also .gz)
                  FWHM     Gaussian full width half-max in degrees
                  prefix   prefix prepended to output files

                  filename1 filename2: same as filename but values in filename1 are divided by values in filename2
                                       for example counts/exposure to give intensity


   restrictions: must be run in same directory as input files since output filenames are based on input filenames
   
   Note Aladin will read .gz files.

   Version 1.2

   A version factored into subroutines will follow.
A. W. Strong, MPE, 31 Oct 2011
History
20111031 Version 1.0
20111102 map2alm_iter used for more accuracy 
20111212 column names: blank replaced by underscore
20111216 executable given name
20111219 generalized build instructions
20120126 Version 1.1: also output smoothed energy bins not energy integrated
20120130 Version 1.2: divides by a second file if present in argument list. Files need not be same order, second file is converted to order of first file before division. Energies taken from third extension for interpolation.
*/

using namespace std;
#include<iostream>

#include "healpix_map.h"
#include "healpix_map_fitsio.h"

#include "alm.h"
#include "alm_healpix_tools.h" //2.20a
#include "alm_fitsio.h"
#include "alm_powspec_tools.h"
#include "powspec.h"

#include "xcomplex.h"
#include "arr.h"
#include "fitshandle.h"

///////////////////////////////////////////////////////////////////////////////////

             
int main(int argc, char*argv[]){

  cout<<endl<< "=== healpix_skymap_convert v1.2 ==="<<endl<<endl;

 int debug=0;

 if(debug==1) cout<<"argc="<<argc<<endl;

 if(argc!=4 && argc!=5) {cout<<" incorrect number of arguments ! usage: "                      <<endl
                             <<"     ./healpix_skymap_convert filename            FWHM prefix" <<endl
			     <<" or  ./healpix_skymap_convert filename1 filename2 FWHM prefix" <<endl;
                exit(0);} //AWS20120127

//cout<<"argv= ";  for(int i=0;i<argc;i++){cout<<argv[i]<<" ";} cout<<endl;
//cout<<"argv[0]="<<argv[0]<<endl;// name of executable

 if(argc==4) //AWS20120127
 {
  cout<<"input file: "    <<argv[1]<<endl;
  cout<<"smoothing FWHM= "<<argv[2]<<endl;
  cout<<"prefix =        "<<argv[3]<<endl;
 }

 if(argc==5) //AWS20120127
 {
  cout<<"input file 1: "  <<argv[1]<<endl;
  cout<<"input file 2: "  <<argv[2]<<endl;
  cout<<"smoothing FWHM= "<<argv[3]<<endl;
  cout<<"prefix =        "<<argv[4]<<endl;
 }


 int         dividefiles;
 if(argc==4) dividefiles=0;
 if(argc==5) dividefiles=1;

 if(dividefiles==1) cout<<"dividing input file 1 by input file 2"<<endl;

  
  
    string infile,outfile;
    string infile2;        //AWS20120127
    string prefix;
    



  
                infile =string(argv[1]);
  if(argc==5)   infile2=string(argv[2]);//AWS20120127

  if(argc==4)    prefix=string(argv[3]);//AWS20120127
  if(argc==5)    prefix=string(argv[4]);//AWS20120127

  // use prefix not suffix to preserve .fits if present in input file

  string outfile_energies           =         prefix+"_healpix_standard_"                    +infile; 
  string outfile_energies_smoothed  =         prefix+"_healpix_standard_smoothed_"           +infile;//AWS20120126
  string outfile_integrated_energies=         prefix+"_healpix_standard_integrated_"         +infile;
  string  infile_integrated_energies=         prefix+"_healpix_standard_integrated_"         +infile;
  string outfile_integrated_energies_smoothed=prefix+"_healpix_standard_integrated_smoothed_"+infile;
  string outfile_Alm                         =prefix+"_healpix_Alm_"                         +infile;
  string outfile_2                           =prefix+"_healpix_standard_"                    +infile2;//AWS20120126 
  string outfile_2_rebinned                  =prefix+"_healpix_standard_rebinned_"           +infile2;//AWS20120126 

  double fwhm_deg=1.0;
        
  if(argc==4)sscanf(argv[2],"%lf",&fwhm_deg ); // l needed for double
  if(argc==5)sscanf(argv[3],"%lf",&fwhm_deg ); // l needed for double

  if (debug==1)  cout<<"smoothing FWHM= "<<fwhm_deg<<" deg"<<endl;

  PDT datatype = PLANCK_FLOAT64; // HealPix defines it this way
      datatype = PLANCK_FLOAT32;

  /////////////////////////////////////////////////////////
  // conversion from galprop healpix to standard healpix //
  /////////////////////////////////////////////////////////

    cout<<" ==== conversion from galprop healpix to energy column format"<<endl;

  // use healpix io rather than cfitsio even at low level to simplify dependence
  // see Healpix_2.20a/src/cxx/cxxsupport/fitshandle.h
  
  //==================== read input file (=file 1 if two are specified on argument list)

  cout<<"reading "<<infile<<endl;
  fitshandle inputhandle;
  inputhandle.open(infile);
  int hdu=2;
  inputhandle.goto_hdu (hdu);
  int colnum=1;
  
  int nenergy = 0;
  int LASTPIX = 0;
  inputhandle.get_key(string("LASTPIX"),LASTPIX); // starting from 0
  inputhandle.get_key(string("NBRBINS"),nenergy); // assuming as Skymap class 

  

  int  npix    = LASTPIX+1;
  long ndata   = nenergy*npix;

  cout<<"npix="<<npix<<" nenergy="<<nenergy<<endl;

  arr<double>data; // HealPix array class
  data.alloc(ndata);

  inputhandle.read_entire_column(colnum, data); // reads complete block of pixels for all energies



  int i,j,k;
  // for( k=0;k<ndata;k++)  cout<<data[k]<<endl;

 
  // data are read energy first then pixel

  double **val;

  val=new double*[npix];


  for( i=0;i<npix;i++)val[i]=new double[nenergy];

  if(debug==1) cout<<"allocated val"<<endl;

  k=0;
  for( i=0;i<npix;i++)  for( j=0;j<nenergy;j++,k++)
    { val[i][j]=data[k];
      //                cout<<i<<" "<<j<<" "<< val[i][j]<<endl;
    }

  // print all pixels for each energy
  if(0)
    {
  for( j=0;j<nenergy;j++)  for( i=0;i<npix;i++)
    {
         cout<<"to compare with fv: energy:"<<j+1<<" pixel:"<<i+1<<" value:"<< val[i][j]<<endl;
    }
    }


  // read the energies assuming in HDU 3 as for galprop healpix (header name ENERGIES) or EBOUNDS (then two columns)
  arr<double> energies_input_1;
  energies_input_1.alloc(nenergy);
  hdu=3;
  inputhandle.goto_hdu(hdu);
  colnum=1;
  inputhandle.read_entire_column(colnum,energies_input_1);
  for( j=0;j<nenergy;j++) {cout<<"energies_input_1["<<j<<"]="<<energies_input_1[j]<<endl;}


  inputhandle.close();

  //==================== read input file 2 if required

  arr<double> energies_input_2;
  int nenergy2 = 0;

 if(dividefiles==1)
 {
  cout<<"reading input file 2 ="<<infile2<<endl;
  fitshandle inputhandle2;
  inputhandle2.open(infile2);
      hdu=2;
  inputhandle2.goto_hdu (hdu);
      colnum=1;
  
  // file 2 may have different dimensions from file 1 !


  int LASTPIX2 = 0;
  inputhandle2.get_key(string("LASTPIX"),LASTPIX2); // starting from 0
  inputhandle2.get_key(string("NBRBINS"),nenergy2); // assuming as Skymap class 

  

  int  npix2   = LASTPIX2+1;
       ndata   = nenergy2*npix2;

  cout<<"npix2="<<npix2<<" nenergy2="<<nenergy2<<endl;

  if(npix2   !=npix   ) cout<<"warning: number of   pixels for file 1 differs from file2: converting order of file 2 "<<endl;
  if(nenergy2!=nenergy) cout<<"warning: number of energies for file 1 differs from file2: interpolating "             <<endl;
  
  data.alloc(ndata);

  inputhandle2.read_entire_column(colnum, data); // reads complete block of pixels for all energies

  energies_input_2.alloc(nenergy2);
  hdu=3;
  inputhandle2.goto_hdu(hdu);
  colnum=1;
  inputhandle2.read_entire_column(colnum,energies_input_2);
  for( j=0;j<nenergy;j++) {cout<<"energies_input_2["<<j<<"]="<<energies_input_2[j]<<endl;}



  inputhandle2.close();
  

 
  // data are read energy first then pixel

  double **val2;

  val2=new double*[npix2];


  for( i=0;i<npix2;i++)val2[i]=new double[nenergy2];

  if(debug==1) cout<<"allocated val2"<<endl;

  k=0;
  for( i=0;i<npix2;i++)  for( j=0;j<nenergy2;j++,k++)
    { val2[i][j]=data[k];
      //                cout<<i<<" "<<j<<" "<< val[i][j]<<endl;
    }

  // print all pixels for each energy
  if(0)
    {
  for( j=0;j<nenergy2;j++)  for( i=0;i<npix2;i++)
    {
         cout<<"to compare with fv: energy:"<<j+1<<" pixel:"<<i+1<<" value:"<< val2[i][j]<<endl;
    }
    }

  




  // put map 2 into standard HealPix

  int order=0; //  determine it from input since it is not present as a FITS keyword

  // npix=12 Nside^2 where Nside=2^order

  order=log(sqrt(npix2/12.))/log(2.)+.1;
  cout<<"HealPix order of map 2  computed from npix2 = "<<order<<endl;
  
  Healpix_Map<double> map2_standard(order,RING);

  data.alloc(npix2);
   cout<<"npix2="<<npix2<<" data.size() before Set = "<<data.size()<<endl;

 outfile=outfile_2;
  cout<<"writing to "<<outfile<<endl;

  //  write_Healpix_map_to_fits(outfile,map_standard, datatype); works but 1 column only,  want column per energy

  fitshandle out;
  out.create("!"+outfile);  // ! to overwrite
  arr<string> colname;

 
  colname.alloc(nenergy2);
  char energy_string[4];
  for( j=0;j<nenergy2;j++)
  { 
   sprintf(energy_string,"%i",j+1);
   colname[j]="energy_"+string(energy_string); //AWS20111212 FITS standard deprecates blanks in column names
   // cout<<"column name = "<<colname[j]<<endl;
  }


  prepare_Healpix_fitsmap(out,map2_standard, datatype, colname);//healpix_map_fitsio.h 
  out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name

  for( j=0;j<nenergy2;j++)
    {
       
  for( i=0;i<npix2;i++)
    {
      data[i]=val2[i][j];
    }
  //  cout<<"npix="<<npix<<" before write data.size() = "<<data.size()<<endl;


    colnum=j+1;
    out.write_column(colnum,data); // fitshandle.h
    }

  out.close();






  // use the standard healpix fits as input to convert order of map 2 and output it

  infile=outfile_2;
  fitshandle in;
  in.open(infile);
  in.goto_hdu(2);

  
  int ncolnum=in.ncols();
  cout<<" number of columns ="<<ncolnum<<endl;

  Healpix_Map<double> map;

  
 // npix=12 Nside^2 where Nside=2^order

  order=log(sqrt(npix/12.))/log(2.)+.1;
  cout<<"HealPix order of input file 1 computed from npix = "<<order<<endl;

  // input map2 to be converted to order of input file 1
  Healpix_Map<double> map_RING_out(order,RING);

  outfile= outfile_2_rebinned;
  cout<<"writing to "<<outfile<<endl;

  out.create("!"+outfile); // ! to overwrite
  prepare_Healpix_fitsmap(out,map_RING_out, datatype, colname); // colname is array, already filled above

  out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name
  debug=1;

  for (colnum=1;colnum<=ncolnum;colnum++)

  {
    if(debug==1) cout<< "reading input file 2 as converted to standard healpix: "<<infile<<",  column number: "<< colnum;

  read_Healpix_map_from_fits(infile,map,colnum,2);// HDU 2 (1=primary)

  if(debug==1)
  {
   cout<<" Npix  = "<<map.Npix(); 
   cout<<" Nside = "<<map.Nside();
   cout<<" Order = "<<map.Order();
   cout<<" Scheme= "<<map.Scheme()<<endl; // 0 = RING, 1 = NESTED



  }

  map_RING_out.Import(map); // converts map order to that of input file 1

  out.write_column(colnum,map_RING_out.Map()); // fitshandle.h

 }// colnum

  out.close();

 }// if dividefiles==1











  //=================== now put input map 1 into standard healpix columns
  
  //  Healpix_Map<double> map_standard(map.Order(),RING);
  int order=0; //  determine it from input since it is not present as a FITS keyword

  // npix=12 Nside^2 where Nside=2^order

  order=log(sqrt(npix/12.))/log(2.)+.1;
  cout<<"HealPix order computed from npix = "<<order<<endl;
  
  Healpix_Map<double> map_standard(order,RING);

  data.alloc(npix);
  if(debug==1)  cout<<"npix="<<npix<<" data.size() before Set = "<<data.size()<<endl;

  j=0;
  for( i=0;i<npix;i++)
    {
      data[i]=val[i][j];
    }

  //  map_standard.Set(data,RING); // data.size()=0 after this!

  outfile=outfile_energies;
  cout<<"writing to "<<outfile<<endl;

  //  write_Healpix_map_to_fits(outfile,map_standard, datatype); works but 1 column only,  want column per energy

  fitshandle out;
  out.create("!"+outfile);  // ! to overwrite
  arr<string> colname;

 
  colname.alloc(nenergy);
  char energy_string[4];
  for( j=0;j<nenergy;j++)
  { 
   sprintf(energy_string,"%i",j+1);
   colname[j]="energy_"+string(energy_string); //AWS20111212 FITS standard deprecates blanks in column names
   // cout<<"column name = "<<colname[j]<<endl;
  }


  prepare_Healpix_fitsmap(out,map_standard, datatype, colname);//healpix_map_fitsio.h 
  out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name

  for( j=0;j<nenergy;j++)
    {
       
  for( i=0;i<npix;i++)
    {
      data[i]=val[i][j];
    }
  //  cout<<"npix="<<npix<<" before write data.size() = "<<data.size()<<endl;



    colnum=j+1;



  if(dividefiles==1)
  {

  infile=outfile_2_rebinned;

  int colnum2;
  colnum2=colnum; // first guess, then need to find the corresponding energy

  double energies_input_1_mean;
  // eventually use second column of EBOUNDS to avoid assumption about ranges
  if(j<=nenergy-2)energies_input_1_mean=sqrt(energies_input_1[j]*energies_input_1[j+1]);
  if(j==nenergy-1)energies_input_1_mean=sqrt(energies_input_1[j]*
                                                                 energies_input_1[j]
                                                                *energies_input_1[1]/energies_input_1[0]);//for log energy scale

  cout<<"mean energy in file 1 for interpolating in file2 = "<< energies_input_1_mean<<endl;


  int jj2;
  for (jj2=0;jj2<nenergy2;jj2++)
  {
    cout<< "search in file 2 for energy "<<energies_input_1_mean<<": "<<energies_input_2[jj2]<<endl;
    if(energies_input_2[jj2]>=energies_input_1_mean) break;
  }

  if(jj2>nenergy2-1) jj2=nenergy2-1;

  cout<< "energy search for "<<energies_input_1_mean<<": first higher energy="<<energies_input_2[jj2]<<endl;

  colnum2=jj2+1;

  cout<< "reading input file 2 as converted to standard healpix and rebinned to order of file 1 for division: "<<infile<<",  column number high: "<< colnum2<<endl;

  Healpix_Map<double>map2_low;
  Healpix_Map<double>map2_high;
  Healpix_Map<double>map2;

  read_Healpix_map_from_fits(infile,map2_high,colnum2,2);// HDU 2 (1=primary)

  map2=map2_high;

  colnum2=jj2;
  if(colnum2<1)colnum2=1;
  cout<< "reading input file 2 as converted to standard healpix and rebinned to order of file 1 for division: "<<infile<<",  column number low: "<< colnum2<<endl;
  read_Healpix_map_from_fits(infile,map2_low,colnum2,2);// HDU 2 (1=primary)

  map2 = map2_low;



  for( i=0;i<npix;i++)
   
    map2[i]=                  map2_low[i]
            + (map2_high[i] - map2_low[i])*(energies_input_1_mean - energies_input_2[jj2-1])
                                          /(energies_input_2[jj2] - energies_input_2[jj2-1] );

 


  debug=1;
  if(debug==1)
  {
   cout<<" Npix  = "<<map2.Npix(); 
   cout<<" Nside = "<<map2.Nside();
   cout<<" Order = "<<map2.Order();
   cout<<" Scheme= "<<map2.Scheme()<<endl; // 0 = RING, 1 = NESTED
  }

  // the maps now have the same order so this will work
  for( i=0;i<npix;i++)
    {
      data[i] /= map2[i];
       val[i][j]=data[i]; //used for energy integrated maps
    }


  }//dividefiles==1






    out.write_column(colnum,data); // fitshandle.h
    }

  out.close();



  //////////////////////////////
  // energy integrated fluxes //
  //////////////////////////////

    cout<<"==== conversion from galprop healpix to energy column format,  values integrated over energy"<<endl;

  
   
  outfile=outfile_integrated_energies;
  cout<<"writing to "<<outfile<<endl;

  
  out.create("!"+outfile);// ! to overwrite

  

  prepare_Healpix_fitsmap(out,map_standard, datatype, colname);//healpix_map_fitsio.h 

  out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name

  for( j=0;j<nenergy;j++)
    {
       
  for( i=0;i<npix;i++)
    {
      data[i]=0.;
      for (int jj=j;jj<nenergy;jj++) data[i]+=val[i][jj];
    }
  


    colnum=j+1;
    out.write_column(colnum,data); // fitshandle.h
    }

  out.close();



  //////////////////////////////////////////////////////////////////////////
  // energy bins  with smoothing, using already converted format as input // //AWS20120126
  //////////////////////////////////////////////////////////////////////////

 {

  cout<<"==== smoothing  values for energy bins       "<<endl;


  cout<<" Gaussian FWHM = "<<fwhm_deg<< " deg"<<endl;

 

  Healpix_Map<double> map;
  Healpix_Map<double> map_RING_out(order,RING);
  

  outfile= outfile_energies_smoothed;
  cout<<"writing to "<<outfile<<endl;

  out.create("!"+outfile); // ! to overwrite
  prepare_Healpix_fitsmap(out,map_RING_out, datatype, colname);

  out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name

  infile=outfile_energies;
  int ncolnum=0;

  cout<< "reading "<<infile<<endl;

  
  fitshandle(in);
  in.open(infile);
  in.goto_hdu(2);

  
  ncolnum=in.ncols();
  cout<<" number of columns ="<<ncolnum<<endl;


  for (colnum=1;colnum<=ncolnum;colnum++)

  {
    if(debug==1) cout<< "reading "<<infile<<",  column number: "<< colnum;

  read_Healpix_map_from_fits(infile,map,colnum,2);// HDU 2 (1=primary)

  if(debug==1)
  {
   cout<<" Npix  = "<<map.Npix(); 
   cout<<" Nside = "<<map.Nside();
   cout<<" Order = "<<map.Order();
   cout<<" Scheme= "<<map.Scheme()<<endl; // 0 = RING, 1 = NESTED
  }

  /*
  // convert to RING (WMAP is NESTED, galprop is RING already)
  Healpix_Map<double> map_RING(map.Order(),RING); // RING needed for Alm
  map_RING.Import_nograde(map);
  */
  
  //  write_Healpix_map_to_fits(outfile,map_RING     , datatype);

  int lmax=3*map.Nside();// too large values give wrong results, this is recommended value
  int mmax=3*map.Nside();

  Alm<xcomplex<double> > alm(lmax,mmax);

  arr<double>  weight;
  weight.alloc(2*map.Nside());// must have at least 2*map.Nside() entries
  weight.fill(1.0);

  bool add_alm=false;
    
  //  map2alm(map,alm,weight,add_alm); // NB consider iter version
  //            (0 equivalent to map2alm)
  int num_iter = 3;                      // AWS20111102  
  map2alm_iter(map,alm,num_iter,weight); // AWS20111102

  double dtr=acos(-1.0)/180.;


  double fwhm = fwhm_deg*dtr;// radian !!
  smoothWithGauss(alm,fwhm);

  // filter
  int filter = 0;
  int lfilter_min=0;
  int lfilter_max=30;
  int mfilter_min=0;
  int mfilter_max=1000;

  //  for m<0  defined as (-1)^m Al|m| so not explicitly required ? Healpix manual A.4 eq. (22)
  if(filter==1)
  for(int l=0;l<=lmax;l++)              
  for(int m=0;m<=l;   m++)if(l<lfilter_min || l>lfilter_max || m<mfilter_min || m>mfilter_max) alm(l,m)=0.;

  alm2map(alm,map_RING_out);

  cout<<"energy bin map average before smoothing = "<< map.average() << " after smoothing = " <<map_RING_out.average()
      <<" after/before = "<<map_RING_out.average()/map.average() <<endl;
  
  out.write_column(colnum,map_RING_out.Map()); // fitshandle.h


  if(colnum==1)// do this only once
  {
   outfile= outfile_Alm;
   cout<<"writing Alm to "<<outfile<<endl;

   fitshandle almout;
   almout.create("!"+outfile); // ! to overwrite
   write_Alm_to_fits (almout, alm,  lmax,  mmax, datatype);
   almout.close();
  }

    }

  out.close();

  }// energy bins smoothed block







  ////////////////////////////////////////////////////////////////////////////////
  // energy integrated  with smoothing, using already converted format as input //
  ////////////////////////////////////////////////////////////////////////////////

  {

  cout<<"==== smoothing  fluxes integrated over energy"<<endl;


  cout<<" Gaussian FWHM = "<<fwhm_deg<< " deg"<<endl;

  Healpix_Map<double> map;
  Healpix_Map<double> map_RING_out(order,RING);

  outfile= outfile_integrated_energies_smoothed;
  cout<<"writing to "<<outfile<<endl;

  out.create("!"+outfile); // ! to overwrite
  prepare_Healpix_fitsmap(out,map_RING_out, datatype, colname);

  out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name

  infile=outfile_integrated_energies;
  int ncolnum=0;

  cout<< "reading "<<infile<<endl;

  
  
  fitshandle in;
  in.open(infile);
  in.goto_hdu(2);

  
  ncolnum=in.ncols();
  cout<<" number of columns ="<<ncolnum<<endl;


  for (colnum=1;colnum<=ncolnum;colnum++)

  {
    if(debug==1) cout<< "reading "<<infile<<",  column number: "<< colnum;

  read_Healpix_map_from_fits(infile,map,colnum,2);// HDU 2 (1=primary)

  if(debug==1)
  {
   cout<<" Npix  = "<<map.Npix(); 
   cout<<" Nside = "<<map.Nside();
   cout<<" Order = "<<map.Order();
   cout<<" Scheme= "<<map.Scheme()<<endl; // 0 = RING, 1 = NESTED
  }

  /*
  // convert to RING (WMAP is NESTED, galprop is RING already)
  Healpix_Map<double> map_RING(map.Order(),RING); // RING needed for Alm
  map_RING.Import_nograde(map);
  */
  
  //  write_Healpix_map_to_fits(outfile,map_RING     , datatype);

  int lmax=3*map.Nside();// too large values give wrong results, this is recommended value
  int mmax=3*map.Nside();

  Alm<xcomplex<double> > alm(lmax,mmax);

  arr<double>  weight;
  weight.alloc(2*map.Nside());// must have at least 2*map.Nside() entries
  weight.fill(1.0);

  bool add_alm=false;
    
  //  map2alm(map,alm,weight,add_alm); // NB consider iter version
  //            (0 equivalent to map2alm)
  int num_iter = 3;                      // AWS20111102  
  map2alm_iter(map,alm,num_iter,weight); // AWS20111102

  double dtr=acos(-1.0)/180.;


  double fwhm = fwhm_deg*dtr;// radian !!
  smoothWithGauss(alm,fwhm);

  // filter
  int filter = 0;
  int lfilter_min=0;
  int lfilter_max=30;
  int mfilter_min=0;
  int mfilter_max=1000;

  //  for m<0  defined as (-1)^m Al|m| so not explicitly required ? Healpix manual A.4 eq. (22)
  if(filter==1)
  for(int l=0;l<=lmax;l++)              
  for(int m=0;m<=l;   m++)if(l<lfilter_min || l>lfilter_max || m<mfilter_min || m>mfilter_max) alm(l,m)=0.;

  alm2map(alm,map_RING_out);

  cout<<"energy integrated map average before smoothing = "<< map.average() << " after smoothing = " <<map_RING_out.average()
      <<" after/before = "<<map_RING_out.average()/map.average() <<endl;
  
  out.write_column(colnum,map_RING_out.Map()); // fitshandle.h




    }

  out.close();

}


  cout<<endl<<" ==== healpix_skymap_convert complete"<<endl;

 return 0;
}