lightcone.c File Reference

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <sys/stat.h>
#include <sys/types.h>
#include "utils.h"
#include "timefac.h"
#include "partmanager.h"
#include "cosmology.h"
#include "physconst.h"

Include dependency graph for lightcone.c:

Go to the source code of this file.

Macros
#define	NENTRY   4096

Functions
static double	lightcone_get_horizon (double a)
static void	lightcone_cross (int p, double ddrift)
static void	lightcone_set_time (double a, const double BoxSize)
static double	kernel (double loga, void *params)
static void	lightcone_init_entry (Cosmology *CP, int i, const double UnitLength_in_cm)
void	lightcone_init (Cosmology *CP, double timeBegin, const double UnitLength_in_cm, const char *OutputDir)
static void	update_replicas (double a, double BoxSize)
void	lightcone_compute (double a, double BoxSize, Cosmology *CP, inttime_t ti_curr, inttime_t ti_next)

Variables
static double	tab_loga [NENTRY]
static double	dloga
static double	tab_Dc [NENTRY]
static int	Nreplica
static int	BoxBoost = 20
static double	Reps [8192][3]
static double	HorizonDistance2
static double	HorizonDistance
static double	HorizonDistancePrev
static double	HorizonDistance2Prev
static double	HorizonDistanceRef
static double	zmin = 0.1
static double	zmax = 80.0
static double	ReferenceRedshift = 2.0
static double	SampleFraction
static FILE *	fd_lightcone

Macro Definition Documentation

NENTRY

#define NENTRY   4096

Definition at line 19 of file lightcone.c.

Function Documentation

kernel()

static double kernel ( double  loga, void *  params  )

static

Definition at line 59 of file lightcone.c.

                                                 {
    double a = exp(loga);
      Cosmology * CP = (Cosmology *) params;
    return 1 / hubble_function(CP, a) * CP->Hubble / a;
}

References CP, Cosmology::Hubble, and hubble_function().

Referenced by density_kernel_desnumngb(), density_kernel_dW(), density_kernel_dwk(), density_kernel_init(), density_kernel_init_with_type(), density_kernel_volume(), density_kernel_wk(), density_ngbiter(), GetDensityKernelType(), lightcone_init_entry(), metal_return_ngbiter(), set_density_params(), and stellar_density_ngbiter().

Here is the call graph for this function:

Here is the caller graph for this function:

lightcone_compute()

void lightcone_compute	(	double	a,
		double	BoxSize,
		Cosmology *	CP,
		inttime_t	ti_curr,
		inttime_t	ti_next
	)

Definition at line 177 of file lightcone.c.

{
    int i;
    lightcone_set_time(a, BoxSize);
    const double ddrift = get_exact_drift_factor(CP, ti_curr, ti_next);
    #pragma omp parallel for
    for(i = 0; i < PartManager->NumPart; i++)
    {
        lightcone_cross(i, ddrift);
    }
}

References CP, get_exact_drift_factor(), lightcone_cross(), lightcone_set_time(), part_manager_type::NumPart, and PartManager.

Referenced by run().

Here is the call graph for this function:

Here is the caller graph for this function:

lightcone_cross()

static void lightcone_cross ( int  p, double  ddrift  )

static

Definition at line 222 of file lightcone.c.

                                                  {
    if(SampleFraction <= 0.0) return;
    int i;
    int k;
    /* DM only */
    if(P[p].Type != 1) return;
 
    for(i = 0; i < Nreplica; i++) {
        double r = get_random_number(P[p].ID + i);
        if(r > SampleFraction) continue;
 
        double pnew[3];
        double pold[3];
        double p3[4];
        double dnew = 0, dold = 0;
        for(k = 0; k < 3; k ++) {
            pold[k] = P[p].Pos[k] + Reps[i][k] - PartManager->CurrentParticleOffset[k];
            pnew[k] = P[p].Pos[k] + P[i].Vel[k] * ddrift - PartManager->CurrentParticleOffset[k];
            dnew += pnew[k] * pnew[k];
            dold += pold[k] * pold[k];
        }
        if(
            (dold <= HorizonDistance2Prev && dnew >= HorizonDistance2)
         ) {
            double u1, u2;
            if(dold != dnew) {
                double cnew, cold;
                dnew = sqrt(dnew);
                dold = sqrt(dold);
                cnew = dnew - HorizonDistance;
                cold = dold - HorizonDistancePrev;
                u1 = -cold / (cnew - cold);
                u2 = cnew / (cnew - cold);
            } else {
                /* really should write all particles along the line:
                 * this partilce is moving along the horizon! */
                u1 = u2 = 0.5;
            }
 
            /* write particle position */
            for(k = 0; k < 3; k ++) {
                p3[k] = pold[k] * u2 + pnew[k] * u1;
            }
            p3[3] = SampleFraction;
            fwrite(p3, sizeof(double), 4, fd_lightcone);
        }
    }
}

References part_manager_type::CurrentParticleOffset, fd_lightcone, get_random_number(), HorizonDistance, HorizonDistance2, HorizonDistancePrev, Nreplica, P, PartManager, Reps, and SampleFraction.

Referenced by lightcone_compute().

Here is the call graph for this function:

Here is the caller graph for this function:

lightcone_get_horizon()

static double lightcone_get_horizon ( double  a )

static

Definition at line 119 of file lightcone.c.

                                              {
    double loga = log(a);
    int bin = (log(a) -tab_loga[0]) / dloga;
    if (bin < 0) {
        return tab_Dc[0];
    }
    if (bin >= NENTRY - 1) {
        return tab_Dc[NENTRY - 1];
    }
    double u1 = loga - tab_loga[bin];
    double u2 = tab_loga[bin + 1] - loga;
    u1 /= (tab_loga[bin + 1] - tab_loga[bin]);
    u2 /= (tab_loga[bin + 1] - tab_loga[bin]);
    return tab_Dc[bin] * u2 + tab_Dc[bin + 1] * u1;
}

References dloga, NENTRY, tab_Dc, and tab_loga.

Referenced by lightcone_init(), and lightcone_set_time().

Here is the caller graph for this function:

lightcone_init()

void lightcone_init	(	Cosmology *	CP,
		double	timeBegin,
		const double	UnitLength_in_cm,
		const char *	OutputDir
	)

Definition at line 91 of file lightcone.c.

{
    int i;
    dloga = (0.0 - log(timeBegin)) / (NENTRY - 1);
    for(i = 0; i < NENTRY; i ++) {
        lightcone_init_entry(CP, i, UnitLength_in_cm);
    };
    char buf[1024];
    int chunk = 100;
    int ThisTask;
    MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
 
    sprintf(buf, "%s/lightcone/", OutputDir);
    mkdir(buf, 02755);
    sprintf(buf, "%s/lightcone/%03d/", OutputDir, (int)(ThisTask / chunk));
    mkdir(buf, 02755);
    sprintf(buf, "%s/lightcone/%03d/lightcone-%05d.raw", OutputDir, (int)(ThisTask / chunk), ThisTask);
 
    fd_lightcone = fopen(buf, "a+");
    if(fd_lightcone == NULL) {
        endrun(1, "failed to open %s\n", buf);
    }
    HorizonDistanceRef = lightcone_get_horizon(1 / (1 + ReferenceRedshift));
    printf("lightcone reference redshift = %g distance = %g\n",
            ReferenceRedshift, HorizonDistanceRef);
}

References CP, dloga, endrun(), fd_lightcone, HorizonDistanceRef, lightcone_get_horizon(), lightcone_init_entry(), NENTRY, ReferenceRedshift, ThisTask, and UnitLength_in_cm.

Referenced by begrun().

Here is the call graph for this function:

Here is the caller graph for this function:

lightcone_init_entry()

static void lightcone_init_entry ( Cosmology *  CP, int  i, const double  UnitLength_in_cm  )

static

Definition at line 65 of file lightcone.c.

                                                                                       {
    tab_loga[i] = - dloga * (NENTRY - i - 1);
 
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);
 
    double result, error;
 
    gsl_function F;
    F.function = &kernel;
    F.params = CP;
    gsl_integration_qags (&F, tab_loga[i], 0, 0, 1e-7, 1000,
            w, &result, &error);
 
    /* result is in DH, hubble distance */
    /* convert to cm / h */
    result *= LIGHTCGS / HUBBLE;
    /* convert to Kpc/h or internal units */
    result /= UnitLength_in_cm;
 
    gsl_integration_workspace_free (w);
    tab_Dc[i] = result;
//    double a = exp(tab_loga[i]);
//    double z = 1 / a - 1;
//    printf("a = %g z = %g Dc = %g\n", a, z, result);
}

References CP, dloga, HUBBLE, kernel(), LIGHTCGS, NENTRY, tab_Dc, tab_loga, and UnitLength_in_cm.

Referenced by lightcone_init().

Here is the call graph for this function:

Here is the caller graph for this function:

lightcone_set_time()

void lightcone_set_time ( double  a, const double  BoxSize  )

static

Definition at line 189 of file lightcone.c.

                                                        {
    double z = 1 / a - 1;
    if(z > zmin && z < zmax) {
        HorizonDistancePrev = HorizonDistance;
        HorizonDistance2Prev = HorizonDistance2;
        HorizonDistance = lightcone_get_horizon(a);
        HorizonDistance2 = HorizonDistance * HorizonDistance;
        update_replicas(a, BoxSize);
        fflush(fd_lightcone);
        if (z < ReferenceRedshift) {
            SampleFraction = 1.0;
        } else {
            /* write a smaller fraction of the points at high redshift
             */
            /* This is the angular resolution rule */
            SampleFraction = HorizonDistanceRef / HorizonDistance;
            SampleFraction *= SampleFraction;
            SampleFraction *= SampleFraction;
            /* This is the luminosity resolution rule */
#if 0
            SampleFraction = HorizonDistanceRef / HorizonDistance;
            SampleFraction *= (1 + ReferenceRedshift) / (1 + z);
            SampleFraction *= SampleFraction;
 
#endif
        }
        message(0,"RefRedeshit=%g, SampleFraction=%g HorizonDistance=%g\n", ReferenceRedshift, SampleFraction, HorizonDistance);
    } else {
        SampleFraction = 0;
    }
}

References fd_lightcone, HorizonDistance, HorizonDistance2, HorizonDistance2Prev, HorizonDistancePrev, HorizonDistanceRef, lightcone_get_horizon(), message(), ReferenceRedshift, SampleFraction, update_replicas(), zmax, and zmin.

Referenced by lightcone_compute().

Here is the call graph for this function:

Here is the caller graph for this function:

update_replicas()

static void update_replicas ( double  a, double  BoxSize  )

static

Definition at line 136 of file lightcone.c.

                                                      {
    int Nmax = BoxBoost * BoxBoost * BoxBoost;
    int i;
    int rx, ry, rz;
    rx = ry = rz = 0;
    Nreplica = 0;
 
    for(i = 0; i < Nmax; i ++) {
        double dx = BoxSize * rx;
        double dy = BoxSize * ry;
        double dz = BoxSize * rz;
        double d1, d2;
        d1 = dx * dx + dy * dy + dz * dz;
        dx += BoxSize;
        dy += BoxSize;
        dz += BoxSize;
        d2 = dx * dx + dy * dy + dz * dz;
        if(d1 <= HorizonDistance2 && d2 >= HorizonDistance2) {
            Reps[Nreplica][0] = rx * BoxSize;
            Reps[Nreplica][1] = ry * BoxSize;
            Reps[Nreplica][2] = rz * BoxSize;
            Nreplica ++;
            if(Nreplica > 1000) {
                endrun(951234, "too many replica");
            }
        }
        rz ++;
        if(rz == BoxBoost) {
            rz = 0;
            ry ++;
        }
        if(ry == BoxBoost) {
            ry = 0;
            rx ++;
        }
    }
}

References BoxBoost, endrun(), HorizonDistance2, Nreplica, and Reps.

Referenced by lightcone_set_time().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

BoxBoost

int BoxBoost = 20

static

Definition at line 35 of file lightcone.c.

Referenced by update_replicas().

dloga

double dloga

static

Definition at line 21 of file lightcone.c.

Referenced by cooling_and_starformation(), cooling_direct(), density(), density_ngbiter(), get_helium_neutral_fraction_sfreff(), get_long_range_timestep_dloga(), get_neutral_fraction_sfreff(), get_PM_timestep_ti(), get_sfr_eeqos(), get_timestep_dloga(), get_timestep_ti(), hydro_ngbiter(), lightcone_get_horizon(), lightcone_init(), lightcone_init_entry(), print_timebin_statistics(), SPH_EntVarPred(), and winds_evolve().

fd_lightcone

FILE* fd_lightcone

static

Definition at line 46 of file lightcone.c.

Referenced by lightcone_cross(), lightcone_init(), and lightcone_set_time().

HorizonDistance

double HorizonDistance

static

Definition at line 38 of file lightcone.c.

Referenced by lightcone_cross(), and lightcone_set_time().

HorizonDistance2

double HorizonDistance2

static

Definition at line 37 of file lightcone.c.

Referenced by lightcone_cross(), lightcone_set_time(), and update_replicas().

HorizonDistance2Prev

double HorizonDistance2Prev

static

Definition at line 40 of file lightcone.c.

Referenced by lightcone_set_time().

HorizonDistancePrev

double HorizonDistancePrev

static

Definition at line 39 of file lightcone.c.

Referenced by lightcone_cross(), and lightcone_set_time().

HorizonDistanceRef

double HorizonDistanceRef

static

Definition at line 41 of file lightcone.c.

Referenced by lightcone_init(), and lightcone_set_time().

Nreplica

int Nreplica

static

Definition at line 34 of file lightcone.c.

Referenced by lightcone_cross(), and update_replicas().

ReferenceRedshift

double ReferenceRedshift = 2.0

static

Definition at line 44 of file lightcone.c.

Referenced by lightcone_init(), and lightcone_set_time().

Reps

double Reps[8192][3]

static

Definition at line 36 of file lightcone.c.

Referenced by lightcone_cross(), and update_replicas().

SampleFraction

double SampleFraction

static

Definition at line 45 of file lightcone.c.

Referenced by lightcone_cross(), and lightcone_set_time().

tab_Dc

double tab_Dc[NENTRY]

static

Definition at line 22 of file lightcone.c.

Referenced by lightcone_get_horizon(), and lightcone_init_entry().

tab_loga

double tab_loga[NENTRY]

static

Definition at line 20 of file lightcone.c.

Referenced by lightcone_get_horizon(), and lightcone_init_entry().

zmax

double zmax = 80.0

static

Definition at line 43 of file lightcone.c.

Referenced by lightcone_set_time().

zmin

double zmin = 0.1

static

Definition at line 42 of file lightcone.c.

Referenced by lightcone_set_time().