def solarAzimuth(_day, _hour = None):
    """Solar azimuth in radians"""  
    def cos_azmth(_hourAngle):
        """Definition of the cosine of the azimuth replacing the altitude by 
        it formula"""
        return (np.sin(np.arcsin(np.cos(latitude)*np.cos(declination(_day,radianes))* \
        np.cos(_hourAngle)+np.sin(latitude)*np.sin(declination(_day,radianes))))*np.sin(latitude)- \
        np.sin(declination(_day,radianes))) \
        / (np.cos(np.arcsin(np.cos(latitude)*np.cos(declination(_day,radianes))* \
        np.cos(_hourAngle)+np.sin(latitude)*np.sin(declination(_day,radianes))))*np.cos(latitude))

    if _hour is None:
        for i in range(np.size(hourAngle)):
            derivada[i] = sc.derivative(cos_azmth,hourAngle[i],dx=1e-6)
            _signo = -np.sign(derivada[i])
            azimuth[i] =_signo * np.arccos(cos_azmth(hourAngle[i]))
            azimuth_deg[i] = azimuth[i]*180/np.pi
            
    else:
         _hourAngle = (12-int(_hour))*np.pi/12
         _derivada = sc.derivative(cos_azmth,_hourAngle,dx=1e-6)
         _signo = - np.sign(_derivada)
         if cos_azmth(_hourAngle) > 1:
             _azimuth = _signo * np.arccos(1) * 180 / np.pi
         else:
             _azimuth = _signo * np.arccos(cos_azmth(_hourAngle)) * 180 / np.pi
         return _azimuth
    
    return

		def fisher(theta, i, j = None):
			"""
			Fisher information using the first order derivative

			:param theta: the theta of the density
			:param i: The ith component of the diagonal of the fisher information matrix will be returned (if j is None)
			:param j: The i,j th component of the fisher information matrix will be returned
			"""

			#Bring it in a form that we can derive
			fh = lambda ti, t0, tn, x: np.log(self.density(x, list(t0) + [ti] + list(tn)))

			# The derivative
			f_d_theta_i = lambda x: derivative(fh, theta[i], dx=1e-5, n=1, args=(theta[0:i], theta[i + 1:], x))

			if j is not None:
				f_d_theta_j = lambda x: derivative(fh, theta[j], dx=1e-5, n=1, args=(theta[0:j], theta[j + 1:], x))
				f = lambda x: np.float128(0) if fabs(self.density(x, theta)) < 1e-5 else f_d_theta_i(x) * f_d_theta_j(x) * self.density(x, theta)
			else:
				# The function to integrate
				f = lambda x: np.float128(0) if fabs(self.density(x, theta)) < 1e-5 else f_d_theta_i(x) ** 2 * self.density(x, theta)


			#First order
			result = integrate(f, self.support)
			return result

def newton(fn=0, fp=0, x=0.1, mn=0):
    t = time()
    test = 0
    if not fn and not fp:
        from scipy.misc import derivative

        func = input("equation(use x as unknown):")
        fn = lambda x: eval(func)
        fp = lambda x: derivative(fn, x)
        print("div")
    elif not fp:
        from scipy.misc import derivative

        fp = lambda x: derivative(fn, x)
    if mn:
        while fn(x, mn) != 0:
            if fp(x, mn) == 0:
                x += 1
                print("error")
                test += 1
                print(test)
                break
            x, lastx = x - (fn(x, mn) / fp(x, mn)), x
    else:
        while fn(x) != 0 or rerr > 10 ** -12:
            x, lastx = x - (fn(x) / fp(x)), x
            rerr = abs(lastx - x) / x
    return x

    def eval_function(self, individual, radius, representatives):
        distanceSum = 0.001+sum([self.share_function(self.distance(individual, r), 1, radius) for r in representatives])
        # print "distanceSum: %s" % distanceSum
        calcWithY = partial(self.function_object.calculateWithXY, individual[0])
        calcWithX = partial(self.function_object.calculateWithYX, individual[1])
        derivY = abs(misc.derivative(calcWithY, individual[1], dx=1e-6))
        derivX = abs(misc.derivative(calcWithX, individual[0], dx=1e-6))

        return self.function_object.calculate(individual)/distanceSum, distanceSum*derivX, distanceSum*derivY

def calc_curvature(func, d_point=0):
    #曲率半径、曲率を求める
    #曲率が大きいほど関数の曲がり具合が大きい

    dx1 = scmisc.derivative(func, d_point, n=1, dx=1e-6) #dxを指定しないとダメ
    dx2 = scmisc.derivative(func, d_point, n=2, dx=1e-6)
    R = (1.0 + (dx1**2.0))**(3.0/2.0) / np.fabs(dx2) #曲率半径
    curvature_of_func = 1.0 / R

    return curvature_of_func

def C_Meng(T, Tc, Pc, D, B):
    r"""Calculate the 3rd virial coefficient using the Meng-Duan-Li correlation

    Parameters
    ----------
    T : float
        Temperature [K]
    Tc : float
        Critical temperature [K]
    Pc : float
        Critical pressure
    D : float
        dipole moment [debye]
    B : list
        Second virial coefficient tuple with B, ∂B/∂T, ∂²B/∂T²

    Returns
    -------
    C : float
        Third virial coefficient [m⁶/mol²]
    C1 : float
        T(∂C/∂T) [m⁶/mol²]
    C2 : float
        T²(∂²C/∂T²) [m⁶/mol²]

    Examples
    --------
    Selected date from Table 2, pag 74 for neon

    >>> from lib.mEoS import Ne
    >>> D = Ne.momentoDipolar.Debye
    >>> B = B_Meng(273.15, Ne.Tc, Ne.Pc, Ne.f_acent, D)
    >>> "%.2f" % (C_Meng(273.15, Ne.Tc, Ne.Pc, D, B)[0]*1e-5)
    '0.24'

    References
    ----------
    [7]_ Meng, L., Duan, Y.Y. Li, L. Correations for Second and Third Virial
        Coefficients of Pure Fluids. Fluid Phase Equilibria 226 (2004) 109-120
    """
    mur = D**2*Pc/1.01325/Tc**2

    def f(T, *args):
        B = args[-1]
        Br = B*Pc/R/Tc
        Tr = T/Tc
        f0 = 1094.051 - 3334.145/Tr**0.1 + 3389.848/Tr**0.2 - 1149.58/Tr**0.3
        f1 = (2.0243-0.85902/Tr)*1e-10
        return 5.476e-3 + (Br-0.0936)**2*(f0+mur**4*f1)

    C = f(T, B[0])
    C1 = derivative(f, T, n=1, args=(T, B[1]))
    C2 = derivative(f, T, n=2, args=(T, B[2]))

    return C, C1, C2

    def gruneisen_parameter(self, temperature, volume):
        """
        Slater-gamma formulation(the default):
            gruneisen paramter = - d log(theta)/ d log(V)
                               = - ( 1/6 + 0.5 d log(B)/ d log(V) )
                               = - (1/6 + 0.5 V/B dB/dV),
                                    where dB/dV = d^2E/dV^2 + V * d^3E/dV^3

        Mie-gruneisen formulation:
            Eq(31) in doi.org/10.1016/j.comphy.2003.12.001
            Eq(7) in Blanco et. al. Joumal of Molecular Structure (Theochem)
                368 (1996) 245-255
            Also se J.P. Poirier, Introduction to the Physics of the Earth’s
                Interior, 2nd ed. (Cambridge University Press, Cambridge,
                2000) Eq(3.53)

        Args:
            temperature (float): temperature in K
            volume (float): in Ang^3

        Returns:
            float: unitless
        """
        if isinstance(self.eos, PolynomialEOS):
            p = np.poly1d(self.eos.eos_params)
            # first derivative of energy at 0K wrt volume evaluated at the
            # given volume, in eV/Ang^3
            dEdV = np.polyder(p, 1)(volume)
            # second derivative of energy at 0K wrt volume evaluated at the
            # given volume, in eV/Ang^6
            d2EdV2 = np.polyder(p, 2)(volume)
            # third derivative of energy at 0K wrt volume evaluated at the
            # given volume, in eV/Ang^9
            d3EdV3 = np.polyder(p, 3)(volume)
        else:
            func = self.ev_eos_fit.func
            dEdV = derivative(func, volume, dx=1e-3)
            d2EdV2 = derivative(func, volume, dx=1e-3, n=2, order=5)
            d3EdV3 = derivative(func, volume, dx=1e-3, n=3, order=7)

        # Mie-gruneisen formulation
        if self.use_mie_gruneisen:
            p0 = dEdV
            return (self.gpa_to_ev_ang * volume *
                    (self.pressure + p0 / self.gpa_to_ev_ang) /
                    self.vibrational_internal_energy(temperature, volume))

        # Slater-gamma formulation
        # first derivative of bulk modulus wrt volume, eV/Ang^6
        dBdV = d2EdV2 + d3EdV3 * volume
        return -(1./6. + 0.5 * volume * dBdV /
                 FloatWithUnit(self.ev_eos_fit.b0_GPa, "GPa").to("eV ang^-3"))

def radius_of_curvature(theta, func_R, *args_for_func_R):
    """Returns R_c = (R^2 + R'^2)^{3/2} / |R^2 + 2 R'^2 - R R''| 

Uses numerical differentiation.  NOT RECOMMENDED SINCE NOT ACCURATE ON
THE AXIS.  Use `axis_Rc` instead.

    """
    R = func_R(theta, *args_for_func_R)
    dR = derivative(func_R, theta,
                    dx=DX_FOR_NUMERICAL_DERIVATIVE, args=args_for_func_R)
    d2R = derivative(func_R, theta, n=2,
                     dx=DX_FOR_NUMERICAL_DERIVATIVE, args=args_for_func_R)
    return (R**2 + dR**2)**1.5 / np.abs(R**2 + 2*dR**2 - R*d2R)

def C_Orbey_Vera(T, Tc, Pc, w):
    """Calculate the third virial coefficient using the Orbey-Vera correlation

    Parameters
    ----------
    T : float
        Temperature [K]
    Tc : float
        Critical temperature [K]
    Pc : float
        Critical pressure
    w : float
        Acentric factor [-]

    Returns
    -------
    C : float
        Third virial coefficient [m⁶/mol²]
    C1 : float
        T(∂C/∂T) [m⁶/mol²]
    C2 : float
        T²(∂²C/∂T²) [m⁶/mol²]

    Examples
    --------
    Selected points from Table 2 of paper

    >>> from lib.mEoS.Benzene import Benzene as Bz
    >>> "%.1f" % (C_Orbey_Vera(0.877*Bz.Tc, Bz.Tc, Bz.Pc, Bz.f_acent)[0]*1e9)
    '41.5'
    >>> "%.1f" % (C_Orbey_Vera(1.019*Bz.Tc, Bz.Tc, Bz.Pc, Bz.f_acent)[0]*1e9)
    '35.8'

    References
    ----------
    [4]_ Orbey, H., Vera, J.H.: Correlation for the third virial coefficient
        using Tc, Pc and ω as parameters, AIChE Journal 29, 107 (1983)
    """
    def f(T):
        Tr = T/Tc
        g0 = 0.01407+0.02432/Tr**2.8-0.00313/Tr**10.5
        g1 = -0.02676+0.0177/Tr**2.8+0.04/Tr**3-0.003/Tr**6-0.00228/Tr**10.5
        g = g0+w*g1
        return g*R**2*Tc**2/Pc**2

    C = f(T)
    C1 = derivative(f, T, n=1)
    C2 = derivative(f, T, n=2)

    return C, C1, C2

def derivative_check(x,popt,left,right):
	if len(left)==0 or len(right)==0:
		return True
	if len(left)==1 or len(right)==1:
		return True
	maxlx=np.max(left[:,0])
	minlx=np.min(left[:,0])
	maxly=np.max(left[:,1])
	minly=np.min(left[:,1])
	maxrx=np.max(right[:,0])
	minrx=np.min(right[:,0])
	maxry=np.max(right[:,1])
	minry=np.min(right[:,1])
	t1=-1
	t2=-1
	pend1=[]
	pend2=[]
	pend1.append(deriv.derivative(fu,minlx,n=1,args=popt))
	pend1.append(deriv.derivative(fu,maxlx,n=1,args=popt))
	pend1.append(deriv.derivative(fu,minrx,n=1,args=popt))
	pend1.append(deriv.derivative(fu,maxrx,n=1,args=popt))
	sign=pend1[0]
	for i in pend1:
		if i<0:
			return False
	pend1=[]
	pend1.append(deriv.derivative(fu,minlx,n=2,args=popt))
	pend1.append(deriv.derivative(fu,maxlx,n=2,args=popt))
	pend1.append(deriv.derivative(fu,minrx,n=2,args=popt))
	pend1.append(deriv.derivative(fu,maxrx,n=2,args=popt))
	sign=pend1[0]
	for i in pend1:
		if (sign>0 and i<0) or (sign<0 and i>0):
			return False
	return True

    def get_corr_pred(self, eps, d_eps, sig, t_n, t_n1, alpha, q, kappa):
        #         g = lambda k: 0.8 - 0.8 * np.exp(-k)
        #         g = lambda k: 1. / (1 + np.exp(-2 * k + 6.))
        n_e, n_ip, n_s = eps.shape
        D = np.zeros((n_e, n_ip, 3, 3))
        D[:, :, 0, 0] = self.E_m
        D[:, :, 2, 2] = self.E_f
        sig_trial = sig[:, :, 1]/(1-self.g(kappa)) + self.E_b * d_eps[:,:, 1]
        xi_trial = sig_trial - q
        f_trial = abs(xi_trial) - (self.sigma_y + self.K_bar * alpha)
        elas = f_trial <= 1e-8
        plas = f_trial > 1e-8
        d_sig = np.einsum('...st,...t->...s', D, d_eps)
        sig += d_sig

        d_gamma = f_trial / (self.E_b + self.K_bar + self.H_bar) * plas
        alpha += d_gamma
        kappa += d_gamma
        q += d_gamma * self.H_bar * np.sign(xi_trial)
        w = self.g(kappa)

        sig_e = sig_trial - d_gamma * self.E_b * np.sign(xi_trial)
        sig[:, :, 1] = (1-w)*sig_e

        E_p = -self.E_b / (self.E_b + self.K_bar + self.H_bar) * derivative(self.g, kappa, dx=1e-6) * sig_e \
            + (1 - w) * self.E_b * (self.K_bar + self.H_bar) / \
            (self.E_b + self.K_bar + self.H_bar)

        D[:, :, 1, 1] = (1-w)*self.E_b*elas + E_p*plas

        return sig, D, alpha, q, kappa

  def _testCompareToExplicitDerivative(self, dtype):
    """Compare to the explicit reparameterization derivative.

    Verifies that the computed derivative satisfies
    dsample / dalpha = d igammainv(alpha, u) / dalpha,
    where u = igamma(alpha, sample).

    Args:
      dtype: TensorFlow dtype to perform the computations in.
    """
    delta = 1e-3
    np_dtype = dtype.as_numpy_dtype
    try:
      from scipy import misc  # pylint: disable=g-import-not-at-top
      from scipy import special  # pylint: disable=g-import-not-at-top

      alpha_val = np.logspace(-2, 3, dtype=np_dtype)
      alpha = constant_op.constant(alpha_val)
      sample = random_ops.random_gamma([], alpha, np_dtype(1.0), dtype=dtype)
      actual = gradients_impl.gradients(sample, alpha)[0]

      (sample_val, actual_val) = self.evaluate((sample, actual))

      u = special.gammainc(alpha_val, sample_val)
      expected_val = misc.derivative(
          lambda alpha_prime: special.gammaincinv(alpha_prime, u),
          alpha_val, dx=delta * alpha_val)

      self.assertAllClose(actual_val, expected_val, rtol=1e-3, atol=1e-3)
    except ImportError as e:
      tf_logging.warn("Cannot use special functions in a test: %s" % str(e))

def func_2nd_derivative_numerical(x):
    value = 1.
### START YOUR CODE HERE ###
    value = misc.derivative(func,x,1e-5,2);

#### END YOUR CODE HERE ####
    return value

def func_first_derivative(x):
    value = 1.
### START YOUR CODE HERE ###
    value = misc.derivative(func,x,1e-5);

#### END YOUR CODE HERE ####
    return value

    def simulateDeltaWithTime(self, timeSteps):
        opt=copy.copy(self.option)
        spotPrice=copy.copy(opt.S0)
        mechanism=self.mechanism
        
        plotDf=pd.DataFrame(np.zeros([len(timeSteps), 3]), columns=['TimeToMaturity', 'DA_Delta', 'BLS_Delta'])
        plotDf['TimeToMaturity']=timeSteps
        
        def optPrice(assetPrice):
            opt.S0=assetPrice

            orders=self.traders.getOrders(opt, self.numOrders)
            orders=mechanism.clearOrders(orders)
       
            return mechanism.getPrice(orders)
        
        for i, p in enumerate(timeSteps):
            opt.T=p            
            plotDf.ix[i]['DA_Delta']=derivative(optPrice, x0=spotPrice, dx=20)
            
            opt.S0=spotPrice
            plotDf.ix[i]['BLS_Delta']=opt.delta()
            
            print i
     
        return plotDf

 def simulateGamma(self, assetPrices):
     opt=copy.copy(self.option)
     mechanism=self.mechanism
     
     plotDf=pd.DataFrame(np.zeros([len(assetPrices), 3]), columns=['AssetPrice', 'DA_Gamma', 'BLS_Gamma'])
     plotDf['AssetPrice']=assetPrices.values
     
     def optPrice(assetPrice):
         opt.S0=assetPrice
         orders=self.traders.getOrders(opt, self.numOrders)
         orders=mechanism.clearOrders(orders)
         
         return mechanism.getPrice(orders)
     
     for i, p in enumerate(assetPrices.values):
         plotDf.ix[i]['DA_Gamma']=derivative(optPrice, x0=p,dx=20, n=2)
         
         opt.S0=p            
         plotDf.ix[i]['BLS_Gamma']=opt.gamma()
         
         print plotDf.ix[i]['DA_Gamma'], plotDf.ix[i]['BLS_Gamma']
         
         print i
  
     return plotDf

    def fitBinCuts(self):
        print("INFO: Fit bin edges iteratively")
        self.splines =[lambda x: 0.5]*21
        fitboundaries = lambda n: [-0.001,-1.5+np.sqrt(n)*0.1]
        converged = True
        for n in range(2,16):
            scan_c = np.logspace(*fitboundaries(n))
            roots = []
            guess = min(0.095 * n,0.95)
            scan_copy = []
            for c in scan_c:

                x_0_p = guess * c

                d_sign_dx = lambda x: derivative(lambda y: self.signN(y,c,n),x,dx=0.0001*c)
                opt = optimize.root(d_sign_dx,x_0_p)
                x_0 = opt.x[0]
                if opt.success and 0 < x_0 < c:
                    guess = x_0/c
                    roots.append(x_0)
                    scan_copy.append(c)
                else:
                    print "root not found" 
                    
            try:
                q = UnivariateSpline(np.log10(scan_copy[::-1]),np.log10(np.array(roots)[::-1]),s=0.0001)
                self.splines[n] = q
                print "parameters for n="+str(n)+" found"
            except:
                print "ERROR: fit failed"
                converged = False
                break

        return converged

def ACM(Cf,i):
	
    q=d.a_3-d.a_4
    L=(Cf/111.)
    e_0=((d.a_7*L**2)/(L**2+d.a_9))
    g_c=((np.abs(float(d.phi_d)))**(d.a_10))/(0.5*float(max(dat.ta_fill[(dat.jd==d.D[i])]-min(dat.ta_fill[(dat.jd==d.D[i])])+d.a_6*float(d.R_tot))))
    p=(((d.a_1*d.N*L)/g_c)*np.exp(d.a_8*float(max(dat.ta_fill[(dat.jd==d.D[i])]))))
    C_i=(0.5*(d.C_a+q-p+np.sqrt((d.C_a+q-p)**2-4*(d.C_a*q-p*d.a_3))))
    delta=(-0.408*np.cos(((360*(float(d.D[i])+10)*np.pi)/(365*180))))
    s=(24*np.arccos((-np.tan(d.bigdelta)*np.tan(delta)))/np.pi)
	
    def GPPfunc(Cf):
		
        #L=(Cf/111.)
        #e_0=((d.a_7*L**2)/(L**2+d.a_9))
        #GPP=((e_0*float(d.I[i])*g_c*(d.C_a-C_i)*(d.a_2*s+d.a_5))/(e_0*float(d.I[i])+g_c*(d.C_a-C_i))) #0.2*float(d.I[i])*(1-np.exp(-0.5*Cf/111.))
        #GPP_diff=((2*float(d.I[i])*d.a_7*d.a_9*((111.0*g_c*(d.C_a-C_i))**2)*(d.a_2*s+d.a_5)*Cf)/(((g_c*(d.C_a-C_i)+d.a_7*float(d.I[i]))*Cf**2+d.a_9*(111.0**2)*g_c*(d.C_a-C_i))**2))  #(0.2*0.5*float(d.I[i])/111.)*np.exp(-0.5*Cf/111.)
        GPP= (d.a_7*(Cf**2)*float(d.I[i])*g_c*(d.C_a-C_i)*(d.a_2*s+d.a_5))/((d.a_7*float(d.I[i])+g_c*(d.C_a-C_i))*(Cf**2)+d.a_9*(111.**2)*g_c*(d.C_a-C_i))
        
        return GPP 

    def GPPdifffunc(Cf):
		
        GPP_diff=(24642.*float(d.I[i])*d.a_7*Cf*(g_c**2)*((d.C_a-C_i)**2)*(d.a_2*s+d.a_5)*d.a_9)/((Cf**2+12321.*d.a_9)*(d.C_a-C_i)*g_c+float(d.I[i])*d.a_7*Cf**2)**2 #(2*Cf*float(d.I[i])*d.a_7*d.a_9*(111.**2)*(g_c**2)*((d.C_a-C_i)**2)*(d.a_2*s+d.a_5))/((d.a_7*float(d.I[i])+g_c*(d.C_a-C_i))*(Cf**2)+d.a_9*(111.**2)*g_c*(d.C_a-C_i))**2
        
        return GPP_diff
		
    GPP = GPPfunc(Cf)
    GPP_diff=  GPPdifffunc(Cf)
    GPP_diff2 = misc.derivative(GPPfunc,Cf) 
    return GPP, GPP_diff

def ACM(Cf,i):
	
    q=d.a_3-d.a_4
    #L=(Cf/111.)
    #e_0=((d.a_7*L**2.)/(L**2.+d.a_9))
    g_c=(((np.abs(float(d.phi_d[i])))**(d.a_10))/(0.5*float(d.T_range[i])+d.a_6*float(d.R_tot[i])))
    #p=np.exp(d.a_8*float(d.T_max[i]))*((d.a_1*d.N*L)/g_c)
    #C_i=(0.5*(d.d.C_a+q-p+np.np.sqrt((d.d.C_a+q-p)**2.-4.*(d.d.C_a*q-p*d.a_3))))
    delta=(-0.408*np.cos(((360.*(float(d.D[i])+10.)*np.pi)/(365.*180.))))
    s=(24.*np.arccos((-np.tan(d.bigdelta)*np.tan(delta)))/np.pi)
    psi=np.exp(d.a_8*float(d.T_max[i]))*((d.a_1*d.N)/(g_c*111.))
    a=np.exp(d.a_8*float(d.T_max[i]))*((d.a_1*d.N)/(g_c*111.))
    phi=d.a_7*float(d.I[i])*g_c
    I=float(d.I[i])
	
    def GPPfunc(Cf):
        GPP=.5*phi*Cf**2*(d.C_a-q+psi*Cf-np.sqrt((d.C_a+q-psi*Cf)**2-4.*(d.C_a*q-d.a_3*psi*Cf)))*(d.a_2*s+d.a_5)/(I*d.a_7*Cf**2+.5*g_c*(Cf**2+d.a_9*111.**2)*(d.C_a-q+psi*Cf-np.sqrt((d.C_a+q-psi*Cf)**2-4.*(d.C_a*q-d.a_3*psi*Cf))))        
        return GPP 

    def GPPdifffunc(Cf):
        GPP_diff=(((2.*(d.a_2*s+d.a_5))*(-49284.*q*d.a_9*psi*Cf+49284.*psi*Cf*d.a_3*d.a_9-49284.*d.C_a*d.a_9*q+24642.*psi**2*Cf**2*d.a_9+24642.*d.a_9*d.C_a**2+24642.*d.a_9*q**2)*phi*Cf*g_c+(2.*I)*(d.a_2*s+d.a_5)*psi*Cf**4*d.a_7*phi)*np.sqrt(psi**2*Cf**2+(-2.*d.C_a-2.*q+4.*d.a_3)*psi*Cf+d.C_a**2-2.*d.C_a*q+q**2)+(-(49284.*(d.a_2*s+d.a_5))*psi**3*d.a_9*phi*Cf**4+(2.*(d.a_2*s+d.a_5))*(24642.*d.C_a-98568.*d.a_3+73926.*q)*phi*psi**2*d.a_9*Cf**3+(2.*(d.a_2*s+d.a_5))*(-98568.*d.a_3*d.C_a+98568.*d.a_3*q+49284.*d.C_a*q-73926.*q**2+24642.*d.C_a**2)*phi*psi*d.a_9*Cf**2+(2.*(d.a_2*s+d.a_5))*(-73926.*d.C_a*q**2-24642.*d.C_a**3+24642.*q**3+73926.*d.C_a**2*q)*phi*d.a_9*Cf)*g_c-(2.*I)*(d.a_2*s+d.a_5)*psi**2*d.a_7*phi*Cf**5+(2.*(d.a_2*s+d.a_5))*((1.*I)*q*d.a_7+(1.*I)*d.C_a*d.a_7-(2.*I)*d.a_3*d.a_7)*phi*psi*Cf**4)/(np.sqrt(psi**2*Cf**2+(-2.*d.C_a-2.*q+4.*d.a_3)*psi*Cf+d.C_a**2-2.*d.C_a*q+q**2)*((2.*I)*d.a_7*Cf**2+g_c*Cf**2*d.C_a-1.*g_c*Cf**2*q+g_c*Cf**3*psi-1.*g_c*Cf**2*np.sqrt(psi**2*Cf**2+(-2.*d.C_a-2.*q+4.*d.a_3)*psi*Cf+d.C_a**2-2.*d.C_a*q+q**2)+12321.*g_c*d.a_9*d.C_a-12321.*g_c*d.a_9*q+12321.*g_c*d.a_9*psi*Cf-12321.*g_c*d.a_9*np.sqrt(psi**2*Cf**2+(-2.*d.C_a-2.*q+4.*d.a_3)*psi*Cf+d.C_a**2-2.*d.C_a*q+q**2))**2)    
        return GPP_diff
		
    GPP = GPPfunc(Cf)
    GPP_diff=  GPPdifffunc(Cf)
    GPP_diff2 = misc.derivative(GPPfunc,Cf) 
    return GPP, GPP_diff #, GPP_diff2

 def core_energy_balance(self, T_cmb, core_flux):
     pc = self.params.core
     core_surface_area = self.outer_surface_area
     if self.params.source == 'Stevenson_1983' :
         inner_core_surface_area = np.power(self.inner_core_radius(T_cmb), 2.0) * 4. * np.pi
         dRi_dTcmb = 0.
         try:
             dRi_dTcmb = derivative( self.inner_core_radius, T_cmb, dx=1.0)
         except ValueError:
             pass
     elif self.params.source == 'Driscoll_2014' :
         dRi_dTcmb = 0.
         inner_core_surface_area = 0.
         # Eqn 29 & 30 from the Driscoll_2014 paper. Note that Eqn 32 in the paper has an error and the form in the code is correct
         sqrt_term_num = (pc.Dn/self.params.R_c0)**2.*np.log(pc.TFe/T_cmb) - 1.
         sqrt_term_den = 2.*(1. - 1./(3.*pc.gamma_c))*(pc.Dn/pc.DFe)**2. - 1.
         if sqrt_term_num/sqrt_term_den > 0 :
             R_ic = self.params.R_c0*np.sqrt(sqrt_term_num/sqrt_term_den)
             print('\n\n R_ic={0:.3f}'.format(R_ic/1e3))
             inner_core_surface_area = np.power(R_ic, 2.0) * 4. * np.pi
             dRi_dTcmb = -1.*((self.params.R_c0/2./T_cmb)*(pc.Dn/self.params.R_c0)**2.)/(sqrt_term_num*sqrt_term_den)
     else :
         raise ValueError('parameter class is not recognized')
     # print('\n\n ratio={0:.1f}'.format(inner_core_surface_area/core_surface_area))
     thermal_energy_change = pc.rho*pc.C*self.volume*pc.mu
     latent_heat = -pc.L_Eg * pc.rho * inner_core_surface_area * dRi_dTcmb
     dTdt = -core_flux * core_surface_area / (thermal_energy_change-latent_heat)
     dTdt = -core_flux * core_surface_area / (thermal_energy_change)
     return dTdt