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

pls need help

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.

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.

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?

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

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