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

#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 "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();
    
    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;
    }
    // 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
    SphericalTransformServer<r_8> ylmserver;
    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 << " Writing clin, map, clout to ecpylm.ppf" << endl;
    POutPersist so("ecpylm.ppf");
    so << PPFNameTag("clin") << clin <<  PPFNameTag("map") << map
       << PPFNameTag("clout") << clout ;
    cout << " Writing mappar, cloutpar to ecpylm.ppf" << endl;
    so << PPFNameTag("mappar") << mappar
       << PPFNameTag("cloutpar") << cloutpar ;
  }
  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;
}