#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)*V(t);

dfac1(t):=diff(fac1(t),t)/k;

simplify(subs(sublist,%));

diff(%,t)/k;
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(H(t),t)=dH};



x;
simplify(%);
simplify(subs(sublist2,%));
ddfac1:=simplify(%);

vterm:=g(t)*(sigma(t)+v(t))*2*(1-1/H(t)/(t0-t));
ddfac1+diff((diff(V(t),t)+g(t))*v(t)+vterm,t)/k;
simplify(subs(sublist2,%));

x:=%;

hprime:=1/3*k*z(t);

term1:= g(t)*phi(t)*(3-1/H(t)/(t0-t))-4/chistar*phi(t)*F(t)+4*(1-1/Hstar/chistar)*diff(phi(t),t)*F(t);

source:=-(x+1/3*fac1(t)-hprime*V(t) - g(t)*clxb(t) );


source;
subs(diff(sigma(t),t)=dsigma,%):
simplify(subs(sublist2,%));
simplify(%);

x:=simplify(subs(sublist2,%));
subs(diff(eta(t),t)=deta,x):
x:=simplify(%);

subs(K=(1-Kfac)*k^2/3,%);
x:=simplify(%);

subs(rhoq(t)=rhoq_sub,x);
x:=simplify(%);
###########
#calculational stuff


simplify(subs(H(t)=adot/S(t),%));
##Approx here - pb(t)=0
simplify(subs(pb(t)=0,%));
subs(diff(rhopi(t),t)=0,%);

subs(rhox(t)=dgrho/S(t)^2/kappa,%);
subs(rhopi(t)=0,%);

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:={clxb(t)=clxb,diff(qg(t),t)=qgdot,diff(pig(t),t)=pigdot,diff(pir(t),t)=pirdot,diff(qr(t),t)=qrdot,diff(v(t),t)=0,rhoq(t)=dgq/kappa/a^2,S(t)=a,sigma(t)=sigma,diff(diff(V(t),t),t)=ddWinV[j],diff(V(t),t)=dwinV[j],V(t)=winV[j],diff(diff(g(t),t),t)=ddwing[j],diff(g(t),t)=dwing[j],g(t)=wing[j],v(t)=0,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};



precision:=double;
with(codegen,fortran);
subs(fortsubs,x);
subs(t0=chi+t,%);

simplify(%);
subs(adot=adotoa*a,%);
#simplify(%,[expmmu(j)*opac(j)=vis(j)]);
tmp:=expand(%);
tmp:=factor(%);

collect(%,[k,adotoa,sigma,pig]);
collect(%,[ddWinV[j],dWinV[j],dWing[j],wing[j],winV[j]]);
collect(%,[Kf[1]]);
collect(%,k);

fortran(%);