f2 = ((1/2)*r_c*c**2+(1/4)*u_c*c**4+(1/6)*v_c*c**6+(1/2)*r_p*p**2+(1/4)*u_p*p**4-E*p-gamma*c*p-e*c**2*p**2/2)
nsimplify(f2)

"""
Explanation: $f_{2}(c,p) = \dfrac{1}{2}r_{c}c^{2}+\dfrac{1}{4}u_{c}c^{4}+\dfrac{1}{6}v_{c}c^{6}+\dfrac{1}{2}r_{p}p^{2}+\dfrac{1}{4}u_{p}p^{4}-\gamma cp-\dfrac{1}{2}ec^{2}p^{2}-Ep$
End of explanation
"""


P, C, w, rho, beta, tau, Epr = symbols('P C w rho beta tau E^{\prime}')

fP = nsimplify(f2.subs(p,w*Pprm))
fP

fP = nsimplify(fP.subs(w,sqrt(r_p/u_p)))
fP

fP = expand(fP/(r_p**2/u_p),r_p)
fP

fP = fP.subs([(E,Epr/sqrt(u_p)*r_p**(3/2)),(c,cprm*sqrt(r_p/e))])
fP

fP = fP.subs([(cprm,thetaprm),(r_c,rtheta*e*r_p/u_p),(u_c,gammaprm*e**2/u_p),(v_c,rhoprm*e**3/(u_p*r_p)),(gamma,betaprm*r_p*sqrt(e/u_p))])
fP

fhelp = nsimplify((1/2)*rtheta*thetaprm**2+(1/4)*gammaprm*thetaprm**4+(1/6)*rhoprm*thetaprm**6+(1/2)*Pprm**2+(1/4)*Pprm**4
        +(1/6)*sigmaprm*Pprm**6-betaprm*thetaprm*Pprm-(1/2)*Pprm**2*thetaprm**2-Eprm*Pprm)
fhelp

Ep = solve(fhelp.diff(Pprm),Eprm)[0]
Ep

fhelp.diff(thetaprm)

P_min = solve(fhelp.diff(thetaprm),Pprm)[0]
P_min

Ep = Ep.subs(Pprm,P_min)
Ep

series(Ep,thetaprm,n=7)

# Series expansion of $E$ out to $\mathcal{O(\theta^{\prime 7})}$:

# $\beta E = \theta^{\prime}\left(-\beta^{\prime 2} + r_{\theta}^{\prime}\right) + \theta^{\prime 3}\left(\gamma^{\prime} - r_{\theta}^{\prime} + \dfrac{r_{\theta}^{\prime 3}}{\beta^{\prime 2}} - \dfrac{r_{\theta}^{\prime 2}}{\beta^{\prime 2}}\right) + \theta^{\prime 5}\left(-\gamma^{\prime} + \rho^{\prime} + \dfrac{3\gamma^{\prime}r_{\theta}^{\prime 2}}{\beta^{\prime 2}} - \dfrac{2\gamma^{\prime}r_{\theta}}{\beta^{\prime 2}} + \dfrac{r_{\theta}^{\prime 2}}{\beta^{\prime 2}} - \dfrac{3r_{\theta}^{\prime 4}}{\beta^{\prime 4}} + \dfrac{2r_{\theta}^{\prime 3}}{\beta^{\prime 4}}\right)$

# The coefficients of $E$ i.t.o. $a$ where $r_{\theta}^{\prime} = \beta^{\prime 2} + a$:

# $B(a) = \dfrac{a^{3}}{\beta^{\prime 2}} + a^{2}\left(3 - \dfrac{1}{\beta^{\prime 2}}\right) + a(3\beta^{\prime 2} - 3) + \beta^{\prime 4} - 2\beta^{\prime 2} + \gamma^{\prime}$

# $C(a) = -\dfrac{3a^{4}}{\beta^{\prime 4}} + a^{3}\left(\dfrac{2}{\beta^{\prime 4}} - \dfrac{12}{\beta^{\prime 2}}\right) + a^{2}\left(\dfrac{3\gamma^{\prime}}{\beta^{\prime 2}} - 18 + \dfrac{7}{\beta^{\prime 2}}\right) + a\left(6\gamma^{\prime} - 12\beta^{\prime 2} - \dfrac{2\gamma^{\prime}}{\beta^{\prime 2}} + 8\right) - 3\beta^{\prime 4} + 3\beta^{\prime 2}\gamma^{\prime} + 3\beta^{\prime 2} - 3\gamma^{\prime} + \rho^{\prime}$

# $R(a) = B(a)^{2} - \dfrac{20aC(a)}{9}$ (not going to write this one out, it's long)

rth = betaprm**2+a*betaprm

B = (gammaprm-rtheta+(rtheta**3-rtheta**2)/betaprm**2)/betaprm

C = (-gammaprm+rhoprm+3*gammaprm*(rtheta/betaprm)**2-2*gammaprm*rtheta/betaprm**2
     +(rtheta/betaprm)**2+rtheta**5*sigmaprm/betaprm**4-3*(rtheta/betaprm)**4+2*rtheta**3/betaprm**4)/betaprm

B = collect(expand(B.subs(rtheta,rth)),a)
B

C = collect(expand(C.subs(rtheta,rth)),a)
C

b0 = betaprm**3-2*betaprm+gammaprm/betaprm
b1 = 3*betaprm**2-3
b2 = 3*betaprm-1/betaprm
b3 = 1

c0 = betaprm**5*sigmaprm-3*betaprm**3+3*betaprm*gammaprm+3*betaprm-3*gammaprm/betaprm+rhoprm/betaprm
c1 = 5*betaprm**4*sigmaprm+6*gammaprm-12*betaprm**2+8-2*gammaprm/betaprm**2
c2 = 10*betaprm**3*sigmaprm+(3*gammaprm+7)/betaprm-18*betaprm
c3 = 10*betaprm**2*sigmaprm-12+2/betaprm**2
c4 = 5*betaprm*sigmaprm-3/betaprm
c5 = sigmaprm

gammap = expand(solve(b0-b_0,gammaprm)[0])
sigmap = expand(solve(c0-c_0,sigmaprm)[0])

replacements = [(gammaprm,gammap),(sigmaprm,sigmap)]

sigmap = simplify(sigmap.subs(gammaprm,gammap))

sigmap

gammap

b0

b1

b2

b3

B_a = b3*a**3+b2*a**2+b1*a+b_0
B_a

# c0 = c0.subs(sigmaprm,sigmap).subs(gammaprm,gammap)
c0

c1 = c1.subs(sigmaprm,sigmap).subs(gammaprm,gammap)
expand(c1)

c2 = c2.subs(sigmaprm,sigmap).subs(gammaprm,gammap)
expand(c2)

c3 = c3.subs(sigmaprm,sigmap).subs(gammaprm,gammap)
expand(c3)

c4 = c4.subs(sigmaprm,sigmap).subs(gammaprm,gammap)
expand(c4)

c5 = c5.subs(sigmaprm,sigmap).subs(gammaprm,gammap)
c5

C_a = c_5*a**5+c4*a**4+c3*a**3+c2*a**2+c1*a+c_0
collect(expand(C_a),a)

series(B_a**2-20*a*C_a/9,a,n=7)

series(B**2-20*a*C/9,a,n=7)

Etrun = a*thetaprm+b*thetaprm**3+c*thetaprm**5
Etrun

collect(Etrun.subs([(b,B),(c,C)]),thetaprm)

theta_L = solve(Etrun.diff(thetaprm),thetaprm)[1]
theta_U = solve(Etrun.diff(thetaprm),thetaprm)[3]
theta_L,theta_U

E_L = simplify(Etrun.subs(thetaprm,theta_U))
E_U = simplify(Etrun.subs(thetaprm,theta_L))
E_L,E_U

"""
Explanation: Rescaling
End of explanation
"""


