# The scalar perturbation equations in any frame
# Assume maple 6 or higher (maple 4 is buggy)

# Use conformal time. 
# Includes cdm, baryons, lambda, massive neutrinos and a scalar field
# We set the CDM (accel=0) frame from the start

restart;
no_quint:={psi(t)=0,V(t)=0,V1(t)=0,V2(t)=0,rhopsi(t)=0,ppsi(t)=0};
no_numassive:={pinu(t)=0,qnu(t)=0,rhonu(t)=0,pnu(t)=0};

assign(no_quint);

H_t:=diff(S(t),t)/S(t);

dS:=S(t)*H(t);
dH:=-1/6*S(t)^2*kappa*(rho(t)+3*p(t));


Friedmann:=H(t)^2=1/3*S(t)^2*kappa*rho(t)-K;

#Eliminate phi(t) [Weyl tensor variable] using constraint equation
cons1:=k^3*Kfac*phi(t)+1/2*kappa*S(t)^2*k*(rhox(t)+rhopi(t)*Kfac)+3/2*kappa*S(t)^2*H(t)*rhoq(t);
solve(cons1,phi(t));
phi_sub:=simplify(%); 

#Eliminate z [perturbation in expansion rate] in favour of the curvature perturbation eta
z_sub:=solve(k^2*eta(t)=kappa*S(t)^2*rhox(t) -2*k*H(t)*z(t),z(t));

#Solve for sigma(t) [shear] using constraint
cons2:=2/3*(k/S(t))^2*(z(t)-Kfac*sigma(t))+kappa*rhoq(t);
solve(cons2,sigma(t));
sigma_sub:=subs(z(t)=z_sub,%);

eta_sub:=solve(sigma(t)=sigma_sub,eta(t));
hdash_sub:=k/3*z(t)-H(t)*A(t);
A_sub:=solve(hdash(t)=hdash_sub,A(t));

solve(cons2,rhoq(t));
rhoq_sub:=%;

varsubs:={z(t)=z_sub,sigma(t)=sigma_sub,phi(t)=phi_sub,hdash(t)=hdash_sub};
#assign(varsubs);

drag_t:=opac(t)*(4/3*v(t)-qg(t));
photbar_t:=rhog(t)/rhob(t);
dgrho_t:=kappa*S(t)^2*rhox(t);


drho:=-3*H(t)*(rho(t)+p(t));

drhob:=-3*H(t)*(rhob(t)+pb(t));
drhoc:=-3*H(t)*rhoc(t);
drhog:=-4*H(t)*rhog(t);
drhor:=-4*H(t)*rhor(t);
drhonu:=-3*H(t)*(rhonu(t)+pnu(t));
drhopsi:=-3*H(t)*(rhopsi(t)+ppsi(t));

dclxr:=-4*hdash(t)-k*qr(t);
dclxg:=-4*hdash(t)-k*qg(t);
dclxb:=-(1+pb(t)/rhob(t))*k*(3*hdash(t)/k+v(t)) +(pb(t)/rhob(t)-c2(t))*3*H(t)*clxb(t);
dclxc:=-3*hdash(t) - k*vc(t);
dclxnu:=-3*hdash(t)*(1+pnu(t)/rhonu(t)) - k*qnu(t) +3*H(t)*(-clxpnu(t) + clxnu(t)*pnu(t)/rhonu(t));

dclvdot:= - 2*H(t)*diff(clv(t),t) - k*z(t)*diff(psi(t),t) - k^2*clv(t) - clv(t)*S(t)^2*V2(t);

dz:=diff(z_sub,t);

dpsidot:=-2*H(t)*diff(psi(t),t)-S(t)^2*V1(t);

dphi:=-H(t)*phi(t)+1/2/k^2*kappa*S(t)^2*(k*(rho(t)+p(t))*sigma(t)+k*rhoq(t)-diff(rhopi(t),t)-H(t)*rhopi(t));
solve(diff(phi(t),t)=dphi,diff(rhopi(t),t));
drhopi:=%;

rhopsi_t:=1/2*diff(psi(t),t)^2/S(t)^2 + V(t);
ppsi_t:=1/2*diff(psi(t),t)^2/S(t)^2 - V(t);


Kfac_sub:= (1-3*K/k^2);
K_sub:=solve(Kfac=Kfac_sub,K);
dqr:=-2/3*k*pir(t)*Kfac+1/3*k*clxr(t)-4*k/3*A(t);
dqg:=-2/3*k*pig(t)*Kfac+1/3*k*clxg(t)+drag(t)-4*k/3*A(t);
dv:=-k*A(t) -(1-3*c2(t))*H(t)*v(t)+1/(1+pb(t)/rhob(t))*(c2(t)*k*clxb(t)-photbar(t)*drag(t));
dqnu:=-H(t)*(1-3*pnu(t)/rhonu(t))*qnu(t) -k/3*(2*Kfac*pinu(t) - 3*clxpnu(t))-(1+pnu(t)/rhonu(t))*k*A(t);
dvc:=-H(t)*vc(t)-k*A(t);

dsigma:=simplify(-H(t)*sigma(t)+k*(phi(t)+A(t))-1/2*kappa*S(t)^2/k*(rhopi(t)));
simplify(subs(rhoq(t)=rhoq_sub,dsigma));
dsigma:=%;

rho_t:=rhob(t)+rhoc(t)+rhor(t)+rhog(t)+rhonu(t)+rhopsi(t)+rhov(t);
p_t:=1/3*(rhor(t)+rhog(t))+pb(t)+pnu(t)+ppsi(t)-rhov(t);
dpb:=c2(t)*drhob;


dpig:=-opac(t)*(pig(t)-polter(t)) - 3/5*k*Kf[2]*J_3(t) + 2/5*k*qg(t) + 8/15*k*sigma(t);

dpir:=- 3/5*k*Kf[2]*G_3(t) + 2/5*k*qr(t) + 8/15*k*sigma(t);
dG_3:=k*(3/7*pir(t)-4/7*Kf[3]*G_4(t));
dG_4:=k*(4/9*G_3(t)-5/9*Kf[4]*G_5(t));

rhopi_t:=rhog(t)*pig(t)+rhor(t)*pir(t) + rhonu(t)*pinu(t);
rhox_t:=rhog(t)*clxg(t)+rhor(t)*clxr(t)+rhob(t)*clxb(t)+rhoc(t)*clxc(t) +rhonu(t)*clxnu(t) + 1/S(t)^2*(diff(psi(t),t)*diff(clv(t),t)+clv(t)*S(t)^2*V1(t));
rhoq_t:=rhog(t)*qg(t)+rhor(t)*qr(t)+(rhob(t)+pb(t))*v(t) + rhonu(t)*qnu(t) +rhoc(t)*vc(t) +k*diff(psi(t),t)*clv(t)/S(t)^2;
#pressure perturbation
rhoxp_t:=rhonu(t)*clxpnu(t) + 1/3*(rhog(t)*clxg(t)+rhor(t)*clxr(t)) + 1/S(t)^2*(diff(psi(t),t)*diff(clv(t),t)-clv(t)*S(t)^2*V1(t))+c2(t)*rhob(t)*clxb(t);

subtots:={rhopi(t)=rhopi_t,rhoq(t)=rhoq_t,rhox(t)=rhox_t,rho(t)=rho_t,p(t)=p_t,rhoxp(t)=rhoxp_t};

subtotsquint := {rhopsi(t)=rhopsi_t, ppsi(t)=ppsi_t};

deta:=simplify(-1/k*(2*K*z(t) + kappa*S(t)^2*rhoq(t) +2*Kfac*k*H(t)*A(t)));

drhoq:=k*rhoxp(t)-4*H(t)*rhoq(t)-2/3*k*Kfac*rhopi(t)-(rho(t)+p(t))*k*A(t);
drhox:=-k*rhoq(t) - 3*H(t)*(rhox(t) + rhoxp(t)) - 3*hdash(t)*(rho(t)+p(t));


polter_t:=2/15*(3*pig(t)/4 + 9*E2(t)/2);

dE2:=-opac(t)*(E2(t) - polter(t)) - k/3*Kf[2]*E3(t);




sublist:={diff(vc(t),t)=dvc,diff(S(t),t)=dS,diff(rhob(t),t)=drhob,diff(rhoc(t),t)=drhoc,diff(rhor(t),t)=drhor,diff(rhog(t),t)=drhog,diff(clxr(t),t)=dclxr,diff(clxg(t),t)=dclxg,diff(clxb(t),t)=dclxb,diff(clxc(t),t)=dclxc,diff(z(t),t)=dz,diff(qr(t),t)=dqr,diff(qg(t),t)=dqg,diff(v(t),t)=dv,diff(sigma(t),t)=dsigma,diff(H(t),t)=dH,diff(exptau(t),t)=g(t),diff(pig(t),t)=dpig,diff(rhonu(t),t)=drhonu,diff(clxnu(t),t)=dclxnu,diff(eta(t),t)=deta,diff(diff(clv(t),t),t)=dclvdot,diff(diff(psi(t),t),t)=dpsidot,diff(V1(t),t)=diff(psi(t),t)*V2(t)};

subtotderivs:={diff(rhox(t),t)=drhox,diff(rhoq(t),t)=drhoq,diff(S(t),t)=dS,diff(H(t),t)=dH,diff(eta(t),t)=deta};

#End of main definitions and equations
#########################################



#first check consistency

diff(cons2,t);
x:=simplify(subs(sublist,%));

subs(subtots,%);
simplify(%);
x:=simplify(subs(sublist,%));

subs(varsubs,%);
subs(subtots,%);
simplify(%);
subs(photbar(t)=photbar_t,%);
subs(diff(pb(t),t)=dpb,%);
subs(diff(qnu(t),t)=dqnu,%);
subs(varsubs,%);
simplify(%);

subs(H(t)=sqrt(rhs(Friedmann)),%);
subs(subtots,%);
simplify(%);
#Should give 0


#Define some frame invariant quantities
GI_vars:={hclxc(t)=clxc(t) + 3/2*eta(t)/Kfac,
hclxb(t)=clxb(t) + 3/2*eta(t)/Kfac,
hclxr(t)=clxr(t) + 2*eta(t)/Kfac,
hclxg(t)=clxg(t) + 2*eta(t)/Kfac,
hvb(t)=v(t)+sigma(t),
hvc(t)=vc(t)+sigma(t),
hqr(t)=qr(t)+ 4/3*sigma(t),
hqg(t)=qg(t)+ 4/3*sigma(t),
Ph(t)=eta(t)/2/Kfac+H(t)*sigma(t)/k};

rhoxbar(t):=rhox(t)+3*H(t)/k*rhoq(t);
etahat(t):=eta(t)+2*Kfac*rhox(t)/3/(rho(t)+p(t));
etabar(t) := eta(t) - 2*H(t)*Kfac*rhoq(t)/(rho(t)+p(t))/k;
sigmabar(t):=sigma(t) + rhoq(t)/(rho(t)+p(t));
dpnad:=rhoxp(t) - diff(p(t),t)/diff(rho(t),t)*rhox(t);
Phi(t):=-phi(t)-1/2*S(t)^2*kappa*rhopi(t)/k^2;
Psi(t):=phi(t)-1/2*S(t)^2*kappa*rhopi(t)/k^2;
A_sub2:=solve(diff(sigma(t),t)=dsigma,A(t));
Newt_subs:={A(t)=-Psi(t),sigma(t)=0,eta(t)=2*Phi(t)};
varsubsN:=solve({z(t)=z_sub,sigma(t)=sigma_sub,phi(t)=phi_sub},{phi(t),eta(t),z(t)});


####################################END