CAMB - Plotting variable evolution

[Paul Anthony Fronteri]

Hey Antony.

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

[Antony Lewis]

The plots are not shown here, but I'd do the CDM case first, forgetting about neutrinos, and make sure that makes sense first.
What exactly are you plotting from the output of the code you quoted?

[Paul Anthony Fronteri]

Hey again.

It was me that misunderstood. I have been plotting delta * rho and not delta. It is now cleared and I have changed the thing to:

write (1,'(10E15.5)') redshift_z, clxnu, grho_nu/a**2/kappa2, y(3), grhoc/a**3/kappa2, y(4), grhob/a**3/kappa2, delta, delta_tot, growth

Sincerely,

Paul