Index: trunk/SophyaProg/PrgMap/cl2map.cc
===================================================================
--- trunk/SophyaProg/PrgMap/cl2map.cc	(revision 1429)
+++ trunk/SophyaProg/PrgMap/cl2map.cc	(revision 1430)
@@ -27,5 +27,5 @@
     cout << " cl2map : Spherical harmonics synthesis - Power spectrum C_l -> HEALPix map \n" 
 	 << " Usage: cl2map [-float/-r4] [-msph pixp] [-beam fwhm] \n" 
-	 << "        [fitsin] [-fitsout] [-autoinirand] InFileName OutFileName \n" 
+	 << "        [-fitsin] [-fitsout] [-autoinirand] InFileName OutFileName \n" 
 	 << "   -float (-r4): single precision C_l and map (default = double) \n"
 	 << "   -msph pix: Spherical map pixelization parameter \n" 
Index: trunk/SophyaProg/PrgMap/map2cl.cc
===================================================================
--- trunk/SophyaProg/PrgMap/map2cl.cc	(revision 1429)
+++ trunk/SophyaProg/PrgMap/map2cl.cc	(revision 1430)
@@ -27,5 +27,5 @@
     cout << " map2cl : Spherical harmonics analysis - HEALPix map -> Power spectrum C_l \n" 
 	 << " Usage: map2cl [-float/-r4] [-lmax lval] [-thetacut dtdeg] \n" 
-	 << "        [fitsin] [-fitsout] InFileName OutFileName \n" 
+	 << "        [-fitsin] [-fitsout] InFileName OutFileName \n" 
 	 << "   -float (-r4): single precision C_l and map (default = double) \n"
 	 << "   -lmax lval: Maximum value for l index (default=255)\n" 
Index: trunk/SophyaProg/PrgMap/prjsmap.cc
===================================================================
--- trunk/SophyaProg/PrgMap/prjsmap.cc	(revision 1430)
+++ trunk/SophyaProg/PrgMap/prjsmap.cc	(revision 1430)
@@ -0,0 +1,312 @@
+#include "machdefs.h"    // Definitions specifiques SOPHYA
+
+#include <math.h>
+#include <iostream.h>
+
+#include "nbmath.h"
+#include "timing.h"
+
+#include "array.h"
+#include "skymap.h"
+#include "samba.h"
+#include "sambainit.h"
+#include "fitsspherehealpix.h"
+#include "fitstarray.h"
+
+
+
+/* Nouvelle-Fonction */
+void Usage(bool fgerr)
+{
+  cout << endl;
+  if (fgerr) {
+    cout << " prjsmap : Argument Error ! prjsmap -h for usage " << endl; 
+    exit(1);
+  }
+  else {
+    cout << " prjsmap : Molleweide and Sinus projection for sky maps \n" 
+	 << " Usage: prjsmap [-float/-r4] [-sinus] [-sx sizeX] [-sy sizeY] [-scale fact] \n" 
+	 << "        [-off val][-ump val] [-invy] [-fitsin] [-fitsout] InFileName OutFileName\n" 
+	 << "   -float (-r4): single precision sky map and projection (default = double) \n"
+	 << "   -sinus: Selects Sinus projection (default Molleweide projection) \n"
+	 << "   -sx sizeX: Projected map size in X \n"
+	 << "   -sy sizeY: Projected map size in Y \n"
+	 << "       default: sizeX = 2*sizeY ; sizeX*sizeY ~ map.NbPixels()/2 \n"
+	 << "   -scale fact: scale factor applied to skymap pixels (def=1.) \n"
+	 << "   -off val: offset value for projected map (def=0.) \n"
+	 << "             prj(ix, jy) = offset + scale*skymap(theta, phi) \n"
+	 << "   -ump val: Unmapped pixel value in projected map (def=0.) \n" 
+	 << "   -invy: reverses the projection direction along Y (Theta) \n"
+	 << "   -fitsout: Select the FITS format for the output map (default PPF format) \n"
+	 << "   -fitsin : Select the FITS format for the input vector (default PPF format) \n"
+	 << "   InFileName : Input file name (HEALPix map) \n" 
+	 << "   OutFileName : Output file name (the C_l vector) \n" << endl; 
+    exit(0);
+  }  
+}
+
+/* Nouvelle-Fonction */
+template <class T> 
+class _PrjMap {
+public :
+/* Methode */
+static void PM_MeanSigma(PixelMap<T> const & map, double& gmoy, double& gsig, T& min, T& max)
+{
+  min = max = map(0);
+  gmoy=0.;
+  gsig = 0.; 
+  double valok;
+  for(int k=0; k<map.NbPixels(); k++) {
+    valok = map(k);
+    gmoy += valok;  gsig += valok*valok; 
+    if (map(k)<min)  min = map(k);
+    else if (map(k)>max)  max = map(k);      
+  }
+  gmoy /= (double)map.NbPixels();
+  gsig = gsig/(double)map.NbPixels() - gmoy*gmoy;
+  if (gsig >= 0.) gsig = sqrt(gsig);
+  return;
+}
+
+/* Methode */
+static TArray<T> Project_Molleweide(SphericalMap<T> const & map, int sx=0, int sy=0, T ump=0, 
+				    T offset=0, T scale=1, bool invy=false)
+{
+  if ( (sx <= 0) && (sy <= 0) ) {
+    sy = sqrt((double)(map.NbPixels()/2));
+    sx = sy*2;    
+  }
+  else {
+    if (sx <= 0) sx = 2*sy;
+    else sy = sx/2;
+  }
+  TArray<T> prj(sx, sy);
+  prj = ump;   // On met tout a ump
+  
+  cout << " Molleweide projection (sizeX= " << sx << " sizeY= " << sy << " )" 
+       << " Scale=" << scale << " Offset=" << offset << endl;
+  r_8 xa, yd, teta,phi, facteur;
+  int_4 ix, jy, k;
+  for(jy=0; jy<sy; jy++) {
+    yd = (r_8)(jy+0.5)/(r_8)sy-0.5;
+    if (invy) teta = (0.5-yd)*M_PI;
+    else teta = (yd+0.5)*M_PI;
+    facteur=2.*M_PI/sin(acos((double)yd*2));
+    for(ix=0; ix<sx; ix++)  {
+      xa = (r_8)(ix+0.5)/(r_8)sx-0.5;
+      phi = xa*facteur+M_PI;
+      if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
+        k = map.PixIndexSph(teta, phi);
+        prj(ix,jy,0) = offset + scale*map(k);
+      }
+    }
+  }
+
+  prj.SetTemp(true);
+  return prj;   
+}
+
+/* Methode */
+static TArray<T> Project_Sinus(SphericalMap<T> const & map, int sx=0, int sy=0, T ump=0, 
+			       T offset=0, T scale=1, bool invy=false)
+{
+  if ( (sx <= 0) && (sy <= 0) ) {
+    sy = sqrt((double)(map.NbPixels()/2));
+    sx = sy*2;    
+  }
+  else {
+    if (sx <= 0) sx = 2*sy;
+    else sy = sx/2;
+  }
+  TArray<T> prj(sx, sy);
+  prj = ump;   // On met tout a ump
+
+  cout << " Sinus projection (sizeX= " << sx << " sizeY= " << sy << " )"
+       << " Scale=" << scale << " Offset=" << offset << endl;
+  r_8 xa, yd, teta,phi, facteur;
+  int_4 ix, jy, k;
+  for(jy=0; jy<sy; jy++) {
+    yd = (r_8)(jy+0.5)/(r_8)sy-0.5;
+    if (invy) teta = (0.5-yd)*M_PI;
+    else teta = (yd+0.5)*M_PI;
+    facteur=1./sin(teta);
+    for(ix=0; ix<sx; ix++)  {
+      xa = (r_8)(ix+0.5)/(r_8)sx-0.5;
+      phi = 2*M_PI*xa*facteur+M_PI;
+      if ( (phi <= 2*M_PI) && (phi >= 0.) ) {
+        k = map.PixIndexSph(teta, phi);
+        prj(ix,jy,0) = offset + scale*map(k);
+      }
+    }
+  }
+  prj.SetTemp(true);
+  return prj;   
+}
+
+/* Methode */
+static void Project(string & infile, string & outfile, int sx, int sy, T ump, T offset, 
+                    T scale, bool fgsin, bool invy, bool fgfitsin, bool fgfitsout) 
+{
+
+  double deg2rad =  M_PI/180.;
+  double minute2rad =  M_PI/(180.*60.);
+
+  SphericalMap<T> * sphp = NULL;
+  SphereHEALPix<T> sphh;
+
+  PPersist * pp = NULL;
+
+  if (fgfitsin) {
+    cout << "--- Reading Input FITS file " << infile << endl; 
+    FitsInFile fii(infile);
+    fii >> sphh;
+    sphp = &sphh;
+  }
+  else {
+    cout << "--- Reading Input PPF file " << infile << endl; 
+    PInPersist ppi(infile);
+    pp = ppi.ReadObject();
+    if (pp) 
+      sphp = dynamic_cast< SphericalMap<T> * > (pp->DataObj());
+    else sphp = NULL;
+    if (sphp == NULL) {
+      cout << " Object in PPF file is not a SphericalMap<T> or wrong data type ! " << endl;
+      throw TypeMismatchExc("_PrjMap::Project(...) Error reading input PPF file !");
+    }
+  }
+
+  SphericalMap<T> & sph = *sphp;
+  T min, max;
+  double mean, sigma;
+  PM_MeanSigma(sph, mean, sigma, min, max);
+  
+  cout << " Input map type= " << sph.TypeOfMap() << " PixParm (=NSide for HEALPix)= " 
+       << sph.SizeIndex() << endl;
+  cout << " Input map : NbPixels= " << sph.NbPixels() << " Resolution= "
+       << sqrt(sph.PixSolAngle(0))/minute2rad << " Arcminutes " << endl;
+  cout << " map.Min= " << min << " map.Max= " << max 
+       << " map.Mean= " << mean << " map.Sigma= " << sigma << endl;
+  
+  cout << "--- Sky map projection " << endl;
+  TArray<T> prj;
+
+  if (fgsin) prj = Project_Sinus(sph, sx, sy, ump, offset, scale, invy); 
+  else prj = Project_Molleweide(sph, sx, sy, ump, offset, scale, invy); 
+
+  cout << "--- Statistics on the projected map :" ; 
+  prj.Show();
+  prj.MinMax(min, max);
+  MeanSigma(prj, mean, sigma);  
+  cout << " prj.Min= " << min << " prj.Max= " << max 
+       << " prj.Mean= " << mean << " prj.Sigma= " << sigma << endl;
+  
+  if (fgfitsout) {
+    cout << "--- Writing projection to Output FITS file " << outfile << endl; 
+    FitsOutFile fio(outfile);
+    fio << prj;
+  }
+  else {
+    cout << "--- Writing projection to Output PPF file " << outfile << endl; 
+    POutPersist ppo(outfile);
+    ppo << prj;
+  }
+
+  if (pp) delete pp;
+}
+
+};
+
+/* Main-Program */
+int main(int narg, char *arg[])
+{
+  if ((narg > 1) && (strcmp(arg[1],"-h") == 0) ) Usage(false);
+  
+  double ump = 0.;
+  double off = 0.;
+  double scale = 1.;
+  int sx = 0;
+  int sy = 0;
+  string infile;
+  string outfile;
+  bool fgfitsin = false;
+  bool fgfitsout = false;
+  bool fgr4 = false;
+  bool fgsinus = false;
+  bool fginvy = false;
+
+  cout << " prjsmap : Decoding command line options ... " << endl;
+
+  int ko=1;
+  for (int k=1; k<narg; k++)   {
+    if (strcmp(arg[k], "-sx") == 0)  {
+      if (k == narg-1) Usage(true);  // -sx est suivi d'un argument 
+      sx = atoi(arg[k+1]);  k++;       // k++ pour sauter au suivant
+    }
+    if (strcmp(arg[k], "-sy") == 0)  {
+      if (k == narg-1) Usage(true);  // -sy est suivi d'un argument 
+      sy = atoi(arg[k+1]);  k++;       // k++ pour sauter au suivant
+    }
+    else if (strcmp(arg[k], "-scale") == 0)  {
+      if (k == narg-1) Usage(true);  
+      scale = atof(arg[k+1]);  k++;   
+    }
+    else if (strcmp(arg[k], "-off") == 0)  {
+      if (k == narg-1) Usage(true);  
+      off = atof(arg[k+1]);  k++;   
+    }
+    else if (strcmp(arg[k], "-ump") == 0)  {
+      if (k == narg-1) Usage(true);  
+      ump = atof(arg[k+1]);  k++;   
+    }
+    else if (strcmp(arg[k], "-sinus") == 0) {
+      fgsinus = true;  
+    }
+    else if (strcmp(arg[k], "-invy") == 0) {
+      fginvy = true;  
+    }
+    else if (strcmp(arg[k], "-fitsin") == 0) {
+      fgfitsin = true;  
+    }
+    else if (strcmp(arg[k], "-fitsout") == 0) {
+      fgfitsout = true;  
+    }
+    else if ((strcmp(arg[k], "-float") == 0) || (strcmp(arg[k], "-r4") == 0) ) {
+      fgr4 = true;  
+    }
+   
+    else { ko = k; break; }  // Debut des noms
+  }
+
+  if ((narg-ko) < 1)  Usage(true); 
+  infile = arg[ko];
+  outfile = arg[ko+1];
+  
+  // Bloc try de gestion des exception 
+  try {
+    InitTim();
+    SophyaInit();
+    if (fgr4) {
+      cout << " Projecting SphericalMap<r_4>  -->  Array<r_4> (float)" << endl;
+      _PrjMap<r_4>::Project(infile, outfile, sx, sy, ump, off, scale, 
+			    fgsinus, fginvy, fgfitsin, fgfitsout);
+    }
+    else {
+      cout << " Projecting SphericalMap<r_8>  -->  Array<r_8> (double)" << endl;
+      _PrjMap<r_8>::Project(infile, outfile, sx, sy, ump, off, scale, 
+			    fgsinus, fginvy, fgfitsin, fgfitsout);
+    }
+  }
+  catch (PThrowable & exc) {   // Exceptions de SOPHYA
+    cerr << " prjsmap: Catched Exception " << (string)typeid(exc).name()
+         << " - Msg= " << exc.Msg() << endl;
+  }
+  catch (...) {    // Autres Exceptions
+    cerr << " prjsmap: some other exception was caught ! " << endl;
+  }
+  
+  PrtTim("End of prjsmap ");
+  return(0);
+}
+
+
+  
Index: trunk/SophyaProg/PrgUtil/programs
===================================================================
--- trunk/SophyaProg/PrgUtil/programs	(revision 1429)
+++ trunk/SophyaProg/PrgUtil/programs	(revision 1430)
@@ -4,3 +4,4 @@
 cl2map cl2map.o
 map2cl map2cl.o
+prjsmap prjsmap.o
 
