#Get the tight coupling equations for the scalar modes pb(t):=0; drag(t):= drag_t; opac(t):=opc(t)/eps; photbar(t):=photbar_t; rhob(t):=4*rhog(t)/3/R(t); polter(t):=polter_t; qg(t):=Delta(t)+4/3*v(t); #without pol: #polter(t):=2/15*3*pig(t)/4; eqs:={simplify(dqg-diff(qg(t),t)),simplify(diff(pig(t),t)-dpig),simplify(diff(v(t),t)-dv),simplify(E_eq(2)),simplify(E_eq(3)),simplify(diff(J_3(t),t)-dJ_3)}; makeseries:=proc(v,x,ord) local Res,i; Res:=0; for i from 0 to ord do Res:=Res + cat('s',v)[i](t)*x^i; end do; v(t)=series(Res,x,ord); end; solvevars := {Delta,pig,J_3,v,E2,E3,E4,J_3,J_4}; series_subs:=map(makeseries,solvevars,eps,4); subs(series_subs,eqs); eqs:=map(series,%,eps,3); map(coeff,eqs,eps,-1); spig[0](t):=0; sDelta[0](t):=0; sJ_3[0](t):=0; sJ_3[1](t):=0; sE2[0](t):=0; sE3[0](t):=0; sE3[1](t):=0; sE4[0](t):=0; sE4[1](t):=0; sE4[2](t):=0; sJ_3[0](t):=0; sJ_3[1](t):=0; sJ_4[0](t):=0; sJ_4[1](t):=0; sJ_4[2](t):=0; eqs0:=map(coeff,eqs,eps,0); solve(eqs0,{spig[1](t),sE2[1](t),diff(sv[0](t),t),sDelta[1](t)}); sol_subs:=%; dsv0:=subs(sol_subs,diff(sv[0](t),t)); dsDelta1:=subs(sol_subs,diff(sDelta[1](t),t)); Delta1_sub:=subs(sol_subs,sDelta[1](t)); submap:=proc(x,y) simplify(subs(y,x)); end; eqs1:=map(simplify,map(coeff,eqs,eps,1)); map(submap,eqs1,diff(sDelta[1](t),t)=dsDelta1); map(simplify,%); map(submap,%,diff(R(t),t)=-H(t)*R(t)); #This is an approx - opac = S(t)*ne(t)*sigma_t #map(submap,%,diff(opc(t),t)=-2*H(t)*opc(t)); map(submap,%,diff(c2(t),t)=-H(t)*c2(t)); map(submap,%,subs(diff(sv[0](t),t) = dsv0)); eqs1:=%; sol_sub2:=solve(eqs1,{sJ_3[2](t),sDelta[2](t),sE2[2](t),sE3[2](t),spig[2](t),diff(sv[1](t),t)}); dsv1:=simplify(subs(sol_sub2,diff(sv[1](t),t))); solve(sDelta[1](t)=Delta1_sub,clxb(t)); clxb_sub:=%; subs(diff(clxb(t),t)=XXX,dsv1); subs(clxb(t)=clxb_sub,%); subs(XXX=diff(clxb(t),t),%); dsv1:=%; #various factors simplify(coeff(dsv1,spig[1](t))); factor(simplify(coeff(dsv1,diff(opc(t),t)))); factor(simplify(coeff(dsv1,sDelta[1](t)))); factor(simplify(coeff(dsv1,sv[0](t)))); factor(simplify(coeff(dsv1,diff(clxg(t),t)))); factor(simplify(coeff(simplify(dsv1-%*diff(clxg(t),t)),clxg(t)))); #higher order sE2[1](t):=spig[1](t)/4;