#include<cmath>
using namespace std;

#define NLAYERS_BCAL 10
double BCAL_tdiff(double geometric_mean, int layer)
{
	double tw_par[NLAYERS_BCAL][7] ={
		{-1287.04, 1371.24, -603.738, 140.589, -18.2519, 1.25202, -0.0354467},	// layer 1
		{201.506, -141.801, 32.9462, -1.38203, -0.564504, 0.0850295, -0.00359389},	// layer 2
		{20.1508, -16.8756, 4.469, -0.164429, -0.111515, 0.0177557, -0.000809127},	// layer 3
		{183.991, -123.778, 24.7619, 0.719147, -0.886138, 0.112469, -0.00459028},	// layer 4
		{-148.837, 105.97, -21.6155, -1.24335, 1.03665, -0.135614, 0.00580366},	// layer 5
		{-519.598, 398.161, -88.4846, -4.97345, 4.83967, -0.695167, 0.0326592},	// layer 6
		{-9736.55, 8762.01, -2405.99, -49.603, 140.351, -24.5959, 1.3708},	// layer 7
		{-15505.4, 14613.7, -4178.08, -92.2772, 262.486, -47.3609, 2.71171},	// layer 8
		{1283.3, -620.292, -8.03484, 30.6972, 2.55503, -1.93823, 0.168936},	// layer 9
		{674.223, -346.91, -0.571775, 17.9206, 1.40735, -1.21613, 0.112509}	// layer 10
		};

	layer = layer<1 ? 1:layer>NLAYERS_BCAL ? NLAYERS_BCAL:layer;
	double *p = tw_par[layer-1];
	double x = log(geometric_mean);
	return p[0]+x*(p[1]+x*(p[2]+x*(p[3]+x*(p[4]+x*(p[5]+x*p[6])))));
}