Blame view

tira/optics_old/esphere.cuh 5.33 KB
ce6381d7   David Mayerich   updating to TIRA
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
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
  #ifndef RTS_ESPHERE
  #define RTS_ESPHERE
  
  #include "../math/complex.h"
  #include "../math/bessel.h"
  #include "../visualization/colormap.h"
  #include "../optics/planewave.h"
  #include "../cuda/devices.h"
  #include "../optics/efield.cuh"
  
  namespace stim{
  
  /*  This class implements a discrete representation of an electromagnetic field
      in 2D scattered by a sphere. This class implements Mie scattering.
  */
  template<typename P>
  class esphere : public efield<P>
  {
  private:
  	stim::complex<P> n;		//sphere refractive index
  	P a;					//sphere radius
  
  	//parameters dependent on wavelength
  	unsigned int Nl;		//number of orders for the calculation
  	rts::complex<P>* A;		//internal scattering coefficients
  	rts::complex<P>* B;		//external scattering coefficients
  
  	void calcNl(P kmag)
  	{
          //return ceil( ((P)6.282 * a) / lambda + 4 * pow( ((P)6.282 * a) / lambda, ((P)1/(P)3)) + 2);
  		Nl = ceil( kmag*a + 4 * pow(kmag * a, (P)1/(P)3) + 2);
  	}
  
  	void calcAB(P k, unsigned int Nl, rts::complex<P>* A, rts::complex<P>* B)
  	{
  		/*  These calculations require double precision, so they are computed
  			using doubles and converted to P at the end.
  
  			Input:
  
  			k	=	magnitude of the k vector (tau/lambda)
  			Nl	=	highest order coefficient ([0 Nl] are computed)
  		*/
  
  		//clear the previous coefficients
  		rts::complex<P>* cpuA = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));
  		rts::complex<P>* cpuB = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));
  
  		//convert to an std complex value
  		complex<double> nc = (rts::complex<double>)n;
  
  		//compute the magnitude of the k vector
  		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)
  		unsigned 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
  			rts::complex<double> An = (2.0 * l + 1.0) * pow(i, l) * (a / b);
  			cpuA[l] = (rts::complex<P>)An;
  
  			//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
  			rts::complex<double> Bn = (2.0 * l + 1.0) * pow(i, l) * (c / d);
  			cpuB[l] = (rts::complex<P>)Bn;
  
  			std::cout<<"A: "<<cpuA[l]<<"     B: "<<cpuB[l]<<std::endl;
  		}
  
  		
  		if(A != NULL)	cudaFree(A);	//free any previous coefficients
  		if(B != NULL)	cudaFree(B);
  		cudaMalloc(&A, sizeof(rts::complex<P>) * (Nl+1));	//allocate memory for new coefficients
  		cudaMalloc(&B, sizeof(rts::complex<P>) * (Nl+1));
  
  		//copy the calculations from the CPU to the GPU
  		cudaMemcpy(A, cpuA, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
  		cudaMemcpy(B, cpuB, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
  	}
  
  public:
  
  	esphere(unsigned int res0, unsigned int res1, P _a, rts::complex<P> _n, bool _scalar = false) :
  		efield<P>(res0, res1, _scalar)
  	{
  		std::cout<<"Sphere scattered field created."<<std::endl;
  		n = _n;		//save the refractive index
  		a = _a;		//save the radius
  	}
  
  	//assignment operator: build an electric field from a plane wave
      efield<P> & operator= (const planewave<P> & rhs)
  	{
  		calcNl(rhs.kmag());		//compute the number of scattering coefficients
  		std::cout<<"Nl: "<<Nl<<std::endl;
  		calcAB(rhs.kmag(), Nl, A, B);	//compute the scattering coefficients
  
  		//determine important parameters for the scattering domain
  		unsigned int sR = ceil(sqrt( (P)(pow((P)esphere::R[0],2) + pow((P)esphere::R[1],2))) );
  		//unsigned int thetaR = 256;
  
  		/////////////////////continue scattering code here/////////////////////////
  
  		esphere::clear();				//clear any previous field data
  		from_planewave(rhs);	//create a field from the planewave
  		return *this;			//return the current object
  	}
  
  	string str()
  	{
  		stringstream ss;
  		ss<<"Mie Scattered Field"<<std::endl;
  		ss<<(*this).efield<P>::str()<<std::endl;
  		ss<<"a = "<<a<<std::endl;
  		ss<<"n = "<<n;
  
  		return ss.str();
  	}
  
  };
  
  }	//end namespace rts
  
  #endif