#include<cmath>
using namespace std;

#define NLAYERS_BCAL 10
double BCAL_tavg(double geometric_mean, int layer)
{
	double tw_par[NLAYERS_BCAL][7] ={
		{4.32461, -2.36229, 0.0425649, 0.106322, 0.00534692, -0.0055173, 0.000451825},	// layer 1
		{136.334, -108.793, 26.7735, -0.0335108, -0.982805, 0.149128, -0.00699155},	// layer 2
		{6.32458, -3.87951, 0.274943, 0.144338, -0.00171109, -0.0064233, 0.000598128},	// layer 3
		{6.65447, -4.36303, 0.345083, 0.166557, -0.00309341, -0.00751594, 0.00071634},	// layer 4
		{0.454535, -1.3315, 0.199717, 0.0724794, -0.00308274, -0.00391478, 0.000407283},	// layer 5
		{-0.665457, -0.712152, 0.146551, 0.0553372, -0.00203001, -0.00336454, 0.000351284},	// layer 6
		{0.272997, -1.0827, 0.139123, 0.0682597, -0.00105811, -0.00390159, 0.000386524},	// layer 7
		{-394.96, 344.185, -93.2602, 0.145909, 4.21215, -0.713165, 0.0373399},	// layer 8
		{-476.981, 405.858, -105.675, -1.38943, 5.15214, -0.84351, 0.0434686},	// layer 9
		{-593.79, 521.22, -140.061, -1.6852, 7.1427, -1.20222, 0.0635351}	// 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])))));
}