Commits

Michael Kuhlen  committed 82b3f1e

Added support for OmegaRadiationNow > 0, which might be important for early (z>100) starting redshifts. This requires a numerical integration of the Friedman equation, which is too slow for on the fly computation. Hence I've introduced a lookup table for t(a), a(t), and da/dt(a) in CosmologyComputeExpansionFactor(), CosmologyComputeTimeFromRedshift(). This lookup table can be turned on using 'gmake use-cosmotable-yes'.

  • Participants
  • Parent commits a4919d8
  • Branches cosmotable

Comments (0)

Files changed (10)

File src/enzo/CosmologyCalculateTable.C

+#ifdef USE_COSMOTABLE
+
+/***********************************************************************
+/
+/ CALCULATES A LOOKUP TABLE FOR 'Time*H0' and '(da/dt)/H0' AS A FUNCTION
+/ OF SCALE FACTOR.
+/
+/  written by: Michael Kuhlen
+/  date:       February 2010
+/
+/
+/  PURPOSE:
+/
+/  NOTE: some routines adopted and adapted from PKDGRAV's cosmo.c.
+/
+************************************************************************/
+ 
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "CosmologyParameters.h"
+
+// needed for dRombergO
+#define EPSCOSMO 1e-5
+#define MAXLEV 13
+#define FLOAT_MAXVAL 1e20
+#define MAX_ITER 100
+
+//double arccosh(double x);
+//double arcsinh(double x);
+
+double csmExp2Time(double ScaleFactor);
+double csmExp2Hub(double ScaleFactor);
+double csmTime2Exp(double Time);
+double csmTime2Hub(double Time);
+double csmCosmoTint(double Y);
+double dRombergO(double (*func)(double),double a,double b,double eps);
+
+
+int CosmologyCalculateTable(void)
+{  
+  int i;
+
+  FLOAT ScaleFactor, MinScaleFactor, MaxScaleFactor;
+
+//   printf("ScaleFactor = %e\n",ScaleFactor);
+//   printf("HubbleConstantNow = %e\n",HubbleConstantNow);
+//   printf("OmegaMatterNow = %e\n",OmegaMatterNow);
+//   printf("OmegaRadiationNow = %e\n",OmegaRadiationNow);
+//   printf("OmegaLambdaNow = %e\n",OmegaLambdaNow);
+
+        
+  // Note: don't know where to free this memory...
+  CosmologyTable_a = new FLOAT[CosmologyTableSize];
+  CosmologyTable_dadt = new FLOAT[CosmologyTableSize];
+  CosmologyTable_time = new FLOAT[CosmologyTableSize];
+  
+  MinScaleFactor = 1.0/1001.0;
+  MaxScaleFactor = 1.0;
+  
+  CosmologyTableDelta = (MaxScaleFactor - MinScaleFactor) / (FLOAT)(CosmologyTableSize - 1);
+  
+  for(i=0, ScaleFactor=MinScaleFactor; i<CosmologyTableSize; i++, ScaleFactor+=CosmologyTableDelta) {
+    CosmologyTable_a[i] = ScaleFactor;
+
+    CosmologyTable_dadt[i] = csmExp2Hub(ScaleFactor) * ScaleFactor / HubbleConstantNow;      
+    
+    CosmologyTable_time[i] = csmExp2Time(ScaleFactor) * HubbleConstantNow;    
+  }
+  
+#ifdef USE_MPI
+  printf("MyProcessorNumber = %d\n",MyProcessorNumber);
+  if (CosmologyTableWriteToFile && MyProcessorNumber == ROOT_PROCESSOR) {
+#else
+  if (CosmologyTableWriteToFile) {
+#endif
+
+    FILE *fptr;
+    if ((fptr = fopen(CosmologyTableFilename, "w")) == NULL) {
+      fprintf(stderr, "Error opening CosmologyTable output file %s\n", CosmologyTableFilename);
+      ENZO_FAIL("");
+    }
+    
+    fprintf(fptr,"#          a           z    (da/dt)/H0       H0*Time\n");
+    fprintf(fptr,"#\n");
+    
+    for(i=0;i<CosmologyTableSize;i++) {
+      FLOAT Redshift = 1.0/CosmologyTable_a[i] - 1.0;
+      fprintf(fptr,"%12.6f %11.6f  %12.6e  %12.6e\n",CosmologyTable_a[i],Redshift,CosmologyTable_dadt[i],CosmologyTable_time[i]);
+    }
+    
+    fclose(fptr);
+  }
+  
+  CosmologyTableCalculated = TRUE;
+
+  return SUCCESS;
+}
+
+
+double csmTime2Exp(double Time) {
+  double al=0,ah=1,a0,a1=1,at,a;
+  double th,f,f1,h,ho;
+  int j;
+  
+  assert(Time > 0);
+  th = csmExp2Time(ah);
+  /*
+  ** Search for upper bracket if needed.
+  */
+  while (Time > th) {
+    a0 = a1;
+    a1 = ah;
+    ah = a1+a0;
+    th = csmExp2Time(ah);
+  }
+  a = 0.5*(al+ah);
+  ho = ah-al;
+  h = ho;
+  f = Time - dRombergO(csmCosmoTint,0.0,POW(a,1.5),EPSCOSMO);
+  f1 = 1/(a*csmExp2Hub(a));
+  for (j=0;j<MAX_ITER;++j) {
+    if (a+f/f1 < al || a+f/f1 > ah || fabs(2*f) > fabs(ho*f1)) {
+      /*
+      ** Bisection Step.
+      */
+      ho = h;
+      h = 0.5*(ah-al);
+      a = al+h;
+      /*
+	printf("bisect al:%.14g ah:%.14g a:%.14g\n",al,ah,a);
+      */
+      if (a == al) return a;
+    }
+    else {
+      /*
+      ** Newton Step.
+      */
+      ho = h;
+      h = f/f1;
+      at = a;
+      a += h;
+      /*
+	printf("newton al:%.14g ah:%.14g a:%.14g\n",al,ah,a);
+      */
+      if (a == at) return a;
+    }
+    if (fabs(h) < EPSCOSMO) {
+      /*
+	printf("converged al:%.14g ah:%.14g a:%.14g t:%.14g == %.14g\n",
+	al,ah,a,dRombergO(csmCosmoTint,0.0,POW(a,1.5),EPSCOSMO*1e-1),
+	Time);
+      */
+      return a;
+    }
+    f = Time - dRombergO(csmCosmoTint,0.0,POW(a,1.5),EPSCOSMO*1e-1);
+    f1 = 1/(a*csmExp2Hub(a));
+    if (f < 0) ah = a;
+    else al = a;
+  }
+  assert(0);
+  
+  return 0.0; /* We never reach here, but this keeps the compiler happy */
+}
+
+double csmExp2Time(double ScaleFactor) {
+  double Time;
+
+  Time = dRombergO(csmCosmoTint, 0.0, POW(ScaleFactor, 1.5), EPSCOSMO);
+  return Time;
+}
+
+double csmCosmoTint(double Y) {
+  double ScaleFactor = POW(Y, 2.0/3.0);
+  
+  assert(ScaleFactor > 0.0);
+  return 2.0/(3.0*Y*csmExp2Hub(ScaleFactor));
+}
+
+double csmExp2Hub(double ScaleFactor) {
+  double OmegaCurvatureNow = 1.0 - OmegaMatterNow - OmegaLambdaNow - OmegaRadiationNow;
+  double a2,a3,a4;
+  double val;
+
+  assert(ScaleFactor > 0.0);
+  
+  a2 = ScaleFactor*ScaleFactor;
+  a3 = a2 * ScaleFactor;
+  a4 = a3 * ScaleFactor;
+
+  //  printf("ScaleFactor = %e\n",ScaleFactor);
+
+  val = HubbleConstantNow * sqrt( OmegaMatterNow/a3 +
+				  OmegaCurvatureNow/a2 +
+				  OmegaRadiationNow/a4 + 
+				  OmegaLambdaNow );
+  //  printf("val = %e\n",val);
+
+  return val;
+
+}
+
+
+double csmTime2Hub(double Time) {
+  double ScaleFactor = csmTime2Exp(Time) * (1 + InitialRedshift);
+  
+  assert(ScaleFactor > 0.0);
+  return csmExp2Hub(ScaleFactor);
+}
+
+
+/*
+** Romberg integrator for an open interval.
+*/
+
+double dRombergO(double (*func)(double),double a,double b,double eps) {
+  double tllnew;
+  double tll;
+  double tlk[MAXLEV+1];
+  Eint32 n = 1;
+  Eint32 nsamples = 1;
+  
+  tlk[0] = tllnew = (b-a)*(*func)(0.5*(b+a));
+  if (a == b) return tllnew;
+  tll = FLOAT_MAXVAL;
+  
+  eps*=0.5;
+  
+  while ((fabs((tllnew-tll)/(tllnew+tll)) > eps) && (n < MAXLEV)) {
+    /*
+     * midpoint rule.
+     */
+    double deltax;
+    double tlktmp;
+    Eint32 i;
+    
+    nsamples *= 3;
+    deltax = (b-a)/nsamples;
+    tlktmp = tlk[0];
+    tlk[0] = tlk[0]/3.0;
+    
+    for (i=0;i<nsamples/3;i++) {
+      tlk[0] += deltax*(*func)(a + (3*i + 0.5)*deltax);
+      tlk[0] += deltax*(*func)(a + (3*i + 2.5)*deltax);
+    }
+    
+    /*
+     * Romberg extrapolation.
+     */
+    
+    for (i=0;i<n;i++) {
+      double tlknew = (POW(9.0, i+1.)*tlk[i] - tlktmp)
+	/(POW(9.0, i+1.) - 1.0);
+      
+      tlktmp = tlk[i+1];
+      tlk[i+1] = tlknew;
+    }
+    tll = tllnew;
+    tllnew = tlk[n];
+    n++;
+  }
+  
+  assert((fabs((tllnew-tll)/(tllnew+tll)) < eps));
+  
+  return tllnew;
+}
+
+
+#endif

File src/enzo/CosmologyComputeExpansionFactor.C

 /
 /  written by: Greg Bryan
 /  date:       April, 1995
-/  modified1:
+/  modified1:  Michael Kuhlen
+/  date:       February, 2010
+/              modified to use a look-up table, #ifdef USE_COSMOTABLE.
 /
 /  PURPOSE:
 /
 #ifdef CONFIG_PFLOAT_16
 #define ETA_TOLERANCE 1.0e-20
 #endif
- 
+
 // function prototypes
  
 double arccosh(double x);
 double arcsinh(double x);
- 
+
 int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt)
 {
+
+  *a = FLOAT_UNDEFINED;
+ 
+  /* Find Omega due to curvature. */
+ 
+  float OmegaCurvatureNow = 1 - OmegaMatterNow - OmegaLambdaNow;
+
+  /* Convert the time from code units to Time * H0 (c.f. CosmologyGetUnits). */
+ 
+  float TimeUnits = 2.52e17/sqrt(OmegaMatterNow)/HubbleConstantNow/
+                    POW(1 + InitialRedshift,FLOAT(1.5));
+ 
+  FLOAT TimeHubble0 = time * TimeUnits * (HubbleConstantNow*3.24e-18);
+
+
+#ifdef USE_COSMOTABLE
+
+  if (CosmologyTableCalculated == FALSE)
+    CosmologyCalculateTable();
+
+  // printf("CosmologyTableSize = %d\n",(Eint32)CosmologyTableSize);
+  // printf("CosmologyTable_time = %e ... %e\n",CosmologyTable_time[0],CosmologyTable_time[CosmologyTableSize-1]);
+  // printf("TimeHubble0 = %e\n",TimeHubble0);
+
+  // Check whether TimeHubble0 lies right on the table boundary
+
+  if ( TimeHubble0 == CosmologyTable_time[0] ) {
+    (*a) = CosmologyTable_a[0] * (1 + InitialRedshift);
+    (*dadt) = CosmologyTable_dadt[0] * sqrt(2.0/(3.0*OmegaMatterNow*(1.0+InitialRedshift)));
+  }
+  if ( TimeHubble0 == CosmologyTable_time[CosmologyTableSize-1] ) {
+    (*a) = CosmologyTable_a[CosmologyTableSize-1] * (1 + InitialRedshift);
+    (*dadt) = CosmologyTable_dadt[CosmologyTableSize-1] * sqrt(2.0/(3.0*OmegaMatterNow*(1.0+InitialRedshift)));
+  }
+
+
+  // Find table index for TimeHubble0. The table is NOT equally spaced
+  // in this direction (time -> a), so we need to use bisection.
+
+  int j, jl, ju, jm;
+  jl = -1;
+  ju = CosmologyTableSize+1;
+  while (ju-jl > 1) {
+    jm = (ju+jl) >> 1;
+    if ( TimeHubble0 >= CosmologyTable_time[jm] )
+      jl = jm;
+    else
+      ju = jm;
+  }
+  j = jl;
+
+  // error check bounds, if j<0 or j>(CosmologyTableSize-1).
+  if( (j<0) || (j>(CosmologyTableSize-1)) ) {
+    fprintf(stderr,"Attempting to use CosmologyTable outside of bounds!\n");
+    ENZO_FAIL("");
+  }
+  
+  FLOAT dt = (TimeHubble0 - CosmologyTable_time[j]) / (CosmologyTable_time[j+1] - CosmologyTable_time[j]);
+
+  // printf("j = %d   (%e, %e)   dt = %e\n",j,CosmologyTable_time[j],CosmologyTable_time[j+1],dt);
+  // printf("%e %e\n",CosmologyTable_a[j],CosmologyTable_a[j+1]);
+
+  // Look up 'a'
+  (*a) = CosmologyTable_a[j] * (1.0 - dt) +
+    CosmologyTable_a[j+1] * dt;
+  (*a) *= (1 + InitialRedshift);  
+
+  // Look up 'da/dt'
+  (*dadt) = CosmologyTable_dadt[j] * (1.0 - dt) +
+    CosmologyTable_dadt[j+1] * dt;
+  (*dadt) *= sqrt(2.0/(3.0*OmegaMatterNow*(1.0+InitialRedshift)));
+  
+
+  return SUCCESS;
+#endif
+
  
   /* Error check. */
  
     ENZO_FAIL("");
   }
  
-  *a = FLOAT_UNDEFINED;
  
-  /* Find Omega due to curvature. */
- 
-  float OmegaCurvatureNow = 1 - OmegaMatterNow - OmegaLambdaNow;
- 
-  /* Convert the time from code units to Time * H0 (c.f. CosmologyGetUnits). */
- 
-  float TimeUnits = 2.52e17/sqrt(OmegaMatterNow)/HubbleConstantNow/
-                    POW(1 + InitialRedshift,FLOAT(1.5));
- 
-  FLOAT TimeHubble0 = time * TimeUnits * (HubbleConstantNow*3.24e-18);
- 
+
   /* 1) For a flat universe with OmegaMatterNow = 1, it's easy. */
  
   if (fabs(OmegaMatterNow-1) < OMEGA_TOLERANCE &&
       OmegaLambdaNow < OMEGA_TOLERANCE)
     *a      = POW(time/InitialTimeInCodeUnits, FLOAT(2.0/3.0));
- 
+
 #define INVERSE_HYPERBOLIC_EXISTS
  
 #ifdef INVERSE_HYPERBOLIC_EXISTS
 		OmegaLambdaNow*TempVal*TempVal*TempVal));
  
   /* Someday, we'll implement the general case... */
- 
+
+
   if ((*a) == FLOAT_UNDEFINED) {
     fprintf(stderr, "Cosmology selected is not implemented.\n");
     ENZO_FAIL("");

File src/enzo/CosmologyComputeTimeFromRedshift.C

 /
 /  written by: Greg Bryan
 /  date:       April, 1995
-/  modified1:
+/  modified1:  Michael Kuhlen
+/  date:       February, 2010
+/              modified to use a look-up table, #ifdef USE_COSMOTABLE.
 /
 /  PURPOSE:
 /
 #include "ErrorExceptions.h"
 #include "macros_and_parameters.h"
 #include "CosmologyParameters.h"
- 
+
 // function prototypes
  
 double arccosh(double x);
 double arcsinh(double x);
+
+int CosmologyComputeTimeFromRedshift(FLOAT Redshift, FLOAT *TimeCodeUnits)
+{ 
+
+  /* Error check [INCOMPLETE]. */
+
+  FLOAT TimeHubble0 = FLOAT_UNDEFINED;
+  FLOAT TimeUnits = 2.52e17/sqrt(OmegaMatterNow)/HubbleConstantNow/
+    POW(1 + InitialRedshift, FLOAT(1.5));  /*  see CosmologyGetUnits */
+
+
+#ifdef USE_COSMOTABLE
+
+  if (CosmologyTableCalculated == FALSE)
+    CosmologyCalculateTable();
+  
+  FLOAT ScaleFactor = 1.0 / (1.0 + Redshift);
+
+  // Find table index for ScaleFactor. The table is equally spaced in
+  // this direction (a -> time), so we can just use direct
+  // positioning.
+
+  int j = (int)( (ScaleFactor - CosmologyTable_a[0]) / CosmologyTableDelta );
+  FLOAT dt = (ScaleFactor - CosmologyTable_a[j]) / CosmologyTableDelta;
+
+  // printf("j = %d   (%e, %e)   dt = %e\n",j,CosmologyTable_a[j],CosmologyTable_a[j+1],dt);
+  // printf("%e %e\n",CosmologyTable_time[j],CosmologyTable_time[j+1]);
+
+  // Look up TimeHubble0
+  TimeHubble0 = CosmologyTable_time[j] * (1.0 - dt) +
+    CosmologyTable_time[j+1] * dt;
+
+  /* Now convert from Time * H0 to code units (see also CosmologyGetUnits). */
+   
+  (*TimeCodeUnits) = TimeHubble0 / (HubbleConstantNow*3.24e-18) / TimeUnits;
+   
+
+  return SUCCESS;
+#endif
  
-int CosmologyComputeTimeFromRedshift(FLOAT Redshift, FLOAT *TimeCodeUnits)
-{
- 
-  /* Error check [INCOMPLETE]. */
- 
-  FLOAT eta, TimeHubble0 = FLOAT_UNDEFINED;
+  FLOAT eta;
  
   /* Find Omega due to curvature. */
  
   float OmegaCurvatureNow = 1 - OmegaMatterNow - OmegaLambdaNow;
- 
+
+
   /* 1) For a flat universe with OmegaMatterNow = 1, things are easy. */
  
   if (OmegaMatterNow == 1 && OmegaLambdaNow == 0)
  
 #endif /* INVERSE_HYPERBOLIC_EXISTS */
  
+
   /* 5) Someday, we'll implement the general case... */
  
   if (TimeHubble0 == FLOAT_UNDEFINED) {
     fprintf(stderr, "Cosmology selected is not implemented.\n");
     ENZO_FAIL("");
   }
- 
+  
+
   /* Now convert from Time * H0 to code units (see also CosmologyGetUnits). */
  
-  float TimeUnits = 2.52e17/sqrt(OmegaMatterNow)/HubbleConstantNow/
-                    POW(1 + InitialRedshift, FLOAT(1.5));
+  // FLOAT TimeUnits = 2.52e17/sqrt(OmegaMatterNow)/HubbleConstantNow/
+  //                   POW(1 + InitialRedshift, FLOAT(1.5));
  
   *TimeCodeUnits = TimeHubble0 / (HubbleConstantNow*3.24e-18) / TimeUnits;
  
  
   return SUCCESS;
 }
+

File src/enzo/CosmologyParameters.h

 
 EXTERN float OmegaMatterNow;
 
+/* The value of Omega due to relativistic particles at z=0. */
+
+EXTERN float OmegaRadiationNow;
+
 /* The value of Omega due to lamba (the cosmological constant) at z=0. */
 
 EXTERN float OmegaLambdaNow;
 EXTERN FLOAT CosmologyOutputRedshift[MAX_NUMBER_OF_OUTPUT_REDSHIFTS];
 EXTERN char *CosmologyOutputRedshiftName[MAX_NUMBER_OF_OUTPUT_REDSHIFTS];
 EXTERN FLOAT CosmologyOutputRedshiftTime[MAX_NUMBER_OF_OUTPUT_REDSHIFTS];
+
+
+#ifdef USE_COSMOTABLE
+
+EXTERN int CosmologyTableCalculated;
+EXTERN int CosmologyTableSize;
+EXTERN int CosmologyTableWriteToFile;
+
+EXTERN char CosmologyTableFilename[200];
+
+EXTERN FLOAT CosmologyTableDelta;
+EXTERN FLOAT *CosmologyTable_a, *CosmologyTable_dadt, *CosmologyTable_time;
+
+// function prototypes
+int CosmologyCalculateTable(void);
+#endif

File src/enzo/CosmologyReadParameters.C

  
   HubbleConstantNow    = 0.701;
   OmegaMatterNow       = 0.279;
+  OmegaRadiationNow    = 0.0;
   OmegaLambdaNow       = 0.721;
   ComovingBoxSize      = 64;
   MaxExpansionRate     = 0.01;
     CosmologyOutputRedshift[i]     = -1;  // Never!!
     CosmologyOutputRedshiftName[i] = NULL;
   }
+
+
+#ifdef USE_COSMOTABLE
+  CosmologyTableCalculated = FALSE;
+  CosmologyTableWriteToFile = TRUE;
+  CosmologyTableSize = 10000;
+  sprintf(CosmologyTableFilename,"CosmologyTable.out");
+#endif
+
  
   /* read input from file */
  
     ret += sscanf(line, "CosmologyHubbleConstantNow = %"FSYM,
 		  &HubbleConstantNow);
     ret += sscanf(line, "CosmologyOmegaMatterNow = %"FSYM, &OmegaMatterNow);
+    ret += sscanf(line, "CosmologyOmegaRadiationNow = %"FSYM, &OmegaRadiationNow);
     ret += sscanf(line, "CosmologyOmegaLambdaNow = %"FSYM, &OmegaLambdaNow);
     ret += sscanf(line, "CosmologyComovingBoxSize = %"FSYM, &ComovingBoxSize);
     ret += sscanf(line, "CosmologyMaxExpansionRate = %"FSYM,
       fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line);
  
   }
+
+  /* #ifndef USE_COSMOTABLE, check that OmegaRadiationNow=0, since
+      it's not supported in the standard version. */
+#ifndef USE_COSMOTABLE
+  if (OmegaRadiationNow > 0.0) {
+    fprintf(stderr, "OmegaRadiationNow > 0 is not supported without use of a cosmological look-up table. Compile with -DUSE_COSMOTABLE.\n");    
+    ENZO_FAIL("");
+  }
+#endif
  
   /* Initialize by finding the time at the initial redshift. */
  

File src/enzo/CosmologyWriteParameters.C

  
   fprintf(fptr, "CosmologyHubbleConstantNow = %"GSYM"\n", HubbleConstantNow);
   fprintf(fptr, "CosmologyOmegaMatterNow    = %"GSYM"\n", OmegaMatterNow);
+  fprintf(fptr, "CosmologyOmegaRadiationNow = %"GSYM"\n", OmegaRadiationNow);
   fprintf(fptr, "CosmologyOmegaLambdaNow    = %"GSYM"\n", OmegaLambdaNow);
   fprintf(fptr, "CosmologyComovingBoxSize   = %"GSYM"\n", ComovingBoxSize);
   fprintf(fptr, "CosmologyMaxExpansionRate  = %"GSYM"\n", MaxExpansionRate);

File src/enzo/Make.config.assemble

 	$(error Illegal value '$(CONFIG_USE_HDF4)' for $$(CONFIG_USE_HDF4))
     endif
 
+#=======================================================================
+# DETERMINE COSMOTABLE USAGE
+#=======================================================================
+
+    ERROR_USE_COSMOTABLE = 1
+
+    # compilers and settings if USE_COSMOTABLE is yes
+
+    ifeq ($(CONFIG_USE_COSMOTABLE),yes)
+        ERROR_USE_COSMOTABLE = 0
+        ASSEMBLE_COSMOTABLE_DEFINES  = -DUSE_COSMOTABLE
+    endif
+
+    # compilers and settings if USE_COSMOTABLE is no
+
+    ifeq ($(CONFIG_USE_COSMOTABLE),no)
+        ERROR_USE_COSMOTABLE = 0
+        ASSEMBLE_COSMOTABLE_DEFINES  =
+    endif
+
+    # error if CONFIG_USE_COSMOTABLE is incorrect
+
+    ifeq ($(ERROR_USE_COSMOTABLE),1)
+       .PHONY: error_compilers
+       error_compilers:
+	$(error Illegal value '$(CONFIG_USE_COSMOTABLE)' for $$(CONFIG_USE_COSMOTABLE))
+    endif
+
 
 #=======================================================================
 # ASSIGN ALL OUTPUT VARIABLES
               $(ASSEMBLE_FLUX_FIX_DEFINES) \
               $(ASSEMBLE_ECUDA_DEFINES) \
               $(ASSEMBLE_CUDADEBUG_DEFINES) \
-              $(ASSEMBLE_HDF4_DEFINES)
+              $(ASSEMBLE_HDF4_DEFINES) \
+              $(ASSEMBLE_COSMOTABLE_DEFINES)
 
     INCLUDES = $(MACH_INCLUDES) \
     	       $(ASSEMBLE_MPI_INCLUDES) \

File src/enzo/Make.config.objects

         CopyOverlappingParticleMassFields.o \
         CopyOverlappingZones.o \
 	CopyZonesFromOldGrids.o \
+	CosmologyCalculateTable.o \
         CosmologyComputeExpansionFactor.o \
         CosmologyComputeExpansionTimestep.o \
         CosmologyComputeTimeFromRedshift.o \

File src/enzo/Make.config.settings

 #----------------------------------------------------------------------- 
  
      CONFIG_ECUDA = no
+
+#======================================================================= 
+# CONFIG_USE_COSMOTABLE
+#======================================================================= 
+#    yes           Use a lookup table for Time(a), da/dt(a).
+#    no            Direct calculation.
+#----------------------------------------------------------------------- 
+ 
+     CONFIG_USE_COSMOTABLE = no

File src/enzo/Make.config.targets

 	@echo "      gmake cuda-yes"
 	@echo "      gmake cuda-no"
 	@echo
+	@echo "   Set whether to use a cosmology lookup table"
+	@echo
+	@echo "      gmake use-cosmotable-yes"
+	@echo "      gmake use-cosmotable-no"
+	@echo
 
 #-----------------------------------------------------------------------
 
 	@echo "   CONFIG_BITWISE_IDENTICALITY:  $(CONFIG_BITWISE_IDENTICALITY)"
 	@echo "   CONFIG_FAST_SIB:              $(CONFIG_FAST_SIB)"
 	@echo "   CONFIG_FLUX_FIX:              $(CONFIG_FLUX_FIX)"
+	@echo "   CONFIG_USE_COSMOTABLE:        $(CONFIG_USE_COSMOTABLE)"
 	@echo
 
 #-----------------------------------------------------------------------
 
 #-----------------------------------------------------------------------
 
+VALID_USE-COSMOTABLE = use-cosmotable-yes use-cosmotable-no
+.PHONY: $(VALID_USE-COSMOTABLE)
+
+use-cosmotable-yes: CONFIG_USE-COSMOTABLE-yes
+use-cosmotable-no: CONFIG_USE-COSMOTABLE-no
+use-cosmotable-%:
+	@printf "\n\tInvalid target: $@\n\n\tValid targets: [$(VALID_USE-COSMOTABLE)]\n\n"
+CONFIG_USE-COSMOTABLE-%: suggest-clean
+	@tmp=.config.temp; \
+        grep -v CONFIG_USE_COSMOTABLE $(MAKE_CONFIG_OVERRIDE) > $${tmp}; \
+        mv $${tmp} $(MAKE_CONFIG_OVERRIDE); \
+        echo "CONFIG_USE_COSMOTABLE = $*" >> $(MAKE_CONFIG_OVERRIDE); \
+	$(MAKE)  show-config | grep CONFIG_USE_COSMOTABLE:
+	@echo
+
+
+#-----------------------------------------------------------------------
+