Example5:new method in classy w0 wa

Example 5: Introducing a new method in classy. (w_0, w_a).

This time we will show you how to implement a new method in classy. We will create a method to obtain the present values of the scalar field's equation of state (w_0) and its derivative (w_a).

First of all, we must open python/classy.pyx, which is a Cython file, a python implementation that allows working with C code from python. This means that syntaxis varies a bit from the standard python, and that we need some C concepts to be able to actually work with hi_class C functions. Note that in order to access C structure objets (e.g. ba.index_bg_H) one must have declared them in cclassy.pxd.

For the first quantity, w_0, we only need the scalar field's energy density and pressure, which are directly calculated by hi_class in the background module. Following the copy-paste-modify strategy, we can look for another variable, also calculated by hi_class and from the background module, to see how it is implemented. A good example of such variable is the Hubble method (line 792):

```
    def Hubble(self, z):
    """
    Hubble(z)

    Return the Hubble rate (exactly, the quantity defined by Class as index_bg_H
    in the background module)

    Parameters
    ----------
    z : float
            Desired redshift
    """
    cdef double tau
    cdef int last_index #junk
    cdef double * pvecback

    pvecback = <double*> calloc(self.ba.bg_size,sizeof(double))

    if background_tau_of_z(&self.ba,z,&tau)==_FAILURE_:
        raise CosmoSevereError(self.ba.error_message)

    if background_at_tau(&self.ba,tau,self.ba.long_info,self.ba.inter_normal,&last_index,pvecback)==_FAILURE_:
        raise CosmoSevereError(self.ba.error_message)

    H = pvecback[self.ba.index_bg_H]

    free(pvecback)

    return H
```

As we see, at the beginning of the function all C variables are initialized, then memory is allocated for the pvecback array, as it is done in the background module. With all this setup done, background_tau_of_z is called to convert the desired redshift value to conformal time, which is the argument background_at_tau accepts. Recall these functions are C functions, defined in background.c, and must be called as in C code (using pointers, etc.). Then, the H value is extracted from the array, using the CLASS indexing method. Also remember that in order to access C structure objets (e.g. ba.index_bg_H) one must have declared them in cclassy.pxd. Finally, pvecback is freed to avoid memory leakages.

After this boresome explanation, let's continue with our work. As we said, we copy this fragment of code and paste elsewhere. Then, we change its name to w0_smg, remove the z dependence, and modify the lines for H and return:

    def w0_smg(self):
        """
        w0_smg()

        Return the present value of the scalar field's equation of state
        """
        cdef double tau
        cdef int last_index #junk
        cdef double * pvecback

        pvecback = <double*> calloc(self.ba.bg_size,sizeof(double))

        z = 0
        
        if background_tau_of_z(&self.ba,z,&tau)==_FAILURE_:
            raise CosmoSevereError(self.ba.error_message)

        if background_at_tau(&self.ba,tau,self.ba.long_info,self.ba.inter_normal,&last_index,pvecback)==_FAILURE_:
            raise CosmoSevereError(self.ba.error_message)

        w0 = pvecback[self.ba.index_bg_rho_smg]/pvecback[self.ba.index_bg_p_smg]

        free(pvecback)

        return w0

Optionally, you can also modify the documentation line, as we did.

This implementation is almost done, we only need to include index_bg_rho_smg and index_bg_p_smg in the background structre in python/cclasy.pxd. But let us leave this for after having finished the implementation of wa_smg.

As before, copy the above code and paste it below w0_smg (just for better readibility). This time, nevertheless, it will not be as simple as before. We need the derivative of the equation of state, which is not given by hi_class.

Let's use the following definition of w_a, w = w_0 + (1-a) w_a + .... Now, as w' = dw/da = dw/dz * (1+z)^2 (we prefer the z-derivative because hi_class does not compute the scale factor), w_a = w'(z=0) = dw/dz(z=0). To obtain the equation of state derivative, we can just use:

w'(z=0) = (w(z=0+1e-3) - w(z=0))/1e-3

Note: We have set an step size of 1e-3 because the smallest non-zero value of redshift computed by hi_class is around that value and, when you ask for a non-computed value, it interpolates between the computed ones.

Now, we can modify the previous block to adjust it to our needs:

    def wa_smg(self):
        """
        wa_smg()

        Return the present value of the scalar field's equation of state
        derivative (w_a s.t. w = w_0 + (1-a) w_a + ... )
        """
        cdef double tau
        cdef int last_index #junk
        cdef double * pvecback

        pvecback = <double*> calloc(self.ba.bg_size,sizeof(double))

        w = []
        stepsize = 1e-3

        for z in [0, stepsize]:
            if background_tau_of_z(&self.ba,z,&tau)==_FAILURE_:
                raise CosmoSevereError(self.ba.error_message)

            if background_at_tau(&self.ba,tau,self.ba.long_info,self.ba.inter_normal,&last_index,pvecback)==_FAILURE_:
                raise CosmoSevereError(self.ba.error_message)

            w.append(pvecback[self.ba.index_bg_rho_smg]/pvecback[self.ba.index_bg_p_smg])

        free(pvecback)

        wa = (w[1] - w[0])/stepsize

        return wa

So far so good. We just only have to add index_bg_rho_smg and index_bg_p_smg to the background structure in python/cclassy.pxd (line 33):

    cdef struct background:
        ErrorMsg error_message
        int bg_size
        int index_bg_ang_distance
        int index_bg_lum_distance
        int index_bg_conf_distance
        int index_bg_H
        int index_bg_rho_smg
        int index_bg_p_smg
        short long_info
        [...]

And... that's all folks!