When observing the reflected light of object surface, the intensity of the 36 images obtained by rotating the linear polarizer from 0 degree to 175 degree with 5 degree interval forms a sinusoidal curve. Here is the program which calculates the sine curve with least square method.
const double PI2 =
6.28318530717958647692528676655901;
const double EPS =
1.0e-15;
double sin2x[36] =
{0.000000000000000,
0.173648177666930,
0.342020143325669,
0.500000000000000,
0.642787609686539,
0.766044443118978,
0.866025403784439,
0.939692620785908,
0.984807753012208,
1.000000000000000,
0.984807753012208,
0.939692620785908,
0.866025403784439,
0.766044443118978,
0.642787609686539,
0.500000000000000,
0.342020143325669,
0.173648177666931,
0.000000000000000,
-0.173648177666930,
-0.342020143325669,
-0.500000000000000,
-0.642787609686539,
-0.766044443118978,
-0.866025403784438,
-0.939692620785908,
-0.984807753012208,
-1.000000000000000,
-0.984807753012208,
-0.939692620785908,
-0.866025403784439,
-0.766044443118978,
-0.642787609686540,
-0.500000000000000,
-0.342020143325669,
-0.173648177666930};
double cos2x[36] =
{1.000000000000000,
0.984807753012208,
0.939692620785908,
0.866025403784439,
0.766044443118978,
0.642787609686539,
0.500000000000000,
0.342020143325669,
0.173648177666930,
0.000000000000000,
-0.173648177666930,
-0.342020143325669,
-0.500000000000000,
-0.642787609686539,
-0.766044443118978,
-0.866025403784439,
-0.939692620785908,
-0.984807753012208,
-1.000000000000000,
-0.984807753012208,
-0.939692620785908,
-0.866025403784439,
-0.766044443118978,
-0.642787609686539,
-0.500000000000000,
-0.342020143325669,
-0.173648177666930,
-0.000000000000000,
0.173648177666930,
0.342020143325669,
0.500000000000000,
0.642787609686539,
0.766044443118978,
0.866025403784438,
0.939692620785908,
0.984807753012208};
/*******************************************************************/
// y = asin2x + bcos2x + c
// S = sum(y - asin2x - bcos2x - c)^2 = min
// y = asin2x + bcos2x + c = Asin(2x - 2B) + C = Asin2(x - B) + C
// A = sqrt(a^2 + b^2)
// sin(-2B) = b / sqrt(a^2 + b^2)
// cos(-2B) = a / sqrt(a^2 + b^2)
// sum1 = sum
// sumy = sum(y)
// sums = sum(sin2x)
// sumc = sum(cos2x)
// sumys = sum(ysin2x)
// sumyc = sum(ycos2x)
// sumss = sum((sin2x)^2)
// sumsc = sum(sin2x cos2x)
// sumcc = sum((cos2x)^2)
/*
numeratora
= sums * sumc * sumyc
+ sumy * sumc * sumsc
+ sum1 * sumys * sumcc
- sum1 * sumyc * sumsc
- sumc * sumc * sumys
- sumy * sums * sumcc;
numeratorb
= sum1 * sumyc * sumss
+ sums * sumc * sumys
+ sumy * sums * sumsc
- sums * sums * sumyc
- sumy * sumc * sumss
- sum1 * sumys * sumsc;
numeratorc
= sums * sumsc * sumyc
+ sumc * sumsc * sumys
+ sumy * sumss * sumcc
- sumc * sumss * sumyc
- sumy * sumsc * sumsc
- sums * sumcc * sumys;
denominator
= 2.0 * sums * sumc * sumsc
+ sum1 * sumss * sumcc
- sums * sums * sumcc
- sumc * sumc * sumss
- sum1 * sumsc * sumsc;
*/
/*
a = numeratora / denominator;
b = numeratorb / denominator;
C = numeratorc / denominator;
*/
// 36: 0, 5, 10, 15, ... 175
// 18: 0, 10, 20, 30, ... 170
// sum1 = 36 18
// sumy = ? ?
// sums = 0 0
// sumc = 0 0
// sumys = ? ?
// sumyc = ? ?
// sumss = 18 9
// sumsc = 0 0
// sumcc = 18 9
void CalcSinusoid(double sinusoidy[36], double *resulta, double *resultb, double *resultc)
{
double sum1;
double sumy;
double sumys;
double sumyc;
double sumss;
double sumcc;
double r;
int i;
double y;
double sx;
double cx;
double a, b, A, B, C;
double asinalpha, acosalpha;
double minus2B;
sum1 = 36.0;
sumy = 0.0;
sumys = 0.0;
sumyc = 0.0;
sumss = 18.0;
sumcc = 18.0;
r = 5.0;
//
Calculate sum
for(i = 0; i < 36; i++)
{
y = sinusoidy[i];
sx = sin2x[i];
cx = cos2x[i];
sumy += y;
sumys += y * sx;
sumyc += y * cx;
}
a = sumys / sumss;
b = sumyc / sumcc;
C = sumy / sum1;
if(C == 0.0)
{
A = B = 0.0;
C = EPS;
}
else if(a == 0.0 && b == 0.0)
{
A = B = 0.0;
}
else
{
A = sqrt(a * a + b * b);
asinalpha = asin(b / A);
acosalpha = acos(a / A);
if(asinalpha < 0.0) minus2B = PI2 -
acosalpha;
else minus2B = acosalpha;
B = -0.5 * minus2B;
}
//
Return A, B, C
*resulta = A;
*resultb = B;
*resultc = C;
}