#include<cmath>
using namespace std;

#define NLAYERS_BCAL 10
double BCAL_tdiff(double geometric_mean, int layer)
{
	double tw_par[NLAYERS_BCAL][7] ={
		{9.99045, -5.03414, 0.117623, 0.202467, 0.00741133, -0.00994784, 0.00085767},	// layer 1
		{-88.5997, 77.7064, -21.9108, 0.713688, 0.736763, -0.130798, 0.00684174},	// layer 2
		{-6.10784, 2.95437, -0.148162, -0.0904252, 0.00032093, 0.00353417, -0.000313743},	// layer 3
		{-15.7033, 8.05419, -0.401011, -0.276547, -0.000258669, 0.011821, -0.00107316},	// layer 4
		{4.5762, -1.96803, 0.0094246, 0.0706581, 0.00323348, -0.0030805, 0.000243647},	// layer 5
		{-1.59688, 0.619361, 0.0106937, -0.0209298, -0.00155328, 0.000883331, -6.05655e-05},	// layer 6
		{2.26819, -1.37097, 0.0946307, 0.0566466, -0.000780712, -0.00279983, 0.000281277},	// layer 7
		{15.0116, -7.52821, 0.172251, 0.308543, 0.0101206, -0.0150841, 0.00131854},	// layer 8
		{11.0285, -5.74203, 0.14336, 0.246885, 0.00794972, -0.0125292, 0.00111506},	// layer 9
		{-7.21862, 3.51847, -0.0602393, -0.145976, -0.00564632, 0.00745769, -0.000663856}	// 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])))));
}