C code for Mie scattering simulation


#include
#include
#include

#define SphereRadius 1.0 // Radius (um) of the sphere
#define SphRefInd_Re 2 // Complex Relative Refractive index of the sphere (real part)
#define SphRefInd_Im 0.1 // Complex Relative Refractive index of the sphere (imag part)
#define MedRefInd_Re 1.3 // Complex Relative Refractive index of the medium (real part)
#define MedRefInd_Im 0.0 // Complex Relative Refractive index of the medium (imag part)
#define Wavelength 0.92 // Wavelength (um)
#define PrintDebug 0 // Set to 1 to print debug info

// Gobal variables

//Form factor
double x = 2 * M_PI * SphereRadius/Wavelength*MedRefInd_Re;

// Max index
int nmax;

// Relative ref. index
complex m = (SphRefInd_Re + SphRefInd_Im * I) / (MedRefInd_Re + MedRefInd_Im * I);

// BesselJ function from SLATEC
extern void zbesj_(double*, double*, double*, int*, int*, double*, double*, int*, int*);

complex besselj(double nu, complex z){
int kode=1;
int n=1;
double zr=creal(z);
double zi=cimag(z);
int nz,ierr;
double cyr[1],cyi[1];
complex res;
zbesj_(&zr,&zi,&nu,&kode,&n,cyr,cyi,&nz,&ierr);
if(ierr!=0){
printf("error!\n");
}
res=cyr[0]+I*cyi[0];
return res;
}

// BesselY function from SLATEC
extern void zbesy_(double*, double*, double*, int*, int*, double*, double*, int*, double*, double*, int*);

complex bessely(double nu, complex z){
int kode=1;
int n=1;
double zr=creal(z);
double zi=cimag(z);
int nz,ierr;
double cyr[1],cyi[1];
double cwkr[1],cwkri[1];
complex res;
zbesy_(&zr,&zi,&nu,&kode,&n,cyr,cyi,&nz,cwkr,cwkri,&ierr);
if(ierr!=0){
printf("error!\n");
}
res=cyr[0]+I*cyi[0];
return res;
}

void Mie_abcd(complex m, double x, complex abcd[nmax][2]) {

// Computes Mie coefficients an and bn

int i;
double n;
complex jn_x;
complex jn_1_x;
complex jn_mx;
complex jn_1_mx;
complex yn_x;
complex yn_1_x;

complex an;
complex bn;
complex cn;
complex dn;

complex d_xjn_x;
complex d_mxjn_mx;
complex hn_x;
complex d_xhn_x;

double sqx = sqrt(M_PI/(2*x));

for (i=0; i