Jump to content
• Advertisement

# Mie Scattering

## Recommended Posts

I implemented the paper, Numerical Methods for Mie Theory of Scattering by a Sphere, Link:https://prints.iiap.res.in/bitstream/2248/72/1/Numerical Methods for Mie Theory of Scattering by a Sphere.pdf

Unfortunately, my implementation is not accurate. I couldn't find any mistakes in my code. Can anyone help me?

struct XMDOUBLE2
{
double x;
double y;

XMDOUBLE2() {}
XMDOUBLE2(double _x, double _y) : x(_x), y(_y) {}
explicit XMDOUBLE2(_In_reads_(2) const double *pArray) : x(pArray[0]), y(pArray[1]) {}

XMDOUBLE2& operator= (const XMFLOAT2& Float2) { x = Float2.x; y = Float2.y; return *this; }
};

XMDOUBLE2 Complex_Add(XMDOUBLE2 z1, XMDOUBLE2 z2)// z1 + z2
{
return XMDOUBLE2((z1.x + z2.x), (z1.y + z2.y));
}

XMDOUBLE2 Complex_Subtract(XMDOUBLE2 z1, XMDOUBLE2 z2)// z1 - z2
{
return XMDOUBLE2((z1.x - z2.x), (z1.y - z2.y));
}

XMDOUBLE2 Complex_Multiply(XMDOUBLE2 z1, XMDOUBLE2 z2)// z1*z2
{
return XMDOUBLE2(((z1.x*z2.x) - (z1.y*z2.y)), ((z1.x*z2.y) + (z1.y*z2.x)));
}

double Complex_Norm(XMDOUBLE2 z)// |z|
{
return sqrt((z.x*z.x) + (z.y*z.y));
}

XMDOUBLE2 Complex_Division(XMDOUBLE2 z1, XMDOUBLE2 z2)// z1/z2
{
XMDOUBLE2 tmp;
if (Complex_Norm(z2) != 0)
{
tmp.x = ((z1.x*z2.x)+(z1.y*z2.y))/((z2.x*z2.x)+(z2.y*z2.y));
tmp.y = ((z1.y*z2.x)-(z1.x*z2.y))/((z2.x*z2.x)+(z2.y*z2.y));
}
return tmp;
}

//m = 1.5+i0.0
void MieCoefficient(XMDOUBLE2 m, float r, float Lamda)//m = (m')+i(-m"): complex index of refraction, r: radius of the spherical particle in um, lamda: wavelength in nm
{
#define TOTAL_MAX_TERMS 2100
double Alpha[TOTAL_MAX_TERMS], Beta[TOTAL_MAX_TERMS], P[TOTAL_MAX_TERMS], Q[TOTAL_MAX_TERMS], R[TOTAL_MAX_TERMS], S[TOTAL_MAX_TERMS], B[TOTAL_MAX_TERMS], C[TOTAL_MAX_TERMS], C_[TOTAL_MAX_TERMS];
XMDOUBLE2 A[TOTAL_MAX_TERMS], G[TOTAL_MAX_TERMS], H[TOTAL_MAX_TERMS], a[TOTAL_MAX_TERMS], b[TOTAL_MAX_TERMS];

for (int n=0; n<TOTAL_MAX_TERMS; n++)
{
Alpha[n] = 0; Beta[n] = 0; P[n] = 0; Q[n] = 0; R[n] = 0; S[n] = 0; B[n] = 0; C[n] = 0; C_[n] = 0;
A[n] = XMDOUBLE2(0, 0);	G[n] = XMDOUBLE2(0, 0);	H[n] = XMDOUBLE2(0, 0);	a[n] = XMDOUBLE2(0, 0);	b[n] = XMDOUBLE2(0, 0);
}

//	float Lamdaum = (float)Lamda*1e-3;
float x = 23.1f;//((2*XM_PI*r) / Lamdaum);

//int MAX_TERMS = (int)(x + 7.5f*pow(x, .34f) + 2);
int MAX_TERMS = (int)(x + 4*pow(x, 0.33f) + 2); //Wiscombe

if (MAX_TERMS<TOTAL_MAX_TERMS)
{
double y1 = x*m.x;
double y2 = x*m.y;
XMDOUBLE2 z = XMDOUBLE2(y1, y2);
double zNorm = Complex_Norm(z);
double y = zNorm*zNorm;

for (int n=MAX_TERMS; n>0; n--)
{
Alpha[n] = ((((2*(double)(n))-1)*y1) / y) - P[n];
Beta[n] = ((((2*(double)(n))-1)*y2) / y) + Q[n];
P[n-1] = (Alpha[n] / ((Alpha[n]*Alpha[n]) + (Beta[n]*Beta[n])));
Q[n-1] = (Beta[n] / ((Alpha[n]*Alpha[n]) + (Beta[n]*Beta[n])));
R[n-1] = (x / (((2*(double)(n))-1) - (x*R[n])));
}

S[0] = sin(x);
for (int n=1; n<MAX_TERMS; n++)
{
S[n] = R[n]*S[n-1];
}

for (int n=0; n<MAX_TERMS; n++)
{
XMDOUBLE2 t1 = Complex_Division(XMDOUBLE2(1.0f, 0), XMDOUBLE2(P[n], Q[n]));
XMDOUBLE2 t2 = Complex_Division(XMDOUBLE2((double)(n), 0), z);
A[n] = Complex_Subtract(t1, t2);

B[n] = (((1/R[n]) - ((double)(n)/x)));
}

C[0] = cos(x);
C[1] = ((cos(x)/x) + sin(x));
for (int n=2; n<MAX_TERMS; n++)
{
C[n] = (((((2*(double)(n))-1) / x) * C[n-1]) - C[n-2]);
}

C_[0] = 0;
for (int n=1; n<MAX_TERMS; n++)
{
C_[n] = (((((double)(-n)) / x) * C[n]) - C[n-1]);
}

for (int n=0; n<MAX_TERMS; n++)
{
G[n] = XMDOUBLE2(1, (C[n]/S[n]));
H[n] = XMDOUBLE2(B[n], (C_[n]/S[n]));
}

for (int n=0; n<MAX_TERMS; n++)
{
XMDOUBLE2 t1 = Complex_Multiply(m, XMDOUBLE2(B[n], 0));
XMDOUBLE2 t2 = Complex_Multiply(A[n], G[n]);
XMDOUBLE2 t3 = Complex_Multiply(m, H[n]);
a[n] = Complex_Division(Complex_Subtract(A[n], t1), Complex_Subtract(t2, t3));

XMDOUBLE2 t4 = Complex_Multiply(m, A[n]);
XMDOUBLE2 t5 = Complex_Multiply(t4, G[n]);
b[n] = Complex_Division(Complex_Subtract(t4, XMDOUBLE2(B[n], 0)), Complex_Subtract(t5, H[n]));
}

double Qext = 0, tmp = 0;
for (int n=0; n<MAX_TERMS; n++)
{
tmp += (((2*(double)(n))+1) * Complex_Add(a[n], b[n]).x);
}
Qext = ((tmp*2) / (x*x));
char s[256];
sprintf_s(s, "%.25f", Qext);
MessageBoxA(0, s, "", MB_OK);
}
}

Advertisement

pls need help

#### Share this post

##### Share on other sites

A few things that might help:

• "My implementation is not accurate" doesn't give me any hints of where to start exploring your code. Explain what evidence you have that the code doesn't do the right thing. Something like "If I plug in these inputs, I get these outputs, while I expected to get these other outputs."
• You are using way too many parentheses.
• You could use std::complex<double> instead of defining your own type. That way you can be pretty sure that any mistakes are not in your implementation of complex numbers.

#### Share this post

##### Share on other sites

When I input m = XMDOUBLE2(1.5, -0.0) (means 1.5 - i0.0) and x = 23.1 (size-to-wavelength param.) it should output (QExt - extinction efficiency) at least 2.46225831, what I get is 1.89330136775

I changed the code to use std::complex<double> still the result is same.

Edited by isu diss

#### Share this post

##### Share on other sites

I have a hard time reading the type on that article. I can't tell if it says "1.5" or "1.6", for instance. Is there a more modern reference?

#### Share this post

##### Share on other sites

Sorry, no modern reference. That's the only one from that author.

#### Share this post

##### Share on other sites

I implemented another paper "Mie Scattering calculation" by Hong Du. It's now working. Thanks for your effort, alvaro.

## Create an account or sign in to comment

You need to be a member in order to leave a comment

## Create an account

Sign up for a new account in our community. It's easy!

Register a new account

## Sign in

Already have an account? Sign in here.

Sign In Now

• Advertisement
• Advertisement

• ### Popular Contributors

1. 1
Rutin
65
2. 2
3. 3
4. 4
5. 5
• Advertisement

• 18
• 10
• 30
• 20
• 9
• ### Forum Statistics

• Total Topics
633416
• Total Posts
3011771
• ### Who's Online (See full list)

There are no registered users currently online

×

## Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

We are the game development community.

Whether you are an indie, hobbyist, AAA developer, or just trying to learn, GameDev.net is the place for you to learn, share, and connect with the games industry. Learn more About Us or sign up!

Sign me up!