> 		lp = (1 + n * c) * l;
> 		r = lp * lp * 0.25 / sqrt(lp * lp - l * l);
> 		do
> 		{
> 			double temp = lp / (2 * r);
> 			rp = r;
> 			r += r * (l - 2 * r * sin(temp)) / (l - lp * cos(temp));
> 		}
> 		while(fabs(r - rp) > r * FLT_EPSILON);
> 		printf("%.3lf\n", r - sqrt(r * r - l * l * 0.25));