Jump to content
  • Advertisement
isu diss

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);
	}
}

 

Share this post


Link to post
Share on other sites
Advertisement

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


Link to 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


Link to 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


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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
×

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!