#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

##Only approx for now

polter_line(t):=F2(t)/10 + 9/15*E_line2(t);

dF2:=k/5*( 2*F1(t) - 3*F3(t))  -2/15*k*sigma(t)*V(t) + opac(t)*(polter_line(t)-F2(t));

dE_line2:=-opac(t)*(E_line2(t) - polter_line(t)) - k/3*Kf[2]*E_line3(t);

######################################
###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(F2(t),t)=dF2,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,diff(E_line2(t),t)=dE_line2};

#Approximate distortion term with cs_2=0
pb(t):=0;
dist_term:= exptau(t)*k/dS*S(t)*(sigma(t)+v(t))*wing(t);
diff(dist_term,t)/k;
simplify(subs(sublist,%));
dist_term:= diff(%,t)/k;
simplify(%);
simplify(subs(sublist2,%));
simplify(subs(sublist,%));
subs(diff(rho(t),t)=-3*H(t)*(rho(t)+p(t)),%);
subs({c2(t)=0,rhopi(t)=0,p(t)=0,photbar(t)=0},%);

dist_term:=simplify(%);


fac1(t):=-exptau(t)*(k*sigma(t)*V(t) + q_fac*g2(t)*15/8*pig(t)) + reion*g(t)*3/4*F2(t);

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

simplify(subs(sublist,%));

diff(%,t)/k;
x:=%;


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

ddfac1+ diff( reion*g(t)*V(t)*v(t)  -exptau(t)*(diff(V(t),t)*v(t)+wing(t)*v(t)-g2(t)*d_fac*(v(t)-qg(t)*3/4)),t)/k;
simplify(subs(sublist2,%));

x:=%;

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

source:=(x+1/3*fac1(t)+reion*g(t)*F0(t)+exptau(t)*hprime*V(t) + exptau(t)*wing(t)*Delta_source(t) + exptau(t)*g2(t)*Delta_source2(t));

source:=source+ distort*dist_term - dist_term;

source;

simplify(subs(polter(t)=polter_t,%)):
simplify(subs(diff(E2(t),t)=dE2,%));

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(%);

subs(diff(qg(t),t)=dqg,x);
subs(drag(t)=drag_t,%);

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


simplify(subs(H(t)=adot/S(t),%));
##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(diff(F(t),t)=wing(t),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)=vbdot,rhoq(t)=dgq/kappa/a^2,S(t)=a,sigma(t)=sigma,diff(diff(V(t),t),t)=xxddWinV[j],diff(V(t),t)=xxdwinV[j],V(t)=xxwinV[j],diff(diff(wing(t),t),t)=xxddwing[j],diff(wing(t),t)=xxdwing[j],wing(t)=xxwing[j],g2(t)=xxwing2[j],diff(g2(t),t)=xxdwing2[j],diff(diff(g2(t),t),t)=xxddwing2[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,Delta_source(t)=Delta_source,Delta_source2(t)=Delta_source2,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],diff(F1(t),t)=yprime[lineoff+1],diff(F3(t),t)=yprime[lineoff+2],F0(t)=y[lineoff],F1(t)=y[lineoff+1],F2(t)=y[lineoff+2],F3(t)=y[lineoff+3],F(t)=xxwinF[j], E_line2(t)=y[lineoffpol+2],E_line3(t)=y[lineoffpol+3]};

precision:=double;
with(codegen,fortran);
subs(fortsubs,x);

x:=%;

makefort:=proc(ain)
local tmp;
simplify(ain);
subs(adot=adotoa*a,%);
#simplify(%,[expmmu(j)*opac(j)=vis(j)]);
tmp:=expand(%);
tmp:=factor(%);
collect(%,[k,adotoa,sigma,pig]);
collect(%,[xxddWinV[j],xxdWinV[j],xxdWing[j],xxwing[j],xxwing2[j],xxdwing2[j],xxwinV[j],d_fac,q_fac]);
collect(%,[Kf[1]]);
collect(%,[k,expmmu[j]]);
fortran(%);
end;

base:=coeff(x,Delta_source)*Delta_source+coeff(x,Delta_source2)*Delta_source2;

dist:=simplify(coeff(x,distort))*distort;


makefort(coeff(x,distort));

makefort(coeff(x,reion));
makefort(coeff(x,d_fac));
makefort(coeff(x,q_fac));
makefort(subs({d_fac=0,q_fac=0,reion=0},x - base-dist));