Main Page | Alphabetical List | Data Structures | File List | Data Fields | Globals | Related Pages

driftfac.c

Go to the documentation of this file.
00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <mpi.h>
00006 #include <gsl/gsl_math.h>
00007 #include <gsl/gsl_integration.h>
00008 
00009 #include "allvars.h"
00010 #include "proto.h"
00011 
00016 static double logTimeBegin;
00017 static double logTimeMax;
00018 
00019 
00026 void init_drift_table(void)
00027 {
00028 #define WORKSIZE 100000
00029   int i;
00030   double result, abserr;
00031   gsl_function F;
00032   gsl_integration_workspace *workspace;
00033 
00034   logTimeBegin = log(All.TimeBegin);
00035   logTimeMax = log(All.TimeMax);
00036 
00037   workspace = gsl_integration_workspace_alloc(WORKSIZE);
00038 
00039   for(i = 0; i < DRIFT_TABLE_LENGTH; i++)
00040     {
00041       F.function = &drift_integ;
00042       gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), All.Hubble,  /* note: absolute error just a dummy */
00043                           1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
00044       DriftTable[i] = result;
00045 
00046 
00047       F.function = &gravkick_integ;
00048       gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), All.Hubble,  /* note: absolute error just a dummy */
00049                           1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
00050       GravKickTable[i] = result;
00051 
00052 
00053       F.function = &hydrokick_integ;
00054       gsl_integration_qag(&F, exp(logTimeBegin), exp(logTimeBegin + ((logTimeMax - logTimeBegin) / DRIFT_TABLE_LENGTH) * (i + 1)), All.Hubble,  /* note: absolute error just a dummy */
00055                           1.0e-8, WORKSIZE, GSL_INTEG_GAUSS41, workspace, &result, &abserr);
00056       HydroKickTable[i] = result;
00057     }
00058 
00059   gsl_integration_workspace_free(workspace);
00060 }
00061 
00062 
00067 double get_drift_factor(int time0, int time1)
00068 {
00069   double a1, a2, df1, df2, u1, u2;
00070   int i1, i2;
00071 
00072   /* note: will only be called for cosmological integration */
00073 
00074   a1 = logTimeBegin + time0 * All.Timebase_interval;
00075   a2 = logTimeBegin + time1 * All.Timebase_interval;
00076 
00077   u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
00078   i1 = (int) u1;
00079   if(i1 >= DRIFT_TABLE_LENGTH)
00080     i1 = DRIFT_TABLE_LENGTH - 1;
00081 
00082   if(i1 <= 1)
00083     df1 = u1 * DriftTable[0];
00084   else
00085     df1 = DriftTable[i1 - 1] + (DriftTable[i1] - DriftTable[i1 - 1]) * (u1 - i1);
00086 
00087 
00088   u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
00089   i2 = (int) u2;
00090   if(i2 >= DRIFT_TABLE_LENGTH)
00091     i2 = DRIFT_TABLE_LENGTH - 1;
00092 
00093   if(i2 <= 1)
00094     df2 = u2 * DriftTable[0];
00095   else
00096     df2 = DriftTable[i2 - 1] + (DriftTable[i2] - DriftTable[i2 - 1]) * (u2 - i2);
00097 
00098   return df2 - df1;
00099 }
00100 
00101 
00105 double get_gravkick_factor(int time0, int time1)
00106 {
00107   double a1, a2, df1, df2, u1, u2;
00108   int i1, i2;
00109 
00110   /* note: will only be called for cosmological integration */
00111 
00112   a1 = logTimeBegin + time0 * All.Timebase_interval;
00113   a2 = logTimeBegin + time1 * All.Timebase_interval;
00114 
00115   u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
00116   i1 = (int) u1;
00117   if(i1 >= DRIFT_TABLE_LENGTH)
00118     i1 = DRIFT_TABLE_LENGTH - 1;
00119 
00120   if(i1 <= 1)
00121     df1 = u1 * GravKickTable[0];
00122   else
00123     df1 = GravKickTable[i1 - 1] + (GravKickTable[i1] - GravKickTable[i1 - 1]) * (u1 - i1);
00124 
00125 
00126   u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
00127   i2 = (int) u2;
00128   if(i2 >= DRIFT_TABLE_LENGTH)
00129     i2 = DRIFT_TABLE_LENGTH - 1;
00130 
00131   if(i2 <= 1)
00132     df2 = u2 * GravKickTable[0];
00133   else
00134     df2 = GravKickTable[i2 - 1] + (GravKickTable[i2] - GravKickTable[i2 - 1]) * (u2 - i2);
00135 
00136   return df2 - df1;
00137 }
00138 
00142 double get_hydrokick_factor(int time0, int time1)
00143 {
00144   double a1, a2, df1, df2, u1, u2;
00145   int i1, i2;
00146 
00147   /* note: will only be called for cosmological integration */
00148 
00149   a1 = logTimeBegin + time0 * All.Timebase_interval;
00150   a2 = logTimeBegin + time1 * All.Timebase_interval;
00151 
00152   u1 = (a1 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
00153   i1 = (int) u1;
00154   if(i1 >= DRIFT_TABLE_LENGTH)
00155     i1 = DRIFT_TABLE_LENGTH - 1;
00156 
00157   if(i1 <= 1)
00158     df1 = u1 * HydroKickTable[0];
00159   else
00160     df1 = HydroKickTable[i1 - 1] + (HydroKickTable[i1] - HydroKickTable[i1 - 1]) * (u1 - i1);
00161 
00162 
00163   u2 = (a2 - logTimeBegin) / (logTimeMax - logTimeBegin) * DRIFT_TABLE_LENGTH;
00164   i2 = (int) u2;
00165   if(i2 >= DRIFT_TABLE_LENGTH)
00166     i2 = DRIFT_TABLE_LENGTH - 1;
00167 
00168   if(i2 <= 1)
00169     df2 = u2 * HydroKickTable[0];
00170   else
00171     df2 = HydroKickTable[i2 - 1] + (HydroKickTable[i2] - HydroKickTable[i2 - 1]) * (u2 - i2);
00172 
00173   return df2 - df1;
00174 }
00175 
00176 
00179 double drift_integ(double a, void *param)
00180 {
00181   double h;
00182 
00183   h = All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda;
00184   h = All.Hubble * sqrt(h);
00185 
00186   return 1 / (h * a * a * a);
00187 }
00188 
00191 double gravkick_integ(double a, void *param)
00192 {
00193   double h;
00194 
00195   h = All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda;
00196   h = All.Hubble * sqrt(h);
00197 
00198   return 1 / (h * a * a);
00199 }
00200 
00201 
00204 double hydrokick_integ(double a, void *param)
00205 {
00206   double h;
00207 
00208   h = All.Omega0 / (a * a * a) + (1 - All.Omega0 - All.OmegaLambda) / (a * a) + All.OmegaLambda;
00209   h = All.Hubble * sqrt(h);
00210 
00211   return 1 / (h * pow(a, 3 * GAMMA_MINUS1) * a);
00212 }
00213 
00214 double growthfactor_integ(double a, void *param)
00215 {
00216   double s;
00217 
00218   s = All.Omega0 + (1 - All.Omega0 - All.OmegaLambda) * a + All.OmegaLambda * a * a * a;
00219   s = sqrt(s);
00220 
00221   return pow(sqrt(a) / s, 3);
00222 }
00223 
00224 

Generated on Sun May 22 17:33:28 2005 for GADGET-2 by  doxygen 1.3.9.1