#Derive the initial series expansions for CAMB
#AL, December 02. 
#Assume neutrinos highly relativistic, baryons cold,
#Don't use tight coupling, assume n_e \propto ne/S^3 (opac propto 1/S^2)

restart;
read "C:/Web/antony/CAMB/maple/vec_eqs.txt";


#vc(t):=vc0/S(t);
#this mode is singular so can do
vc(t):=0;

opac(t):= ne0/S(t)^2;
ne0:=omm^2*H0^4/omega^3/delta_tight;

pb(t):=0;
drag(t):= drag_t;
photbar(t):=photbar_t;
polter(t):=polter_t;


#B0:=0;

#include terms up to t^Ord
Ord:=6;
DoQuint = "F";

Order:=Ord+3;

makeseries:=proc(v,x,ord)
 local Res,i;
 Res:=0;
 for i from 0 to ord do
   Res:=Res + cat('s',v)[i]*x^i;
 end do;
 v(x)=series(Res,x,ord); 
end;
    

expandfully:=proc(x)
 eval(x);
 %;
 %;
end;


Rb:=1-Rc;

#v(t):=3/4*qg(t);

omb:=Rb*omm;
omc:=Rc*omm;
omr:=Rv*omm^2*H0^2/omega^2;
omg:=omm^2*H0^2/omega^2-omr;

pb(t):=0;c2(t):=0;

#pig(t):=0;


rhog(t):=3*omg*H0^2/kappa/S(t)^4;
rhor(t):=3*omr*H0^2/kappa/S(t)^4;
rhob(t):=3*omb*H0^2/kappa/S(t)^3;
rhoc(t):=3*omc*H0^2/kappa/S(t)^3;
rhov(t):=3*omv*H0^2/kappa;
H(t):=H_t;

assign(subtots);

K:=0;
Kfac:=beta[2];



#Start series solving
dsolve({expandfully(Friedmann),S(0)=0,D(S)(0)=H0^2*omm/omega},S(t),type=series);
assign(%);

# S(t):=omm*H0^2/omega^2*(omega*t+(omega*t)^2/4-K/6*omega*t^3-K/48*omega^2*t^4);



Rb:=1-Rc;

#v(t):=3/4*qg(t);

omb:=Rb*omm;
omc:=Rc*omm;
omr:=Rv*omm^2*H0^2/omega^2;
omg:=omm^2*H0^2/omega^2-omr;

pb(t):=0;c2(t):=0;
#pig(t):=0;

#Use tight coupling result
#R:=4/3*rhog(t)/rhob(t);
#drag(t):=1/(1+R)*(-4/3*H(t)*v(t)  - k/2*B0*R);

#Some frame invariant variables

solvevars:= {B2,J_3,E2,pig,v,qr,qg,pir,G_3,G_4,G_5};

eqs:={simplify(diff(J_3(t),t)-dJ_3),simplify(diff(B2(t),t)-dB2),simplify(diff(E2(t),t)-dE2),simplify(diff(pig(t),t)-dpig),simplify(diff(v(t),t)-dv),simplify(dqr-diff(qr(t),t)),simplify(dqg-diff(qg(t),t)),G_eq(2),G_eq(3),G_eq(4),G_eq(5)};

knownvars:={H0,omm,omv,k,beta[2],Kf[2],Kf[3],Kf[4],omega,Rv,Rc,B0,delta_tight,sqg[0],sG_3[0]};

series_subs:=map(makeseries,solvevars,t,Ord+3);
assign(%);

#spir[0]:=0;   
#sG_3[1]:=0;
#sG_4[2]:=0;

E3(t):=O(t^2);

#sG_3[0]:=0;
sG_4[0]:=0;
G_6(t):=O(t^2);
sG_5[0]:=0;
sG_5[1]:=0;

#sqg[0]:=0;
B0:=0;

series_eqs:=map(expandfully,eqs);

map(series,series_eqs,t):
map(simplify,%):
series_eqs:=%;

dosolve:=proc(ineqs)
 local tmp,vars;
 map(series,ineqs,t);
 map(simplify,%);
 map(convert,%,polynom);
 map(tcoeff,%,t):
 tmp:=map(simplify,%);
 vars:=indets(tmp) minus knownvars;
 print(vars);
 solve(tmp,vars);
 assign(%);
end;



for i from 1 to Ord+1 do
 dosolve(series_eqs);
end do;


sqg[0]:=0;

sig_ser:=simplify(series(expandfully(sigma(t)),t,7));



sG_3[0]:=solve(sig0=coeff(sig_ser,t,0),sG_3[0]);


delta_tight:=0;

simplify(sig_ser);

simplify(series(expandfully(qg(t)),t,4));
simplify(subs(delta_tight=0,subs(B0=0,%)));
simplify(coeff(%,t,3)/coeff(%,t,2));