// Generate input spectra  a + b* l + c * gaussienne(l, 50, 20)
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;
int m = 128;  // HealPix pixelisation parameter
SphereHEALPix<r_8> map(m);
ylmserver.GenerateFromCl(map, m,  clin, 0.);
// Compute power spectrum from map
Vector clout = ylmserver.DecomposeToCl(map, lmax,  0.);

// Copy the synthetised map
SphereHEALPix<r_8> mapdip;
mapdip = map;

// Define the dipole direction
double thetadipole = 1.2;
double phidipole = 0.3;
UnitVector vd(thetadipole, phidipole);
UnitVector vc(thetadipole, phidipole);

// Add a dipole component to the map
for(int i=0; i<mapdip.NbPixels(); i++) {  // Boucle sur les pixels
  double Thetar,Phir;
  mapdip.PixThetaPhi(i,Thetar,Phir);  // direction du pixel      
  vc.SetThetaPhi(Thetar,  Phir);      
  mapdip(i) += 5.*vd.Psc(vc);  // Produit scalaire 
}

// Compute the new map power spectrum 
Vector cldip = ylmserver.DecomposeToCl(mapdip, lmax,  0.);

// Display objects (kept automatically through Display)
ExecuteCommand("zone 1 2");
DisplayObj(map, " ");
DisplayObj(mapdip, " ");
DisplayObj(clin, "win");
DisplayObj(clout, "same,red,fcirclemarker5,defline");
KeepObj(cldip);