#include<cmath>
using namespace std;

#define NLAYERS_BCAL 10
double BCAL_tdiff(double geometric_mean, int layer)
{
	double tw_par[NLAYERS_BCAL][7] ={
		{808.586, -693.44, 188.134, -4.26566, -6.52423, 1.0974, -0.0552376},	// layer 1
		{-101.003, 67.9936, -13.5937, -0.158227, 0.371202, -0.0432308, 0.00153581},	// layer 2
		{-235.26, 170.446, -39.4841, 1.09227, 0.870144, -0.125967, 0.00535562},	// layer 3
		{-120.163, 61.2668, -7.13265, -0.841886, 0.178191, -0.000729844, -0.000731371},	// layer 4
		{-266.726, 193.64, -44.5708, 0.96511, 1.07154, -0.153535, 0.00656072},	// layer 5
		{-286.231, 213.037, -50.0289, 0.964621, 1.30463, -0.189916, 0.00828705},	// layer 6
		{2225.2, -2709.76, 1342.64, -347.675, 49.7571, -3.73922, 0.115465},	// layer 7
		{-412.733, 311.077, -72.2212, 0.149959, 2.32796, -0.332787, 0.0146726},	// layer 8
		{40.2953, -68.302, 26.5327, -1.59454, -0.916355, 0.179147, -0.00957573},	// layer 9
		{-298.351, 216.884, -45.1379, -2.14004, 1.99339, -0.256942, 0.010692}	// 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])))));
}