Index: /trunk/SophyaProg/PrgMap/extrap2sph.cc
===================================================================
--- /trunk/SophyaProg/PrgMap/extrap2sph.cc	(revision 1526)
+++ /trunk/SophyaProg/PrgMap/extrap2sph.cc	(revision 1526)
@@ -0,0 +1,113 @@
+// Pour boucher les trous d'une sphere HEALPIX utilisant une autre sphere
+//             cmv 13/6/01
+#include "machdefs.h"
+#include <unistd.h>
+#include <stdexcept>
+#include <iostream.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "skymapinit.h"
+#include "skymap.h"
+#include "fitsspherehealpix.h" 
+
+void usage();
+void usage()
+{
+cout<<"extrap2sph [-w sphinw.fits -0 nulval]"<<endl
+    <<"           sphin.fits sphred.fits sphout.fits"<<endl
+    <<" -w sphinw.fits : input sphere  weight (def=NULL)"<<endl
+    <<" -0 nulval : value to identify empty pixel in sphin (def=0.)"<<endl
+    <<"             (used if sphinw.fits==NULL)"<<endl
+    <<" sphin.fits : input sphere"<<endl
+    <<" sphred.fits  : input reducted sphere"<<endl
+    <<" sphout.fits  : output sphere"<<endl;
+}
+
+int main(int narg, char* arg[])
+{
+double nulval=0.;
+char * fsphinw = NULL;
+int c;
+while((c = getopt(narg,arg,"h0:w:")) != -1) {
+  switch (c) {
+  case '0' :
+    sscanf(optarg,"%lf",&nulval);
+    break;
+  case 'w' :
+    fsphinw = optarg;
+    break;
+  case 'h' :
+  default:
+    usage(); exit(1);
+  }
+}
+
+if(optind+2>=narg) {usage(); exit(1);}
+char * fsphin   = arg[optind];
+char * fsphred = arg[optind+1];
+char * fsphout  = arg[optind+2];
+
+cout<<"Sphere Input : "<<fsphin<<endl
+    <<"Weight Sphere Input : "<<fsphinw<<endl
+    <<"Reducted Sphere Input : "<<fsphred<<endl
+    <<"Sphere Output : "<<fsphout<<endl
+    <<"nulval : "<<nulval<<endl;
+
+// Lecture des spheres
+SphereHEALPix<r_8> sphin;
+{FitsInFile sfits(fsphin); sfits >> sphin;}
+int nlat = sphin.SizeIndex();
+cout<<"Opening Sphere Input :"<<endl
+    <<"          Type of map : "<<sphin.TypeOfMap()<<endl
+    <<"     Number of pixels : "<<sphin.NbPixels()<<endl
+    <<"                 Nlat : "<<sphin.SizeIndex()<<endl;
+
+SphereHEALPix<r_8> sphinw;
+if(fsphinw) {
+  FitsInFile sfits(fsphinw);
+  sfits >> sphinw;
+  cout<<"Opening Weight Sphere Input :"<<endl
+      <<"          Type of map : "<<sphinw.TypeOfMap()<<endl
+      <<"     Number of pixels : "<<sphinw.NbPixels()<<endl
+      <<"                 Nlat : "<<sphinw.SizeIndex()<<endl;
+  if(sphinw.SizeIndex()!=nlat) {
+    cout<<"Incompatible sphin sphinw"<<endl;
+    exit(2);
+  }
+}
+
+SphereHEALPix<r_8> sphred;
+{FitsInFile sfits(fsphred); sfits >> sphred;}
+cout<<"Opening Reducted Sphere Input :"<<endl
+    <<"          Type of map : "<<sphred.TypeOfMap()<<endl
+    <<"     Number of pixels : "<<sphred.NbPixels()<<endl
+    <<"                 Nlat : "<<sphred.SizeIndex()<<endl;
+
+// Creation de la sphere output
+cout<<"Creating Sphere for output"<<endl;
+SphereHEALPix<r_8> sphout(sphin,false);
+
+// Filling hole for Output Sphere
+cout<<"...Filling hole for Output Sphere"<<endl;
+uint_4 n=0;
+for(int_4 i=0;i<sphout.NbPixels();i++) {
+  if(fsphinw) {if(sphinw(i)>0.) continue;}
+     else     {if(sphout(i)!=nulval) continue;}
+  double theta,phi;
+  sphout.PixThetaPhi(i,theta,phi);
+  int_4 ir = sphred.PixIndexSph(theta,phi);
+  sphout(i) = sphred(ir);
+  n++;
+}
+cout<<"      Number of filled pixels "<<n<<endl;
+
+// Ecriture de la sphere
+{
+string dum = "rm -f "; dum += fsphout; system(dum.c_str());
+FitsOutFile sfits(fsphout);
+cout<<"Writing Output Sphere Fits file"<<endl;
+sfits<<sphout;
+}
+
+exit(0);
+}
Index: /trunk/SophyaProg/PrgMap/extrapsph.cc
===================================================================
--- /trunk/SophyaProg/PrgMap/extrapsph.cc	(revision 1526)
+++ /trunk/SophyaProg/PrgMap/extrapsph.cc	(revision 1526)
@@ -0,0 +1,156 @@
+// Pour boucher les trous d'une sphere HEALPIX en clusterisant dans une sphere
+// plus petite             cmv 13/6/01
+#include "machdefs.h"
+#include <unistd.h>
+#include <stdexcept>
+#include <iostream.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "skymapinit.h"
+#include "skymap.h"
+#include "fitsspherehealpix.h" 
+
+void usage();
+void usage()
+{
+cout<<"extrapsph [-r reduc] sphin.fits sphinw.fits"<<endl
+    <<"                     sphout.fits [sphout_w.fits"
+    <<" sphred.fits sphred_w.fits]"<<endl
+    <<" -r reduc : reduction factor for nlat in clustering must be 2^n (def=2)"<<endl
+    <<" sphin.fits : input sphere"<<endl
+    <<" sphin_w.fits : input sphere filling weight"<<endl
+    <<" sphout.fits  : output sphere"<<endl
+    <<" sphout_w.fits  : output sphere filling weight"<<endl
+    <<" sphred.fits  : output reducted sphere"<<endl
+    <<" sphred_w.fits  : output reducted sphere filling weight"<<endl;
+}
+
+int main(int narg, char* arg[])
+{
+int red=2;
+int c;
+while((c = getopt(narg,arg,"hr:")) != -1) {
+  switch (c) {
+  case 'r' :
+    sscanf(optarg,"%d",&red);
+    break;
+  case 'h' :
+  default:
+    usage(); exit(1);
+  }
+}
+
+if(optind+2>=narg) {usage(); exit(1);}
+char * fsphin   = arg[optind];
+char * fsphinw  = arg[optind+1];
+char * fsphout  = arg[optind+2];
+char * fsphoutw = NULL;
+char * fsphred = NULL;
+char * fsphredw = NULL;
+if(optind+3<narg) fsphoutw = arg[optind+3];
+if(optind+4<narg) fsphred = arg[optind+4];
+if(optind+5<narg) fsphredw = arg[optind+5];
+
+cout<<"Sphere Input : "<<fsphin<<endl
+    <<"Weight Sphere Input : "<<fsphinw<<endl
+    <<"Sphere Output : "<<fsphout<<endl
+    <<"Weight Sphere Output : "<<fsphoutw<<endl
+    <<"Reducted Sphere Output : "<<fsphred<<endl
+    <<"Reducted Weight Sphere Output : "<<fsphredw<<endl
+    <<"Reduction : "<<red<<endl;
+
+// Lecture des spheres
+SphereHEALPix<r_8> sphin;
+{FitsInFile sfits(fsphin); sfits >> sphin;}
+int nlat = sphin.SizeIndex();
+int nlatr = nlat/red;
+cout<<"Opening Sphere Input :"<<endl
+    <<"          Type of map : "<<sphin.TypeOfMap()<<endl
+    <<"     Number of pixels : "<<sphin.NbPixels()<<endl
+    <<"                 Nlat : "<<sphin.SizeIndex()<<endl;
+
+SphereHEALPix<r_8> sphinw;
+{FitsInFile sfits(fsphinw); sfits >> sphinw;}
+cout<<"Opening Weight Sphere Input :"<<endl
+    <<"          Type of map : "<<sphinw.TypeOfMap()<<endl
+    <<"     Number of pixels : "<<sphinw.NbPixels()<<endl
+    <<"                 Nlat : "<<sphinw.SizeIndex()<<endl;
+if(sphinw.SizeIndex()!=nlat) {
+  cout<<"Incompatible sphin sphinw"<<endl;
+  exit(2);
+}
+
+// Creation des spheres output
+cout<<"Creating Spheres for output"<<endl;
+SphereHEALPix<r_8> sphout(sphin,false);
+SphereHEALPix<r_8> sphoutw(sphinw,false);
+
+// Creation des spheres reduites
+cout<<"Creating Reducted Spheres : nlatr = "<<nlatr<<endl;
+SphereHEALPix<r_8> sphred(nlatr);
+  sphred.SetPixels(0.);
+SphereHEALPix<r_8> sphredw(nlatr);
+  sphredw.SetPixels(0.);
+cout<<"Creating Reducted Sphere :"<<endl
+    <<"     Number of pixels : "<<sphred.NbPixels()<<endl
+    <<"                 Nlat : "<<sphred.SizeIndex()<<endl;
+
+// Filling reducted spheres
+cout<<"...Filling Reducted Spheres"<<endl;
+for(int_4 i=0;i<sphin.NbPixels();i++) {
+  if(sphinw(i)<=0.) continue;
+  double theta,phi;
+  sphin.PixThetaPhi(i,theta,phi);
+  int_4 ir = sphred.PixIndexSph(theta,phi);
+  sphred(ir)  += sphin(i)*sphinw(i);
+  sphredw(ir) += sphinw(i);
+}
+cout<<"...Computing Reducted Spheres"<<endl;
+uint_4 n=0;
+for(int_4 ir=0;ir<sphred.NbPixels();ir++)
+  if(sphredw(ir) > 0.) {sphred(ir) /= sphredw(ir); n++;}
+cout<<"      Number of filled pixels "<<n<<endl;
+
+// Filling hole for Output Spheres
+cout<<"...Filling hole for Output Spheres"<<endl;
+n=0;
+for(int_4 i=0;i<sphout.NbPixels();i++) {
+  if(sphoutw(i)>0.) continue;
+  double theta,phi;
+  sphout.PixThetaPhi(i,theta,phi);
+  int_4 ir = sphred.PixIndexSph(theta,phi);
+  if(sphredw(ir)<=0.) continue;
+  sphout(i)  = sphred(ir);
+  sphoutw(i) = -sphredw(ir);
+  n++;
+}
+cout<<"      Number of filled pixels "<<n<<endl;
+
+// Ecriture des spheres
+{
+string dum = "rm -f "; dum += fsphout; system(dum.c_str());
+FitsOutFile sfits(fsphout);
+cout<<"Writing Output Sphere Fits file"<<endl;
+sfits<<sphout;
+}
+if(fsphoutw) {
+string dum = "rm -f "; dum += fsphoutw; system(dum.c_str());
+FitsOutFile sfits(fsphoutw);
+cout<<"Writing Output Sphere Weight Fits file"<<endl;
+sfits<<sphoutw;
+}
+if(fsphred) {
+string dum = "rm -f "; dum += fsphred; system(dum.c_str());
+FitsOutFile sfits(fsphred);
+cout<<"Writing Reducted Sphere Fits file"<<endl;
+sfits<<sphred;
+}
+if(fsphredw) {
+string dum = "rm -f "; dum += fsphredw; system(dum.c_str());
+FitsOutFile sfits(fsphredw);
+cout<<"Writing Reducted Sphere Weight Fits file"<<endl;
+sfits<<sphredw;
+}
+
+exit(0);
+}
Index: /trunk/SophyaProg/PrgUtil/Makefile
===================================================================
--- /trunk/SophyaProg/PrgUtil/Makefile	(revision 1525)
+++ /trunk/SophyaProg/PrgUtil/Makefile	(revision 1526)
@@ -16,8 +16,8 @@
 LDLIBS := $(LIBS) $(LDLIBS)
 PROGS = $(EXE)runcxx $(EXE)scanppf $(EXE)scanfits $(EXE)cl2map $(EXE)map2cl $(EXE)prjsmap \
-        $(EXE)cremskfrsph
+        $(EXE)cremskfrsph $(EXE)extrapsph $(EXE)extrap2sph
 
 PROBJS = $(OBJ)runcxx.o $(OBJ)scanppf.o $(OBJ)scanfits.o $(OBJ)cl2map.o $(OBJ)map2cl.o \
-         $(OBJ)prjsmap.o $(OBJ)cremskfrsph.o
+         $(OBJ)prjsmap.o $(OBJ)cremskfrsph.o $(OBJ)extrapsph.o $(OBJ)extrap2sph.o
 
 all: $(PROGS)
@@ -61,4 +61,14 @@
 	@ echo cremskfrsph done
 $(EXE)cremskfrsph: $(OBJ)cremskfrsph.o
+	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
+
+extrapsph: $(EXE)extrapsph
+	@ echo extrapsph done
+$(EXE)extrapsph: $(OBJ)extrapsph.o
+	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
+
+extrap2sph: $(EXE)extrap2sph
+	@ echo extrap2sph done
+$(EXE)extrap2sph: $(OBJ)extrap2sph.o
 	$(LINK.cc) $^ $(LOADLIBES) $(LDLIBS) -o $@ 
 
@@ -292,2 +302,4 @@
  $(INC)tarrinit.h
 $(OBJ)cremskfrsph.o: cremskfrsph.cc
+$(OBJ)extrapsph.o: extrapsph.cc
+$(OBJ)extrap2sph.o: extrap2sph.cc
Index: /trunk/SophyaProg/PrgUtil/programs
===================================================================
--- /trunk/SophyaProg/PrgUtil/programs	(revision 1525)
+++ /trunk/SophyaProg/PrgUtil/programs	(revision 1526)
@@ -6,2 +6,4 @@
 prjsmap prjsmap.o
 cremskfrsph cremskfrsph.o
+extrapsph extrapsph.o
+extrap2sph extrap2sph.o
