/*
 *  Methods of class Sym_tensor_tt related to eta and mu
 *
 *   (see file sym_tensor.h for documentation)
 *
 */

/*
 *   Copyright (c) 2003-2004 Eric Gourgoulhon & Jerome Novak
 *
 *   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
 *
 */


char sym_tensor_tt_etamu_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $" ;

/*
 * $Id: sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
 * $Log: sym_tensor_tt_etamu.C,v $
 * Revision 1.18  2014/10/13 08:53:44  j_novak
 * Lorene classes and functions now belong to the namespace Lorene.
 *
 * Revision 1.17  2014/10/06 15:13:19  j_novak
 * Modified #include directives to use c++ syntax.
 *
 * Revision 1.16  2006/10/24 13:03:19  j_novak
 * New methods for the solution of the tensor wave equation. Perhaps, first
 * operational version...
 *
 * Revision 1.15  2005/04/01 14:28:32  j_novak
 * Members p_eta and p_mu are now defined in class Sym_tensor.
 *
 * Revision 1.14  2004/06/04 09:25:58  e_gourgoulhon
 * Method eta(): eta is no longer computed from h^rr but from khi (in the
 *   case where khi is known).
 *
 * Revision 1.13  2004/05/25 15:08:44  f_limousin
 * Add parameters in argument of the functions update, eta and mu for
 * the case of a Map_et.
 *
 * Revision 1.12  2004/05/24 13:45:29  e_gourgoulhon
 * Added parameter dzp to method Sym_tensor_tt::update.
 *
 * Revision 1.11  2004/05/05 14:24:54  e_gourgoulhon
 * Corrected a bug in method set_khi_mu: the division of khi by r^2
 * was ommitted in the case dzp=2 !!!
 *
 * Revision 1.10  2004/04/08 16:38:43  e_gourgoulhon
 * Sym_tensor_tt::set_khi_mu: added argument dzp (dzpuis of resulting h^{ij}).
 *
 * Revision 1.9  2004/03/04 09:53:04  e_gourgoulhon
 * Methods eta(), mu() and upate(): use of Scalar::mult_r_dzpuis and
 * change of dzpuis behavior of eta and mu.
 *
 * Revision 1.8  2004/03/03 13:16:21  j_novak
 * New potential khi (p_khi) and the functions manipulating it.
 *
 * Revision 1.7  2004/02/05 13:44:50  e_gourgoulhon
 * Major modif. of methods eta(), mu() and update() to treat
 * any value of dzpuis, thanks to the new definitions of
 * Scalar::mult_r(), Scalar::dsdr(), etc...
 *
 * Revision 1.6  2004/01/28 13:25:41  j_novak
 * The ced_mult_r arguments have been suppressed from the Scalar::*dsd* methods.
 * In the div/mult _r_dzpuis, there is no more default value.
 *
 * Revision 1.5  2003/11/05 15:28:31  e_gourgoulhon
 * Corrected error in update.
 *
 * Revision 1.4  2003/11/04 23:03:34  e_gourgoulhon
 * First full version of method update().
 * Add method set_rr_mu.
 * Method set_eta_mu ---> set_rr_eta_mu.
 *
 * Revision 1.3  2003/11/04 09:35:27  e_gourgoulhon
 * First operational version of update_tp().
 *
 * Revision 1.2  2003/11/03 22:33:36  e_gourgoulhon
 * Added methods update_tp and set_eta_mu.
 *
 * Revision 1.1  2003/11/03 17:08:37  e_gourgoulhon
 * First version
 *
 *
 * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_tt_etamu.C,v 1.18 2014/10/13 08:53:44 j_novak Exp $
 *
 */

// Headers C
#include <cstdlib>
#include <cassert>

// Headers Lorene
#include "tensor.h"

			//--------------//
			//     khi      //
			//--------------//

namespace Lorene {
const Scalar& Sym_tensor_tt::khi() const {

  if (p_khi == 0x0) {   // a new computation is necessary
		
    // All this has a meaning only for spherical components:
    assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ; 

    // khi is computed from $h^{rr}$ component

    p_khi = new Scalar(operator()(1,1)) ;
    p_khi->mult_r() ;
    p_khi->mult_r() ;
  }

  return *p_khi ; 

}


			//--------------//
			//     eta      //
			//--------------//
			
			
const Scalar& Sym_tensor_tt::eta(Param* par) const {


    if (p_eta == 0x0) {   // a new computation is necessary
	

	    // All this has a meaning only for spherical components:
	    assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ; 

        // eta is computed from the divergence-free condition:

        int dzp = operator()(1,1).get_dzpuis() ; 
        int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;

        Scalar source_eta(*mp) ; 
        
        if (p_khi == 0x0) { //  eta is computed from h^rr
                            //  -------------------------
        
	        source_eta = - operator()(1,1).dsdr() ; 	
        
        // dhrr contains - dh^{rr}/dr in all domains but the CED,                                           
        // in the CED:   - r^2 dh^{rr}/dr        if dzp = 0          (1)
        //               - r^(dzp+1) dh^{rr}/dr  if dzp > 0          (2)
                                                    
		        
	        // Multiplication by r of (-d h^{rr}/dr) (with the same dzpuis as h^{rr})
	        source_eta.mult_r_dzpuis( dzp ) ;                           

            // Substraction of the h^rr part and multiplication by r :
            source_eta -= 3. * operator()(1,1) ;                          

            source_eta.mult_r_dzpuis(dzp_resu) ;
        }
        else {      // eta is computed from khi
                    // ------------------------

            source_eta = - p_khi->dsdr() ;
            int diff_dzp = source_eta.get_dzpuis() - dzp_resu ; 
            assert( diff_dzp >= 0 ) ;
            source_eta.dec_dzpuis(diff_dzp) ;
            
            Scalar tmp(*p_khi) ; 
            tmp.div_r_dzpuis(dzp_resu) ; 
            
            source_eta -= tmp ;  
            
        }

        
	    // Resolution of the angular Poisson equation for eta
	    // --------------------------------------------------
	    if (dynamic_cast<const Map_af*>(mp) != 0x0) {
	        p_eta = new Scalar( source_eta.poisson_angu() ) ; 
	    }
	    else {
	        Scalar resu (*mp) ;
	        resu = 0. ;
	        mp->poisson_angu(source_eta, *par, resu) ;
	        p_eta = new Scalar( resu ) ;  	    
	    }
	
    }

    return *p_eta ; 

}

			

			//-------------------//
			//  set_rr_eta_mu    //
			//-------------------//
			

void Sym_tensor_tt::set_rr_eta_mu(const Scalar& hrr, const Scalar& eta_i, 
		const Scalar& mu_i) {

		// All this has a meaning only for spherical components:
		assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ; 
						
		set(1,1) = hrr ; 	// h^{rr}
							// calls del_deriv() and therefore delete previous
							// p_eta and p_mu
		
		p_eta = new Scalar( eta_i ) ; 	// eta

		p_mu = new Scalar( mu_i ) ; 	// mu 
		
		update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
		
}
			
			//---------------//
			//  set_rr_mu    //
			//---------------//
			

void Sym_tensor_tt::set_rr_mu(const Scalar& hrr, const Scalar& mu_i) {

		// All this has a meaning only for spherical components:
		assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ; 
						
		set(1,1) = hrr ; 	// h^{rr}
							// calls del_deriv() and therefore delete previous
							// p_eta and p_mu
		
		p_mu = new Scalar( mu_i ) ; 	// mu 
		
		eta() ; // computes eta form the divergence-free condition
		
		update( hrr.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
		
}
			
			//-------------------//
			//  set_khi_eta_mu    //
			//-------------------//
			

void Sym_tensor_tt::set_khi_eta_mu(const Scalar& khi_i, const Scalar& eta_i, 
		const Scalar& mu_i) {

  // All this has a meaning only for spherical components:
  assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ; 
			
  set(1,1) = khi_i ;
  set(1,1).div_r() ;
  set(1,1).div_r() ;     // h^{rr}

  // calls del_deriv() and therefore delete previous
  // p_khi, p_eta and p_mu
		
  p_khi = new Scalar( khi_i ) ;        // khi

  p_eta = new Scalar( eta_i ) ; 	// eta
  
  p_mu = new Scalar( mu_i ) ; 	// mu 
  
  update( khi_i.get_dzpuis() ) ; // all h^{ij}, except for h^{rr}
		
}
			
			//---------------//
			//  set_khi_mu    //
			//---------------//
			

void Sym_tensor_tt::set_khi_mu(const Scalar& khi_i, const Scalar& mu_i, 
                              int dzp, Param* par1, Param* par2, Param* par3) {

  // All this has a meaning only for spherical components:
  assert( dynamic_cast<const Base_vect_spher*>(triad) != 0x0 ) ; 
						
  set(1,1) = khi_i ; 
                        // calls del_deriv() and therefore delete previous
                        // p_eta and p_mu

  assert( khi_i.check_dzpuis(0) ) ; 
  if (dzp == 0) {
    set(1,1).div_r() ;
    set(1,1).div_r() ;	// h^{rr}
  }
  else {
    assert(dzp == 2) ; //## temporary: the other cases are not treated yet
    set(1,1).div_r_dzpuis(1) ;
    set(1,1).div_r_dzpuis(2) ;
  }		
  
  p_khi = new Scalar ( khi_i ) ;  // khi

  p_mu = new Scalar( mu_i ) ; 	// mu 

  if (dynamic_cast<const Map_af*>(mp) != 0x0) {		
      eta() ; // computes eta form the divergence-free condition
      // dzp = 0 ==> eta.dzpuis = 0
      // dzp = 2 ==> eta.dzpuis = 1
      
      update(dzp) ; // all h^{ij}, except for h^{rr}
  }
  else {
      eta(par1) ; // computes eta form the divergence-free condition
      // dzp = 0 ==> eta.dzpuis = 0
      // dzp = 2 ==> eta.dzpuis = 1
      
      update(dzp, par2, par3) ; // all h^{ij}, except for h^{rr}
  }

}

			

			//-------------//
			//   update    //
			//-------------//
			


void Sym_tensor_tt::update(int dzp, Param* par1, Param* par2) {

    // All this has a meaning only for spherical components:
    assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ; 

    assert( (p_eta != 0x0) && (p_mu != 0x0) ) ; 
	        
    Itbl idx(2) ;
    idx.set(0) = 1 ;	// r index
	
	// h^{r theta} : 
	// ------------
	idx.set(1) = 2 ;	// theta index
	*cmp[position(idx)] = p_eta->srdsdt() - p_mu->srstdsdp() ; 
    
    if (dzp == 0) {
        assert( cmp[position(idx)]->check_dzpuis(2) ) ; 
        cmp[position(idx)]->dec_dzpuis(2) ;
    }
    
    assert( cmp[position(idx)]->check_dzpuis(dzp) ) ; 
    
	// h^{r phi} :
	// ------------
	idx.set(1) = 3 ;	// phi index
	*cmp[position(idx)] = p_eta->srstdsdp() + p_mu->srdsdt() ; 

    if (dzp == 0) {
        assert( cmp[position(idx)]->check_dzpuis(2) ) ; 
        cmp[position(idx)]->dec_dzpuis(2) ;
    }
	
    assert( cmp[position(idx)]->check_dzpuis(dzp) ) ; 
	
	// h^{theta phi} and h^{phi phi}
	// -----------------------------
	
	//--------------  Computation of T^theta   --> taut : 
    
    Scalar tautst = operator()(1,2).dsdr() ; 

    // dhrr contains  dh^{rt}/dr in all domains but the CED,                                           
    // in the CED:    r^2 dh^{rt}/dr        if dzp = 0          (1)
    //                r^(dzp+1) dh^{rt}/dr  if dzp > 0          (2)
                                                    
	// Multiplication by r of dh^{rt}/dr (with the same dzpuis than h^{rt})
    tautst.mult_r_dzpuis( operator()(1,2).get_dzpuis() ) ; 	
    
    	        
    // Addition of the remaining parts :	
	tautst += 3 * operator()(1,2) - operator()(1,1).dsdt() ; 
	tautst.mult_sint() ; 
	
	Scalar tmp = operator()(1,1) ;
	tmp.mult_cost() ; 		// h^{rr} cos(th)
	
	tautst -= tmp ; 	// T^th / sin(th)
	
	Scalar taut = tautst ; 
	taut.mult_sint() ; 	// T^th
	

	//----------- Computation of T^phi   --> taup : 
    
	Scalar taupst = - operator()(1,3).dsdr() ; 

    // dhrr contains  - dh^{rp}/dr in all domains but the CED,                                           
    // in the CED:    - r^2 dh^{rp}/dr        if dzp = 0          (3)
    //                - r^(dzp+1) dh^{rp}/dr  if dzp > 0          (4)
                                                    	        
	// Multiplication by r of -dh^{rp}/dr  (with the same dzpuis than h^{rp})
    taupst.mult_r_dzpuis( operator()(1,3).get_dzpuis() ) ; 	

                          
    // Addition of the remaining part :	
	
	taupst -= 3 * operator()(1,3) ; 
	taupst.mult_sint() ; 	// T^ph / sin(th)
	
	Scalar taup = taupst ; 
	taup.mult_sint() ; 		// T^ph 
	
    
    //------------------- Computation of F and h^[ph ph}
    
	tmp = tautst ; 
	tmp.mult_cost() ; 	// T^th / tan(th)
	
	// dT^th/dth + T^th / tan(th) + 1/sin(th) dT^ph/dph :
	tmp = taut.dsdt() + tmp + taup.stdsdp() ;

	Scalar tmp2 (*mp) ;		
	if (dynamic_cast<const Map_af*>(mp) != 0x0) {		
	    tmp2 = tmp.poisson_angu() ;  // F
	}
	else {
	    tmp2 = 0. ;
	    mp->poisson_angu(tmp, *par1, tmp2) ; // F
	}
	    

	tmp2.div_sint() ; 
	tmp2.div_sint() ; // h^{ph ph}
	
	idx.set(0) = 3 ;	// phi index
	idx.set(1) = 3 ;	// phi index
	*cmp[position(idx)] = tmp2 ; 		// h^{ph ph} is updated
	
    
    //------------------- Computation of G and h^[th ph}
    
	tmp = taupst ; 
	tmp.mult_cost() ; // T^ph / tan(th)
	
	// - 1/sin(th) dT^th/dph + dT^ph/dth + T^ph / tan(th) :
	tmp = - taut.stdsdp() + taup.dsdt() + tmp ; 
	
	if (dynamic_cast<const Map_af*>(mp) != 0x0) {		
	    tmp2 = tmp.poisson_angu() ;  // G
	}
	else {
	    tmp2 = 0. ;
	    mp->poisson_angu(tmp, *par2, tmp2) ;  // G
	}
	
	tmp2.div_sint() ; 
	tmp2.div_sint() ; // h^{th ph}
	
	idx.set(0) = 2 ;	// theta index
	idx.set(1) = 3 ;	// phi index
	*cmp[position(idx)] = tmp2 ; 		// h^{th ph} is updated
	
	// h^{th th}  (from the trace-free condition)
	// ---------
	idx.set(1) = 2 ;	// theta index
	*cmp[position(idx)] = - operator()(1,1) - operator()(3,3) ; 
	

	Sym_tensor_trans::del_deriv() ; //## in order not to delete p_eta and p_mu
	


}			


			//-----------------//
			//  set_A_tilde_B  //
			//-----------------//

void Sym_tensor_tt::set_A_tildeB(const Scalar& a_in, const Scalar& tb_in, 
				 Param* par_bc, Param* par_mat) {

    Scalar zero(*mp) ;
    zero.set_etat_zero() ;
    set_AtB_trace(a_in, tb_in, zero, par_bc, par_mat) ;
    return ;
}






}