Index: trunk/Cosmo/RadioBeam/calcpk2.cc
===================================================================
--- trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3784)
+++ trunk/Cosmo/RadioBeam/calcpk2.cc	(revision 3784)
@@ -0,0 +1,162 @@
+/*  ------------------------ Projet BAORadio -------------------- 
+  Programme de calcul du spectre de puissance (3D) a partir d'un 
+  cube de  delta T/T LSS, d'un cube  delta T/T LSS synchrotron
+  ou radio-sources, apres ajustement / soustraction d'une loi de 
+  puissance en frequence (l'axe Z du tableau doit etre en frequence)
+
+    R. Ansari , C. Magneville - Juin 2010 
+
+Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk 
+---------------------------------------------------------------  */
+
+#include "machdefs.h"
+#include "sopnamsp.h"
+#include <iostream>
+#include <string>
+#include <math.h>
+
+#include <typeinfo>
+
+#include "specpk.h"
+#include "histats.h"
+
+#include "qhist.h" 
+#include "histinit.h"
+#include "fftwserver.h"
+#include "randr48.h"      
+
+#include "ctimer.h"
+
+typedef ThSDR48RandGen RandomGenerator ;
+//--- Declaration des fonctions 
+TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx);
+
+//-------------------------------------------------------------------------
+//      ------------------ MAIN PROGRAM ------------------------------
+//-------------------------------------------------------------------------
+/* --Fonction-- */
+int main(int narg, const char* arg[])
+{
+  if (narg<6) {
+    cout << " Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync OutPk  " << endl;
+    return 1;
+  }
+  Timer tm("calcpk2");
+  int rc = 0;
+  try {
+    string inppflss = arg[1];
+    r_4 rfaclss = atof(arg[2]);
+    string inppfsync = arg[3];
+    r_4 rfacsync = atof(arg[4]);
+    string outname = arg[5];
+    TArray<r_4> inmaplss, inmapsync;
+    const char * tits[2]={"LSS", "Sync"};
+    for(int ks=0; ks<2; ks++)  {
+      string& ppfname=inppflss;
+      r_4 rfac=rfaclss;
+      TArray<r_4>* inmap=&inmaplss; 
+      if (ks==1) {  ppfname=inppfsync;  rfac=rfacsync;  inmap=&inmapsync; }
+      cout << "calcpk2[" << ks+1 << "] : reading 3D map " << tits[ks] << " from file " << ppfname 
+	   << "  RenormFactor=" << rfac << endl;
+      PInPersist pin(ppfname);
+      pin >> (*inmap);
+      (*inmap) *= rfac;
+      double mean, sigma;
+      MeanSigma(*inmap, mean, sigma);
+      cout << " ...InMap sizes " << inmap->InfoString() << endl;
+      inmap->Show(); 
+      cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
+    }
+    tm.Split(" After read ");
+    TArray<r_8> synctemp, specidx; 
+    cout << "calcpk2[3] : calling CleanForeground(...) " << endl;
+    TArray<r_4> inmap = CleanForeground(inmaplss, inmapsync, synctemp, specidx);
+    double mean, sigma;
+    MeanSigma(inmap, mean, sigma);
+    cout << " After cleaning: Mean=" << mean << " Sigma=" << sigma << endl;      
+    tm.Split(" After CleanForeground");
+
+    cout << "calcpk2[4] : computing 3D Fourier coefficients ... " << endl;
+    FFTWServer ffts;
+    TArray< complex<r_4> > four3d;
+    ffts.FFTForward(inmap, four3d);
+    tm.Split(" After FFTForward ");
+
+    cout << "calcpk2[5] : computing power spectrum ... " << endl;
+    RandomGenerator rg;
+    Four3DPk pkc(four3d, rg);
+    
+    HProf hp = pkc.ComputePk(0.,256);
+
+    tm.Split(" Done ComputePk ");    
+    
+    cout << "calcpk2[4] : writing profile P(k) and foreground maps to  " << outname << endl;
+    POutPersist po(outname);
+    po << PPFNameTag("Pk") << hp;
+    po << PPFNameTag("Tsync") << synctemp;
+    po << PPFNameTag("async") << specidx;
+
+  }  // End of try bloc 
+  catch (PThrowable & exc) {  // catching SOPHYA exceptions
+    cerr << " calcpk2.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
+         << "\n...exc.Msg= " << exc.Msg() << endl;
+    rc = 99;
+  }
+  catch (std::exception & e) {  // catching standard C++ exceptions
+    cerr << " calcpk2.cc: Catched std::exception "  << " - what()= " << e.what() << endl;
+    rc = 98;
+  }
+  catch (...) {  // catching other exceptions
+    cerr << " calcpk2.cc: some other exception (...) was caught ! " << endl;
+    rc = 97;
+  }
+  cout << " ==== End of calcpk.cc program  Rc= " << rc << endl;
+  return rc;    
+}
+
+
+/* --Fonction-- */
+TArray<r_4> CleanForeground(TArray<r_4>& maplss, TArray<r_4>& mapsync, TArray<r_8>& synctemp, TArray<r_8>& specidx)
+{
+  bool smo;
+  if (!maplss.CompareSizes(mapsync,smo) ) {
+    cout << " CleanForeground/ERROR sizes " << endl;
+    maplss.Show();  mapsync.Show(); 
+    throw ParmError("CleanForeground- maplss , mapsync not the same size !");
+  }
+  sa_size_t sz[5];   sz[0]=maplss.SizeX();  sz[1]=maplss.SizeY();
+  synctemp.SetSize(2, sz); 
+  specidx.SetSize(2, sz);
+  TArray<r_4>& omap=maplss;
+  Vector vlnf(maplss.SizeZ());
+  int nprt = 0;
+  double freq0=840.;  // Frequence premier index en k (MHz)
+  double dfreq=1.;   // largeur en frequence de chaque plan (Mhz)  
+  for(sa_size_t i=0; i<maplss.SizeX(); i++) 
+    for(sa_size_t j=0; j<maplss.SizeY(); j++)  {
+      r_8 s1, sx, sx2, sy, sxy;
+      s1=sx=sx2=sy=sxy=0.;
+      for(sa_size_t k=0; k<maplss.SizeZ(); k++)  {
+	double lnf=log((double)k*dfreq+freq0);
+	vlnf(k)=lnf;
+	double ttot=(r_8)(mapsync(i,j,k)+maplss(i,j,k));
+	if (ttot < 1.e-5) continue;
+	double lntt=log(ttot);
+	s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
+	sy+=lntt;    sxy+=(lnf*lntt);
+      }
+      double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
+      double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
+      double T0 = exp(beta+alpha*vlnf(0));
+      if ((i%6==0)&&(j%7==0)) 
+	cout << "CleanForeground[" << i << "," << j << "]: T0=" << T0 << " alpha=" << alpha 
+	     << " (mapsync=" << mapsync(i,j,0) << " ... " << mapsync(i,j,125) << ")" << endl; 
+      synctemp(i,j) = T0;
+      specidx(i,j) = alpha;
+      for(sa_size_t k=0; k<maplss.SizeZ(); k++) {
+	r_4 fittedtemp = (r_4)(exp(beta+alpha*vlnf(k)));
+	omap(i,j,k) = mapsync(i,j,k)+maplss(i,j,k)-fittedtemp;
+      }
+    }
+  return omap;
+}
