### Blame view

sphere.cpp 8.5 KB
 `1` `````` #include "sphere.h" `````` `2` `````` #include "defaults.h" `````` `3` `````` `````` `4` `````` #include "rts/math/complex.h" `````` ```5 6``` `````` #include #include `````` `7` `````` #include `````` ```8 9 10 11 12 13 14 15 16 17``` `````` using namespace rts; using namespace std; int cbessjyva(double v,complex z,double &vm,complex*cjv, complex*cyv,complex*cjvp,complex*cyvp); int cbessjyva_sph(int v,complex z,double &vm,complex*cjv, complex*cyv,complex*cjvp,complex*cyvp); `````` ```18 19 20``` `````` int bessjyv_sph(int v, double z, double &vm, double* cjv, double* cyv, double* cjvp, double* cyvp); `````` ```21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66``` `````` void sphere::calcCoeff(ptype lambda, rtsComplex 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 nc(ri.real(), ri.imag()); n = ri; //compute the magnitude of the k vector double k = 2 * PI / lambda; complex kna = nc * k * (double)a; //compute the arguments k*a and k*n*a complex ka(k * a, 0.0); //allocate space for the Bessel functions of the first and second kind (and derivatives) int bytes = sizeof(complex) * (Nl + 1); complex* cjv_ka = (complex*)malloc(bytes); complex* cyv_ka = (complex*)malloc(bytes); complex* cjvp_ka = (complex*)malloc(bytes); complex* cyvp_ka = (complex*)malloc(bytes); complex* cjv_kna = (complex*)malloc(bytes); complex* cyv_kna = (complex*)malloc(bytes); complex* cjvp_kna = (complex*)malloc(bytes); complex* cyvp_kna = (complex*)malloc(bytes); //allocate space for the spherical Hankel functions and derivative complex* chv_ka = (complex*)malloc(bytes); complex* chvp_ka = (complex*)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); `````` ```67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84``` `````` //compute A for each order complex i(0, 1); complex a, b, c, d; complex 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())); `````` `85` `````` `````` ```86 87 88 89 90 91 92 93``` `````` //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())); `````` `94` `````` `````` ```95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189``` `````` } } 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) * (Nl + 1); complex* cjv_knr = (complex*)malloc(bytes); complex* cyv_knr = (complex*)malloc(bytes); complex* cjvp_knr = (complex*)malloc(bytes); complex* cyvp_knr = (complex*)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 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 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; `````` `296` `````` } ``````