#include<cmath>
using namespace std;

#define NLAYERS_BCAL 10
double BCAL_tdiff(double geometric_mean, int layer)
{
	double tw_par[NLAYERS_BCAL][7] ={
		{125.58, -90.4879, 21.2074, -0.92698, -0.340234, 0.0503592, -0.00206276},	// layer 1
		{82.6629, -67.0763, 17.0693, -0.652289, -0.367534, 0.0560427, -0.0024239},	// layer 2
		{138.512, -133.418, 36.8925, -0.367645, -1.36581, 0.210593, -0.0097389},	// layer 3
		{880.116, -714.243, 168.36, 7.09186, -8.72353, 1.27494, -0.0599132},	// layer 4
		{187.51, -95.6733, 1.34752, 4.27537, 0.233895, -0.243178, 0.0215259},	// layer 5
		{268.245, -138.18, 0.746405, 6.5842, 0.442014, -0.399537, 0.0354304},	// layer 6
		{463.635, -222.942, -1.97731, 10.0412, 0.851, -0.578803, 0.0462016},	// layer 7
		{204.808, -105.725, -0.121178, 5.39634, 0.398771, -0.358147, 0.0335638},	// layer 8
		{1.68181, 0.44492, 0.117702, 0.0311379, 0.00823744, 0.00217919, 0.0005765},	// layer 9
		{-25.0387, 20.2092, 3.00509, -0.604143, -0.385751, -0.0645002, 0.025439}	// 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])))));
}