sphere.cpp 8.5 KB
``````#include "sphere.h"
#include "defaults.h"

#include "rts/math/complex.h"
#include <complex>
#include <stdlib.h>
#include <fstream>

//using namespace rts;
using namespace std;

int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);

int cbessjyva_sph(int v,complex<double> z,double &vm,complex<double>*cjv,
complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);

int bessjyv_sph(int v, double z, double &vm, double* cjv,
double* cyv, double* cjvp, double* cyvp);

void sphere::calcCoeff(ptype lambda, bsComplex ri)
{
/*  These calculations are done at high-precision on the CPU
since they are only required once for each wavelength.

Input:

lambda  =   wavelength of the incident field
n       =   complex refractive index of the sphere
*/

//clear the previous coefficients
A.clear();
B.clear();

//convert to an std complex value
complex<double> nc(ri.real(), ri.imag());
n = ri;

//compute the magnitude of the k vector
double k = 2 * PI / lambda;
complex<double> kna = nc * k * (double)a;

//compute the arguments k*a and k*n*a
complex<double> ka(k * a, 0.0);

//allocate space for the Bessel functions of the first and second kind (and derivatives)
int bytes = sizeof(complex<double>) * (Nl + 1);
complex<double>* cjv_ka = (complex<double>*)malloc(bytes);
complex<double>* cyv_ka = (complex<double>*)malloc(bytes);
complex<double>* cjvp_ka = (complex<double>*)malloc(bytes);
complex<double>* cyvp_ka = (complex<double>*)malloc(bytes);
complex<double>* cjv_kna = (complex<double>*)malloc(bytes);
complex<double>* cyv_kna = (complex<double>*)malloc(bytes);
complex<double>* cjvp_kna = (complex<double>*)malloc(bytes);
complex<double>* cyvp_kna = (complex<double>*)malloc(bytes);

//allocate space for the spherical Hankel functions and derivative
complex<double>* chv_ka = (complex<double>*)malloc(bytes);
complex<double>* chvp_ka = (complex<double>*)malloc(bytes);

//compute the bessel functions using the CPU-based algorithm
double vm;
cbessjyva_sph(Nl, ka, vm, cjv_ka, cyv_ka, cjvp_ka, cyvp_ka);
cbessjyva_sph(Nl, kna, vm, cjv_kna, cyv_kna, cjvp_kna, cyvp_kna);

//compute A for each order
complex<double> i(0, 1);
complex<double> a, b, c, d;
complex<double> An, Bn;
for(int l=0; l<=Nl; l++)
{
//compute the Hankel functions from j and y
chv_ka[l] = cjv_ka[l] + i * cyv_ka[l];
chvp_ka[l] = cjvp_ka[l] + i * cyvp_ka[l];

//Compute A (internal scattering coefficient)
//compute the numerator and denominator for A
a = cjv_ka[l] * chvp_ka[l] - cjvp_ka[l] * chv_ka[l];
b = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;

//calculate A and add it to the list
An = (2.0 * l + 1.0) * pow(i, l) * (a / b);
A.push_back(bsComplex(An.real(), An.imag()));

//Compute B (external scattering coefficient)
c = cjv_ka[l] * cjvp_kna[l] * nc - cjv_kna[l] * cjvp_ka[l];
d = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;

//calculate B and add it to the list
Bn = (2.0 * l + 1.0) * pow(i, l) * (c / d);
B.push_back(bsComplex(Bn.real(), Bn.imag()));

}
}

void sphere::calcBesselLut(bsComplex* j, ptype k, bsComplex n, int aR)
{
/*Compute the look-up-table for spherical bessel functions used inside of the sphere
j    =   (Nl + 1) x aR array of values
aR      =   resolution of j
*/

//allocate space for the Bessel functions of the first and second kind (and derivatives -- which will be ignored)
int bytes = sizeof(complex<double>) * (Nl + 1);
complex<double>* cjv_knr = (complex<double>*)malloc(bytes);
complex<double>* cyv_knr = (complex<double>*)malloc(bytes);
complex<double>* cjvp_knr = (complex<double>*)malloc(bytes);
complex<double>* cyvp_knr = (complex<double>*)malloc(bytes);

//compute the bessel functions using the CPU-based algorithm
double vm;

//for each sample along r
ptype dr = a / (aR - 1);
ptype r;
for(int ir = 0; ir < aR; ir++)
{
r = ir * dr;
complex<double> knr( (k*n*r).real(), (k*n*r).imag() );
cbessjyva_sph(Nl, knr, vm, cjv_knr, cyv_knr, cjvp_knr, cyvp_knr);

//copy the double data to the bsComplex array
for(int l=0; l<=Nl; l++)
{
//deal with the NaN case at the origin
if(ir == 0)
{
if(l == 0)
j[ir * (Nl+1)] = 1;
else
j[ir * (Nl+1) + l] = 0;
}
else
j[ir * (Nl+1) + l] = bsComplex(cjv_knr[l].real(), cjv_knr[l].imag());
}
}

/*ofstream outfile("besselout.txt");
for(int ir = 0; ir < aR; ir++)
{
for(int l = 0; l<Nl+1; l++)
{
outfile<<j[ir * (Nl+1) + l].real()<<"     ";
}
outfile<<endl;
}
outfile.close();*/

}

void sphere::calcHankelLut(bsComplex* h, ptype k, int rR)
{
/*Compute the look-up-table for spherical bessel functions used inside of the sphere
h_out   =   (Nl + 1) x aR array of values
rmin	=	minimum value of r
d_max	=	maximum value of r
rR      =   resolution of h_out
*/

//allocate space for the Bessel functions of the first and second kind (and derivatives -- which will be ignored)
int bytes = sizeof(double) * (Nl + 1);
double* cjv_kr = (double*)malloc(bytes);
double* cyv_kr = (double*)malloc(bytes);
double* cjvp_kr = (double*)malloc(bytes);
double* cyvp_kr = (double*)malloc(bytes);

//compute the bessel functions using the CPU-based algorithm
double vm;

//for each sample along r
ptype dr = (d_max - max(a, d_min)) / (rR - 1);
ptype r;
for(int ir = 0; ir < rR; ir++)
{
r = ir * dr + max(a, d_min);
double kr = k*r;
bessjyv_sph(Nl, kr, vm, cjv_kr, cyv_kr, cjvp_kr, cyvp_kr);

//copy the double data to the bsComplex array
for(int l=0; l<=Nl; l++)
{
//h[ir * (Nl+1) + l] = bsComplex(cjv_kr[l].real(), cyv_kr[l].real());
h[ir * (Nl+1) + l] = bsComplex(cjv_kr[l], cyv_kr[l]);
}
}

/*ofstream outfile("hankelout.txt");
for(int ir = 0; ir < rR; ir++)
{
outfile<<ir*dr + max(a, d_min)<<"     ";
for(int l = 0; l<=0; l++)
{
outfile<<h[ir * (Nl+1) + l].real()<<"     "<<h[ir * (Nl+1) + l].imag()<<"     ";
}
outfile<<endl;
}
outfile.close();*/
}

void sphere::calcLut(bsComplex* j, bsComplex* h, ptype lambda, bsComplex n, int aR, int rR)
{
/*Compute the look-up-tables for spherical bessel functions used both inside and outside of the sphere.
j       =   (Nl + 1) x aR array of values
j       =   (Nl + 1) x rR array of values
d_max    =   maximum distance for the LUT
aR      =   resolution of j_in
rR      =   resolution of j_out
*/

//compute the magnitude of the k vector
double k = 2 * PI / lambda;

calcBesselLut(j, k, n, aR);
calcHankelLut(h, k, rR);
}

void sphere::calcUp(ptype lambda, bsComplex n, rts::quad<ptype, 3> nfPlane, unsigned int R)
{
//calculate the parameters of the lookup table

//first find the distance to the closest and furthest points on the nearfield plane
d_min = nfPlane.dist(p);
d_max = nfPlane.dist_max(p);

//compute the radius of the cross-section of the sphere with the plane
ptype a_inter = 0;
if(d_min < a)
a_inter = sqrt(a - d_min);

//calculate the resolution of the Usp and Uip lookup tables
int aR = 1 + 2 * R * a_inter / (nfPlane(0, 0) - nfPlane(1, 1)).len();
int dR = 2 * R;
int thetaR = DEFAULT_SPHERE_THETA_R;

//allocate space for the bessel function LUTs
bsComplex* j = (bsComplex*)malloc(sizeof(bsComplex) * (Nl + 1) * aR);
bsComplex* h = (bsComplex*)malloc(sizeof(bsComplex) * (Nl + 1) * dR);

calcLut(j, h, lambda, n, aR, dR);

//allocate space for the Usp lookup texture
Usp.R[0] = dR;
Usp.R[1] = thetaR;
Usp.init_gpu();

//allocate space for the Uip lookup texture
Uip.R[0] = aR;
Uip.R[1] = thetaR;
Uip.init_gpu();

scalarUsp(h, dR, thetaR);
scalarUip(j, aR, thetaR);

scalarslice UspMag = Usp.Mag();
UspMag.toImage("Usp.bmp", true);

scalarslice UipMag = Uip.Mag();
UipMag.toImage("Uip.bmp", true);

//free memory
free(j);
free(h);

}

sphere& sphere::operator=(const sphere &rhs)
{
p = rhs.p;
a = rhs.a;
iMaterial = rhs.iMaterial;
Nl = rhs.Nl;
n = rhs.n;
B = rhs.B;
A = rhs.A;

return *this;
}

sphere::sphere(const sphere &rhs)
{
p = rhs.p;
a = rhs.a;
iMaterial = rhs.iMaterial;
Nl = rhs.Nl;
n = rhs.n;
B = rhs.B;
A = rhs.A;
}``````