double f4 ( double x ) 
{
  double alpha = 0.138E-06; 
  double t = 60 * 24 * 60 * 60;
  double temp_init = 20.0; 
  double temp_cold = -15.0; 
  double value;

  value = temp_cold + ( temp_init - temp_cold ) * erf ( 0.5 * x / sqrt ( alpha * t ) );

  return value;
}