bug fixes related to the development of bimsym

David Mayerich
1 parent ca99f951
Showing 8 changed files with 338 additions and 46 deletions Show diff stats
stim/envi/bsq.h
stim/math/fft.h
stim/optics/scalarbeam.h
stim/optics/scalarfield.h
stim/optics/mie.h → stim/optics/scalarmie.h
stim/optics/scalarwave.h
stim/parser/arguments.h
stim/parser/filename.h
@@ -1021,12 +1021,13 @@ public:
 		//for each band
 		for (unsigned long long z = b0; z < b1; z++)
 		{
+			//std::cout<<z<<std::endl;
 			for (unsigned long long y = 0; y < lines; y++)
 			{
 				file.read((char *)(temp + y * samples), sizeof(T) * samples);
 				file.seekg(jumpl, std::ios::cur);    //go to the next band
-				if(PROGRESS) progress = (double)((z+1) * lines + y + 1) / ((b1 - b0) * lines) * 100;
+				if(PROGRESS) progress = (double)((z - b0 + 1) * lines + y + 1) / ((b1 - b0) * lines) * 100;
 			}
 			out.write(reinterpret_cast<const char*>(temp), L);   //write slice data into target file
 			file.seekg(jumpb, std::ios::cur);
@@ -3,7 +3,7 @@
 namespace stim{
-	template<class T>
+	/*template<class T>
 	void circshift(T *out, const T *in, size_t xdim, size_t ydim, size_t xshift, size_t yshift){
 		size_t i, j, ii, jj;
 		for (i =0; i < xdim; i++) {
@@ -13,16 +13,31 @@ namespace stim{
 				out[ii * ydim + jj] = in[i * ydim + j];
 			}
 		}
+	}*/
+
+	template<typename T>
+	void circshift(T *out, const T *in, int xdim, int ydim, int xshift, int yshift)
+	{
+	 for (int i =0; i < xdim; i++) {
+	   int ii = (i + xshift) % xdim;
+	   if (ii<0) ii = xdim + ii;
+	   for (int j = 0; j < ydim; j++) {
+	     int jj = (j + yshift) % ydim;
+	     if (jj<0) jj = ydim + jj;
+	     //out[ii * ydim + jj] = in[i * ydim + j];
+	     out[jj * xdim + ii] = in[j * xdim + i];
+	   }
+	 }
 	}
 	template<typename T>
 	void cpu_fftshift(T* out, T* in, size_t xdim, size_t ydim){
-		circshift(out, in, xdim, ydim, xdim/2, ydim/2);
+		circshift(out, in, xdim, ydim, std::floor(xdim/2), std::floor(ydim/2));
 	}
 	template<typename T>
 	void cpu_ifftshift(T* out, T* in, size_t xdim, size_t ydim){
-		circshift(out, in, xdim, ydim, xdim/2, ydim/2);
+		circshift(out, in, xdim, ydim, std::ceil(xdim/2), std::ceil(ydim/2));
 	}
@@ -213,13 +213,13 @@ void gpu_scalar_psf_local(stim::complex&lt;T&gt;* E, size_t N, T* r, T* phi, T lambda,
 	T dr = (r_max - r_min) / (Nlut_j-1);												//distance between values in the LUT
 	T jl;
 	for(size_t ri = 0; ri < Nlut_j; ri++){													//for each value in the LUT
-		for(size_t l = 0; l <= Nl; l++){													//for each order
+		for(unsigned l = 0; l <= (unsigned)Nl; l++){													//for each order
 			jl = boost::math::sph_bessel<T>(l, k*(r_min + ri * dr));					//use boost to calculate the spherical bessel function
 			j_lut[ri * (Nl + 1) + l] = jl;												//store the bessel function result
 		}
 	}
-	stim::cpu2image<T>(j_lut, "j_lut.bmp", Nl+1, Nlut_j, stim::cmBrewer);
+	//stim::cpu2image<T>(j_lut, "j_lut.bmp", Nl+1, Nlut_j, stim::cmBrewer);
 	//Allocate device memory and copy everything to the GPU
 	T* gpu_C;
@@ -316,6 +316,7 @@ void gpu_scalar_psf_cart(stim::complex&lt;T&gt;* E, size_t N, T* x, T* y, T* z, T lamb
 	stim::quaternion<T> q;												//create a quaternion
 	q.CreateRotation(d, stim::vec3<T>(0, 0, 1));						//create a mapping from the propagation direction to the PSF space
+	stim::matrix<T, 3> rot = q.toMatrix3();
 	int threads = stim::maxThreadsPerBlock();							//get the maximum number of threads per block for the CUDA device
 	dim3 blocks( (unsigned)(N / threads + 1));							//calculate the optimal number of blocks
 	cuda_cart2psf<T> <<< blocks, threads >>> (gpu_r, gpu_phi, N, x, y, z, f, q);	//call the CUDA kernel to move the cartesian coordinates to PSF space
@@ -442,15 +443,19 @@ public:
 		return samples;
 	}
+	void eval(stim::scalarfield<T>& E, T* X, T* Y, T* Z, int order = 500){
+		cpu_scalar_psf_cart<T>(E.ptr(), E.size(), X, Y, Z, lambda, A, f, d, NA[0], NA[1], order, E.spacing());
+	}
+
 	/// Evaluate the beam to a scalar field using Debye focusing
-	void eval(stim::scalarfield<T>& E, size_t order = 500){
+	void eval(stim::scalarfield<T>& E, int order = 500){
 		size_t array_size = E.grid_bytes();
 		T* X = (T*) malloc( array_size );			//allocate space for the coordinate meshes
 		T* Y = (T*) malloc( array_size );
 		T* Z = (T*) malloc( array_size );
 		E.meshgrid(X, Y, Z, stim::CPUmem);			//calculate the coordinate meshes
-		cpu_scalar_psf_cart<T>(E.ptr(), E.size(), X, Y, Z, lambda, A, f, d, NA[0], NA[1], order, E.spacing());
+		eval(E, X, Y, Z, order);
 		free(X);									//free the coordinate meshes
 		free(Y);
 #ifndef STIM_SCALARFIELD_H
 #define STIM_SCALARFIELD_H
+
 #include "../math/rect.h"
 #include "../math/complex.h"
+#include "../math/fft.h"
 namespace stim{
+	/// Perform a k-space transform of a scalar field (FFT). The given field has a width of x and the calculated momentum space has a
+	///		width of kx (in radians).
+	/// @param K is a pointer to the output array of all plane waves in the field
+	/// @param kx is the width of the frame in momentum space
+	/// @param ky is the height of the frame in momentum space
+	/// @param E is the field to be transformed
+	/// @param x is the width of the field in the spatial domain
+	/// @param y is the height of the field in the spatial domain
+	/// @param nx is the number of pixels representing the field in the x (and kx) direction
+	/// @param ny is the number of pixels representing the field in the y (and ky) direction
+	template<typename T>
+	void cpu_scalar_to_kspace(stim::complex<T>* K, T& kx, T& ky, stim::complex<T>* E, T x, T y, size_t nx, size_t ny){
+
+		kx = stim::TAU * nx / x;			//calculate the width of the momentum space
+		ky = stim::TAU * ny / y;
+
+		stim::complex<T>* dev_FFT;
+		HANDLE_ERROR( cudaMalloc(&dev_FFT, sizeof(stim::complex<T>) * nx * ny) );		//allocate space on the CUDA device for the output array
+
+		stim::complex<T>* dev_E;
+		HANDLE_ERROR( cudaMalloc(&dev_E, sizeof(stim::complex<T>) * nx * ny) );		//allocate space for the field
+		HANDLE_ERROR( cudaMemcpy(dev_E, E, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyHostToDevice) );	//copy the field to GPU memory
+
+		cufftResult result;
+		cufftHandle plan;
+		result = cufftPlan2d(&plan, nx, ny, CUFFT_C2C);
+		if(result != CUFFT_SUCCESS){
+			std::cout<<"Error creating cuFFT plan."<<std::endl;
+			exit(1);
+		}
+
+		result = cufftExecC2C(plan, (cufftComplex*)dev_E, (cufftComplex*)dev_FFT, CUFFT_FORWARD);
+		if(result != CUFFT_SUCCESS){
+			std::cout<<"Error using cuFFT to perform a forward Fourier transform of the field."<<std::endl;
+			exit(1);
+		}
+
+		stim::complex<T>* fft = (stim::complex<T>*) malloc(sizeof(stim::complex<T>) * nx * ny);
+		HANDLE_ERROR( cudaMemcpy(fft, dev_FFT, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyDeviceToHost) );
+
+		stim::cpu_fftshift(K, fft, nx, ny);
+		//memcpy(K, fft, sizeof(stim::complex<T>) * nx * ny);
+	}
+
+	template<typename T>
+	void cpu_scalar_from_kspace(stim::complex<T>* E, T& x, T& y, stim::complex<T>* K, T kx, T ky, size_t nx, size_t ny){
+
+		x = stim::TAU * nx / kx;			//calculate the width of the momentum space
+		y = stim::TAU * ny / ky;
+		
+		stim::complex<T>* fft = (stim::complex<T>*) malloc(sizeof(stim::complex<T>) * nx * ny);
+		stim::cpu_ifftshift(fft, K, nx, ny);
+		//memcpy(fft, K, sizeof(stim::complex<T>) * nx * ny);
+
+		stim::complex<T>* dev_FFT;
+		HANDLE_ERROR( cudaMalloc(&dev_FFT, sizeof(stim::complex<T>) * nx * ny) );		//allocate space on the CUDA device for the output array
+		HANDLE_ERROR( cudaMemcpy(dev_FFT, fft, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyHostToDevice) );	//copy the field to GPU memory
+
+		stim::complex<T>* dev_E;
+		HANDLE_ERROR( cudaMalloc(&dev_E, sizeof(stim::complex<T>) * nx * ny) );		//allocate space for the field
+
+		cufftResult result;
+		cufftHandle plan;
+		result = cufftPlan2d(&plan, nx, ny, CUFFT_C2C);
+		if(result != CUFFT_SUCCESS){
+			std::cout<<"Error creating cuFFT plan."<<std::endl;
+			exit(1);
+		}
+
+		result = cufftExecC2C(plan, (cufftComplex*)dev_FFT, (cufftComplex*)dev_E, CUFFT_INVERSE);
+		if(result != CUFFT_SUCCESS){
+			std::cout<<"Error using cuFFT to perform a forward Fourier transform of the field."<<std::endl;
+			exit(1);
+		}
+
+		HANDLE_ERROR( cudaMemcpy(E, dev_E, sizeof(stim::complex<T>) * nx * ny, cudaMemcpyDeviceToHost) );
+
+		
+	}
+
+	/// Propagate a field slice along its orthogonal direction by a given distance z
+	/// @param Enew is the resulting propogated field
+	/// @param E is the field to be propogated
+	/// @param sx is the size of the field in the lateral x direction
+	/// @param sy is the size of the field in the lateral y direction
+	/// @param z is the distance to be propagated
+	/// @param k is the wavenumber 2*pi/lambda
+	/// @param nx is the number of samples in the field along the lateral x direction
+	/// @param ny is the number of samples in the field along the lateral y direction
+	template<typename T>
+	void cpu_scalar_propagate(stim::complex<T>* Enew, stim::complex<T>* E, T sx, T sy, T z, T k, size_t nx, size_t ny){
+		
+		stim::complex<T>* K = (stim::complex<T>*) malloc( sizeof(stim::complex<T>) * nx * ny );
+
+		T Kx, Ky;											//width and height in k space
+		cpu_scalar_to_kspace(K, Kx, Ky, E ,sx, sy, nx, ny);
+
+		T* mag = (T*) malloc( sizeof(T) * nx * ny );
+		stim::abs(mag, K, nx * ny);
+		stim::cpu2image<float>(mag, "kspace_pre_shift.bmp", nx, ny, stim::cmBrewer);
+		
+		size_t kxi, kyi;
+		size_t i;
+		T kx, kx_sq, ky, ky_sq, k_sq;
+		T kz;
+		stim::complex<T> shift;
+		T min_kx = -Kx / 2;
+		T dkx = Kx / (nx);
+
+		T min_ky = -Ky / 2;
+		T dky = Ky / (ny);
+
+		for(kyi = 0; kyi < ny; kyi++){						//for each plane wave in the ky direction
+			for(kxi = 0; kxi < nx; kxi++){					//for each plane wave in the ky direction
+				i = kyi * nx + kxi;
+
+				kx = min_kx + kxi * dkx;					//calculate the position of the current plane wave
+				ky = min_ky + kyi * dky;
+
+				kx_sq = kx * kx;
+				ky_sq = ky * ky;
+				k_sq = k*k;
+				
+				if(kx_sq + ky_sq < k_sq){
+					kz = sqrt(k_sq - kx_sq - ky_sq);			//estimate kz using the Fresnel approximation				
+					shift = -exp(stim::complex<T>(0, kz * z));
+					K[i] *= shift;
+					K[i] /= (nx*ny);							//normalize the DFT
+				}
+				else{
+					K[i] = 0;
+				}
+			}
+		}
+		
+		stim::abs(mag, K, nx * ny);
+		stim::cpu2image<float>(mag, "kspace_post_shift.bmp", nx, ny, stim::cmBrewer);
+		
+		cpu_scalar_from_kspace(Enew, sx, sy, K, Kx, Ky, nx, ny);
+	}
+
+	/// Apply a lowpass filter to a field slice
+	/// @param Enew is the resulting propogated field
+	/// @param E is the field to be propogated
+	/// @param sx is the size of the field in the lateral x direction
+	/// @param sy is the size of the field in the lateral y direction
+	/// @param highest is the highest spatial frequency that can pass through the filter
+	/// @param nx is the number of samples in the field along the lateral x direction
+	/// @param ny is the number of samples in the field along the lateral y direction
+	template<typename T>
+	void cpu_scalar_lowpass(stim::complex<T>* Enew, stim::complex<T>* E, T sx, T sy, T highest, size_t nx, size_t ny){
+		
+		stim::complex<T>* K = (stim::complex<T>*) malloc( sizeof(stim::complex<T>) * nx * ny );
+
+		T Kx, Ky;											//width and height in k space
+		cpu_scalar_to_kspace(K, Kx, Ky, E ,sx, sy, nx, ny);
+
+		T* mag = (T*) malloc( sizeof(T) * nx * ny );
+		stim::abs(mag, K, nx * ny);
+		stim::cpu2image<float>(mag, "kspace_pre_lowpass.bmp", nx, ny, stim::cmBrewer);
+		
+		size_t kxi, kyi;
+		size_t i;
+		T kx, kx_sq, ky, ky_sq, k_sq;
+		T kz;
+		stim::complex<T> shift;
+		T min_kx = -Kx / 2;
+		T dkx = Kx / (nx);
+
+		T min_ky = -Ky / 2;
+		T dky = Ky / (ny);
+
+		T highest_sq = highest * highest;
+
+		for(kyi = 0; kyi < ny; kyi++){						//for each plane wave in the ky direction
+			for(kxi = 0; kxi < nx; kxi++){					//for each plane wave in the ky direction
+				i = kyi * nx + kxi;
+
+				kx = min_kx + kxi * dkx;					//calculate the position of the current plane wave
+				ky = min_ky + kyi * dky;
+
+				kx_sq = kx * kx;
+				ky_sq = ky * ky;
+				
+				if(kx_sq + ky_sq > highest_sq){
+					K[i] = 0;
+				}
+				else
+					K[i] /= nx * ny;						//normalize the DFT
+			}
+		}
+		
+		stim::abs(mag, K, nx * ny);
+		stim::cpu2image<float>(mag, "kspace_post_lowpass.bmp", nx, ny, stim::cmBrewer);
+		
+		cpu_scalar_from_kspace(Enew, sx, sy, K, Kx, Ky, nx, ny);
+	}
+
 	enum locationType {CPUmem, GPUmem};
 	/// Class represents a scalar optical field.
@@ -22,23 +222,20 @@ protected:
 	size_t R[2];
 	locationType loc;
-	
-
-public:
-
-	CUDA_CALLABLE scalarfield(size_t X, size_t Y, T size = 1, T z_pos = 0) : rect<T>::rect(size, z_pos){
-		R[0] = X;											//set the field resolution
-		R[1] = Y;
-
-		E = (stim::complex<T>*) malloc(sizeof(stim::complex<T>) * R[0] * R[1]);		//allocate in CPU memory
-		loc = CPUmem;
+	/// Convert the field to a k-space representation (do an FFT)
+	void to_kspace(T& kx, T& ky){
+		cpu_scalar_to_kspace(E, kx, ky, E, X.len(), Y.len(), R[0], R[1]);
 	}
-	CUDA_CALLABLE ~scalarfield(){
-		if(loc == CPUmem) free(E);
-		else cudaFree(E);
+	void from_kspace(){
+		kx = stim::TAU * R[0] / X.len();			//calculate the width of the momentum space
+		ky = stim::TAU * R[1] / Y.len();
+		T x, y;
+		cpu_scalar_from_kspace(E, x, y, E, kx, ky, R[0], R[1]);
 	}
+public:
+
 	/// Returns the number of values in the field
 	CUDA_CALLABLE size_t size(){
 		return R[0] * R[1];
@@ -48,6 +245,20 @@ public:
 		return sizeof(stim::complex<T>) * R[0] * R[1];
 	}
+	scalarfield(size_t X, size_t Y, T size = 1, T z_pos = 0) : rect<T>::rect(size, z_pos){
+		R[0] = X;											//set the field resolution
+		R[1] = Y;
+
+		E = (stim::complex<T>*) malloc(grid_bytes());		//allocate in CPU memory
+		memset(E, 0, grid_bytes());
+		loc = CPUmem;
+	}
+
+	~scalarfield(){
+		if(loc == CPUmem) free(E);
+		else cudaFree(E);
+	}	
+
 	/// Calculates the distance between points on the grid
 	T spacing(){
 		T du = rect<T>::X.len() / R[0];
@@ -78,6 +289,16 @@ public:
 		}
 	}
+	/// Propagate the field along its orthogonal direction by a distance d
+	void propagate(T d, T k){
+		cpu_scalar_propagate(E, E, X.len(), Y.len(), d, k, R[0], R[1]);
+	}
+
+	/// Propagate the field along its orthogonal direction by a distance d
+	void lowpass(T highest){
+		cpu_scalar_lowpass(E, E, X.len(), Y.len(), highest, R[0], R[1]);
+	}
+
 	std::string str(){
 		std::stringstream ss;
 		ss<<rect<T>::str()<<std::endl;
@@ -96,11 +317,11 @@ public:
 	/// Evaluate the cartesian coordinates of each point in the field. The resulting arrays are allocated in the same memory where the field is stored.
 	void meshgrid(T* X, T* Y, T* Z, locationType location){
-		size_t array_size = sizeof(T) * R[0] * R[1];
+		//size_t array_size = sizeof(T) * R[0] * R[1];
 		if(location == CPUmem){
-			T du = 1.0 / (R[0] - 1);					//calculate the spacing between points in the grid
-			T dv = 1.0 / (R[1] - 1);
+			T du = (T)1.0 / (R[0] - 1);					//calculate the spacing between points in the grid
+			T dv = (T)1.0 / (R[1] - 1);
 			size_t ui, vi, i;
 			stim::vec3<T> p;
@@ -114,9 +335,9 @@ public:
 					i++;					
 				}
 			}
-			stim::cpu2image(X, "X.bmp", R[0], R[1], stim::cmBrewer);
-			stim::cpu2image(Y, "Y.bmp", R[0], R[1], stim::cmBrewer);
-			stim::cpu2image(Z, "Z.bmp", R[0], R[1], stim::cmBrewer);
+			//stim::cpu2image(X, "X.bmp", R[0], R[1], stim::cmBrewer);
+			//stim::cpu2image(Y, "Y.bmp", R[0], R[1], stim::cmBrewer);
+			//stim::cpu2image(Z, "Z.bmp", R[0], R[1], stim::cmBrewer);
 		}
 		else{
 			std::cout<<"GPU allocation of a meshgrid isn't supported yet. You'll have to write kernels to do the calculation.";
@@ -124,6 +345,14 @@ public:
 		}
 	}
+	//clear the field, setting all values to zero
+	void clear(){
+		if(loc == GPUmem)
+			HANDLE_ERROR(cudaMemset(E, 0, grid_bytes()));
+		else
+			memset(E, 0, grid_bytes());
+	}
+
 	void image(std::string filename, stim::complexComponentType type = complexMag, stim::colormapType cmap = stim::cmBrewer){
 		if(loc == GPUmem) to_cpu();									//if the field is in the GPU, move it to the CPU
@@ -366,7 +366,7 @@ __global__ void cuda_scalar_mie_internal(stim::complex&lt;T&gt;* E, size_t N, T* x, T*
 	(z == NULL) ? p[2] = 0 : p[2] = z[i];
 	T r = p.len();														//calculate the distance from the sphere
-	if(r > a) return;													//exit if the point is inside the sphere (we only calculate the internal field)
+	if(r >= a) return;													//exit if the point is inside the sphere (we only calculate the internal field)
 	T fij = (r - r_min)/dr;											//FP index into the spherical bessel LUT
 	size_t ij = (size_t) fij;											//convert to an integral index
 	T alpha = fij - ij;													//calculate the fractional portion of the index
@@ -608,6 +608,39 @@ void cpu_scalar_mie_internal(stim::complex&lt;T&gt;* E, size_t N, T* x, T* y, T* z, st
 	cpu_scalar_mie_internal(E, N, x, y, z, W, a, n, r_spacing);
 }
-}
+
+/// Class stim::scalarmie represents a scalar Mie scattering model that can be used to calculate the fields produced by a scattering sphere.
+template<typename T>
+class scalarmie
+{
+private:
+	T radius;					//radius of the scattering sphere
+	stim::complex<T> n;			//refractive index of the scattering sphere
+	
+public:
+
+	scalarmie(T r, stim::complex<T> ri){
+		radius = r;
+		n = ri;
+	}
+
+	void eval(stim::scalarfield<T>& E, stim::scalarbeam<T> b, int order = 500, int samples = 1000){
+
+		size_t array_size = E.grid_bytes();											//calculate the number of bytes in the scalar grid
+		float* X = (float*) malloc( array_size );									//allocate space for the coordinate meshes
+		float* Y = (float*) malloc( array_size );
+		float* Z = (float*) malloc( array_size );
+		E.meshgrid(X, Y, Z, stim::CPUmem);											//calculate the coordinate meshes
+
+		b.eval(E, X, Y, Z, order);													//evaluate the incident field using a plane wave expansion
+
+		std::vector< stim::scalarwave<float> > wave_array = b.mc(samples);			//decompose the beam into an array of plane waves
+		stim::cpu_scalar_mie_scatter<float>(E.ptr(), E.size(), X, Y, Z, wave_array, radius, n, E.spacing());
+		stim::cpu_scalar_mie_internal<float>(E.ptr(), E.size(), X, Y, Z, wave_array, radius, n, E.spacing());
+	}
+
+};			//end stim::scalarmie
+
+}			//end namespace stim
 #endif
 \ No newline at end of file
@@ -210,7 +210,7 @@ template&lt;typename T&gt;
 __global__ void cuda_scalarwave(stim::complex<T>* F, size_t N, T* x, T* y, T* z, stim::scalarwave<T>* W, size_t n_waves){
 	extern __shared__ stim::scalarwave<T> shared_W[];		//declare the list of waves in shared memory
-	stim::cuda::sharedMemcpy(shared_W, W, n_waves, threadIdx.x, blockDim.x);	//copy the plane waves into shared memory for faster access
+	stim::cuda::threadedMemcpy(shared_W, W, n_waves, threadIdx.x, blockDim.x);	//copy the plane waves into shared memory for faster access
 	__syncthreads();															//synchronize threads to insure all data is copied
 	size_t i = blockIdx.x * blockDim.x + threadIdx.x;				//get the index into the array
@@ -42,11 +42,13 @@ namespace stim{
 		bool char_numerical(char c){
 			if( (c >= '0' && c <= '9') || c == '-' || c == '.')
 				return true;
+			return false;
 		}
 		bool char_integral(char c){
 			if( (c >= '0' && c <= '9') || c == '-')
 				return true;
+			return false;
 		}
 		/// Test if a given string contains a numerical (real) value
@@ -475,31 +477,31 @@ namespace stim{
         ///Determines of a parameter has been set and returns true if it has
         /// @param _name is the name of the argument
-        bool operator()(std::string _name)
-        {
-            size_t i = find(opts.begin(), opts.end(), _name) - opts.begin();
+        bool operator()(std::string _name){
+			std::vector<cmd_option>::iterator it;
+            it = std::find(opts.begin(), opts.end(), _name);// - opts.begin();
-            if(i < 0){
+            if(it == opts.end()){
                 std::cout<<"ERROR - Unspecified parameter name: "<<_name<<std::endl;
                 exit(1);
             }
-            return opts[i].is_set();
+            return it->is_set();
         }
         ///Returns the number of parameters for a specified argument
         /// @param _name is the name of the argument whose parameter number will be returned
-        size_t nargs(std::string _name)
-        {
-            size_t i = find(opts.begin(), opts.end(), _name) - opts.begin();
+        size_t nargs(std::string _name){
+			std::vector<cmd_option>::iterator it;
+            it = find(opts.begin(), opts.end(), _name);// - opts.begin();
-            if(i < 0){
+            if(it == opts.end()){
                 std::cout<<"ERROR - Unspecified parameter name: "<<_name<<std::endl;
                 exit(1);
             }
-            return opts[i].nargs();
+            return it->nargs();
         }
         ///Returns the number of arguments that have been set
@@ -525,17 +527,16 @@ namespace stim{
         ///Returns an object describing the argument
         /// @param _name is the name of the requested argument
-        cmd_option operator[](std::string _name)
-        {
-            size_t i = find(opts.begin(), opts.end(), _name) - opts.begin();
+        cmd_option operator[](std::string _name){
+			std::vector<cmd_option>::iterator it;
+            it = find(opts.begin(), opts.end(), _name);// - opts.begin();
-            if(i < 0 || i >= opts.size())
-            {
+            if(it == opts.end()){
                 std::cout<<"ERROR - Unspecified parameter name: "<<_name<<std::endl;
                 exit(1);
             }
-            return opts[i];
+            return *it;
         }
@@ -68,6 +68,7 @@ protected:
 	//parse a file locator string
 	void parse(std::string loc){
+		if(loc.size() == 0) return;
 		loc = unix_div(loc);
@@ -307,6 +308,13 @@ public:
 		parse_name(fname);
 	}
+	//append a string to the filename and return a new filename
+	stim::filename append(std::string s){
+		stim::filename result = *this;		//create a new filename, copy the current filename
+		result.prefix = prefix + s;			//append the string to the filename
+		return result;
+	}
+
 	//get a path relative to the current one
 	stim::filename get_relative(std::string rel){