source: Sophya/trunk/SophyaProg/Tests/tsphereecp.cc@ 4038

#include <stdio.h>
#include <stdlib.h>
#include "sopnamsp.h"
#include "math.h"
#include <iostream>
#include <string>

#include "sambainit.h"
#include "pexceptions.h"
#include <typeinfo>
#include "timing.h"
#include "resusage.h"

#include "randr48.h"
#include "sphereecp.h"
#include "fiosphereecp.h"
#include "array.h"
#include "samba.h"

int main(int narg, char* arg[])
{
  cout << " ---- Programme tsphereecp.cc - Test SphereECP --- " << endl;
  int rc = 0;
  try {
    SophyaInit();
    ResourceUsage ru;
    
    InitTim();
    SphereECP<r_8> ecp(M_PI/2.-0.05, M_PI/2.+0.05, 10, 0., 0.2, 10);
    ecp.SetPixels(15.);
    cout << ecp ;
 
    SphereECP<r_8> ecp2;
    ecp2 = ecp;
    ecp2 += 4.;
    cout << " ecp2 = ecp+4. = \n " << ecp2 ;
    
    {
      cout << " Writing ecp ecp2 to POutPersist sphecp.ppf " << endl;
      POutPersist po("sphecp.ppf");
      po << ecp << ecp2;
    }
    {
      SphereECP<r_8> recp, recp2;
      cout << " Reading from sphecp.ppf " << endl;
      PInPersist pi("sphecp.ppf");
      pi >> recp >> recp2;
      cout << " recp recp2 from PInPersist " << endl; 
      cout << recp;
      cout << recp2;
    }
    
    // On cree une carte complete et on la remplit 
    SphereECP<r_4> ecpf(200);
    for(int_4 kk=0; kk<ecpf.NbPixels(); kk++) { 
      double tet, phi;
      ecpf.PixThetaPhi(kk, tet, phi);
      ecpf(kk) = sin(tet)*cos(phi);
    }
    // Extraction d'une carte partielle de ecpf
    SphereECP<r_4> ecpex = ecpf.ExtractPartial(M_PI/2.-M_PI/12., M_PI/2.+M_PI/12, 0.23, 0.23+M_PI/2.5);
    cout << " ecpf (full) and ecpex (extracted partiel) created ... " << endl;
    cout << ecpf << ecpex ;
    {
      cout << " Writing ecpf ecpex to POutPersist sphecpfp.ppf " << endl;
      POutPersist po("sphecpfp.ppf");
      po << PPFNameTag("ecpf") << ecpf << PPFNameTag("ecpex") << ecpex;
    }
         
    // Test avec transforme Ylm 
    int lmax = 92;
    Vector clin(lmax);
    for(int l=0; l<lmax; l++) {
      double xx = (l-50.)/10.;
      clin(l) = 1.e-2 -1.e-4*l + 0.1*exp(-xx*xx);
    }
    // Compute map from spectra
    ThSDR48RandGen rg; 
    SphericalTransformServer<r_8> ylmserver(rg);
    SphereECP<r_8> map(256);
    cout << " Generating map from Cl " << endl;
    ylmserver.GenerateFromCl(map, -1,  clin, 0.);
    cout << " Computing Cl from map " << endl;
    Vector clout = ylmserver.DecomposeToCl(map, lmax,  0.);
    
    SphereECP<r_8> mappar(0.25*M_PI, 0.75*M_PI, 128, 0.25*M_PI, 1.25*M_PI, 256);
    cout << " Generating partial map from Cl " << endl;
    ylmserver.GenerateFromCl(mappar, -1,  clin, 0.);
    cout << " Computing Cl from partialmap " << endl;
    Vector cloutpar = ylmserver.DecomposeToCl(mappar, lmax,  0.);

    cout << " Extracting partial map and Computing Cl from mapex " << endl;
    SphereECP<r_8> mapex = map.ExtractPartial(0.25*M_PI, 0.75*M_PI, 0.25*M_PI, 1.25*M_PI);
    Vector clmex = ylmserver.DecomposeToCl(mapex, lmax,  0.);

    cout << " Writing clin, map, clout to ecpylm.ppf" << endl;
    POutPersist so("ecpylm.ppf");
    so << PPFNameTag("clin") << clin <<  PPFNameTag("map") << map
       << PPFNameTag("clout") << clout ;
    cout << " Writing mapex, clmex to ecpylm.ppf" << endl;
    so << PPFNameTag("clmex") << clmex <<  PPFNameTag("mapex") << mapex;
    cout << " Writing mappar, cloutpar to ecpylm.ppf" << endl;
    so << PPFNameTag("mappar") << mappar
       << PPFNameTag("cloutpar") << cloutpar ;
    
    cout << " --> End of tsphereecp computations \n" << ru; 
   
  }
  catch (PThrowable & exc) {
    cerr << " tsphereecp.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name() 
         << " - Msg= " << exc.Msg() << endl;
    rc = 99;
  }
  catch (std::exception & e) {
    cerr << " tsphereecp.cc: Catched std::xception "  
         << " - what()= " << e.what() << endl;
    rc = 98;
  }
  catch (...) {
    cerr << " ssphereecp.cc: some other exception (...) was caught ! " << endl;
    rc = 97;
  }
  PrtTim("End tsphereecp " );
  cout << " ---- Programme tsphereecp.cc -  FIN  (Rc=" << rc << ") --- " << endl;
  return rc;
}