#include <Metastability.h>

Inheritance diagram for FindAction:

Detailed Description

Definition at line 24 of file Metastability.h.

Public Member Functions
double	computeThValue ()

	FindAction (const StandardModel &SM_i)
	FindAction constructor. More...

gslpp::vector< double >	rkqs (double y01, double y02, double dydr0, double r0, double dr, double epsfrac[2], double epsabs[2])

	~FindAction ()
	FindAction destructor. More...

Public Member Functions inherited from ThObservable
double	getBinMax ()
	A get method to get the maximum value of the bin. More...

double	getBinMin ()
	A get method to get the minimum value of the bin. More...

const StandardModel &	getModel ()
	A get method to get the model. More...

const std::vector< std::string >	getParametersForObservable ()
	A get method to get the parameters for the specific observable. More...

void	setBinMax (double max)
	A set method to set the maximum value of the bin. More...

void	setBinMin (double min)
	A set method to set the minimum value of the bin. More...

void	setParametersForObservable (std::vector< std::string > parametersForObservable_i)
	A set method to get the parameters for the specific observable. More...

	ThObservable (const StandardModel &SM_i)
	Constructor. More...

	ThObservable (const ThObservable &orig)
	The copy constructor. More...

virtual	~ThObservable ()
	The default destructor. More...

Public Attributes
double	d2Vpar

double	delta_phi_cutoffpar

double	dVpar

double	phi0par

gslpp::vector< double >	potentialcoefficientspar

double	rpar

double	x1par

double	x2par

double	x3par

Protected Attributes
SUSYScalarPotential *	mySUSYScalarPotential

Protected Attributes inherited from ThObservable
double	max
	the bin maximum. More...

double	min
	The bin minimum. More...

std::vector< std::string >	parametersForObservable
	a vector of parameter namesfor the specific observable More...

const StandardModel &	SM
	A reference to an object of StandardMode class. More...

Private Member Functions
template<class F >
gsl_function	convertToGslFunctionS (const F &f)

double	deformedV (double phi)

gslpp::vector< double >	dY (double y1, double y2, double r)

int	dYfunc (double r, const double y[], double ODE[], void *flags)

int	dYJac (double r, const double y[], double dfdy, double dfdt[], void order)

gslpp::vector< double >	ExactSolution (double r, double phi0, double dV, double d2V)

double	fPS (double x)

double	func (double rpar)

gslpp::vector< double >	InitialConditions (double delta_phi0, double rmin, double delta_phi_cutoff, double distance, double dV_at_delta_phi0, double d2V_at_phi0)

gslpp::vector< double >	integrateProfile (double r0, double y01, double y02, double dr0, double epsfrac[2], double epsabs[2], double drmin, double rmax, double distance)

double	invertedpotential (double x)

double	Simpsonintegrand (double r, double phi, double dphi, double VphiMin_i)

Static Private Member Functions
template<class F >
static double	gslFunctionAdapterS (double x, void *p)

Private Attributes
const SUSY &	mySUSY

Constructor & Destructor Documentation

◆ FindAction()

FindAction::FindAction ( const StandardModel & SM_i )

FindAction constructor.

Definition at line 19 of file Metastability.cpp.

                                                :
    ThObservable(SM_i),
    potentialcoefficientspar(35,0.),
    mySUSY(static_cast<const SUSY&> (SM_i))
{
    mySUSYScalarPotential=new SUSYScalarPotential(SM_i);
}

◆ ~FindAction()

FindAction::~FindAction ( )

FindAction destructor.

Definition at line 27 of file Metastability.cpp.

{
  delete mySUSYScalarPotential;
}

Member Function Documentation

◆ computeThValue()

double FindAction::computeThValue ( )

virtual

Returns: \(FindAction\)

Implements ThObservable.

Definition at line 32 of file Metastability.cpp.

{
 
    /* 1. Definition of all potential minima */
 
    // define scalar potential
    gslpp::vector<double> potentialcoefficients=mySUSYScalarPotential->coefficients();
 
    double Vmin0=mySUSYScalarPotential->potential(potentialcoefficients, 0.0, 0.0, 0.0);
 
//    std::cout << "Vmin0 = " << Vmin0 << std::endl;
 
    gslpp::vector<double> dV0=mySUSYScalarPotential->potentialderivative(potentialcoefficients, 0.0, 0.0, 0.0);
 
//    std::cout << "dV0 = " << dV0 << std::endl;
 
    // calculate all minima (Ayan)
    // calculate V and dV at origin (V0 and dV0)
    // check whether there is a charge-breaking minimum very close to the origin
//    double S_0=0.0;
    gslpp::vector<double> S(10,0.);
    //define toy minima (remove this as soon as the true minimum finder is available)
    gslpp::vector<double> minimavector(3,0.);
    minimavector(0)=-807192.888;
    minimavector(1)=807260.056;
    minimavector(2)=-803413.309;
//    minimavector(3)=4.0;
//    minimavector(4)=5.0;
//    minimavector(5)=6.0;
//    minimavector(6)=7.0;
//    minimavector(7)=8.0;
//    minimavector(8)=9.0;
    //if mod3 length of minima is not zero, throw error
    int lengthofminima=3;
    int NofMinima=lengthofminima/3;
 
    gslpp::vector<double> deeperminima(lengthofminima+NofMinima,0.), dV(3,0.);
//    gslpp::vector<double> deeperminima(3,0.), dV(3,0.);
    int i,n=0;
    double x1, x2, x3, Vmin;
 
    for(i=0;i<NofMinima;i++)
    {
        x1=minimavector(3*i);
        x2=minimavector(3*i+1);
        x3=minimavector(3*i+2);
        Vmin=mySUSYScalarPotential->potential(potentialcoefficients, x1, x2, x3);
    std::cout << "Vmin = " << Vmin << std::endl;
        if(Vmin>=Vmin0)
        {
            continue;
        }
        else
        {
            deeperminima(4*n)=minimavector(3*i);
            deeperminima(4*n+1)=minimavector(3*i+1);
            deeperminima(4*n+2)=minimavector(3*i+2);
            deeperminima(4*n+3)=Vmin;
            n++;
        }
    }
 
    std::cout << "2." << std::endl;
 
    /* 2. Calculation of the tunneling rate for each of the minima */
 
//    gslpp::vector<double> path;
    //n is now the number of deeper minima
    for(i=0;i<n;i++)
    {
        x1=deeperminima(4*i);
        x2=deeperminima(4*i+1);
        x3=deeperminima(4*i+2);
        Vmin=deeperminima(4*i+3);
        dV=mySUSYScalarPotential->potentialderivative(potentialcoefficients, x1, x2, x3);
//        path=splinepath(x1,x2,x3,Vmin0,dV0,Vmin,dV);
        //the following function is supposed to find the cubic spline of the potential along a straight line from one to the other minimum
                int steps=100;
                gslpp::vector<double> linearpath(steps,0.);
                gslpp::vector<double> V(steps,0.);
//                gslpp::vector<double> spath(steps,0.);
                //maybe the following will not work in extreme cases
                double distance=x1*x1+x2*x2+x3*x3;
                double stepsize=distance/double(steps-1);
                int k;
                for(k=0;k<steps;k++)
                {
                    linearpath(k)=stepsize*k;
                    V(k)=mySUSYScalarPotential->potential(potentialcoefficients, x1*linearpath(k), x2*linearpath(k), x3*linearpath(k));
//    gsl_interp_accel *acc 
//      = gsl_interp_accel_alloc ();
//    gsl_spline *spline 
//      = gsl_spline_alloc (gsl_interp_cspline, 10);
//
//    gsl_spline_init (spline, x, y, 10);
//
//    for (xi = x[0]; xi < x[9]; xi += 0.01)
//      {
//        yi = gsl_spline_eval (spline, xi, acc);
//        printf ("%g %g\n", xi, yi);
//      }
//    gsl_spline_free (spline);
//        gsl_interp_accel_free (acc);
 
                }
 
 
    /* 2.1 Calculation of the tunneling rate along the undeformed (straight) path */
 
                    std::cout << "2.1.1" << std::endl;
 
//        double xmin = 0.001;
//        double xmax = std::numeric_limits<double>::infinity();
        
    /* 2.1.1 Finding the edge of the tunneling barrier */
 
        double min = 1.0e-12;
        double pos1=0.0;  //relative position of the metastable minimum
        double pos2=1.0;  //relative position of the stable minimum
        double pos = 0.5;  //relative position between metastable and stable minimum (first guess)
        while(fabs(pos1-pos2) > min)
        {
            double Vpos = mySUSYScalarPotential->potential(potentialcoefficients, pos*x1, pos*x2, pos*x3);
            if(Vpos > 0.0)  //maybe this should be ">"
            {
                pos1 = pos;
            }
            else
            {
                pos2 = pos;
            }
            pos = 0.5*(pos1+pos2);
        }
        //pos is now the edge of the barrier
 
    /* 2.1.2 Finding the typical scale of the barrier (using gsl 1D minimization) */
 
                    std::cout << "2.1.2" << std::endl;
 
        int status;
        int iter = 0, max_iter = 100;
        const gsl_min_fminimizer_type *T;
        gsl_min_fminimizer *s;
        double m=0.5*pos;
        double a=0.0, b=pos;
                    std::cout << "pos =" << pos << std::endl;
        potentialcoefficientspar=potentialcoefficients;
        x1par=x1;
        x2par=x2;
        x3par=x3;
        gsl_function F = convertToGslFunctionS(bind(&FindAction::invertedpotential, &(*this), _1));
//        F.function = -&invertedpotential; //pointer?
        T = gsl_min_fminimizer_brent;
        s = gsl_min_fminimizer_alloc(T);
        gsl_min_fminimizer_set(s, &F, m, a, b);
        do
        {
            iter++;
            status = gsl_min_fminimizer_iterate(s);
            m = gsl_min_fminimizer_x_minimum(s);
                    std::cout << "m =" << m << std::endl;
            a = gsl_min_fminimizer_x_lower(s);
            b = gsl_min_fminimizer_x_upper(s);
            status = gsl_min_test_interval(a, b, 0.0, 1.0e-3);
        } while(status == GSL_CONTINUE && iter<max_iter);
        gsl_min_fminimizer_free(s);
        double barriertop=m;
        if(barriertop<=0.0 || barriertop>=pos)
        {
            throw std::runtime_error("Error in Metastability.cpp: Potential barrier top outside the barrier range!");
        }
        double barrierheight=mySUSYScalarPotential->potential(potentialcoefficients, barriertop*x1, barriertop*x2, barriertop*x3)-Vmin0;
 
        if(barrierheight<=0.0)
        {
            throw std::runtime_error("Error in Metastability.cpp: No potential barrier!");
        }
        double rscale = barriertop/sqrt(6.0*barrierheight);
 
    /* 2.1.3 Finding a phi solution for the tunneling using an overshooting/undershooting algorithm  */
 
        double x = -log(pos);
//        double xincrease = 5.0;
 
        double rmin = 1.e-4*rscale;
        double rmax = 1.e4*rscale;
        double dr0 = rmin;
        double drmin = 0.01*rmin;
 
        double delta_phi = distance;
        double delta_phi_cutoff = 1.e-2 * delta_phi;
        double epsabs[2] = {fabs(delta_phi)*1.e-4 , fabs(delta_phi/rscale)*1.e-4};
        double epsfrac[2] = {1.e-4 , 1.e-4};
 
        double eps = distance*1.e-3;
        gslpp::vector<double> inconds(3,0.);
        do
        {
            double delta_phi0 = distance - exp(-x)*delta_phi;
            double sdp = delta_phi0/distance; //scaled_delta_phi0
 
            double dV_at_delta_phi0 = (mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0-2.0*eps)*x1/distance, (delta_phi0-2.0*eps)*x2/distance, (delta_phi0-2.0*eps)*x3/distance) 
                                       -8.0*mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0-eps)*x1/distance, (delta_phi0-eps)*x2/distance, (delta_phi0-eps)*x3/distance)
                                       +8.0*mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0+eps)*x1/distance, (delta_phi0+eps)*x2/distance, (delta_phi0+eps)*x3/distance) 
                                       -mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0+2.0*eps)*x1/distance, (delta_phi0+2.0*eps)*x2/distance, (delta_phi0+2.0*eps)*x3/distance) ) / (12.0*eps);
 
            double d2V_at_phi0 = (-mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0-2.0*eps)*x1/distance, (delta_phi0-2.0*eps)*x2/distance, (delta_phi0-2.0*eps)*x3/distance)
                                  +16.0*mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0-eps)*x1/distance, (delta_phi0-eps)*x2/distance, (delta_phi0-eps)*x3/distance)
                                  -30.0*mySUSYScalarPotential->potential(potentialcoefficients, sdp*x1, sdp*x2, sdp*x3)
                                  +16.0*mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0+eps)*x1/distance, (delta_phi0+eps)*x2/distance, (delta_phi0+eps)*x3/distance)
                                  -mySUSYScalarPotential->potential(potentialcoefficients, (delta_phi0+2.0*eps)*x1/distance, (delta_phi0+2.0*eps)*x2/distance, (delta_phi0+2.0*eps)*x3/distance) ) / (12.0*eps*eps);
 
            inconds = InitialConditions(delta_phi0, rmin, delta_phi_cutoff, distance, dV_at_delta_phi0, d2V_at_phi0);
            if(!std::isfinite(inconds(0)) || !std::isfinite(x))
            {
                break;                
            }
 
            //r_y = rf, yf1, yf2, convergence=1(converged), 2(undershoot), 3(overshoot)
            gslpp::vector<double> r_y(4,0.);
            r_y = integrateProfile(inconds(0), inconds(1), inconds(2), dr0, epsfrac, epsabs, drmin, rmax, distance);
//
//                        # Check for overshoot, undershoot
//                        if ctype == "converged":
//                            break
//                        elif ctype == "undershoot": # x is too low
//                            xmin = x
//                            x = x*xincrease if xmax == np.inf else .5*(xmin+xmax)
//                        elif ctype == "overshoot": # x is too high
//                            xmax = x
//                            x = .5*(xmin+xmax)
//                        # Check if we've reached xtol
//                        if (xmax-xmin) < xtol:
//                            break
 
        }
        while(true);
 
            /* 2.1.4 Integrate the barrier to get the tunneling action */
 
    /* 2.2 Calculation of the tunneling rate along the deformed path */
 
        /* 2.2.1 Deforming the path */
        /* 2.2.2 Calculating the tunneling rate */
 
    }//loop over the minima
 
 
 
    //          deform path between origin and i (old path is a line, new path is p_i)
        //define splines between origin and i
            //take origin and i
            //for i_spline < ?
            //add one knot, calculate its V and dV
        //split into parallel and orthogonal components
        //solve parallel tunneling
        //calculate N and minimize it wrt orthogonal directions
    //          calculate action integrating along p_i
        int rlength = 10;
    gslpp::vector<double> r(rlength,0.);
    gslpp::vector<double> phi(rlength,0.);
    gslpp::vector<double> dphi(rlength,0.);
//        double rlength = length of r;
        double VphiMin_i = deformedV(phi(rlength));
        //integrate.simps(integrand,r)
        double integral = 0.0;
        for(int j=1;j<rlength;j++)
        {
            integral += (r(j)-r(j-1))*(Simpsonintegrand(r(j-1),phi(j-1),dphi(j-1),VphiMin_i)
                                       +4.0*Simpsonintegrand((r(j)+r(j-1))/2.0,(phi(j)+phi(j-1))/2.0,(dphi(j)+dphi(j-1))/2.0,VphiMin_i)
                                       +Simpsonintegrand(r(j),phi(j),dphi(j),VphiMin_i))/(3.0*(double)rlength);
        }
 
//        double volume = r(0)*r(0)*r(0)*r(0)*M_PI*M_PI/2.0;
        
//        S.assign(j, integral+volume*(deformedV(phi(0))-VphiMin_i));
    //          check whether action is larger than S_(i-1)
 
    return 0.0;
}

◆ convertToGslFunctionS()

template<class F >

gsl_function FindAction::convertToGslFunctionS ( const F & f )

inlineprivate

Definition at line 89 of file Metastability.h.

{
    gsl_function gslFunction;
    
    const void* p = &f;
    assert (p != 0);
    
    gslFunction.function = &gslFunctionAdapterS<F>;
    // Just to eliminate the const.
    gslFunction.params = const_cast<void*>( p );
    
    return gslFunction;
}

◆ deformedV()

double FindAction::deformedV ( double phi )

private

Definition at line 318 of file Metastability.cpp.

{
    return 0.0;
}

◆ dY()

gslpp::vector< double > FindAction::dY	(	double	y1,
		double	y2,
		double	r
	)

private

Definition at line 549 of file Metastability.cpp.

{
    gslpp::vector<double> values(2,0.);
//    double alpha=3.0; /*spatial dimensions*/
    values(0)=y2;
//    values(1)=(mySUSYScalarPotential->potential(potentialcoefficients, (y1-2.0*eps)*x1/distance, (y1-2.0*eps)*x2/distance, (y1-2.0*eps)*x3/distance) 
//                                       -8.0*mySUSYScalarPotential->potential(potentialcoefficients, (y1-eps)*x1/distance, (y1-eps)*x2/distance, (y1-eps)*x3/distance)
//                                       +8.0*mySUSYScalarPotential->potential(potentialcoefficients, (y1+eps)*x1/distance, (y1+eps)*x2/distance, (y1+eps)*x3/distance) 
//                                       -mySUSYScalarPotential->potential(potentialcoefficients, (y1+2.0*eps)*x1/distance, (y1+2.0*eps)*x2/distance, (y1+2.0*eps)*x3/distance) ) / (12.0*eps)
//              -alpha*y2/r;
    return values;
}

◆ dYfunc()

int FindAction::dYfunc	(	double	r,
		const double	y[],
		double	ODE[],
		void *	flags
	)

private

◆ dYJac()

int FindAction::dYJac	(	double	r,
		const double	y[],
		double *	dfdy,
		double	dfdt[],
		void *	order
	)

private

◆ ExactSolution()

gslpp::vector< double > FindAction::ExactSolution	(	double	r,
		double	phi0,
		double	dV,
		double	d2V
	)

private

Definition at line 410 of file Metastability.cpp.

{
        gslpp::vector<double> phi_dphi(2,0.);
        double phi = 0.0;
        double dphi = 0.0;
 
        double curv = sqrt(fabs(d2V));
        double curv_r = curv*r;
        double nu = 1.0;
        if(curv_r<1.e-2)
        {
            double s = (d2V>0) ? 1.0 : -1.0;
            double add_to_phi;
            for(int k=1;k<=4;k++)
            {
                add_to_phi = pow(0.5*curv_r,2.0*k-2.0) * pow(s,k) / (gsl_sf_gamma(k+1)*gsl_sf_gamma(k+1+nu));
                phi += add_to_phi;
                dphi += add_to_phi * (2.0*k);
            }
            phi *= 0.25 * gsl_sf_gamma(nu+1) * r*r * dV * s;
            dphi *= 0.25 * gsl_sf_gamma(nu+1) * r * dV * s;
            phi += phi0;
        }
        else if(d2V>0)
        {
            phi = (gsl_sf_gamma(nu+1)*pow(0.5*curv_r,-nu) *gsl_sf_bessel_In(int(nu),curv_r)-1.0) * dV/d2V;
            dphi = -nu*(pow(0.5*curv_r,-nu) / r) * gsl_sf_bessel_In(int(nu), curv_r);
            dphi += pow(0.5*curv_r,-nu) * 0.5*curv * (gsl_sf_bessel_In(int(nu)-1, curv_r)+gsl_sf_bessel_In(int(nu)+1, curv_r));
            dphi *= gsl_sf_gamma(nu+1) * dV/d2V;
            phi += phi0;
        }
        else
        {
            phi = (gsl_sf_gamma(nu+1)*pow(0.5*curv_r,-nu) * gsl_sf_bessel_Jn(int(nu), curv_r)-1.0) * dV/d2V;
            dphi = -nu*(pow(0.5*curv_r,-nu) / r) * gsl_sf_bessel_Jn(int(nu), curv_r);
            dphi += pow(0.5*curv_r,-nu) * 0.5*curv * (gsl_sf_bessel_Jn(int(nu)-1, curv_r)-gsl_sf_bessel_Jn(int(nu)+1, curv_r));
            dphi *= gsl_sf_gamma(nu+1) * dV/d2V;
            phi += phi0;
        }
        phi_dphi(0)=phi;
        phi_dphi(1)=dphi;
        return phi_dphi;
}

◆ fPS()

double FindAction::fPS ( double x )

private

◆ func()

double FindAction::func ( double rpar )

inlineprivate

Definition at line 53 of file Metastability.h.

    {
        double exsol=ExactSolution(rpar, phi0par, dVpar, d2Vpar)(0);
        return fabs(exsol)-fabs(delta_phi_cutoffpar);
    };

◆ gslFunctionAdapterS()

template<class F >

static double FindAction::gslFunctionAdapterS	(	double	x,
		void *	p
	)

inlinestaticprivate

Definition at line 80 of file Metastability.h.

{
    // Here I do recover the "right" pointer, safer to use static_cast
    // than reinterpret_cast.
    F* function = static_cast<F*>( p );
    return (*function)( x );
}

◆ InitialConditions()

gslpp::vector< double > FindAction::InitialConditions	(	double	delta_phi0,
		double	rmin,
		double	delta_phi_cutoff,
		double	distance,
		double	dV_at_delta_phi0,
		double	d2V_at_phi0
	)

private

Definition at line 341 of file Metastability.cpp.

{
    double phi0 = distance + delta_phi0;
    double dV = dV_at_delta_phi0;
    double d2V = d2V_at_phi0;
    gslpp::vector<double> exsol(2,0.);
//switch this on!    exsol=ExactSolution(rmin, phi0, dV, d2V);
    double phi_r0=exsol(0);
    double dphi_r0=exsol(1);
 
    gslpp::vector<double> inc(3,0.);
    inc(0)=rmin;
    inc(1)=phi_r0;
    inc(2)=dphi_r0;
 
    if(fabs(phi_r0)<fabs(delta_phi_cutoff) || dphi_r0*delta_phi0>0)
    {
        double r = rmin;
        double rlast=r;
        while(std::isfinite(r))
        {
            rlast = r;
            r *= 10;
//switch this on!            exsol=ExactSolution(r, phi0, dV, d2V);
            if(fabs(exsol(0))>fabs(delta_phi_cutoff))
            {
                break;
            }
        }
 
        rpar=r;
        phi0par=phi0;
        dVpar=dV;
        d2Vpar=d2V;
        delta_phi_cutoffpar=delta_phi_cutoff;
 
        gsl_function F = convertToGslFunctionS(bind(&FindAction::func, &(*this), _1));
        //gsl root finding algorithm
                      int status;
                      int iter = 0, max_iter = 100;
                      const gsl_root_fsolver_type *T;
                      gsl_root_fsolver *s;
                      double r2 = 0;
                      double x_lo = rlast, x_hi = r;
                      T = gsl_root_fsolver_brent;
                      s = gsl_root_fsolver_alloc (T);
                      gsl_root_fsolver_set (s, &F, x_lo, x_hi);
                      do
                        {
                          iter++;
                          status = gsl_root_fsolver_iterate (s);
                          r2 = gsl_root_fsolver_root (s);
                          x_lo = gsl_root_fsolver_x_lower (s);
                          x_hi = gsl_root_fsolver_x_upper (s);
                          status = gsl_root_test_interval (x_lo, x_hi,
                                                           0, 0.001);
                        }
                      while (status == GSL_CONTINUE && iter < max_iter);
                      gsl_root_fsolver_free (s);
        //now r2 is the root
        exsol = ExactSolution(r2, phi0, dV, d2V);
        inc(1)=exsol(0);
        inc(2)=exsol(1);        
    }
 
    return inc;
}

◆ integrateProfile()

gslpp::vector< double > FindAction::integrateProfile	(	double	r0,
		double	y01,
		double	y02,
		double	dr0,
		double	epsfrac[2],
		double	epsabs[2],
		double	drmin,
		double	rmax,
		double	distance
	)

private

Definition at line 458 of file Metastability.cpp.

{
    double dr = dr0;
    gslpp::vector<double> dydr0(2,0.);
    dydr0 = dY(y01, y02, r0);
    rmax += r0;
 
//    int i=1;
//    int convergence_type = 0;
 
    gslpp::vector<double> ret_values(4,0.);
 
    do
    {
        gslpp::vector<double> r_y_drnext(4,0.);
//        r_y_drnext = rkqs(y01, y02, dydr0, r0, dr, epsfrac, epsabs);
 
 
        double r1 = r_y_drnext(0);
        double y11 = y01+r_y_drnext(1);
        double y12 = y02+r_y_drnext(2);
 
        gslpp::vector<double> dydr1 = dY(y11,y12,r1);
 
        if (r1 > rmax)
        {
            throw std::runtime_error("r > rmax");
        }
        else if (dr < drmin)
        {
            throw std::runtime_error("dr < drmin");
        }
//        else if ( fabs(y11 - distance) < 3.0*epsabs[0] && fabs(y12) < 3.0*epsabs[1] )
//        {
//            ret_values(0)=r1;
//            ret_values(1)=y11;
//            ret_values(2)=y12;
//            ret_values(3)=1;    //converged
//            break;
//        }
        else if ( y12*(y01-distance) > 0 || (y11-distance)*(y01-distance) < 0 )
        {
//            double f=y0*(1-t)**3 + 3*y1*(1-t)*(1-t)*t + 3*y2*(1-t)*t*t + y1*t**3;
//(self, y0, dy0, y1, dy1)=            y0*mt**3 + 3*y1*mt*mt*t + 3*y2*mt*t*t + y3*t**3
//            ret_values(0)=r1;
//            ret_values(1)=y11;
//            ret_values(2)=y12;
            ret_values(3)=2;    //undershoot
            break;
            
        }
//        else if ( fabs(y11 - distance) < 3.0*epsabs && fabs(y12) < 3.0*epsabs )
//        {
            
//        }
        
//          elif( y1[1]*(y0[0]-self.phi_metaMin) > 0 or (y1[0]-self.phi_metaMin)*(y0[0]-self.phi_metaMin) < 0 ):
//                f = cubicInterpFunction(y0, dr*dydr0, y1, dr*dydr1)
//                if(y1[1]*(y0[0]-self.phi_metaMin) > 0):
//                    # Extrapolate to where dphi(r) = 0
//                    x = optimize.brentq(lambda x: f(x)[1], 0, 1 )
//                    convergence_type = "undershoot"
//                else:
//                    # Extrapolate to where phi(r) = phi_metaMin
//                    x = optimize.brentq(lambda x: f(x)[0]-self.phi_metaMin, 0,1)
//                    convergence_type = "overshoot"
//                r = r0 + dr*x
//                y = f(x)
//                break
//            # Advance the integration variables
//            r0,y0,dydr0 = r1,y1,dydr1
//            dr = drnext
    }
    while(true);
//        # Check convergence for a second time. 
//        # The extrapolation in overshoot/undershoot might have gotten us within
//        # the acceptable error.
//        if (abs(y - np.array([self.phi_metaMin,0])) < 3*epsabs).all():
//            convergence_type = "converged"
//        return self._integrateProfile_rval(r, y, convergence_type)
 
    return ret_values;
}

◆ invertedpotential()

double FindAction::invertedpotential ( double x )

inlineprivate

Definition at line 59 of file Metastability.h.

    {
        return -mySUSYScalarPotential->potential(potentialcoefficientspar, x*x1par, x*x2par, x*x3par);
    }

◆ rkqs()

gslpp::vector< double > FindAction::rkqs	(	double	y01,
		double	y02,
		double	dydr0,
		double	r0,
		double	dr,
		double	epsfrac[2],
		double	epsabs[2]
	)

◆ Simpsonintegrand()

double FindAction::Simpsonintegrand	(	double	r,
		double	phi,
		double	dphi,
		double	VphiMin_i
	)

private

Definition at line 313 of file Metastability.cpp.

{
    return (0.5*dphi*dphi + deformedV(phi)-VphiMin_i) * 2.0*M_PI*M_PI*r*r*r;
}

Member Data Documentation

◆ d2Vpar

double FindAction::d2Vpar

Definition at line 42 of file Metastability.h.

◆ delta_phi_cutoffpar

double FindAction::delta_phi_cutoffpar

Definition at line 42 of file Metastability.h.

◆ dVpar

double FindAction::dVpar

Definition at line 42 of file Metastability.h.

◆ mySUSY

const SUSY& FindAction::mySUSY

private

Definition at line 51 of file Metastability.h.

◆ mySUSYScalarPotential

SUSYScalarPotential* FindAction::mySUSYScalarPotential

protected

Definition at line 48 of file Metastability.h.

◆ phi0par

double FindAction::phi0par

Definition at line 42 of file Metastability.h.

◆ potentialcoefficientspar

gslpp::vector<double> FindAction::potentialcoefficientspar

Definition at line 43 of file Metastability.h.

◆ rpar

double FindAction::rpar

Definition at line 42 of file Metastability.h.

◆ x1par

double FindAction::x1par

Definition at line 42 of file Metastability.h.

◆ x2par

double FindAction::x2par

Definition at line 42 of file Metastability.h.

◆ x3par

double FindAction::x3par

Definition at line 42 of file Metastability.h.

The documentation for this class was generated from the following files:

Detailed Description

Public Member Functions

Public Attributes

Protected Attributes

Private Member Functions

Static Private Member Functions

Private Attributes

Constructor & Destructor Documentation

◆ FindAction()

◆ ~FindAction()

Member Function Documentation

◆ computeThValue()

◆ convertToGslFunctionS()

◆ deformedV()

◆ dY()

◆ dYfunc()

◆ dYJac()

◆ ExactSolution()

◆ fPS()

◆ func()

◆ gslFunctionAdapterS()

◆ InitialConditions()

◆ integrateProfile()

◆ invertedpotential()

◆ rkqs()

◆ Simpsonintegrand()

Member Data Documentation

◆ d2Vpar

◆ delta_phi_cutoffpar

◆ dVpar

◆ mySUSY

◆ mySUSYScalarPotential

◆ phi0par

◆ potentialcoefficientspar

◆ rpar

◆ x1par

◆ x2par

◆ x3par