 		std::cout << "Warning: origin read from file not equal to origin from constants.h" << std::endl;
 	}
 
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
+
 	//Allocate array on device
 	cudaArray* dataArray;
 	cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32,32,32,32,cudaChannelFormatKindFloat);

	free(hostvfield[0][0]);
	free(hostvfield[0]);
	free(hostvfield);
	free(h_lines);
	free(h_lengths);
	cudaFree(d_origins);
//	free(h_radii);
	cudaFree(d_alpha);
	cudaFree(d_beta);
	free(h_windingdata);
	cudaFree(d_radii);
	cudaFree(d_lines);
	cudaFree(d_lengths);
	cudaFreeArray(dataArray);
	return 0;
}

	reduceSum<<<dataCount/(2*blockSize),blockSize,blockSize*sizeof(float)>>>(d_lengths,d_lengths);
	reduceSum<<<dataCount/steps,steps/(4*blockSize),steps/(4*blockSize)*sizeof(float)>>>(d_lengths,d_lengths);

	//Copy lengths from device to host
	checkCudaErrors(cudaMemcpy(h_lengths, d_lengths, (dataCount/steps)*sizeof(float), cudaMemcpyDeviceToHost));

	cudaFree(d_lengths);

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

	//Allocating the array to store the radius data, both for host and device
	float *d_radii, *h_radii;
	checkCudaErrors(cudaMalloc(&d_radii, dataCount*sizeof(float)));
	h_radii = (float*) malloc((dataCount/steps)*sizeof(float));


	//Compute the distance from the origin in the xy plane of each point
	rxy<<<dataCount/blockSize,blockSize>>>(d_lines, d_radii, (float)steps, d_origins, steps);

	//Average these distances to find the torus radius
	reduceSum<<<dataCount/(2*blockSize),blockSize,blockSize*sizeof(float)>>>(d_radii,d_radii);
	reduceSum<<<dataCount/steps,steps/(4*blockSize),steps/(4*blockSize)*sizeof(float)>>>(d_radii,d_radii);

	//Copy radii from device to host
	checkCudaErrors(cudaMemcpy(h_radii, d_radii, (dataCount/steps)*sizeof(float), cudaMemcpyDeviceToHost));

	//r_radii are still needed, so no memory free just yet.

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//
	
	//Allocating arrays to store the alpha and beta data to compute winding numbers.
	float* h_windingdata = (float*) malloc((dataCount/steps)*sizeof(float));
	float *d_alpha, *d_beta;
	checkCudaErrors(cudaMalloc(&d_alpha, dataCount*sizeof(float)));
	checkCudaErrors(cudaMalloc(&d_beta, dataCount*sizeof(float)));

	winding<<<dataCount/blockSize,blockSize>>>(d_lines, d_alpha, d_beta, d_origins, d_radii, steps);

	//A check to evaluate the alpha/beta steps found
	float* h_alphacheck = (float*) malloc(steps*sizeof(float));
	checkCudaErrors(cudaMemcpy(h_alphacheck, &d_beta[0], steps*sizeof(float), cudaMemcpyDeviceToHost));

	for(unsigned int i=0; i < 500; i++) {
		std::cout << h_alphacheck[i] << ", ";
	}
	std::cout << std::endl;

	

	//Adding the steps Deltaalpha and Deltabeta to find overall windings
	reduceSum<<<dataCount/(2*blockSize),blockSize,blockSize*sizeof(float)>>>(d_alpha, d_alpha);
	reduceSum<<<dataCount/steps,steps/(4*blockSize),steps/(4*blockSize)*sizeof(float)>>>(d_alpha, d_alpha);

	reduceSum<<<dataCount/(2*blockSize),blockSize,blockSize*sizeof(float)>>>(d_beta, d_beta);
	reduceSum<<<dataCount/steps,steps/(4*blockSize),steps/(4*blockSize)*sizeof(float)>>>(d_beta, d_beta);

	//Dividing these windings to compute the winding numbers and store them in d_alpha
	divide<<<gridSizeRK4,blockSizeRK4>>>(d_alpha, d_beta, d_alpha);

	//Copy winding numbers from from device to host
	checkCudaErrors(cudaMemcpy(h_windingdata, d_alpha, (dataCount/steps)*sizeof(float), cudaMemcpyDeviceToHost));
	
	cudaFree(d_alpha);
	cudaFree(d_beta);
	cudaFree(d_radii);

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%//

	//Print some data to screen
	dataprint(h_windingdata,blockSizeRK4);

	//Write all the lines
	datawrite("../datadir/data.bin", dataCount, h_lines);
