How to use the pint.models.timing_model.Cache.use_cache function in Pint

    @Cache.use_cache
    def d_omega_d_OM(self):
        """dOmega/dOM = 1
        """
        return np.longdouble(np.ones((len(self.t))))*u.Unit('')

    @Cache.use_cache
    def d_Pobs_d_P1(self):
        geom = (-np.sin(self.nu())*np.sin(self.om())+(np.cos(self.nu())+self.ecc())*np.cos(self.om()))
        ds = self.Doppler() * geom
        return self.t0() * (1 + ds)

    @Cache.use_cache
    def delayL2(self):
        """Second left-hand term of delay equation"""
        return (self.a1()*np.cos(self.om())*\
                np.sqrt(1-self.ecc()**2)+self.GAMMA)*np.sin(self.E())

    @Cache.use_cache
    def d_M_d_T0(self):
        """dM/dT0  this could be a generic function
        """
        return ((self.PBDOT-self.XPBDOT)*self.tt0/self.PB-1.0)*2*np.pi*u.rad/self.PB

    @Cache.use_cache
    def nhat(self):
        """nhat defined as
           T. Damour and N. Deruelle(1986)equation [51]
           nhat = n/(1-ecc*cos(E))
           n = 2*pi/PB # should here be M()
        """
        return 2.0*np.pi/self.PB.to('second')/(1-self.ecc()*self.cosEcc_A)
    @Cache.use_cache

    @Cache.use_cache
    def d_M_d_PB(self):
        """dM/dPB this could be a generic function
        """
        return 2*np.pi*u.rad*((self.PBDOT+self.XPBDOT)*self.tt0**2/self.PB**3 - \
                        self.tt0/self.PB**2)

    @Cache.use_cache
    def d_nu_d_ecc(self):
        return np.sin(self.E())**2/(self.ecc()*np.cos(self.E())-1)**2/np.sin(self.nu())

    @Cache.use_cache
    def d_E_d_PB(self):
        brack1 = self.PBDOT*self.tt0()**2/self.pb()**3 - self.tt0()/self.pb()**2
        brack2 = np.pi * 2 / (1-self.ecc()*np.cos(self.E()))
        return self.SECS_PER_DAY * brack1 * brack2

    @Cache.use_cache
    def d_E_d_PBDOT(self):
        RHS = self.d_M_d_PBDOT()
        return RHS/(1.0-np.cos(self.E())*self.ecc())
    @Cache.use_cache

    @Cache.use_cache
    def d_delayI_d_par(self,par):
        """ddelayI/dPar = dDre/dPar*delayI/Dre       [1] half
                          +Dre*d(delayI/Dre)/dPar    [2] half
           d(delayI/Dre)/dPar = -nhat*dDrep/dPar-Drep*dnhat/dPar   part (1)
                                +2*Drep*nhat*(nhat*dDrep/dPar+Drep*dnhat/dPar  part (2)
                                +1/2*nhat*(Drepp*nhat*dDre/dPar+Dre*nhat*dDrepp/dPar+2*Dre*Drepp*dnhat/dPar) part (3)
                                +Part4    part (4)
           Define x:= -1.0/2*ecct*sin(E)/(1-ecct*cos(E))
           part4 = nhat*(Dre*Drep*nhat*dx/dPar+x*(Drep*nHat*dDre/dPar+Dre*nHat*dDrep/dPar+
                   2*Dre*Drep*dnhat/dPar))
           dx/dPar = (-ecc*cos(E)*dE/dPar-sin(E)*decc/dPar
                     +ecc*sin(E)*(-cos(E)*decc/dPar+ecc*sin(E)*dE/dPar)/(1-ecc*cos(E)))/(2*(1-ecc*cos(E)))
        """
        e = self.ecc()
        sE = self.sinEcc_A
        cE = self.cosEcc_A