 
 
 int main(int argc, char *argv[]) {
+
+	//Check if the input is sensible
+
 	if (argc == 1) {
 		std::cout << "Please specify as an argument the path to the .vtk file to use"  << std::endl;
 		return 1;

	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;
	checkCudaErrors(cudaMalloc(&d_radii, dataCount*sizeof(float)));
//	float *h_radii;
//	h_radii = (float*) malloc((dataCount/steps)*sizeof(float));

	//Compute the length of each line (locally)
	rxy<<<dataCount/blockSize,blockSize>>>(d_lines, d_radii, (float)steps, d_origins, steps);

	//Add shit together per line
	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);

	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);

	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;

	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);
	divide<<<gridSizeRK4,blockSizeRK4>>>(d_alpha, d_beta, d_alpha);
	float* h_windingdata = (float*) malloc((dataCount/steps)*sizeof(float));
	//Copy lengths from device to host
	checkCudaErrors(cudaMemcpy(h_windingdata, d_alpha, (dataCount/steps)*sizeof(float), cudaMemcpyDeviceToHost));
	
	//Print some data to screen
	dataprint(h_windingdata,blockSizeRK4);
		

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

	//Copy our texture properties (linear interpolation, texture access) to data array on device
	checkCudaErrors(cudaBindTextureToArray(dataTex, dataArray, channelDesc));

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

	//Declare pointers to arrays with line data (output of integration), one each on device and host
	float4 *d_lines, *h_lines;

