#The vector perturbation equations, in A=0 gauge
#See astro-ph/0403583
#AL March 04

# Use conformal time. 
# Includes cdm, baryons, lambda

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);
assign(no_numassive);

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));

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

phi(t):=-kappa*S(t)^2/2/k^3*(k*rhopi(t) + 4 * H(t) * rhoq(t));

drag_t:=opac(t)*(4/3*v(t)-qg(t));
photbar_t:=rhog(t)/rhob(t);

dqr:=-1/2*k*pir(t);
dqg:=-1/2*k*pig(t)+drag(t);
dv:=-(1-3*c2(t))*H(t)*v(t)-photbar(t)*drag(t) - k/2*B0*rhog(t)/rhob(t);
dvc=-H(t)*vc(t);


#dsigma:=simplify(-H(t)*sigma(t)+k*phi(t)-1/2*kappa*S(t)^2/k*(rhopi(t)));
dsigma:=simplify(-2*H(t)*sigma(t)-kappa*S(t)^2/k*(rhopi(t)));


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)) - 8/15*k*J_3(t) + 2/5*k*qg(t) + 8/15*k*sigma(t);
dJ_3:=k*(3/7*pig(t)-15/28*J_4(t))-opac(t)*J_3(t);

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

G_eq:=proc(l)
 local Gl,Eq;
 Gl:=cat('G_',l)(t);
 Eq:=diff(Gl,t)+k*l/(2*l+1)*( (l+2)/(l+1)*cat('G_',l+1)(t) - cat('G_',l-1)(t));
 if (l = 2) then
   Eq:=Eq - 8/15*k*sigma(t);
 fi;
 simplify(subs({G_1(t)=qr(t),G_2(t)=pir(t)},Eq));
end;

rhopi_t:=rhog(t)*pig(t)+rhor(t)*pir(t) + rhonu(t)*pinu(t) + rhog(t)*B0;
rhoq_t:=rhoc(t)*vc(t) + rhog(t)*qg(t)+rhor(t)*qr(t)+(rhob(t)+pb(t))*v(t) + rhonu(t)*qnu(t) +k*diff(psi(t),t)*clv(t)/S(t)^2;

sigma(t):=2*kappa*S(t)^2*rhoq_t/k^2;
#rhoq(t):=k^2*sigma(t)/2/kappa/S(t)^2;


subtots:={rhopi(t)=rhopi_t,rho(t)=rho_t,p(t)=p_t};

drhoq:=-4*H(t)*rhoq(t)-1/2*k*rhopi(t);

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

dE2:=-opac(t)*(E2(t) - polter(t)) - 8/27*k*E3(t) + 1/3*k*B2(t);
dB2:=-opac(t)*B2(t)  - 8/27*k*B3(t) - 1/3*k*E2(t);

sublist:={diff(S(t),t)=dS,diff(qr(t),t)=dqr,diff(qg(t),t)=dqg,diff(v(t),t)=dv,diff(H(t),t)=dH,diff(exptau(t),t)=g(t),diff(pig(t),t)=dpig};

subtotderivs:={diff(S(t),t)=dS,diff(H(t),t)=dH};

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