#Get second order action for scalar perturbations
#in terms of user-defined variable, following astro-ph/9902265
#I paste into Maple what I want to run
#AML 13 Nov 00

#Note Delta here = box+3K = -k^2+3K

restart;
phidot:=diff(phi(t),t);
phiddot:=diff(diff(phi(t),t),t);
phidddot:=diff(phiddot,t);
phidot_sub:=sqrt(2)*z(t)*H(t)/sqrt(kappa);
zdot:=diff(z(t),t);
zddot:=diff(diff(z(t),t),t);
nudot:=diff(nu(t),t);

#Newtonian potential
PsiN_sub:=psi(t)+H(t)*kappa/2/a(t)^2*Pipsi(t)/Delta;
#Curvature perturbation

chi_sub:=(H(t)*dphi(t)/phidot + psi(t));
chi_sub:=(dphi(t) + psi(t)*phidot/H(t));

Gamma_sub:=1+(K-diff(H(t),t))/H(t)^2;

#Define variable nu(t) we want to use

NuDef:=nu(t)=chi_sub;
#nu(t)=(dphi(t) + phidot/H(t)*psi(t))/sqrt(1-K/Delta*kappa/2*phidot^2/H(t)^2);
#nu(t)=a(t)*(dphi(t) + phidot/H(t)*psi(t));
#nu(t)=((H(t)*dphi(t)/phidot + psi(t))-2*K/kappa/phidot^2*PsiN_sub));
#nu(t)=Delta*dphi(t)-kappa*phidot*Pipsi(t)/2/a(t)^2; 
#nu(t)=alpha*H(t)/phidot*(dphi(t)-kappa*phidot*Pipsi(t)/2/a(t)^2/Delta)+(psi(t)+H(t)*kappa/2/a(t)^2*Pipsi(t)/Delta);
#nu(t)=a(t)*(dphi(t)-kappa*phidot*Pipsi(t)/2/a(t)^2/Delta);
#nu(t)=simplify(-(K*PsiN_sub+H(t)^2*chi_sub)/(H(t)^2+K)+Delta/3/Gamma_sub/H(t)^2*PsiN_sub);


#V3 comaptible coeff
coef:= proc(X,Y)
	collect(X,Y);
	coeff(%,Y);
	RETURN(%);
end;


#background eqs:
background:={diff(H(t),t)=-kappa*diff(phi(t),t)^2/2+H(t)^2 + K,diff(a(t),t)=a(t)*H(t)};
Ebackground:={diff(H(t),t)=-kappa*diff(phi(t),t)^2/2+H(t)^2 - K,diff(a(t),t)=a(t)*H(t)};


##End useful definitions

NuDef:=simplify(subs(background,NuDef));

#Do full thing in maple

S2:=Pipsi(t)*diff(psi(t),t)+Piphi(t)*diff(dphi(t),t) - kappa/4/a(t)^2/Delta*(-K*Pipsi(t)^2+2*Delta*Piphi(t)^2/kappa) -kappa/2*phidot*Pipsi(t)*dphi(t)- a(t)^2/kappa*(Delta*psi(t)^2 - kappa*dphi(t)^2/2*(Delta-H(t)^2-diff(H(t),t) + phidddot/phidot));

#A integral:
-H(t)*Pipsi(t) + phidot*Piphi(t) + 2*a(t)^2/kappa*(-Delta*psi(t) + kappa/2*(H(t)*phidot - phiddot)*dphi(t));
solve(%,psi(t));
psi(t):=simplify(%);


solve(NuDef,dphi(t));
dphi(t):=simplify(%);

simplify(S2);
subs(background,%);
subs(background,%);
simplify(%);
S2:=%;

collect(S2,diff(Pipsi(t),t));
coeff(%,diff(Pipsi(t),t));
tmp:=simplify(%);
collect(%,Pipsi(t));
coeff(%,Pipsi(t));
T1:=%;
subs(Pipsi(t)=0,tmp);
-diff(%,t)*Pipsi(t)-1/2*diff(T1,t)*Pipsi(t)^2;

subs(diff(Pipsi(t),t)=0,S2)+%;
subs(background,%);
subs(background,%);
simplify(%);
S2:=%;


collect(S2,diff(Piphi(t),t));
simplify(coeff(%,diff(Piphi(t),t)));
subs(Piphi(t)=1,%);
-1/2*diff(%,t)*Piphi(t)^2+subs(diff(Piphi(t),t)=0,S2);
simplify(%);
subs(background,%);
subs(background,%);
simplify(%);

S2:=%;


collect(S2,diff(nu(t),t));
coeff(%,diff(nu(t),t));
Pinu(t)=%;
solve(%,Piphi(t));
Piphi(t):=simplify(%);


simplify(S2);
subs(background,%);
subs(background,%);
simplify(%);
S2:=%;
#Got rid of Pipsi OK!

collect(S2,Pinu(t));
coeff(%,Pinu(t)^2);
Pinu2:=simplify(%);
collect(S2,Pinu(t));
coeff(%,Pinu(t));
Pinu1:=simplify(%);

tmp:=S2;


#Do gaussian integral
subs(Pinu(t)=0,tmp) - Pinu1^2/4/Pinu2;

tmp:=simplify(%);
collect(tmp,diff(nu(t),t));
coeff(%,diff(nu(t),t));
subs(nu(t)=1,%);
nudotnusurface:=simplify(%);
nudotterm:=-1/2*diff(%,t);

collect(tmp,diff(nu(t),t));
coeff(%,diff(nu(t),t)^2);
%*diff(nu(t),t)^2 + nudotterm*nu(t)^2 + subs(diff(nu(t),t)=0,tmp);
subs(background,%);
subs(background,%);
simplify(%);
res:=%;

#####res is the result for the action#####

z_sub:= proc(X)
    subs(phidot=phidot_sub,X):
    simplify(%):
    subs(background,%):
    subs(background,%):
    simplify(%):
    subs(background,%):
    subs(background,%):
    simplify(%):
    subs(phidot=phidot_sub,%):
    simplify(%):
    subs(background,%):
    subs(phidot=phidot_sub,%):
    simplify(%):
    RETURN(%):
end;

zres:=z_sub(res);

#put in terms of z


#See what individual terms are
collect(zres,nudot);
coeff(%,nudot^2);
simplify(%);

subs(nu(t)=1,zres);
zres2:=simplify(%);

collect(%,zddot);
coeff(%,zddot);
simplify(%);

subs(zddot=0,zres2);
collect(%,zdot);
simplify(coeff(%,zdot));

subs(zdot=0,zres2);
simplify(%);


#collect(zres,nudot);
#coeff(%,nudot^2);
#nudotnusurface*nu(t)^2/2+%*nudot*nu;


############End of general stuff#########################

#Put in schodinger form, for Mukhanov nu
subs(nu(t)=sqrt(1-K*z(t)^2/Delta)*Q(t),zres);

Qres:=simplify(%);

collect(%,diff(Q(t),t));
coeff(%,diff(Q(t),t));
subs(Q(t)=1,%);
-1/2*diff(%,t);
term:=%;

collect(Qres,diff(Q(t),t));
coeff(%,diff(Q(t),t)^2);
%*diff(Q(t),t)^2 + subs(diff(Q(t),t)=0,Qres) + term*Q(t)^2;
Qres:=simplify(%);

collect(Qres,diff(Q(t),t));
coeff(%,diff(Q(t),t)^2);
simplify(%);

subs(diff(Q(t),t)=0,Qres);
subs(Q(t)=1,%);
Q2terms:=%;

collect(Q2terms,diff(diff(z(t),t),t));
simplify(coeff(%,diff(diff(z(t),t),t)));
Z2:=%;

subs(diff(diff(z(t),t),t)=0,Q2terms);
collect(%,diff(z(t),t));
coeff(%,diff(z(t),t));
simplify(%);
Z1:=%;

subs(diff(diff(z(t),t),t)=0,Q2terms);
collect(%,diff(z(t),t));
coeff(%,diff(z(t),t)^2);
simplify(%);
Z12:=%;

subs(diff(z(t),t)=0,Q2terms);
Z0:=simplify(%);

#Got final answer in terms of z

#Get in terms of phi for computation
#U(t):=simplify(subs(z(t)=sqrt(kappa/2)*phidot/H(t),Q2terms)*2-diff(H(t),t)-H(t)^2);
U(t)=simplify(subs(z(t)=sqrt(kappa/2)*phidot/H(t),Q2terms)*2);
subs(background,%);
subs(background,%);
simplify(%);
U(t):=%;

###########





U(t);
#Go Euclidean
subs(K=-K,%);  
subs(Delta=-Delta,%);

UEuc(t):=%;
subs(phiddot=a(t)^2*Vdash -2*H(t)*phidot,UEuc(t));
subs(Ebackground,%);
simplify(%);


subs(H(t)=H,%);
subs(phiddot=dotdotphi,%);
subs(phidot=dotphi,%);
subs(a(t)=a,%);
subs(phi(t)=phi,%);
collect(%,{dotphi,Delta,K,kappa,H});
fortran(%);

#Get value at t=0


exp_subs:={phi(t)=phi0 + phi2*t^2/2 + phi4*t^4 + phi6*t^6,H(t)=-K*t + H3*t^3+H5*t^5,a(t)=1 - K*t^2/2 + a4*t^4};

V := V0+V2*t^2+V4*t^4;

phiddot-a(t)^2*diff(V,t)/phidot + 2*H(t)*phidot;
subs(exp_subs,%);
collect(%,t);
series(%,t,5);
ser1:=%;
coeff(%,t,0);
solve(%,phi2);
phi2:=%;
phi2:=sqrt(2*V2);
ser1:=collect(simplify(ser1),t);

K-H(t)^2 - kappa/3*(-phidot^2/2 + a(t)^2*V);
subs(exp_subs,%);
collect(%,t);
series(%,t,5);
ser2:=%;



coeff(collect(ser1,t),t^2);
solve(%,phi4);
phi4:=%;

diff(a(t),t)/a(t)-H(t);
subs(exp_subs,%);
series(%,t);
coeff(%,t^3);
solve(%,a4);
a4:=%;

collect(ser2,t);
coeff(%,t^4);
simplify(%);
solve(%,H3);
H3:=%;


diff(V,t)/phidot;
subs(exp_subs,%);
simplify(%);
series(%,t);
coeff(%,t,0);
Vdash0=%;
solve(%,V2);
V2:=%;

diff(diff(V,t)/phidot,t)/phidot;
simplify(%);
subs(exp_subs,%);
simplify(%);
series(%,t);
coeff(%,t,0);
Vddash0=%;
solve(%,V4);
V4:=simplify(%);

#Get z

sqrt(kappa/2)*phidot/H(t);
subs(exp_subs,%);
z:=collect(simplify(%),t);
series(z,t);

subs(K=kappa/3*V0,%);
simplify(%);
coeff(%,t,2);


diff(z,t)/z/H(t);
subs(exp_subs,%);
subs(K=kappa/3*V0,%);
series(%,t);
tmp2:=%;

coeff(%,t,0);
simplify(%);

#Get classical solution near origin
p:=subs(exp_subs,-z^2*a(t)^2/(1-z^2*K/Delta));
U:=p*(2*K/(1-z^2*K/Delta)*tmp2+Delta+K);
G:=G0+G1*t+G2*t^2;

diff(p*diff(G,t),t) + U*G;
series(%,t);
Eq:=%;
coeff(%,t,0);
subs(K=kappa/3*V0,%);
simplify(%);
solve(%,G2);
G2:=simplify(%);


coeff(Eq,t,1);
subs(K=kappa/3*V0,%);
simplify(%);


subs(exp_subs,UEuc(t));

simplify(%);
series(%,t,5);
coeff(%,t,0);
simplify(%);
tmp:=%;





simplify(tmp);
simplify(%,{K=kappa*m*phi0^2/3});
subs(phi0^2=3*K/kappa/m,%);
simplify(%);
fortran(%);

######################
##relation between nu and psi_N (nu=bigQ)
#Restart to do this...

nu(t):=a(t)*(dphi(t) + phidot/H(t)*psi(t));
subs(dphi(t)=2/kappa*(diff(psi(t),t)+H(t)*psi(t))/phidot,%);
nu_sub=%;
solve(%,diff(psi(t),t));
psidot_sub:=%; 

diff(nu(t)/a(t)/phidot*H(t),t)*phidot^2/H(t);
subs(dphi(t)=2/kappa*(diff(psi(t),t)+H(t)*psi(t))/phidot,%);
subs(background,%);
subs(diff(diff(psi(t),t),t)=-2*(H(t)-phiddot/phidot)*diff(psi(t),t)-2*(diff(H(t),t)-H(t)*phiddot/phidot)*psi(t) - (-Delta - K)*psi(t),%);
subs(background,%);
simplify(%);
subs(diff(psi(t),t)=psidot_sub,%);
simplify(%);
X=%;
solve(%,psi(t));

subs(X=diff(n(t)*H(t)/a(t)/phidot,t)*phidot^2/H(t),%);
subs(nu_sub=n(t),%);
subs(background,%);
subs(n(t)=Q(t)*sqrt(1-K*kappa/2*phidot^2/H(t)^2/Delta),%);
psi:=simplify(%);

phidot_sub:=sqrt(2)*z(t)*H(t)/sqrt(kappa);

subs(phidot=phidot_sub,psi);
subs(background,%);
subs(phidot=phidot_sub,%);
subs(background,%);
simplify(%);
tmp:=%;
collect(%,diff(Q(t),t));
coeff(%,diff(Q(t),t));
simplify(%);

subs(diff(Q(t),t)=0,tmp);
simplify(%)
#check
sqrt(kappa/2)*H(t)/a(t)/Delta/sqrt(1-K*z(t)^2/Delta)*(z(t)*diff(Q(t),t) - diff(z(t),t)*Q(t)/(1-K*z(t)^2/Delta) - z(t)*Q(t)*(H(t)+K/H(t)))-tmp;
simplify(%);
#OK

##For N...

nu(t):=sqrt(2/kappa)*(H(t)*dphi(t)/phidot + psi(t));
subs(dphi(t)=2/kappa*(diff(psi(t),t)+H(t)*psi(t))/phidot,%);
nu_sub=%;
solve(%,diff(psi(t),t));
psidot_sub:=%; 

diff(nu(t),t);
subs(dphi(t)=2/kappa*(diff(psi(t),t)+H(t)*psi(t))/phidot,%);
subs(background,%);
subs(diff(diff(psi(t),t),t)=-2*(H(t)-phiddot/phidot)*diff(psi(t),t)-2*(diff(H(t),t)-H(t)*phiddot/phidot)*psi(t) - (-Delta - K)*psi(t),%);
subs(background,%);
simplify(%);
subs(diff(psi(t),t)=psidot_sub,%);
simplify(%);
dnu(t)=%;
solve(%,psi(t));
simplify(%);
subs(phidot=phidot_sub,%);
simplify(%);
psi_sub:=%;

#Get Chi
1/(1+K/H(t)^2)*(psidot_sub/H(t) + psi_sub) + psi_sub*(1/(1+K/H(t)^2) - 1);
simplify(%);
subs(psi(t)=psi_sub,%);
subs(phidot=phidot_sub,%);
simplify(%);
tmp:=%;

psi_sub+2/3*tmp;
simplify(%);
subs(H(t)=H,%);
subs(dnu(t)=yprime[5],%);
subs(nu_sub=y[5],%);
subs(z(t)=z,%);
fortran(%);



###################End

subs(phidot=2*z(t)*H(t),phidddot/phidot);

subs(phidot=sqrt(2*Y(t)*H(t)^2),%);


subs(diff(phi(t),t)=H(t)*z(t)/a(t),%);
subs(background,%);
subs(background,%);
simplify(%);
subs(diff(phi(t),t)=H(t)*z(t)/a(t),%);
subs(background,%);
simplify(%);

%check K=0;
subs(phidddot=phidddot_sub,%);
subs(K=0,%):
simplify(%);

###############################################
#tests

#See if get correct term for X


bterm:=H(t)*diff(phi(t),t) - diff(diff(phi(t),t),t);

-diff(beta(t),t)/2*a(t)^2*Gamma(t) + a(t)^2/2*beta(t)*diff(Gamma(t),t) - Gamma(t)^2/2 - a(t)^2/4/Delta*( bterm^2*beta(t)^2*a(t)^2 - 2*bterm*beta(t)*Gamma(t)*diff(phi(t),t) + Gamma(t)^2*diff(phi(t),t)^2/a(t)^2) + a(t)^4/2*beta(t)^2*(Delta-H(t)^2-diff(H(t),t)+diff(diff(diff(phi(t),t),t),t)/diff(phi(t),t));

simplify(%);
subs(diff(a(t),t)=H(t)*a(t),%);
subs(diff(H(t),t)=-diff(phi(t),t)^2/2+H(t)^2 + K,%);
simplify(%);


#Yes


#coeff of v after int

-H(t) - diff(Gamma(t),t)/Gamma(t) - a(t)^2/4/Delta*(-2*(H(t)*diff(phi(t),t) - diff(diff(phi(t),t),t))^2/Gamma(t)*beta(t) + 2*(H(t)*diff(phi(t),t)-diff(diff(phi(t),t),t))*diff(phi(t),t)/a(t)^2);

##Old stuff (junk)



#Work

#coeff of v before int:
bterm:=H(t)*diff(phi(t),t) - diff(diff(phi(t),t),t);

-H(t) - diff(Gamma(t),t)/Gamma(t) - a(t)^2/4/Delta*(-2*bterm^2/Gamma(t)*beta(t) + 2*bterm*diff(phi(t),t)/a(t)^2) -beta(t)/Gamma(t)*a(t)^2*(Delta-H(t)^2-diff(H(t),t)+phidddot/phidot);

simplify(%);
subs(background,%);
simplify(%);

A(t):=%;

#Do integral. nu^2 term is

diff(A(t),t)-A(t)*diff(X(t),t)/X(t) - A(t)^2;
subs(background,%);
simplify(%);

%/4/X(t) - bterm^2/Gamma(t)^2/4/Delta + 1/Gamma(t)^2/2*(Delta-H(t)^2 - diff(H(t),t)+phidddot/phidot);
subs(background,%);
simplify(%);
tmp:=%;



subs(phidddot=phidddot_sub,tmp);
simplify(%);

$check K=0
subs(K=0,%):
simplify(%);     






subs(diff(Pinu(t),t)=0,tmp);
collect(%,Pinu(t));
coeff(%,Pinu(t));
subs(diff(nu(t),t)=0,%);
subs(nu(t)=1,%):
simplify(%);
A2:=%;
 
tmp;
subs(Pinu(t)=0,%);
subs(diff(nu(t),t)=0,%);
collect(%,nu(t));
coeff(%,nu(t)^2);
nu2:=%;

-(diff(A2,t)-A2*diff(X(t),t)/X(t) - A2^2)/4/X(t) + nu2:
subs(background,%);
subs(background,%);
simplify(%);
res:=%;


#Change vars to q = nu/2/sqrt(X)
-res* 4*X(t) - diff(diff(X(t),t)/2/X(t),t) + 1/4*diff(X(t),t)^2/X(t)^2;


simplify(%);
subs(background,%);
subs(background,%);
simplify(%);

U(t):=%;





subs(phidot=phidot_sub,U(t));
simplify(%):
subs(background,%);
subs(background,%);
simplify(%);
subs(background,%);
subs(background,%);
simplify(%);
subs(phidot=phidot_sub,%);
simplify(%);
subs(background,%);
subs(phidot=phidot_sub,%);
simplify(%);
Uz(t):=%;

collect(Uz(t),diff(diff(z(t),t),t));
coeff(%,diff(diff(z(t),t),t));
Z2:=simplify(%);


subs(diff(z(t),t)=0,Uz(t));
simplify(%);
Z0:=%;

collect(Uz(t),diff(z(t),t));
subs(diff(diff(z(t),t),t)=0,%);
coeff(%,diff(z(t),t));
Z1:=simplify(%);

collect(Uz(t),diff(z(t),t));
subs(diff(diff(z(t),t),t)=0,%);
coeff(%,diff(z(t),t)^2);
Z12:=simplify(%);

Uz(t) - Z2*diff(diff(z(t),t),t)- Z1*diff(z(t),t) - Z0 - Z12*diff(z(t),t)^2;
simplify(%);
#so we can write it down moderately neatly