#Get results for scalar source term for the line of sight integral read(scal_eqs.txt); #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)=rhog(t)/rhob(t),%); subs(diff(pb(t),t)=dpb,%); subs(diff(qnu(t),t)=dqnu,%); simplify(%); # Consistent if now have zero ###################################### ###Do line of sight stuff (diff by parts, etc). phi_sub:=simplify(subs(rhoq(t)=rhoq_sub,phi_sub)); z(t):=z_sub; phi(t):=phi_sub; fac1(t):=k*sigma(t)*exptau(t) +g(t)*15/8*polter(t); dfac1(t):=-diff(fac1(t),t)/k; simplify(subs(polter(t)=polter_t,%)): simplify(subs(diff(E2(t),t)=dE2,%)); simplify(subs(sublist,%)); -diff(%,t)/k; simplify(subs(polter(t)=polter_t,%)): x:=%; sublist2:={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(rhonu(t),t)=drhonu,diff(exptau(t),t)=g(t),diff(H(t),t)=dH}; x; #subs(diff(rhopi(t),t)=diff(rhopi_t,t),%); simplify(%); simplify(subs(sublist2,%)); #subs(rhopi(t)=rhopi_t,%); ddfac1:=simplify(%); %+diff(g(t)*v(t),t)/k; simplify(subs(sublist2,%)); x:=%; source:=x+1/3*fac1(t)-1/3*exptau(t)*k*z(t) +1/4*g(t)*clxg(t); source; subs(diff(sigma(t),t)=dsigma,%): simplify(subs(sublist2,%)); simplify(%); #simplify(subs(rhopi(t)=rhopi_t,%)); x:=simplify(subs(sublist2,%)); subs(diff(eta(t),t)=deta,x): x:=simplify(%); #collect(%,theta(t)); #x:=subs(Friedmann,%); #subs(K=(1-Kfac)*k^2/3,%); #x:=simplify(%); #subs(eta(t)=eta_sub,%): #simplify(%): #x:=%; subs(K=(1-Kfac)*k^2/3,%); #x:=simplify(subs(diff(polter(t),t)=diff(polter_t,t),%)); x:=simplify(subs(polter(t)=polter_t,%)); ## subs(Friedmann,expand(x)); subs(K=(1-Kfac)*k^2/3,%); x:=simplify(%); subs(rhoq(t)=rhoq_sub,x); x:=simplify(%); ISW:=coeff(x,exptau(t))*exptau(t); subs(Friedmann,simplify(ISW)); subs(K=(1-Kfac)*k^2/3,%); ########### #calculational stuff simplify(subs(H(t)=adot/S(t),x)); ##Approx here - pb(t)=0 simplify(subs(pb(t)=0,%)); subs(diff(rhopi(t),t)=diff_rhopi/kappa/S(t)^2,%); subs(rhox(t)=dgrho/S(t)^2/kappa,%); subs(rhopi(t)=dgpi/S(t)^2/kappa,%); #subs(z(t)=zs(t)/S(t),%); subs(rhopi(t)=rhopi_t/kappa/S(t)^2,%); subs(rhog(t)=grhog_t/kappa/S(t)^2,rhor(t)=grhor_t/kappa/S(t)^2,rhob(t)=grhob_t/kappa/S(t)^2,rhoc(t)=grhoc_t/kappa/S(t)^2,rhonu(t)=grhonu_t/kappa/S(t)^2,pnu(t)=gpnu_t/kappa/S(t)^2,rhopsi(t)=grhopsi_t/kappa/S(t)^2,ppsi(t)=gppsi_t/kappa/S(t)^2,%); subs(diff(S(t),t)=adot,%); x:=simplify(%); ##Go into fderivs's variables fortsubs:={diff(qg(t),t)=qgdot,diff(pig(t),t)=pigdot,diff(pir(t),t)=pirdot,diff(qr(t),t)=qrdot,diff(v(t),t)=vbdot,rhoq(t)=dgq/kappa/a^2,S(t)=a,sigma(t)=sigma,exptau(t)=expmmu[j],diff(diff(g(t),t),t)=ddvis[j],diff(g(t),t)=dvis[j],g(t)=vis[j],diff(opac(t),t)=dopac[j],opac(t)=opac[j],v(t)=vb,pig(t)=pig,pir(t)=pir,qg(t)=qg,diff(J_3(t),t)=yprime[9],z(t)=z,clxg(t)=clxg,diff(E2(t),t)=ypolprime[2],E2(t)=ypol[2],diff(E3(t),t)=ypolprime[3],E3(t)=ypol[3],polter(t)=polter,J_3(t)=y[9],Kfac=Kf[1],diff(pinu(t),t)=pinudot,pinu(t)=pinu,p(t)=gpres/kappa/a^2,rho(t)=grho/kappa/a^2,eta(t)=-2*etak/k,diff(eta(t),t)=-2*etakdot/k}; ##simplify(subs(subtots,%)); precision:=double; with(codegen,fortran); subs(fortsubs,x); fortresult:=%; ISW:=coeff(fortresult,expmmu[j])*expmmu[j]; #ISW bit simplify(ISW); subs(adot=adotoa*a,%); tmp:=expand(%); tmp:=factor(%); collect(%,[Kf[1],k,adotoa,sigma,pig]); collect(%,[expmmu[j]]); fortran(%); #the rest simplify(fortresult-ISW); subs(adot=adotoa*a,%); tmp:=expand(%); tmp:=factor(%); collect(%,[k,adotoa,sigma,pig]); collect(%,[ddvis(j),exppmu[j],dvis[j],vis[j],opac[j]]); collect(%,[Kf[1]]); fortran(%);