/*
 * Method of class Isol_Hole to compute metric data associated to a quasistationary single 
 * black hole spacetime slice.
 *
 * (see file isol_hole.h for documentation).
 *
 */

/*
 *   Copyright (c) 2009 Nicolas Vasset
 *
 *   This file is part of LORENE.
 *
 *   LORENE is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 *   LORENE is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with LORENE; if not, write to the Free Software
 *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

// Headers Lorene
#include "param_elliptic.h"
#include "proto.h"
#include "excised_slice.h"
#include "unites.h"	    


namespace Lorene {
void Excised_slice::compute_stat_metric(double precis,  double Omega, bool NorKappa, 
					Scalar boundNoK, bool isCF,double relax, int mer_max, 
					int mer_max2, bool isvoid) {
    
    // Fundamental constants and units
    // -------------------------------
    using namespace Unites ;

    // Verifying that we are in the right case 
    assert(type_hor==1);

    // Noticing the user that the iteration has started

    cout << "================================================================" << endl;
    cout << "STARTING THE MAIN ITERATION FOR COMPUTING METRIC FIELDS" << endl;
    cout << "        iteration parameters are the following:        " << endl;
    cout << "        convergence precision required:" << precis << endl;
    cout << "        max number of global steps    :" << mer_max << endl;
    cout << "        relaxation parameter          :" << relax   << endl;
    cout << "================================================================" << endl;


  // Construction of a multi-grid (Mg3d) and an affine mapping from the class mapping
  // --------------------------------------------------------------------------------
  const Map_af* map = dynamic_cast<const Map_af*>(&mp) ;
  const Mg3d* mgrid = (*map).get_mg();
	
  // Construct angular grid for h(theta,phi) 
  const Mg3d* g_angu = (*mgrid).get_angu_1dom() ;
  
  const int nz = (*mgrid).get_nzone(); 	// Number of domains
  int nt = (*mgrid).get_nt(1); 	// Number of collocation points in theta in each domain
  const int np = (*mgrid).get_np(1); 
  const Coord& rr = (*map).r;
   Scalar rrr (*map) ; 
  rrr = rr ; 
  rrr.std_spectral_base();  

  // For now the code handles only horizons at r=1, corresponding to the first shell 
  // inner boundary. This test assures this is the case with our mapping.
  assert((rrr.val_grid_point(1,0,0,0) - 1.) <= 1.e-9); 
 
  // Angular mapping defined as well
  //--------------------------------
  double r_limits2[] = {rrr.val_grid_point(1,0,0,0), rrr.val_grid_point(2,0,0,0)} ; 
  const Map_af map_2(*g_angu, r_limits2); //2D mapping; check if this is useful.
  const Metric_flat& mets = (*map).flat_met_spher() ;
 

  //----------------
  // Initializations
  // ---------------
  Scalar logn (*map) ; logn = log(lapse) ; 
  logn.std_spectral_base();


  Scalar logpsi(*map) ; logpsi = log(conf_fact) ;
  logpsi.std_spectral_base();
  Scalar psi4 (*map) ; psi4 = conf_fact*conf_fact*conf_fact*conf_fact ;
  Scalar npsi (*map) ;  npsi =conf_fact*lapse ;
  

  Scalar conf_fact_new(*map) ; conf_fact_new.annule_hard(); conf_fact_new.std_spectral_base(); 
  Scalar npsi_new(*map); npsi_new.annule_hard(); npsi_new.std_spectral_base();

  Vector shift_new (*map, CON, (*map).get_bvect_spher()); 
  for(int i=1; i<=3; i++){
    shift_new.set(i)=0;
  }
  shift_new.std_spectral_base();
 
  // Non conformally flat variables
  //-------------------------------
  Sym_tensor gamtuu = mets.con() + hij; 
  Metric gamt(gamtuu);
  Metric gam(gamt.cov()*psi4) ;
  Sym_tensor gamma = gam.cov();
  

  // Extrinsic curvature variables
  //------------------------------
  Sym_tensor aa(*map, CON, (*map).get_bvect_spher());
  for (int iii= 1; iii<=3; iii++){ 
    for(int j=1; j<=3; j++){
      aa.set(iii,j)= 0;
    }
  }
  aa.std_spectral_base(); 
  Scalar aa_quad_scal(*map) ; aa_quad_scal = 0. ;

  Sym_tensor aa_hat(*map, CON, (*map).get_bvect_spher());
  for (int iii= 1; iii<=3; iii++){ 
    for(int j=1; j<=3; j++){
      aa_hat.set(iii,j)= 0;
    }
  }
 
  Sym_tensor kuu = aa/psi4 ;
  Sym_tensor kuu2 = aa_hat/(psi4*psi4*sqrt(psi4));
  Sym_tensor kdd = contract (gamma, 0, contract(gamma, 1, kuu, 0),1);
 

  // (2,1)-rank delta tensor: difference between ricci rotation coefficients. 
  //-------------------------------------------------------------------------
   Tensor delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
   Scalar tmp(*map);

   for (int i=1; i<=3; i++) {
     for (int j=1; j<=3; j++) {
       for (int k=1; k<=3; k++) {
	 tmp = 0.;
	 tmp =  -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
	 for (int l=1; l<=3; l++) {
	   tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j) 
			 + gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
	 }
       	 delta.set(k,i,j) += tmp ; 
       }
     }
   }
   

   // Conformal Rstar scalar(eq 61, Bonazzola et al. 2003) 
   Scalar Rstar = 
    0.25 * contract(gamt.con(), 0, 1,
		    contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 1), 0, 1 ) 
     - 0.5  * contract(gamt.con(), 0, 1,
		       contract(hij.derive_cov(mets), 0, 1, gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ; 
  
   Scalar norm(*map);
   Scalar norm3(*map);

   
   if (isvoid == false){
     cout <<"FAIL: case of non-void spacetime not treated yet" << endl;
   }
   else {

     // Parameters for the iteration 
     //-----------------------------
     
     double diff_ent = 1 ; // initialization of difference marker between two iterations; 

     int util = 0; // Tool used to stop tensorial iteration at any wished step "util"
  
  //////////////////////////////////////////////////////////////////////////////////////////////////
  ///////////// ITERATION   ////////////////////////////////////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////////////////////////////
  
  
  for(int mer=0 ;(diff_ent > precis) && (mer<mer_max) ; mer++) {    

    //Global relaxation coefficient

    // Scalar variables linked to the norm of normal vector to horizon.
    norm = sqrt(1. + hij(1,1)); norm.std_spectral_base();
    norm3 = sqrt(1. + hij(3,3)); norm3.std_spectral_base();
    
    ///////////////////////////
    // Solving for (Psi-1) //
    ///////////////////////////

  
    // Setting of the boundary
    //------------------------
    double   diric= 0.; 
    double   neum = 1.; 

    Vector ssalt = rrr.derive_cov(gam);
    Vector ssaltcon = ssalt.up_down(gam);
    Scalar ssnormalt = sqrt(contract (ssalt,0, ssaltcon, 0));
    ssnormalt.std_spectral_base();
      
    ssalt.annule_domain(nz-1);
    ssalt.annule_domain(0);
    ssaltcon.annule_domain(nz-1);
    ssaltcon.annule_domain(0);
    
    ssalt = ssalt/ssnormalt;
    ssaltcon = ssaltcon/ssnormalt;
    // \tilde{s} in the notations of Gourgoulhon and Jaramillo, 2006
    Vector ssconalt = ssaltcon*conf_fact*conf_fact; 
    ssconalt.std_spectral_base();
    ssconalt.annule_domain(nz-1);
    Scalar bound3bis =   -((1./conf_fact)*contract((contract(kdd,1,ssconalt,0)),0, ssconalt,0));
    
    bound3bis.annule_domain(nz-1);
    bound3bis += -conf_fact*ssconalt.divergence(gamt);
    bound3bis.annule_domain(nz-1);     
    bound3bis = 0.25*bound3bis;
    bound3bis += -contract(conf_fact.derive_cov(gamt),0,ssconalt,0) + conf_fact.dsdr();
    bound3bis.annule_domain(nz-1);
    bound3bis.std_spectral_base();
    bound3bis.set_spectral_va().ylm();    
    
    Mtbl_cf *boundd3bis = bound3bis.set_spectral_va().c_cf;
            
    // Computing the source
    //---------------------
    Scalar source_conf_fact(*map) ; source_conf_fact=3. ; // Pour le fun... 
    source_conf_fact.std_spectral_base();                            

    Scalar d2logpsi = contract(conf_fact.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);  
    d2logpsi.inc_dzpuis(1);
    
    source_conf_fact = -(0.125* aa_quad_scal )/(psi4*conf_fact*conf_fact*conf_fact) 
      +  conf_fact* 0.125* Rstar - d2logpsi; 
    
    source_conf_fact.std_spectral_base(); 
    
    if (source_conf_fact.get_etat() == ETATZERO) {
      source_conf_fact.annule_hard() ;
      source_conf_fact.set_dzpuis(4) ;
      source_conf_fact.std_spectral_base() ;
    }
    source_conf_fact.set_spectral_va().ylm();
      
    // System inversion   
    //-----------------
    Param_elliptic source11(source_conf_fact);
    // Resolution has been done for quantity Q-1, 
    // because our solver gives a vanishing solution at infinity! 
    conf_fact_new = 
      source_conf_fact.sol_elliptic_boundary(source11, *boundd3bis, diric , neum) + 1 ; 
      
    // tests for resolution
    //---------------------
    Scalar baba2 = (conf_fact_new-1).laplacian();
//     cout << "psi+1-resolution" << endl;
//     maxabs (baba2 - source_conf_fact);
   
    Scalar psinewbis = conf_fact_new -1. ; psinewbis.annule_domain(nz -1);
    psinewbis.std_spectral_base();
    psinewbis = psinewbis.dsdr();
    Scalar psinewfin2 (map_2) ;
    psinewfin2.allocate_all(); 
    psinewfin2.set_etat_qcq();  
    psinewfin2.std_spectral_base();
    
    for (int k=0; k<np; k++)
      for (int j=0; j<nt; j++) {
	psinewfin2.set_grid_point(0, k, j, 0) = 
	  psinewbis.val_grid_point(1, k,j,0) - bound3bis.val_grid_point(1, k, j, 0);
      }
    //  maxabs (psinewfin2);
    

    // Update during the loop
    //-----------------------
    conf_fact = conf_fact_new* (1-relax) + conf_fact* relax ;
    psi4 = conf_fact*conf_fact*conf_fact*conf_fact;
    logpsi = log(conf_fact) ; 
    logpsi.std_spectral_base();	   



    //////////////////////////    
    // Solving for (N*Psi -1)/ 
    //////////////////////////


    // Setting of the boundary
    //------------------------
    assert (NorKappa == false) ; 
    Scalar bound(*map);
    bound = (boundNoK)*conf_fact -1;
    bound.annule_domain(nz -1);
    bound.std_spectral_base();
    bound.set_spectral_va().ylm();
    Mtbl_cf *boundd = bound.get_spectral_va().c_cf;

    diric =1; 
    neum = 0 ; 

    // Computing the source ...      
    //-------------------------
    Scalar d2lognpsi = contract(npsi.derive_cov(mets).derive_cov(mets), 0, 1, hij, 0,1);
    d2lognpsi.inc_dzpuis(1); //  dzpuis correction.
    
    Scalar source_npsi = npsi*(aa_quad_scal*(7./8.)/(psi4*psi4) + Rstar/8.) - d2lognpsi; 
    source_npsi.std_spectral_base();
    if (source_npsi.get_etat() == ETATZERO) {
      source_npsi.annule_hard() ;
      source_npsi.set_dzpuis(4) ;
      source_npsi.std_spectral_base() ;
    }


    // Inversion of the operator
    //--------------------------
    Param_elliptic source1 (source_npsi); 
    npsi_new = source_npsi.sol_elliptic_boundary(source1, *boundd, diric, neum) ;

    npsi_new = npsi_new +1;
  

    // Resolution tests in npsi
    //-------------------------
    Scalar baba = npsi_new.laplacian();
    //       cout << "resolution_npsi" << endl;
    //      maxabs (baba - source_npsi);
    //  cout << "bound_npsi" << endl;
    Scalar npsibound2 (map_2) ;
    npsibound2.allocate_all(); 
    npsibound2.set_etat_qcq();  
    npsibound2.std_spectral_base();		 
    for (int k=0; k<np; k++)
      for (int j=0; j<nt; j++) {
	npsibound2.set_grid_point(0, k, j, 0) = 
	  npsi_new.val_grid_point(1, k,j,0) - bound.val_grid_point(1, k, j, 0) -1.;	     
      }
    //   maxabs (npsibound2);
    

    // Update during the loop
    //-----------------------
    npsi = npsi_new*(1-relax) + npsi* relax; 
    lapse = npsi/conf_fact; 
    logn = log(lapse);
    logn.std_spectral_base(); 


 
      ///////////////////////
      //Resolution in Beta //
     ////////////////////////


    // Setting of the boundary  
    //------------------------
    bound = (boundNoK)/(conf_fact*conf_fact) ;
    bound.annule_domain(nz -1);
 
    // Rotation parameter for spacetime
    Scalar hor_rot(*map); hor_rot.annule_hard(); hor_rot = Omega; 
    hor_rot.std_spectral_base() ; hor_rot.mult_rsint();
    hor_rot.annule_domain(nz -1);
    
    Vector limit = shift_new;
    Vector ephi(*map, CON, (*map).get_bvect_spher());
    ephi.set(1).annule_hard();
    ephi.set(2).annule_hard();
    ephi.set(3) = 1;
    ephi.std_spectral_base();
    ephi.annule_domain(nz -1);
    
    limit = bound*ssconalt + hor_rot*ephi;
    // Boundary is fixed by value of 3 components of a vector (rather than value of potentials)   
    limit.std_spectral_base(); 
    
    Scalar Vrb = limit(1);  Vrb.set_spectral_va().ylm();
    Scalar mmub = limit.mu(); mmub.set_spectral_va().ylm(); 
    Scalar etab = limit.eta(); etab.set_spectral_va().ylm();
        
    // Computing the source
    //---------------------
    Vector deltaA =  - 2*lapse*contract(delta, 1,2, aa, 0,1);
    Vector hijddb =  - contract (shift.derive_cov(mets).derive_cov(mets), 1,2, hij, 0,1) ;
    Vector hijddivb =  
      - 0.3333333333333* contract (shift.divergence(mets).derive_cov(mets),0, hij,1);
    hijddb.inc_dzpuis(); // dzpuis fixing patch... 
    hijddivb.inc_dzpuis(); 
       
    Vector sourcevect2(*map,CON, (*map).get_bvect_spher()); 
    sourcevect2 = 2.* contract(aa, 1, lapse.derive_cov(mets),0) 
      - 12*lapse*contract(aa, 1, logpsi.derive_cov(mets), 0)  
      + deltaA + hijddb + hijddivb ; 
       
    sourcevect2.set(1).set_dzpuis(4);
    sourcevect2.set(2).set_dzpuis(4);
    sourcevect2.set(3).set_dzpuis(4);
    sourcevect2.std_spectral_base(); 
    if(sourcevect2.eta().get_etat() == ETATZERO)
      { sourcevect2.set(2).annule_hard();}
    
    double lam = (1./3.);    
       
    // System inversion
    //-----------------
    sourcevect2.poisson_boundary2(lam, shift_new, Vrb, etab, mmub, 1., 0., 1. ,0. ,1. ,0.) ;   


    // resolution tests
    //-----------------
    Vector source2 = contract(shift_new.derive_con(mets).derive_cov(mets), 1,2) 
      + lam* contract(shift_new.derive_cov(mets), 0,1).derive_con(mets);
    source2.inc_dzpuis(1);
    //  maxabs (source2 - sourcevect2);   

    Scalar mufin = shift_new.mu();
    mufin.set_spectral_va().coef();
    
    Scalar mufin2 (map_2) ;
    mufin2.allocate_all(); 
    mufin2.set_etat_qcq();  
    mufin2.std_spectral_base();
    
    for (int k=0; k<np; k++)
      for (int j=0; j<nt; j++) {
	mufin2.set_grid_point(0, k, j, 0) = 
	  mufin.val_grid_point(1, k,j,0) - mmub.val_grid_point(1, k, j, 0);
      }
    //  maxabs (mufin2);

    Scalar brfin = shift_new(1);
    brfin.set_spectral_va().coef();
		 
    Scalar brfin2 (map_2) ;
    brfin2.allocate_all(); 
    brfin2.set_etat_qcq();  
    brfin2.std_spectral_base();
		 
    for (int k=0; k<np; k++)
      for (int j=0; j<nt; j++) {		     
	brfin2.set_grid_point(0, k, j, 0) = 
	  brfin.val_grid_point(1, k,j,0) - Vrb.val_grid_point(1, k, j, 0);	     
      }
    //  maxabs (brfin2);
    

    // Update during the loop
    //-----------------------
    for (int ii=1; ii <=3; ii++){
      shift.set(ii) = shift_new(ii)*(1-relax) + shift(ii)* relax;
    }

    diff_ent = max(maxabs(npsi_new - npsi )); // Convergence parameter (discutable relevance...) 


    ////////////////////////////////////
    // Tensor hij resolution        ///
    ///////////////////////////////////

    if (isCF == false){
   
      if (diff_ent <=5.e-3) { // No resolution until we are close to the result.
	
	//WARNING; parameter maybe to be changed according to convergence parameters 
	//in Poisson-Hole/Kerr1.9/pplncp.C
	util = util+1; // Loop marker for NCF equation.
	
	
	///////////////////////////////////////////////////////////////////////////////////////////
	///////////// ITERATION   /////////////////////////////////////////////////////////////////
	///////////////////////////////////////////////////////////////////////////////////////////

	// Local convergence can be asked for hij equation. 
	// Allows to satisfy integrability conditions.
	
	// Here we ask for a local convergence; this can be improved later. 
	// WARNING; parameter maybe to be changed according to convergence parameters 
	// in Poisson-Hole/Kerr1.9/pplncp.C
	double precis2 =  1.e5*precis ; 
	
	double diff_ent2 = 1 ; // Local convergence marker
	// Local relaxation parameter. If not 1, the determinant condition won't be satisfied 
	// on a particular iteration.
	double relax2 = 0.5; 

	Sym_tensor sourcehij = hij; // Random initialization...
	// cuts off high spherical harmonics with threshold being the last argument.
	//  coupe_l_tous (hij, aa, nn, ppsi, bb, nt, 6); 


	for(int mer2=0 ;(diff_ent2 > precis2) && (mer2<mer_max2) ; mer2++) {    

	  // Calculation of the source
	  //--------------------------

	  // The double Lie derivative term is taken care of in the subroutine. 
	  secmembre_kerr(sourcehij);
       
	  //System inversion (note that no boundary condition is imposed)
	  //-------------------------------------------------------------
	  Sym_tensor hij_new = hij;
	  
	  hij_new = boundfree_tensBC (sourcehij, shift , conf_fact, lapse, hij, precis2);

	  cout << "maximum of convergence marker for the subiteration" << endl;

	  diff_ent2 = max(maxabs(hij - hij_new));
	  hij = relax2*hij_new + (1 - relax2)*hij;

	  cout << "mer2, diffent2" << endl;
	  cout << mer2 << endl;
	  cout << diff_ent2 << endl;
	}

	// Resolution tests
	//-----------------

	Sym_tensor gammatilde = mets.con() + hij;   
	Metric gammatilde2(gammatilde); Scalar detgam = gammatilde2.determinant();
	//       cout << "determinant of result" << endl;
	//       maxabs (detgam-1.);
	//      cout << "comment l'equation en hij est elle verifiee?" << endl;
        
	Sym_tensor test =contract (hij.derive_cov(mets).derive_con(mets), 2,3);
	test.annule(nz-1, nz-1);
	test = test 
	  - hij.derive_lie(shift).derive_lie(shift)
	  / ((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));        
	test.annule(nz-1, nz-1);
	Sym_tensor youps = test 
	  - sourcehij
	  / ((lapse/(conf_fact*conf_fact))*(lapse/(conf_fact*conf_fact)));
	//       maxabs (youps); 
	//       maxabs((youps).trace(mets));
	//       cout << " AAABBB" << endl;
	//       maxabs((youps).compute_A());
	//       maxabs((youps).compute_tilde_B());
      }
    }
 

    ////////////////////////////////////////////////////
    // Global variable update after an entire loop ////
    ///////////////////////////////////////////////////

    gamtuu = mets.con() + hij;
    gamt = gamtuu; // Metric
    gam = gamt.cov()*psi4;
    gamma = gam.cov();
    
    for (int i=1; i<=3; i++) {
      for (int j=1; j<=3; j++) {	   
	tmp = 0;       
	tmp = ((shift.derive_con(mets))(i,j) + (shift.derive_con(mets))(j,i) 
	       - (2./3.)*shift.divergence(mets) * mets.con()(i,j))*(1./(2.*lapse));	   
	aa.set(i,j) = tmp ; 
      }
    }

    //Non conformally flat correction; we assume here dhij/dt = 0.
    aa = aa -  (hij.derive_lie(shift) + (2./3.)*shift.divergence(mets)*hij)*(1./(2.*lapse)); 
    
    aa_hat = aa*psi4*sqrt(psi4); // Rescaling of traceless exrinsic curvature.
    aa_hat.std_spectral_base();

    hatA = aa_hat; // Probably obsolete very soon... replace aa_hat by hatA.

    Sym_tensor aaud = aa.up_down(gamt);
    Sym_tensor aaud_hat = aa_hat.up_down(gamt);
    aa_quad_scal =  contract(contract (aa_hat, 0, aaud_hat, 0), 0,1);
    
    delta = -0.5*contract( hij, 1, gamt.cov().derive_cov(mets), 2);
 
    for (int i=1; i<=3; i++) {
      for (int j=1; j<=3; j++) {
	for (int k=1; k<=3; k++) {
	  tmp = 0.;
	  tmp =  -0.5 *(gamt.cov().derive_con(mets))(i,j,k);
	  for (int l=1; l<=3; l++) {
	    tmp += -0.5*( gamt.cov()(i,l)*(hij.derive_cov(mets))(k,l,j) 
			  + gamt.cov()(l,j)*(hij.derive_cov(mets))(k,l,i));
	  }
	  delta.set(k,i,j) += tmp ; 
	}
      }
    } 

    Rstar = 
      0.25 * contract(gamt.con(), 0, 1,
		      contract(gamt.con().derive_cov(mets), 0, 1, 
			       gamt.cov().derive_cov(mets), 0, 1), 0, 1 ) 
      - 0.5  * contract(gamt.con(), 0, 1,
			contract(gamt.con().derive_cov(mets), 0, 1, 
				 gamt.cov().derive_cov(mets), 0, 2), 0, 1 ) ;  
    kuu = aa/(psi4); 
    kdd =  contract (gamma, 0, contract(gamma, 1, kuu, 0),1);  
 
    // Convergence markers
    //--------------------
    cout << "diffent" << endl; 
    cout<< diff_ent << endl;   
    cout <<"mer" << mer << endl;   
   

    //------------------------------------------------
    //     Check of Einstein equations (3+1 form)
    //------------------------------------------------
    

    // Lapse 
    //------
    Scalar lapse2(*map) ;
    lapse2 = lapse ;
    lapse2.std_spectral_base();
    
    // 3-metric
    //---------
    //   const Metric gam(mets.cov()*psi4) ;
   
    Sym_tensor gamt2(*map, COV, (*map).get_bvect_spher()); 
    for (int i=1; i<=3; i++)
      for (int j=1; j<=3; j++) 
	{ gamt2.set(i,j)=gam.cov()(i,j); 
	}

    //Shift
    //-----
    Vector beta = shift ;
  
    // Extrinsic curvature
    //--------------------
    Scalar TrK3(*map);
    Sym_tensor k_uu = aa/(psi4) ;
    k_uu.dec_dzpuis(k_uu(1,1).get_dzpuis()); 
    // Another way of computing the same thing, just to be sure... 
    Sym_tensor k_dd = k_uu.up_down(gam); 
    TrK3 = k_uu.trace(gam);
    //  TrK3.spectral_display("TraceKvraie", 1.e-10);
 
    // Hamiltonian constraint
    //-----------------------
    Scalar ham_constr = gam.ricci_scal() ;
    ham_constr.dec_dzpuis(3) ;
    ham_constr +=  TrK3*TrK3 - contract(k_uu, 0, 1, k_dd, 0, 1) ;
    // maxabs(ham_constr, "Hamiltonian constraint: ") ;

    ham_constr.set_spectral_va().ylm();
    //  ham_constr.spectral_display("ham_constr", 1.e-9);
  
    // Momentum constraint
    //-------------------
    Vector mom_constr = k_uu.divergence(gam)  - TrK3.derive_con(gam) ;
    mom_constr.dec_dzpuis(2) ;
    //   maxabs(mom_constr, "Momentum constraint: ") ;
    //   mom_constr(1).spectral_display("mom1", 1.e-9) ;
    //  mom_constr(2).spectral_display("mom2", 1.e-9) ;
    //  mom_constr(3).spectral_display("mom3", 1.e-9) ;
    
    // Evolution equations
    //--------------------
    Sym_tensor evol_eq = lapse2*gam.ricci() - lapse2.derive_cov(gam).derive_cov(gam);
    evol_eq.dec_dzpuis() ;
    evol_eq += k_dd.derive_lie(beta) ;
    evol_eq.dec_dzpuis(2) ;
    evol_eq += lapse2*(TrK3*k_dd - 2*contract(k_dd, 1, k_dd.up(0, gam), 0) ) ;
    //   maxabs(evol_eq, "Evolution equations: ") ;
    //   evol_eq.trace(gam).spectral_display("evoltrace", 1.e-10);
    //   maxabs (evol_eq.trace(gam));
    //  evol_eq.spectral_display("evol", 1.e-10);
  
  }

  cout << "================================================================" << endl;
  cout << "                THE ITERATION HAS NOW CONVERGED" << endl;
  cout << "================================================================" << endl;
   }
   return; 
}
}