added look-up tables

dmayerich
1 parent b6179de6
Showing 27 changed files with 1518 additions and 588 deletions Show diff stats
bessjy.cpp
cbessjy.cpp
colormap.h
dataTypes.h
defaults.h
fieldslice.cpp
fieldslice.cu
fieldslice.h
fileout.cu
fileout.h
main.cpp
microscope.cu
microscope.h
montecarlo.cpp
nearfield.cpp
nearfield.h
nfScalarUf.cu
nfScalarUfLut.cu
nfScalarUpLut.cu
nfScalarUs.cu
@@ -13,7 +13,9 @@
 //
 #define _USE_MATH_DEFINES
 #include <math.h>
-#include "bessel.h"
+#include "bessel.h"
+
+#define PI  3.14159
 double gamma(double x);
 //
@@ -426,7 +428,7 @@ int bessjynb(int n,double x,int &amp;nm,double *jn,double *yn,
         0.2775764465332031,
        -1.993531733751297,
         2.724882731126854e1};
-        
+
     int i,k,m;
     nm = n;
     if ((x < 0.0) || (n < 0)) return 1;
@@ -702,5 +704,26 @@ int bessjyv(double v,double x,double &amp;vm,double *jv,double *yv,
     }
     vm = n + v0;
     return 0;
+}
+
+int bessjyv_sph(int v, double z, double &vm, double* cjv,
+    double* cyv, double* cjvp, double* cyvp)
+{
+    //first, compute the bessel functions of fractional order
+    bessjyv(v + 0.5, z, vm, cjv, cyv, cjvp, cyvp);
+
+    //iterate through each and scale
+    for(int n = 0; n<=v; n++)
+    {
+
+        cjv[n] = cjv[n] * sqrt(PI/(z * 2.0));
+        cyv[n] = cyv[n] * sqrt(PI/(z * 2.0));
+
+        cjvp[n] = -1.0 / (z * 2.0) * cjv[n] + cjvp[n] * sqrt(PI / (z * 2.0));
+        cyvp[n] = -1.0 / (z * 2.0) * cyv[n] + cyvp[n] * sqrt(PI / (z * 2.0));
+    }
+
+	return 0;
+
 }
- 
+
@@ -724,6 +724,7 @@ int cbessjyva_sph(int v,complex&lt;double&gt; z,double &amp;vm,complex&lt;double&gt;*cjv,
     //iterate through each and scale
     for(int n = 0; n<=v; n++)
     {
+
         cjv[n] = cjv[n] * sqrt(PI/(z * 2.0));
         cyv[n] = cyv[n] * sqrt(PI/(z * 2.0));
-#ifndef RTS_COLORMAP_H
-#define RTS_COLORMAP_H
-
-#include <string>
-#include <qimage.h>
-#include <qcolor.h>
-#include "rts/cuda/error.h"
-
-
-#define BREWER_CTRL_PTS 11
-
-#ifdef __CUDACC__
-texture<float4, cudaTextureType1D> cudaTexBrewer;
-static cudaArray* gpuBrewer;
-#endif
-
-
-
-namespace rts{
-	namespace colormap{
-
-enum colormapType {cmBrewer, cmGrayscale};
-
-static void buffer2image(unsigned char* buffer, std::string filename, unsigned int x_size, unsigned int y_size)
-{
-	//create an image object
-	QImage image(x_size, y_size, QImage::Format_RGB32);
-
-	int i;
-	unsigned char r, g, b;
-	unsigned int x, y;
-	for(y=0; y<y_size; y++)
-		for(x=0; x<x_size; x++)
-		{
-			//calculate the 1D index
-			i = y * x_size + x;
-
-			r = buffer[i * 3 + 0];
-			g = buffer[i * 3 + 1];
-			b = buffer[i * 3 + 2];
-
-			//set the image pixel
-			QColor color(r, g, b);
-			image.setPixel(x, y, color.rgb());
-		}
-
-	image.save(filename.c_str());
-}
-
-#ifdef __CUDACC__
-static void initBrewer()
-{
-	//initialize the Brewer colormap
-
-	//allocate CPU space
-	float4 cpuColorMap[BREWER_CTRL_PTS];
-
-	//define control rtsPoints
-	cpuColorMap[0] = make_float4(0.192157f, 0.211765f, 0.584314f, 1.0f);
-	cpuColorMap[1] = make_float4(0.270588f, 0.458824f, 0.705882f, 1.0f);
-	cpuColorMap[2] = make_float4(0.454902f, 0.678431f, 0.819608f, 1.0f);
-	cpuColorMap[3] = make_float4(0.670588f, 0.85098f, 0.913725f, 1.0f);
-	cpuColorMap[4] = make_float4(0.878431f, 0.952941f, 0.972549f, 1.0f);
-	cpuColorMap[5] = make_float4(1.0f, 1.0f, 0.74902f, 1.0f);
-	cpuColorMap[6] = make_float4(0.996078f, 0.878431f, 0.564706f, 1.0f);
-	cpuColorMap[7] = make_float4(0.992157f, 0.682353f, 0.380392f, 1.0f);
-	cpuColorMap[8] = make_float4(0.956863f, 0.427451f, 0.262745f, 1.0f);
-	cpuColorMap[9] = make_float4(0.843137f, 0.188235f, 0.152941f, 1.0f);
-	cpuColorMap[10] = make_float4(0.647059f, 0.0f, 0.14902f, 1.0f);
-
-
-	int width = BREWER_CTRL_PTS;
-	int height = 0;
-
-
-	// allocate array and copy colormap data
-	cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 32, 32, 32, cudaChannelFormatKindFloat);
-
-	HANDLE_ERROR(cudaMallocArray(&gpuBrewer, &channelDesc, width, height));
-
-	HANDLE_ERROR(cudaMemcpyToArray(gpuBrewer, 0, 0, cpuColorMap, sizeof(float4)*width, cudaMemcpyHostToDevice));
-
-	// set texture parameters
-    cudaTexBrewer.addressMode[0] = cudaAddressModeClamp;
-	//texBrewer.addressMode[1] = cudaAddressModeClamp;
-    cudaTexBrewer.filterMode = cudaFilterModeLinear;
-    cudaTexBrewer.normalized = true;  // access with normalized texture coordinates
-
-	// Bind the array to the texture
-    HANDLE_ERROR(cudaBindTextureToArray( cudaTexBrewer, gpuBrewer, channelDesc));
-
-}
-
-static void destroyBrewer()
-{
-    HANDLE_ERROR(cudaFreeArray(gpuBrewer));
-
-}
-
-template<class T>
-__global__ static void applyBrewer(T* gpuSource, unsigned char* gpuDest, unsigned int N, T minVal = 0, T maxVal = 1)
-{
-	int i = blockIdx.x * blockDim.x + threadIdx.x;
-    if(i >= N) return;
-
-	//compute the normalized value on [minVal maxVal]
-	float a = (gpuSource[i] - minVal) / (maxVal - minVal);
-
-	//lookup the color
-	float shift = 1.0/BREWER_CTRL_PTS;
-	float4 color = tex1D(cudaTexBrewer, a+shift);
-
-	gpuDest[i * 3 + 0] = 255 * color.x;
-	gpuDest[i * 3 + 1] = 255 * color.y;
-	gpuDest[i * 3 + 2] = 255 * color.z;
-}
-
-template<class T>
-__global__ static void applyGrayscale(T* gpuSource, unsigned char* gpuDest, unsigned int N, T minVal = 0, T maxVal = 1)
-{
-    int i = blockIdx.x * blockDim.x + threadIdx.x;
-    if(i >= N) return;
-
-	//compute the normalized value on [minVal maxVal]
-	float a = (gpuSource[i] - minVal) / (maxVal - minVal);
-
-	gpuDest[i * 3 + 0] = 255 * a;
-	gpuDest[i * 3 + 1] = 255 * a;
-	gpuDest[i * 3 + 2] = 255 * a;
-}
-
-template<class T>
-static void gpu2gpu(T* gpuSource, unsigned char* gpuDest, unsigned int nVals, T minVal = 0, T maxVal = 1, colormapType cm = cmGrayscale, int blockDim = 128)
-{
-	//This function converts a scalar field on the GPU to a color image on the GPU
-	int gridDim = (nVals + blockDim - 1)/blockDim;
-	if(cm == cmGrayscale)
-		applyGrayscale<<<gridDim, blockDim>>>(gpuSource, gpuDest, nVals, minVal, maxVal);
-	else if(cm == cmBrewer)
-	{
-		initBrewer();
-		applyBrewer<<<gridDim, blockDim>>>(gpuSource, gpuDest, nVals, minVal, maxVal);
-		destroyBrewer();
-	}
-
-}
-
-template<class T>
-static void gpu2cpu(T* gpuSource, unsigned char* cpuDest, unsigned int nVals, T minVal, T maxVal, colormapType cm = cmGrayscale)
-{
-    //this function converts a scalar field on the GPU to a color image on the CPU
-
-    //first create the color image on the GPU
-
-    //allocate GPU memory for the color image
-    unsigned char* gpuDest;
-    HANDLE_ERROR(cudaMalloc( (void**)&gpuDest, sizeof(unsigned char) * nVals * 3 ));
-
-	//HANDLE_ERROR(cudaMemset(gpuSource, 0, sizeof(T) * nVals));
-
-    //create the image on the gpu
-    gpu2gpu(gpuSource, gpuDest, nVals, minVal, maxVal, cm);
-
-	//HANDLE_ERROR(cudaMemset(gpuDest, 0, sizeof(unsigned char) * nVals * 3));
-
-    //copy the image from the GPU to the CPU
-    HANDLE_ERROR(cudaMemcpy(cpuDest, gpuDest, sizeof(unsigned char) * nVals * 3, cudaMemcpyDeviceToHost));
-
-	HANDLE_ERROR(cudaFree( gpuDest ));
-
-}
-
-template<typename T>
-static void gpu2image(T* gpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, T valMin, T valMax, colormapType cm = cmGrayscale)
-{
-	//allocate a color buffer
-	unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size);
-
-	//do the mapping
-	gpu2cpu<T>(gpuSource, cpuBuffer, x_size * y_size, valMin, valMax, cm);
-
-	//copy the buffer to an image
-	buffer2image(cpuBuffer, fileDest, x_size, y_size);
-
-	free(cpuBuffer);
-}
-
-#endif
-
-template<class T>
-static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T valMin, T valMax, colormapType cm = cmGrayscale)
-{
-	int i;
-	float a;
-	float range = valMax - valMin;
-	for(i = 0; i<nVals; i++)
-	{
-		//normalize to the range [valMin valMax]
-		a = (cpuSource[i] - valMin) / range;
-
-		cpuDest[i * 3 + 0] = 255 * a;
-		cpuDest[i * 3 + 1] = 255 * a;
-		cpuDest[i * 3 + 2] = 255 * a;
-	}
-
-}
-
-
-
-template<typename T>
-static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, T valMin, T valMax, colormapType cm = cmGrayscale)
-{
-    //allocate a color buffer
-	unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size);
-
-	//do the mapping
-	cpu2cpu<T>(cpuSource, cpuBuffer, x_size * y_size, valMin, valMax, cm);
-
-	//copy the buffer to an image
-	buffer2image(cpuBuffer, fileDest, x_size, y_size);
-
-	free(cpuBuffer);
-
-}
-
-}}	//end namespace colormap and rts
-
-#endif
-
@@ -24,6 +24,8 @@ typedef double ptype;
 typedef ptype fieldPoint;
+extern bool verbose;
+
 //hybrid GPU/CPU complex data typ
 #include "rts/math/complex.h"
 #include "rts/math/vector.h"
@@ -15,14 +15,14 @@
 #define DEFAULT_FOCUS_X         0
 #define DEFAULT_FOCUS_Y         0
 #define DEFAULT_FOCUS_Z         0
-#define DEFAULT_INCIDENT_ORDER	100
+//#define DEFAULT_INCIDENT_ORDER	20
 #define DEFAULT_STABILITY_PARM	1.4
 //optics
-#define DEFAULT_CONDENSER_MIN   0.0
+#define DEFAULT_CONDENSER_MIN   0
 #define DEFAULT_CONDENSER_MAX   1
-#define DEFAULT_OBJECTIVE_MIN   0.0
+#define DEFAULT_OBJECTIVE_MIN   0
 #define DEFAULT_OBJECTIVE_MAX   1
 //incident light direction
@@ -36,17 +36,20 @@
 //#define DEFAULT_OUTPUT_POINT	fileoutStruct::imageObjective
-#define DEFAULT_SLICE_MIN_X     -5
-#define DEFAULT_SLICE_MIN_Y     0
-#define DEFAULT_SLICE_MIN_Z     -5
+#define DEFAULT_PLANE_MIN_X     -5
+#define DEFAULT_PLANE_MIN_Y     0
+#define DEFAULT_PLANE_MIN_Z     -5
-#define DEFAULT_SLICE_MAX_X     5
-#define DEFAULT_SLICE_MAX_Y     0
-#define DEFAULT_SLICE_MAX_Z     5
+#define DEFAULT_PLANE_MAX_X     5
+#define DEFAULT_PLANE_MAX_Y     0
+#define DEFAULT_PLANE_MAX_Z     5
-#define DEFAULT_SLICE_NORM_X    0
-#define DEFAULT_SLICE_NORM_Y    1
-#define DEFAULT_SLICE_NORM_Z    0
+#define DEFAULT_PLANE_NORM_X    0
+#define DEFAULT_PLANE_NORM_Y    1
+#define DEFAULT_PLANE_NORM_Z    0
+
+#define DEFAULT_PLANE_SIZE		40
+#define DEFAULT_PLANE_POSITION	0
 /*
@@ -64,21 +67,23 @@
 */
-#define DEFAULT_FIELD_ORDER     200
+#define DEFAULT_FIELD_ORDER     10
-#define DEFAULT_SAMPLES         200
+#define DEFAULT_SAMPLES         400
 #define DEFAULT_SLICE_RES		256
+#define DEFAULT_SPHERE_THETA_R  1000
+
 #define DEFAULT_PADDING			1
 #define DEFAULT_SUPERSAMPLE		1
-#define DEFAULT_INTENSITY_FILE	    "testappend"
+#define DEFAULT_INTENSITY_FILE	    "out_i.bmp"
 #define DEFAULT_TRANSMITTANCE_FILE	""
-#define DEFAULT_ABSORBANCE_FILE	    "out_a"
+#define DEFAULT_ABSORBANCE_FILE	    "out_a.bmp"
 #define DEFAULT_NEAR_FILE		    "out_n.bmp"
 #define DEFAULT_FAR_FILE		    "out_f.bmp"
-#define DEFAULT_EXTENDED_SOURCE     "einstein_small.jpg"
+#define DEFAULT_EXTENDED_SOURCE     ""
 #define DEFAULT_FIELD_TYPE		    "magnitude"
 #define DEFAULT_FORMAT			    fileoutStruct::formatImage
 #define DEFAULT_COLORMAP		    "brewer"
@@ -8,14 +8,16 @@
 using namespace std;
 fieldslice::fieldslice(unsigned int x_size, unsigned int y_size)
-{
+{
+    x_hat = y_hat = z_hat = NULL;
+
 	//save the slice resolution
 	R[0] = x_size;
 	R[1] = x_size;
 	scalarField = true;
-	//init_gpu();
+	init_gpu();
 }
@@ -101,5 +103,5 @@ fieldslice::fieldslice()
 fieldslice::~fieldslice()
 {
-	//kill_gpu();
+	kill_gpu();
 }
 #include "fieldslice.h"
 #include "dataTypes.h"
-#include "rts/cuda/error.h"
+#include "rts/cuda/error.h"
+#include "rts/cuda/threads.h"
 __global__ void field_intensity(bsComplex* x, bsComplex* y, bsComplex* z, ptype* I, unsigned int N)
 {
     //compute the index for this thread
-	int i = blockIdx.x * blockDim.x + threadIdx.x;
+	//int i = blockIdx.x * blockDim.x + threadIdx.x;
+	int i = ThreadIndex1D();
+
 	if(i >= N) return;
 	ptype xm = x[i].abs();
@@ -66,7 +69,8 @@ __global__ void resample_intensity(bsComplex* x, bsComplex* y, bsComplex* z, pty
 __global__ void field_real(bsComplex* field_component, ptype* V, unsigned int N)
 {
     //compute the index for this thread
-	int i = blockIdx.x * blockDim.x + threadIdx.x;
+	//int i = blockIdx.x * blockDim.x + threadIdx.x;
+	int i = ThreadIndex1D();
 	if(i >= N) return;
 	V[i] = field_component[i].real();
@@ -75,7 +79,8 @@ __global__ void field_real(bsComplex* field_component, ptype* V, unsigned int N)
 __global__ void field_imaginary(bsComplex* field_component, ptype* V, unsigned int N)
 {
     //compute the index for this thread
-	int i = blockIdx.x * blockDim.x + threadIdx.x;
+	//int i = blockIdx.x * blockDim.x + threadIdx.x;
+	int i = ThreadIndex1D();
 	if(i >= N) return;
 	V[i] = field_component[i].imag();
@@ -84,7 +89,8 @@ __global__ void field_imaginary(bsComplex* field_component, ptype* V, unsigned i
 __global__ void field_sqrt(ptype* input, ptype* output, unsigned int N)
 {
 	//compute the index for this thread
-	int i = blockIdx.x * blockDim.x + threadIdx.x;
+	//int i = blockIdx.x * blockDim.x + threadIdx.x;
+	int i = ThreadIndex1D();
 	if(i >= N) return;
 	output[i] = sqrt(input[i]);
@@ -115,7 +121,8 @@ scalarslice fieldslice::Mag()
 	//compute the total number of values in the slice
 	unsigned int N = R[0] * R[1];
-	int gridDim = (N+BLOCK-1)/BLOCK;
+	//int gridDim = (N+BLOCK-1)/BLOCK;
+	dim3 gridDim = GenGrid1D(N, BLOCK);
 	field_intensity<<<gridDim, BLOCK>>>(x_hat, y_hat, z_hat, result->S, N);
 	field_sqrt<<<gridDim, BLOCK>>>(result->S, result->S, N);
@@ -132,7 +139,8 @@ scalarslice fieldslice::Real()
 	//compute the total number of values in the slice
 	unsigned int N = R[0] * R[1];
-	int gridDim = (N+BLOCK-1)/BLOCK;
+	//int gridDim = (N+BLOCK-1)/BLOCK;
+	dim3 gridDim = GenGrid1D(N, BLOCK);
 	field_real<<<gridDim, BLOCK>>>(x_hat, result->S, N);
@@ -148,7 +156,8 @@ scalarslice fieldslice::Imag()
 	//compute the total number of values in the slice
 	unsigned int N = R[0] * R[1];
-	int gridDim = (N+BLOCK-1)/BLOCK;
+	//int gridDim = (N+BLOCK-1)/BLOCK;
+	dim3 gridDim = GenGrid1D(N, BLOCK);
 	field_imaginary<<<gridDim, BLOCK>>>(x_hat, result->S, N);
@@ -192,7 +201,6 @@ void fieldslice::ScaleField(ptype v)
     //compute the total number of values in the slice
 	unsigned int N = R[0] * R[1];
-	//cout<<"Size of mag field: "<<N<<endl;
 	int gridDim = (N+BLOCK-1)/BLOCK;
 	field_scale<<<gridDim, BLOCK>>>(x_hat, y_hat, z_hat, N, v);
@@ -200,19 +208,23 @@ void fieldslice::ScaleField(ptype v)
 }
 void fieldslice::init_gpu()
-{
+{
+	//if the field has no size, return
+	if(R[0] == 0 || R[1] == 0)
+		return;
+
+    //free any previous memory allocations
+    if(x_hat)
+        HANDLE_ERROR(cudaFree(x_hat));
+    if(y_hat)
+        HANDLE_ERROR(cudaFree(y_hat));
+    if(z_hat)
+        HANDLE_ERROR(cudaFree(z_hat));
+
     //allocate space on the GPU for the field slice
 	HANDLE_ERROR(cudaMalloc((void**)&x_hat, R[0] * R[1] * sizeof(bsComplex)));
-	//HANDLE_ERROR(cudaMemset(x_hat, 0, R[0] * R[1] * sizeof(bsComplex)));
-	//if the field is scalar, y_hat and z_hat are unused
-	if(scalarField)
-	{
-		y_hat = NULL;
-		z_hat = NULL;
-
-	}
-	else
+	if(!scalarField)
 	{
 		HANDLE_ERROR(cudaMalloc((void**)&y_hat, R[0] * R[1] * sizeof(bsComplex)));
 		//HANDLE_ERROR(cudaMemset(y_hat, 0, R[0] * R[1] * sizeof(bsComplex)));
@@ -233,6 +245,8 @@ void fieldslice::kill_gpu()
     if(z_hat != NULL)
         HANDLE_ERROR(cudaFree(z_hat));
+	x_hat = y_hat = z_hat = NULL;
+
 }
 void fieldslice::clear_gpu()
@@ -275,7 +289,7 @@ fieldslice fieldslice::crop(int u, int v, int su, int sv)
 	result.scalarField = scalarField;
 	//allocate space for the new field
-	result.init_gpu();
+	//result.init_gpu();
 	//create one thread for each pixel of the field slice
 	dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
@@ -291,3 +305,57 @@ fieldslice fieldslice::crop(int u, int v, int su, int sv)
 	return result;
 }
+
+fieldslice::fieldslice(const fieldslice& rhs)
+{
+	R[0] = rhs.R[0];
+	R[1] = rhs.R[1];
+	scalarField = rhs.scalarField;
+
+	x_hat = y_hat = z_hat = NULL;
+
+	unsigned int bytes = sizeof(bsComplex) * R[0] * R[1];
+	if(rhs.x_hat != NULL)
+	{
+		HANDLE_ERROR(cudaMalloc( (void**)&x_hat, bytes));
+		HANDLE_ERROR(cudaMemcpy( x_hat, rhs.x_hat, bytes, cudaMemcpyDeviceToDevice));
+	}
+	if(rhs.y_hat != NULL)
+	{
+		HANDLE_ERROR(cudaMalloc( (void**)&y_hat, bytes));
+		HANDLE_ERROR(cudaMemcpy( y_hat, rhs.y_hat, bytes, cudaMemcpyDeviceToDevice));
+	}
+	if(rhs.z_hat != NULL)
+	{
+		HANDLE_ERROR(cudaMalloc( (void**)&z_hat, bytes));
+		HANDLE_ERROR(cudaMemcpy( z_hat, rhs.z_hat, bytes, cudaMemcpyDeviceToDevice));
+	}
+
+}
+
+fieldslice& fieldslice::operator=(const fieldslice& rhs)
+{
+	//make sure this isn't a self-allocation
+	if(this != &rhs)
+	{
+		//make a shallow copy
+		R[0] = rhs.R[0];
+		R[1] = rhs.R[1];
+		scalarField = rhs.scalarField;
+
+		//initialize to new parameters
+		init_gpu();
+
+		//make a deep copy
+		unsigned int bytes = sizeof(bsComplex) * R[0] * R[1];
+		if(x_hat != NULL)
+			HANDLE_ERROR(cudaMemcpy(x_hat, rhs.x_hat, bytes, cudaMemcpyDeviceToDevice));
+		if(y_hat != NULL)
+			HANDLE_ERROR(cudaMemcpy(y_hat, rhs.y_hat, bytes, cudaMemcpyDeviceToDevice));
+		if(z_hat != NULL)
+			HANDLE_ERROR(cudaMemcpy(z_hat, rhs.z_hat, bytes, cudaMemcpyDeviceToDevice));
+	}
+
+	return *this;
+
+}
@@ -31,6 +31,9 @@ struct fieldslice
 	~fieldslice();
+	//copy constructor
+	fieldslice(const fieldslice& rhs);
+
 	//void setPos(bsPoint pMin, bsPoint pMax, bsVector N);
 	scalarslice Mag();
@@ -47,6 +50,7 @@ struct fieldslice
 	//crop a region from the field
 	fieldslice crop(int u, int v, int su, int sv);
+	fieldslice& operator=(const fieldslice& rhs);
 	void init_gpu();
 	void kill_gpu();
@@ -186,11 +186,21 @@ void fileoutStruct::Save(microscopeStruct* scope)
 	//save images of the fields in the microscope
 	//if the user specifies an extended source
-	if(scope->focalPoints.size() > 1)
+	if(scope->focalPoints.size() > 0)
 	{
 		//simulate the extended source and output the detector image
 		scope->SimulateExtendedSource();
+		//saveNearField(&scope->nf);
+		saveFarField(scope);
+
+		//save the detector images
+		saveDetector(scope);
+
+		//simulate scattering for the last point (so that you have a near field image)
+		scope->SimulateScattering();
+        saveNearField(&scope->nf);
+
 	}
 	else
 	{
@@ -203,12 +213,15 @@ void fileoutStruct::Save(microscopeStruct* scope)
 		//run the far-field simulation
 		scope->SimulateImaging();
+		//saveNearField(&scope->nf);
 		saveFarField(scope);
+		//save the detector images
+		saveDetector(scope);
+
 	}
-	//save the detector images
-    saveDetector(scope);
+	
 }
@@ -5,7 +5,7 @@
 //#include "defaults.h"
 #include "dataTypes.h"
-#include "colormap.h"
+#include "rts/graphics/colormap.h"
 #include "fieldslice.h"
 #include "nearfield.h"
 #include "microscope.h"
@@ -34,7 +34,7 @@ struct fileoutStruct{
 	//image_source source;
 	//color map info
-	rts::colormap::colormapType colormap;
+	rts::colormapType colormap;
 	ptype colorMax;
 	void Save(microscopeStruct* scope);
@@ -24,6 +24,7 @@ microscopeStruct* SCOPE;
 #include "warnings.h"
 fileoutStruct gFileOut;
+bool verbose = false;
 using namespace std;
 int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
@@ -31,32 +32,19 @@ int cbessjyva(double v,complex&lt;double&gt; z,double &amp;vm,complex&lt;double&gt;*cjv,
 int main(int argc, char *argv[])
 {
-	//test Envi loading and saving
-	//EnviFile envi("testenvi", "w");
-
-	//float* data = (float*)malloc(sizeof(float) * 100 * 100);
-	//envi.addBand(data, 100, 100, 100);
-
-	//envi.close();
-
-	//return 0;
 	SCOPE = new microscopeStruct();
-    cout<<SCOPE->nf.Uf.R[0]<<endl;
-
     LoadParameters(argc, argv);
-    //TestSimulation(NF, SCOPE, &gFileOut);
-
 	//initialize GPU memory for fields
 	SCOPE->init();
-    OutputOptions();
-
 	gFileOut.Save(SCOPE);
-	//NF->destroy();
+	if(verbose)
+        OutputOptions();
+
 	SCOPE->destroy();
@@ -4,7 +4,7 @@
 #include "rts/tools/progressbar.h"
 #include "rts/cuda/timer.h"
 #include "dataTypes.h"
-#include "colormap.h"
+#include "rts/graphics/colormap.h"
 #include <QImage>
@@ -112,8 +112,8 @@ void microscopeStruct::getFarField()
     //Compute the Far Field image of the focal plane
     //clear the memory from previous detector fields
-    Ud.kill_gpu();
-    Ufd.kill_gpu();
+    //Ud.kill_gpu();
+    //Ufd.kill_gpu();
 	//first crop the filtered near-field image of the source and scattered fields
 	Ud = nf.U.crop(padding * Ud.R[0], padding * Ud.R[1], Ud.R[0], Ud.R[1]);
@@ -261,9 +261,14 @@ void microscopeStruct::SimulateExtendedSource()
 		t += gpuStopTimer();
 		rtsProgressBar((double)(i+1)/(double)npts * 100);
+		//unsigned char c;
+		//cin>>c;
 	}
-	cout<<endl;
-	cout<<"Time per source: "<<t/npts<<"ms"<<endl;
+	if(verbose)
+	{
+        cout<<endl;
+        cout<<"Time per source: "<<t/npts<<"ms"<<endl;
+    }
 }
@@ -304,3 +309,15 @@ void microscopeStruct::LoadExtendedSource(std::string filename)
 			}
         }
 }
+
+std::string microscopeStruct::toStr()
+{
+	stringstream ss;
+	ss<<nf.toStr();
+
+	ss<<"----------Optics--------------"<<endl<<endl;
+	ss<<"Objective NA: "<<objective[0]<<" to "<<objective[1]<<endl;
+	return ss.str();
+
+
+}
@@ -63,6 +63,8 @@ struct microscopeStruct
 	scalarslice getTransmittance();
 	scalarslice getIntensity();
+	string toStr();
+
 };
@@ -35,18 +35,12 @@ void mcSampleNA(bsVector* samples, int N, bsVector k, ptype NAin, ptype NAout)
     ptype inPhi = asin(NAin);
     ptype outPhi = asin(NAout);
-    //cout<<"inPhi: "<<inPhi<<endl;
-    //cout<<"outPhi: "<<outPhi<<endl;
-
     //calculate the z-values associated with these angles
     ptype inZ = cos(inPhi);
     ptype outZ = cos(outPhi);
     ptype rangeZ = inZ - outZ;
-    //cout<<"inZ: "<<inZ<<endl;
-    //cout<<"outZ: "<<outZ<<endl;
-
     //draw a distribution of random phi, z values
     ptype z, phi, theta;
     for(int i=0; i<N; i++)
@@ -58,7 +52,6 @@ void mcSampleNA(bsVector* samples, int N, bsVector k, ptype NAin, ptype NAout)
         phi = acos(z);
         //compute and store cartesian coordinates
-        //bsVector spherical(1, theta + kSph[1], phi + kSph[2]);
 		bsVector spherical(1, theta, phi);
 		bsVector cart = spherical.sph2cart();
         samples[i] = rotation * cart;
 #include "nearfield.h"
+#include <time.h>
+#include <math.h>
+
+#ifdef _WIN32
+#define isnan(x) _isnan(x)
+#define isinf(x) (!_finite(x))
+#endif
+
+int bessjyv_sph(int v, double z, double &vm, double* cjv,
+    double* cyv, double* cjvp, double* cyvp);
 nearfieldStruct::nearfieldStruct()
 {
     scalarSim = true;
 	planeWave = false;
+	lut_us = true;
+	lut_uf = false;
 	nWaves = 0;
 }
@@ -46,6 +58,8 @@ std::string nearfieldStruct::toStr()
 	ss<<"Condenser NA: "<<condenser[0]<<" to "<<condenser[1]<<std::endl;
 	ss<<"Focal Point: "<<focus[0]<<", "<<focus[1]<<", "<<focus[2]<<std::endl;
 	ss<<"Field Slice: "<<std::endl;
+	if(lut_us)
+        ss<<"LUT Parameters --- min: "<<d_min<<"   max: "<<d_max<<std::endl;
 	ss<<pos<<std::endl;
 	ss<<std::endl<<"---------Materials-----------"<<std::endl;
@@ -61,6 +75,10 @@ std::string nearfieldStruct::toStr()
 	for(unsigned int s=0; s<sVector.size(); s++)
 		ss<<sVector[s].toStr()<<std::endl;
+    ss<<"---------Timings-------------"<<std::endl;
+    ss<<"Uf = "<<t_Uf<<"ms"<<std::endl;
+    ss<<"Us = "<<t_Us<<"ms"<<std::endl;
+
 	return ss.str();
 }
@@ -70,7 +88,8 @@ void nearfieldStruct::calcWaves()
     inWaves.resize(nWaves);
     //re-seed the random number generator
-    //srand(seed);
+    //srand(time(NULL));
+	srand(NULL);
     //calculate the monte-carlo samples
     mcSampleNA(&inWaves[0], nWaves, k, condenser[0], condenser[1]);
@@ -84,6 +103,8 @@ void nearfieldStruct::calcSpheres()
     //calculate all of the constants necessary to evaluate the scattered field
 	//estimate the order required to represent the scattered field for each sphere
+
+
 	//for each sphere
 	for(int i=0; i<sVector.size(); i++)
 	{
@@ -91,12 +112,10 @@ void nearfieldStruct::calcSpheres()
 		//calculate the required order
 		sVector[i].calcNl(lambda);
-		//std::cout<<sVector[i].Nl<<std::endl;
 		//set the refractive index for the sphere
 		int imat = sVector[i].iMaterial;
         rts::rtsComplex<ptype> n = mVector[imat](lambda);
-		//std::cout<<"Sphere refractive index: "<<n<<std::endl;
 		//calculate the scattering coefficients
 		sVector[i].calcCoeff(lambda, n);
@@ -104,18 +123,109 @@ void nearfieldStruct::calcSpheres()
 		//save the refractive index
 		sVector[i].n = n;
+		//if the LUT is used, calculate Usp(theta, r)
+		if(lut_us)
+		{
+			sVector[i].calcUp(lambda, n, pos, max(U.R[0], U.R[1]));
+        }
+
+
 	}
 }
+void nearfieldStruct::calcUs()
+{
+
+
+    if(lut_us)
+        scalarUpLut();
+    else
+        scalarUs();
+}
+
+void nearfieldStruct::calcUf()
+{
+	if(lut_uf)
+		scalarUfLut();
+	else
+		scalarUf();
+}
+
 void nearfieldStruct::Simulate()
 {
+    //initialize timings
+    t_Uf = 0;
+    t_Us = 0;
+
 	//compute a set of plane waves for Monte-Carlo simulation
 	calcWaves();
     //the near field has to be simulated no matter what the output rtsPoint is
-    scalarUf();
+    calcUf();
     calcSpheres();
-    scalarUs();
+    calcUs();
     sumUf();
+
+	//U.Mag().toImage("testU.bmp");
+}
+
+void nearfieldStruct::calcBesselLut(ptype* j, ptype d_min, ptype d_max, int dR)
+{
+    /*Compute the look-up-table for spherical bessel functions used for the incident field
+        j    =   (Nl + 1) x aR array of values
+        aR      =   resolution of j
+    */
+
+	//compute the wavenumber
+	ptype k = 2 * PI / lambda;
+	unsigned int Nl = m;
+
+    //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_kd = (double*)malloc(bytes);
+    double* cyv_kd = (double*)malloc(bytes);
+    double* cjvp_kd = (double*)malloc(bytes);
+    double* cyvp_kd = (double*)malloc(bytes);
+
+    //compute the bessel functions using the CPU-based algorithm
+    double vm;
+
+    //for each sample along r
+    ptype dr = (d_max - d_min) / (dR - 1);
+    ptype d;
+    ptype jv;
+    for(int id = 0; id < dR; id++)
+    {
+        d = id * dr + d_min;
+        double kd = k*d;
+        bessjyv_sph(Nl, kd, vm, cjv_kd, cyv_kd, cjvp_kd, cyvp_kd);
+
+        //copy the double data to the bsComplex array
+        for(int l=0; l<=Nl; l++)
+		{
+            jv = cjv_kd[l];
+			if(isnan(jv) || isinf(jv))
+			{
+                if(kd == 0 && l == 0)
+                    jv = 1;
+                else
+                    jv = 0;
+            }
+            j[id * (Nl+1) + l] = jv;
+		}
+    }
+
+	/*ofstream outfile("uf_besselout.txt");
+    for(int ir = 0; ir < dR; ir++)
+    {
+		outfile<<ir*dr + d_min<<endl;
+        for(int l = 0; l<=Nl; l++)
+        {
+            outfile<<j[ir * (Nl+1) + l]<<" --";
+        }
+        outfile<<endl;
+    }
+	outfile.close();*/
+
 }
@@ -31,6 +31,8 @@ struct nearfieldStruct
 	//slices for the focused field
 	fieldslice Uf;
+	ptype d_min, d_max;
+
 	//	and total field: Uf + sum(Us)
 	fieldslice U;
@@ -43,6 +45,14 @@ struct nearfieldStruct
 	//flag for a plane wave
 	bool planeWave;
+	//flag for using a LUT
+	bool lut_uf;
+	bool lut_us;
+
+	//timings
+	float t_Uf;
+	float t_Us;
+
 	//---------Scatterers------------
@@ -78,10 +88,17 @@ struct nearfieldStruct
 	void setPos(bsPoint pMin, bsPoint pMax, bsVector normal);
 	//this function re-computes the focused field
+	void calcUf();
 	void scalarUf();
+	void scalarUfLut();
+
+	void calcBesselLut(ptype* j, ptype d_min, ptype d_max, int dR);
 	//compute the field scattered by all of the materials
+	void calcUs();
 	void scalarUs();
+	void scalarUpLut();
+
 	//add the incident field to the sum of scattered fields
 	void sumUf();
@@ -5,7 +5,7 @@
 #include "rts/cuda/error.h"
 #include "rts/cuda/timer.h"
-
+//Incident field for a single plane wave
 __global__ void gpuScalarUfp(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR)
 {
 	/*Compute the scalar focused field using Debye focusing
@@ -41,7 +41,8 @@ __global__ void gpuScalarUfp(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, p
 	Uf[i] = exp(d) * A;
 }
-
+
+//Incident field for a focused point source
 __global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR, ptype cosAlpha, ptype cosBeta, int nl, ptype j_conv = 1.4)
 {
 	/*Compute the scalar focused field using Debye focusing
@@ -151,7 +152,6 @@ __global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, pt
 		}
 		sumUf += il * jl * Pl * (Palpha[1] - Palpha[2] - Pbeta[1] + Pbeta[2]);
-		//sumUf += il * Pl * (Palpha[1] - Palpha[2] - Pbeta[1] + Pbeta[2]);
 		il *= im;
 	}
@@ -162,21 +162,12 @@ __global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, pt
 void nearfieldStruct::scalarUf()
 {
-	//Compute the incident field via a scalar simulation
-	//This method uses Debye focusing to approximate the field analytically
-
-	//time the calculation of the focused field
-	//gpuStartTimer();
-
-	//set the field slice to a scalar field
-	//Uf.scalarField = true;
-
-	//initialize the GPU arrays
-	//Uf.init_gpu();
+
+    gpuStartTimer();
 	//create one thread for each pixel of the field slice
 	dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
-	dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
+	dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
 	//if we are computing a plane wave, call the gpuScalarUfp function
 	if(planeWave)
@@ -191,10 +182,7 @@ void nearfieldStruct::scalarUf()
 		ptype cosBeta = cos(asin(condenser[1]));
 		//compute the scalar Uf field (this will be in the x_hat channel of Uf)
 		gpuScalarUf<<<dimGrid, dimBlock>>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1], cosAlpha, cosBeta, m);
-	}
-
-	//float t = gpuStopTimer();
-	//std::cout<<"Scalar Uf Time: "<<t<<"ms"<<std::endl;
-	//std::cout<<focus<<std::endl;
-
+	}
+
+	t_Uf = gpuStopTimer();
 }
+#include "nearfield.h"
+
+#include "rts/math/legendre.h"
+#include "rts/cuda/error.h"
+#include "rts/cuda/timer.h"
+
+texture<float, cudaTextureType2D> texJ;
+
+__global__ void gpuScalarUfp(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR);
+
+__global__ void gpuScalarUfLut(bsComplex* Uf, bsRect ABCD, int uR, int vR, bsPoint f, bsVector k, ptype A, ptype cosAlpha, ptype cosBeta, int nl, ptype dmin, ptype dmax, int dR)
+{
+    /*This function computes the focused field for a 2D slice
+
+    Uf      =   destination field slice
+    ABCD    =   plane representing the field slice in world space
+    uR, vR  =   resolution of the Uf field
+    f       =   focal point of the condenser
+    k       =   direction of the incident light
+    A       =   amplitude of the incident field
+    cosAlpha=   cosine of the solid angle subtended by the condenser obscuration
+    cosBeta =   cosine of the solid angle subtended by the condenser aperature
+    nl      =   number of orders used to compute the field
+    dR      =   number of Bessel function values in the look-up texture
+
+    */
+
+    //get the current coordinate in the plane slice
+	int iu = blockIdx.x * blockDim.x + threadIdx.x;
+	int iv = blockIdx.y * blockDim.y + threadIdx.y;
+
+	//make sure that the thread indices are in-bounds
+	if(iu >= uR || iv >= vR) return;
+
+	//compute the index (easier access to the scalar field array)
+	int i = iv*uR + iu;
+
+	//compute the parameters for u and v
+	ptype u = (ptype)iu / (uR);
+	ptype v = (ptype)iv / (vR);
+
+	
+
+	//get the rtsPoint in world space and then the r vector
+	bsPoint p = ABCD(u, v);
+	bsVector r = p - f;
+	ptype d = r.len();
+
+	if(d == 0)
+	{
+        Uf[i] = A * 2 * PI * (cosAlpha - cosBeta);
+        return;
+    }
+
+	//get info for the light direction and frequency
+	r = r.norm();
+
+	//compute the imaginary factor i^l
+	bsComplex im = bsComplex(0, 1);
+	bsComplex il = bsComplex(1, 0);
+
+	//Legendre functions are computed dynamically to save memory
+	//initialize the Legendre functions
+
+	ptype P[2];
+	//get the angle between k and r (light direction and position vector)
+	ptype cosTheta;
+	cosTheta = k.dot(r);
+
+	rts::init_legendre<ptype>(cosTheta, P[0], P[1]);
+
+	//initialize legendre functions for the cassegrain angles
+	ptype Palpha[3];
+	rts::init_legendre<ptype>(cosAlpha, Palpha[0], Palpha[1]);
+	Palpha[2] = 1;
+
+	ptype Pbeta[3];
+	rts::init_legendre<ptype>(cosBeta, Pbeta[0], Pbeta[1]);
+	Pbeta[2] = 1;
+
+	//for each order l
+	bsComplex sumUf(0.0, 0.0);
+	ptype jl = 0.0;
+	ptype Pl;
+	ptype di = ( (d - dmin)/(dmax - dmin) ) * (dR - 1);
+	for(int l = 0; l<=nl; l++)
+	{
+        jl = tex2D(texJ, l + 0.5, di + 0.5);
+		if(l==0)
+			Pl = P[0];
+		else if(l==1)
+		{
+			Pl = P[1];
+
+			//adjust the cassegrain Legendre function
+			Palpha[2] = Palpha[0];
+			rts::shift_legendre<ptype>(l+1, cosAlpha, Palpha[0], Palpha[1]);
+			Pbeta[2] = Pbeta[0];
+			rts::shift_legendre<ptype>(l+1, cosBeta, Pbeta[0], Pbeta[1]);
+		}
+		else
+		{
+			rts::shift_legendre<ptype>(l, cosTheta, P[0], P[1]);
+
+			Pl = P[1];
+
+			//adjust the cassegrain outer Legendre function
+			Palpha[2] = Palpha[0];
+			rts::shift_legendre<ptype>(l+1, cosAlpha, Palpha[0], Palpha[1]);
+			Pbeta[2] = Pbeta[0];
+			rts::shift_legendre<ptype>(l+1, cosBeta, Pbeta[0], Pbeta[1]);
+		}
+
+		sumUf += il * jl * Pl * (Palpha[1] - Palpha[2] - Pbeta[1] + Pbeta[2]);
+		//sumUf += jl;
+
+		il *= im;
+	}
+
+	Uf[i] = sumUf * 2 * PI * A;
+	//Uf[i] = u;
+	//return;
+}
+
+void nearfieldStruct::scalarUfLut()
+{
+    gpuStartTimer();
+	
+    //calculate the minimum and maximum points in the focused field
+    d_min = pos.dist(focus);
+    d_max = pos.dist_max(focus);
+
+    //allocate space for the Bessel function
+    int dR = 2 * max(Uf.R[0], Uf.R[1]);
+    ptype* j = NULL;
+	j = (ptype*) malloc(sizeof(ptype) * dR * (m+1));
+
+	//calculate Bessel function LUT
+	calcBesselLut(j, d_min, d_max, dR);
+	
+    //create a CUDA array structure and specify the format description
+	cudaArray* arrayJ;
+    cudaChannelFormatDesc channelDesc =
+        cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
+	
+    //allocate memory
+    HANDLE_ERROR(cudaMallocArray(&arrayJ, &channelDesc, m+1, dR));
+	
+    //specify texture properties
+    texJ.addressMode[0] = cudaAddressModeMirror;
+    texJ.addressMode[1] = cudaAddressModeMirror;
+    texJ.filterMode     = cudaFilterModeLinear;
+    texJ.normalized     = false;
+
+    //bind the texture to the array
+    HANDLE_ERROR(cudaBindTextureToArray(texJ, arrayJ, channelDesc));
+
+    //copy the CPU Bessel LUT to the GPU-based array
+    HANDLE_ERROR( cudaMemcpy2DToArray(arrayJ, 0, 0, j, (m+1)*sizeof(float), (m+1)*sizeof(float), dR, cudaMemcpyHostToDevice));
+
+    //----------------Compute the focused field
+    //create one thread for each pixel of the field slice
+	dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
+	dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
+
+	//if we are computing a plane wave, call the gpuScalarUfp function
+	if(planeWave)
+	{
+		gpuScalarUfp<<<dimGrid, dimBlock>>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1]);
+	}
+	//otherwise compute the condenser info and create a focused field
+	else
+	{
+		//pre-compute the cosine of the obscuration and objective angles
+		ptype cosAlpha = cos(asin(condenser[0]));
+		ptype cosBeta = cos(asin(condenser[1]));
+		//compute the scalar Uf field (this will be in the x_hat channel of Uf)
+		gpuScalarUfLut<<<dimGrid, dimBlock>>>(Uf.x_hat, pos, Uf.R[0], Uf.R[1], focus, k, A, cosAlpha, cosBeta, m, d_min, d_max, dR);
+	}
+
+	
+    //free everything
+	free(j);
+	
+	HANDLE_ERROR(cudaFreeArray(arrayJ));
+
+	t_Uf = gpuStopTimer();
+}
+#include "nearfield.h"
+#include "rts/math/spherical_bessel.h"
+#include "rts/math/legendre.h"
+#include <stdlib.h>
+#include "rts/cuda/error.h"
+#include "rts/cuda/timer.h"
+
+texture<float2, cudaTextureType2D> texUsp;
+texture<float2, cudaTextureType2D> texUip;
+
+__global__ void gpuScalarUpLut(bsComplex* Us, bsVector* k, int nk, ptype kmag, ptype a, ptype dmin, ptype dmax, bsPoint f, bsPoint ps, ptype A, bsRect ABCD, int uR, int vR, int dR, int aR, int thetaR)
+{
+    /*This function uses Monte-Carlo integration to sample a texture-based LUT describing the scattered field
+        produced by a plane wave through a sphere.  The MC sampling is used to approximate a focused field.
+
+        Us  =   final scattered field
+        k   =   list of incoming plane waves (Monte-Carlo samples)
+        nk  =   number of incoming MC samples
+        kmag=   magnitude of the incoming field 2pi/lambda
+        dmin=   minimum distance of the Usp texture
+        dmax=   maximum distance of the Usp texture
+        f   =   position of the focus
+        ps  =   position of the sphere
+        A   =   total amplitude of the incident field arriving at the focal spot
+        ABCD=   rectangle representing the field slice
+        uR  =   resolution of the field slice in the u direction
+        vR  =   resolution of the field slice in the v direction
+        dR  =   resolution of the Usp texture in the d direction
+        thetaR= resolution of the Usp texture in the theta direction
+    */
+
+	//get the current coordinate in the plane slice
+	int iu = blockIdx.x * blockDim.x + threadIdx.x;
+	int iv = blockIdx.y * blockDim.y + threadIdx.y;
+
+	//make sure that the thread indices are in-bounds
+	if(iu >= uR || iv >= vR) return;
+
+	//compute the index (easier access to the scalar field array)
+	int i = iv*uR + iu;