Sorry for bothering you. But i have a problem with the plotting out variable evolution in cmbmain.
When I try to plot the perturbation of cdm (divided by the perturbation today), I get the plot LCDM_cmd.png.
Here the perturbation start high and goes down before going back up again. While I taugh that the perturbation was growing with time (decreasing redshift).
The neutrino density, LCDM_neutrino.png, follow the same path, which does not make sense since neutrinos are relativistic at first, then the perturbation grows as they become non-relativistic.
I have also tried dividing by the factor a^2 kappa, but this not change the shape of the perturbation. This I did since they all should come with this factor.
For the neutrino density have I also tried to use
call MassiveNuVars(EV,y,a,grho_nu,gpres_nu,dgrho_nu,dgq_nu)
but realized that this will give me the total density perturbation and not only for neutrinos.
Can you help me, or point out what is wrong. There only a few line,and should be "easy" to stop the mistake or the misunderstanding.
Under you find the code:
Sincerly,
Paul
Code: Select all
if (fixq/=0._dl) then
write(*,*) 'k', EV%q
tol1=tol/exp(AccuracyBoost-1)
call CreateTxtFile(CP%output_name,1)
do j=1,1000
tauend = taustart+(j-1)*(CP%tau0-taustart)/1000
call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
a = y(1)
redshift_z = 1._dl/a - 1._dl
do nu_i = 1, CP%Nu_mass_eigenstates
grhormass_t=grhormass(nu_i)/a**2
call Nu_background(a*nu_masses(nu_i),rhonu,pnu)
!Integrate over q
call Nu_Integrate_L012(EV, y, a, nu_i, clxnu,qnu) !Before FFF
qnu=qnu/rhonu
clxnu = clxnu/rhonu
grhonu_t=grhormass_t*rhonu
gpnu_t=grhormass_t*pnu
grho_nu = grhonu_t
gpres_nu = gpnu_t
dgrho_nu = grhonu_t*clxnu
dgq_nu = grhonu_t*qnu
end do
yprime = 0
call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
adotoa = 1/(y(1)*dtauda(y(1)))
ddelta = (yprime(3)*grhoc+yprime(4)*grhob)/(grhob+grhoc)
delta = (grhoc*y(3)+grhob*y(4))/(grhob+grhoc)
delta_tot = (grhoc*y(3)/a+grhob*y(4)/a+dgrho_nu)/(grhob/a+grhoc/a+grho_nu)
growth = ddelta/delta/adotoa
write (1,'(9E15.5)') redshift_z, dgrho_nu, grho_nu, grhoc*y(3)/a, grhoc/a, grhob*y(4)/a, grhob/a, delta, delta_tot
end do
close(1)
stop
end if