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


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,diff(p(t),t)=diff(p_t,t),diff(rho(t),t)=diff(rho_t,t),diff(rhov(t),t)=0};

simpeq:=proc(x)
 subs(sublist,x);
 simplify(subs(sublist2,simplify(%)));
 subs(diff(sigma(t),t)=dsigma,%):
 simplify(subs(sublist2,%));
 subs(diff(eta(t),t)=deta,x):
 simplify(%);
 subs(K=(1-Kfac)*k^2/3,%);
 simplify(subs(sublist2,%));
 subs(rhoq(t)=rhoq_sub,x);
 subs(sublist,%);
 simplify(%);
 simplify(subs(sublist2,%));
 simplify(%);
 simplify(subs(sublist2,%));
end;


#note g2 includes 1/H factor, g2= g/H

#evolution when g is the actual total source distribution
evterm:=-diff(g2(t),t)*sigma(t);
simpeq(%);
diff(%,t)/k;
evterm:=simpeq(%);

#total redshif tterm for g actual source distribution
fac1(t):=simplify(diff(sigma(t),t)*g2(t));
simpeq(%);
diff(%,t)/k;
redterm:=simpeq(%);

#normal redshift term
fac1(t):=simplify(k*sigma(t)*g2(t));
simpeq(%);
diff(%,t)/k;
simpeq(%);
diff(%,t)/k;
redterm:=simpeq(%);


#radial displacement
vterm:=g2(t)*sigma(t)*2/(t0-t);

vterm:=simpeq(diff(vterm,t)/k);
vterm:=simplify(subs(t0=chi+t,%));

#other
otherterm:=winH(t)*sigma(t);
otherterm:=simpeq(diff(%,t)/k);


###########
#calculational stuff

##Go into fderivs's variables


adot:=adotoa*a;

fortsubs:={H(t)=adotoa,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)=xxddwing[j],diff(g(t),t)=xxdwing[j],g(t)=xxwing[j],diff(diff(winH(t),t),t)=xxddwingtau[j],diff(winH(t),t)=xxdwingtau[j],winH(t)=xxwingtau[j],diff(diff(g2(t),t),t)=xxddwing2[j],diff(g2(t),t)=xxdwing2[j],g2(t)=xxwing2[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};



calcproc:=proc(x)
#Approx here - pb(t)=0
simplify(subs(pb(t)=0,x));
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,%);
simplify(%);
subs(diff(rho(t),t)=-3*H(t)*(rho(t)+p(t)),%);
subs({c2(t)=0,rhopi(t)=0,photbar(t)=0},%);
 subs(K=(1-Kfac)*k^2/3,%);
simplify(%);
subs(fortsubs,%);
simplify(%);
collect(%,[k,adotoa,sigma,pig]);
collect(%,[Kf[1],k,adotoa,sigma,pig]);
collect(%,[sigma,xxdWing2[j],xxwing2[j],xxddWing2[j]]);

end;


precision:=double;
with(codegen,fortran);

calcproc(redterm);
fortran(%);


calcproc(evterm);
fortran(%);


calcproc(vterm);
fortran(%);


calcproc(redterm);
fortran(%);

calcproc(Hterm);
fortran(%);