diff --git a/bioem.cpp b/bioem.cpp
index 9720b80922020f4298c834d8e4bfb0f41d9ad0ae..aac4177eeaa4c192cac2f932de4d54250569d99e 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -54,14 +54,14 @@ int bioem::configure(int ac, char* av[])
 	//*************************************************************************************
 
 	/*** Inizialzing default variables ***/
-	std::string infile,modelfile,mapfile;
-	Model.readPDB=false;
-	param.writeAngles=false;
+	std::string infile, modelfile, mapfile;
+	Model.readPDB = false;
+	param.writeAngles = false;
 	RefMap.dumpMap = false;
 	RefMap.loadMap = false;
-        RefMap.readMRC = false;
+	RefMap.readMRC = false;
 	RefMap.readMultMRC = false;
-	
+
 
 	//*************************************************************************************
 	cout << "+++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
@@ -76,8 +76,8 @@ int bioem::configure(int ac, char* av[])
 		("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file")
 		("Particlesfile", po::value< std::string>(), "(Mandatory) Name of paricles file")
 		("ReadPDB", "(Optional) If reading model file in PDB format")
-                ("ReadMRC", "(Optional) If reading particle file in MRC format")
-		 ("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
+		("ReadMRC", "(Optional) If reading particle file in MRC format")
+		("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
 		("DumpMaps", "(Optional) Dump maps after they were red from maps file")
 		("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
 		("help", "(Optional) Produce help message")
@@ -88,8 +88,8 @@ int bioem::configure(int ac, char* av[])
 		p.add("Modelfile", -1);
 		p.add("Particlesfile", -1);
 		p.add("ReadPDB", -1);
-                p.add("ReadMRC", -1);
-                p.add("ReadMultipleMRC", -1);
+		p.add("ReadMRC", -1);
+		p.add("ReadMultipleMRC", -1);
 		p.add("DumpMaps", -1);
 		p.add("LoadMapDump", -1);
 
@@ -99,9 +99,9 @@ int bioem::configure(int ac, char* av[])
 		po::notify(vm);
 
 		if((ac < 6)) {
-			std::cout << desc << std::endl;			
+			std::cout << desc << std::endl;
 			exit(1);
-			return 0;			
+			return 0;
 		}
 		if (vm.count("help")) {
 			cout << "Usage: options_description [options]\n";
@@ -112,32 +112,32 @@ int bioem::configure(int ac, char* av[])
 		if (vm.count("Inputfile"))
 		{
 			cout << "Input file is: ";
-			cout << vm["Inputfile"].as< std::string >()<< "\n";
-			infile=vm["Inputfile"].as< std::string >();			
+			cout << vm["Inputfile"].as< std::string >() << "\n";
+			infile = vm["Inputfile"].as< std::string >();
 		}
 		if (vm.count("Modelfile"))
 		{
 			cout << "Model file is: "
 				 << vm["Modelfile"].as<  std::string  >() << "\n";
-			modelfile=vm["Modelfile"].as<  std::string  >();
+			modelfile = vm["Modelfile"].as<  std::string  >();
 		}
 
 		if (vm.count("ReadPDB"))
 		{
 			cout << "Reading model file in PDB format.\n";
-			Model.readPDB=true;
+			Model.readPDB = true;
 		}
 		if (vm.count("ReadMRC"))
-                {
-                        cout << "Reading particle file in MRC format.\n";
-                        RefMap.readMRC=true;
-                }
-                
+		{
+			cout << "Reading particle file in MRC format.\n";
+			RefMap.readMRC = true;
+		}
+
 		if (vm.count("ReadMultipleMRC"))
-                {
-                        cout << "Reading Multiple MRCs.\n";
-                        RefMap.readMultMRC=true;
-                }
+		{
+			cout << "Reading Multiple MRCs.\n";
+			RefMap.readMultMRC = true;
+		}
 
 		if (vm.count("DumpMaps"))
 		{
@@ -155,7 +155,7 @@ int bioem::configure(int ac, char* av[])
 		{
 			cout << "Paricle file is: "
 				 << vm["Particlesfile"].as< std::string >() << "\n";
-			mapfile=vm["Particlesfile"].as< std::string >();
+			mapfile = vm["Particlesfile"].as< std::string >();
 		}
 	}
 	catch(std::exception& e)
@@ -203,15 +203,15 @@ int bioem::precalculate()
 	// Generating Grids of orientations
 	param.CalculateGridsParam();
 
-	myfloat_t sum,sumsquare;
+	myfloat_t sum, sumsquare;
 
 	//Precalculating cross-correlations of maps
 	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap ; iRefMap++)
 	{
-		calcross_cor(RefMap.Ref[iRefMap],sum,sumsquare);
+		calcross_cor(RefMap.Ref[iRefMap], sum, sumsquare);
 		//Storing Crosscorrelations in Map class
-		RefMap.sum_RefMap[iRefMap]=sum;
-		RefMap.sumsquare_RefMap[iRefMap]=sumsquare;
+		RefMap.sum_RefMap[iRefMap] = sum;
+		RefMap.sumsquare_RefMap[iRefMap] = sumsquare;
 	}
 
 	// Precalculating CTF Kernels stored in class Param
@@ -238,13 +238,13 @@ int bioem::run()
 	// Inizialzing Probabilites to zero and constant to -Infinity
 	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
 	{
-		pProb[iRefMap].Total=0.0;
-		pProb[iRefMap].Constoadd=-9999999;
-		pProb[iRefMap].max_prob=-9999999;
+		pProb[iRefMap].Total = 0.0;
+		pProb[iRefMap].Constoadd = -9999999;
+		pProb[iRefMap].max_prob = -9999999;
 		for (int iOrient = 0; iOrient < param.nTotGridAngles; iOrient ++)
 		{
-			pProb[iRefMap].forAngles[iOrient]=0.0;
-			pProb[iRefMap].ConstAngle[iOrient]=-99999999;
+			pProb[iRefMap].forAngles[iOrient] = 0.0;
+			pProb[iRefMap].ConstAngle[iOrient] = -99999999;
 		}
 	}
 	//**************************************************************************************
@@ -256,11 +256,11 @@ int bioem::run()
 	mycomplex_t* proj_mapFFT;
 	bioem_map conv_map;
 	mycomplex_t* conv_mapFFT;
-	myfloat_t sumCONV,sumsquareCONV;
+	myfloat_t sumCONV, sumsquareCONV;
 
 	//allocating fftw_complex vector
-	proj_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
-	conv_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
+	proj_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
+	conv_mapFFT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
 
 	HighResTimer timer;
 
@@ -322,7 +322,7 @@ int bioem::run()
 	for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
 	{
 		//**** Total Probability ***
-		outputProbFile << "RefMap " << iRefMap << " Probability  "  << log(pProb[iRefMap].Total)+pProb[iRefMap].Constoadd+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd  << "\n";
+		outputProbFile << "RefMap " << iRefMap << " Probability  "  << log(pProb[iRefMap].Total) + pProb[iRefMap].Constoadd + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " Constant " << pProb[iRefMap].Constoadd  << "\n";
 
 		outputProbFile << "RefMap " << iRefMap << " Maximizing Param: ";
 
@@ -344,7 +344,7 @@ int bioem::run()
 		{
 			for (int iProjectionOut = 0; iProjectionOut < param.nTotGridAngles; iProjectionOut++)
 			{
-				angProbfile << " " << iRefMap << " " << param.angles[iProjectionOut].pos[0] << " " << param.angles[iProjectionOut].pos[1] << " " << param.angles[iProjectionOut].pos[2] << " " << log(pProb[iRefMap].forAngles[iProjectionOut])+pProb[iRefMap].ConstAngle[iProjectionOut]+0.5*log(M_PI)+(1-param.param_device.Ntotpi*0.5)*(log(2*M_PI)+1)+log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
+				angProbfile << " " << iRefMap << " " << param.angles[iProjectionOut].pos[0] << " " << param.angles[iProjectionOut].pos[1] << " " << param.angles[iProjectionOut].pos[2] << " " << log(pProb[iRefMap].forAngles[iProjectionOut]) + pProb[iRefMap].ConstAngle[iProjectionOut] + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
 
 			}
 		}
@@ -364,7 +364,7 @@ int bioem::run()
 	if (param.refCTF)
 	{
 		delete[] param.refCTF;
-		param.refCTF =NULL;
+		param.refCTF = NULL;
 	}
 
 	if(RefMap.RefMapsFFT)
@@ -378,17 +378,17 @@ int bioem::run()
 int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
 {
 
-		//***************************************************************************************
-                //***** BioEM routine for comparing reference maps to convoluted maps *****
+	//***************************************************************************************
+	//***** BioEM routine for comparing reference maps to convoluted maps *****
 
 	if (FFTAlgo) //IF using the fft algorithm
 	{
-#pragma omp parallel
+		#pragma omp parallel
 		{
 			mycomplex_t *localCCT;
 			myfloat_t *lCC;
-			localCCT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
-			lCC= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
+			localCCT = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
+			lCC = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
 
 			const int num_threads = omp_get_num_threads();
 			const int thread_id = omp_get_thread_num();
@@ -398,7 +398,7 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
 
 			for (int iRefMap = iStart; iRefMap < iEnd; iRefMap ++)
 			{
-				calculateCCFFT(iRefMap,iProjectionOut, iConv, sumC, sumsquareC, localmultFFT, localCCT,lCC);
+				calculateCCFFT(iRefMap, iProjectionOut, iConv, sumC, sumsquareC, localmultFFT, localCCT, lCC);
 			}
 			myfftw_free(localCCT);
 			myfftw_free(lCC);
@@ -406,34 +406,34 @@ int bioem::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_m
 	}
 	else
 	{
-#pragma omp parallel for
+		#pragma omp parallel for
 		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
 		{
-			compareRefMapShifted<-1>(iRefMap,iProjectionOut,iConv,conv_map, pProb, param.param_device, RefMap);
+			compareRefMapShifted < -1 > (iRefMap, iProjectionOut, iConv, conv_map, pProb, param.param_device, RefMap);
 		}
 	}
 	return(0);
 }
 
-inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,myfloat_t* lCC)
+inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
 {
 
-		//***************************************************************************************
-                //***** Calculating cross correlation in FFTALGOrithm *****
+	//***************************************************************************************
+	//***** Calculating cross correlation in FFTALGOrithm *****
 
 	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.RefMapSize];
-	for(int i = 0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
+	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
 	{
 		localCCT[i][0] = localConvFFT[i][0] * RefMapFFT[i][0] + localConvFFT[i][1] * RefMapFFT[i][1];
 		localCCT[i][1] = localConvFFT[i][1] * RefMapFFT[i][0] - localConvFFT[i][0] * RefMapFFT[i][1];
 	}
 
-	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,localCCT,lCC);
+	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
 
 	doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
 }
 
-int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
+int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
 {
 	//**************************************************************************************
 	//****  BioEM Create Projection routine in Euler angle predefined grid****************
@@ -442,43 +442,43 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
 
 	myfloat3_t RotatedPointsModel[Model.nPointsModel];
 	myfloat_t rotmat[3][3];
-	myfloat_t alpha, gam,beta;
+	myfloat_t alpha, gam, beta;
 	myfloat_t* localproj;
 
-	localproj= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
-	memset(localproj,0,param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj));
+	localproj = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
+	memset(localproj, 0, param.param_device.NumberPixels * param.param_device.NumberPixels * sizeof(*localproj));
 
-	alpha=param.angles[iMap].pos[0];
-	beta=param.angles[iMap].pos[1];
-	gam=param.angles[iMap].pos[2];
+	alpha = param.angles[iMap].pos[0];
+	beta = param.angles[iMap].pos[1];
+	gam = param.angles[iMap].pos[2];
 
 	//**** To see how things are going: cout << "Id " << omp_get_thread_num() <<  " Angs: " << alpha << " " << beta << " " << gam << "\n"; ***
 
 	//********** Creat Rotation with pre-defiend grid of orientations**********
 
-	rotmat[0][0]=cos(gam)*cos(alpha)-cos(beta)*sin(alpha)*sin(gam);
-	rotmat[0][1]=cos(gam)*sin(alpha)+cos(beta)*cos(alpha)*sin(gam);
-	rotmat[0][2]=sin(gam)*sin(beta);
-	rotmat[1][0]=-sin(gam)*cos(alpha)-cos(beta)*sin(alpha)*cos(gam);
-	rotmat[1][1]=-sin(gam)*sin(alpha)+cos(beta)*cos(alpha)*cos(gam);
-	rotmat[1][2]=cos(gam)*sin(beta);
-	rotmat[2][0]=sin(beta)*sin(alpha);
-	rotmat[2][1]=-sin(beta)*cos(alpha);
-	rotmat[2][2]=cos(beta);
-
-	for(int n=0; n< Model.nPointsModel; n++)
+	rotmat[0][0] = cos(gam) * cos(alpha) - cos(beta) * sin(alpha) * sin(gam);
+	rotmat[0][1] = cos(gam) * sin(alpha) + cos(beta) * cos(alpha) * sin(gam);
+	rotmat[0][2] = sin(gam) * sin(beta);
+	rotmat[1][0] = -sin(gam) * cos(alpha) - cos(beta) * sin(alpha) * cos(gam);
+	rotmat[1][1] = -sin(gam) * sin(alpha) + cos(beta) * cos(alpha) * cos(gam);
+	rotmat[1][2] = cos(gam) * sin(beta);
+	rotmat[2][0] = sin(beta) * sin(alpha);
+	rotmat[2][1] = -sin(beta) * cos(alpha);
+	rotmat[2][2] = cos(beta);
+
+	for(int n = 0; n < Model.nPointsModel; n++)
 	{
-		RotatedPointsModel[n].pos[0]=0.0;
-		RotatedPointsModel[n].pos[1]=0.0;
-		RotatedPointsModel[n].pos[2]=0.0;
+		RotatedPointsModel[n].pos[0] = 0.0;
+		RotatedPointsModel[n].pos[1] = 0.0;
+		RotatedPointsModel[n].pos[2] = 0.0;
 	}
-	for(int n=0; n< Model.nPointsModel; n++)
+	for(int n = 0; n < Model.nPointsModel; n++)
 	{
-		for(int k=0; k< 3; k++)
+		for(int k = 0; k < 3; k++)
 		{
-			for(int j=0; j< 3; j++)
+			for(int j = 0; j < 3; j++)
 			{
-				RotatedPointsModel[n].pos[k]+=rotmat[k][j]*Model.PointsModel[n].pos[j];
+				RotatedPointsModel[n].pos[k] += rotmat[k][j] * Model.PointsModel[n].pos[j];
 			}
 		}
 	}
@@ -486,41 +486,41 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
 	int i, j;
 
 	//************ Projection over the Z axis********************
-	for(int n=0; n< Model.nPointsModel; n++)
+	for(int n = 0; n < Model.nPointsModel; n++)
 	{
 		//Getting pixel that represents coordinates & shifting the start at to Numpix/2,Numpix/2 )
-		i=floor(RotatedPointsModel[n].pos[0]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
-		j=floor(RotatedPointsModel[n].pos[1]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
+		i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
+		j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
 
-		localproj[i*param.param_device.NumberPixels+j]+=Model.densityPointsModel[n]/Model.NormDen;
+		localproj[i * param.param_device.NumberPixels + j] += Model.densityPointsModel[n] / Model.NormDen;
 	}
 
 	//**** Output Just to check****
-	if(iMap==10)
+	if(iMap == 10)
 	{
 		ofstream myexamplemap;
 		ofstream myexampleRot;
 		myexamplemap.open ("MAP_i10");
 		myexampleRot.open ("Rot_i10");
 		myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
-		for(int k=0; k<param.param_device.NumberPixels; k++)
+		for(int k = 0; k < param.param_device.NumberPixels; k++)
 		{
-			for(int j=0; j<param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*param.param_device.NumberPixels+j];
+			for(int j = 0; j < param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j << " " << localproj[k * param.param_device.NumberPixels + j];
 		}
 		myexamplemap << " \n";
-		for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
+		for(int n = 0; n < Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
 		myexamplemap.close();
 		myexampleRot.close();
 	}
 
 	//***** Converting projection to Fourier Space for Convolution later with kernel****
 	//********** Omp Critical is necessary with FFTW*******
-	myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localproj,mapFFT);
+	myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localproj, mapFFT);
 
 	return(0);
 }
 
-int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
+int bioem::createConvolutedProjectionMap(int iMap, int iConv, mycomplex_t* lproj, bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
 {
 	//**************************************************************************************
 	//****  BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
@@ -529,14 +529,14 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
 	//**************************************************************************************
 
 	myfloat_t* localconvFFT;
-	localconvFFT= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t)*param.param_device.NumberPixels*param.param_device.NumberPixels);
+	localconvFFT = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
 	mycomplex_t* tmp;
 	tmp = (mycomplex_t*) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
 
 	//**** Multiplying FFTmap with corresponding kernel ****
 
 	const mycomplex_t* refCTF = &param.refCTF[iConv * param.RefMapSize];
-	for(int i=0;i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D;i++)
+	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
 	{
 		localmultFFT[i][0] = lproj[i][0] * refCTF[i][0] + lproj[i][1] * refCTF[i][1];
 		localmultFFT[i][1] = lproj[i][1] * refCTF[i][0] - lproj[i][0] * refCTF[i][1];
@@ -547,20 +547,20 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
 	memcpy(tmp, localmultFFT, sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
 
 	//**** Bringing convoluted Map to real Space ****
-	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,tmp,localconvFFT);
+	myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, tmp, localconvFFT);
 
 	//****Asigning convolution fftw_complex to bioem_map ****
-	for(int i=0; i < param.param_device.NumberPixels ; i++ )
+	for(int i = 0; i < param.param_device.NumberPixels ; i++ )
 	{
-		for(int j=0; j < param.param_device.NumberPixels ; j++ )
+		for(int j = 0; j < param.param_device.NumberPixels ; j++ )
 		{
-			Mapconv.points[i][j]=localconvFFT[i*param.param_device.NumberPixels+j];
+			Mapconv.points[i][j] = localconvFFT[i * param.param_device.NumberPixels + j];
 		}
 	}
 
 	//*** Calculating Cross-correlations of cal-convoluted map with its self *****
-	sumC=0;
-	sumsquareC=0;
+	sumC = 0;
+	sumsquareC = 0;
 	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberPixels; i++)
 	{
 		sumC += localconvFFT[i];
@@ -580,12 +580,12 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
 	return(0);
 }
 
-int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
+int bioem::calcross_cor(bioem_map& localmap, myfloat_t& sum, myfloat_t& sumsquare)
 {
 	//*********************** Routine to calculate Cross correlations***********************
 
-	sum=0.0;
-	sumsquare=0.0;
+	sum = 0.0;
+	sumsquare = 0.0;
 	for (int i = 0; i < param.param_device.NumberPixels; i++)
 	{
 		for (int j = 0; j < param.param_device.NumberPixels; j++)
@@ -593,7 +593,7 @@ int bioem::calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare)
 			// Calculate Sum of pixels
 			sum += localmap.points[i][j];
 			// Calculate Sum of pixels squared
-			sumsquare += localmap.points[i][j]*localmap.points[i][j];
+			sumsquare += localmap.points[i][j] * localmap.points[i][j];
 		}
 	}
 	return(0);
diff --git a/bioem_algorithm.h b/bioem_algorithm.h
index c80876a7dd7cf6f9228f964c8ce7f715d8b4e6fb..da35df44ff686d0041f1962681b99c7d1e8ca41c 100644
--- a/bioem_algorithm.h
+++ b/bioem_algorithm.h
@@ -61,8 +61,8 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
 	const myfloat_t ForLogProb = (sumsquare * param.Ntotpi - sum * sum);
 
 	// Products of different cross-correlations (first element in formula)
-	const myfloat_t firstele = param.Ntotpi * (sumsquareref * sumsquare-crossproMapConv * crossproMapConv) +
-							2 * sumref * sum * crossproMapConv - sumsquareref * sum * sum - sumref * sumref * sumsquare;
+	const myfloat_t firstele = param.Ntotpi * (sumsquareref * sumsquare - crossproMapConv * crossproMapConv) +
+							   2 * sumref * sum * crossproMapConv - sumsquareref * sum * sum - sumref * sumref * sumsquare;
 
 	//******* Calculating log of Prob*********/
 	// As in fortran code: logpro=(3-Ntotpi)*0.5*log(firstele/pConvMap[iOrient].ForLogProbfromConv[iConv])+(Ntotpi*0.5-2)*log(Ntotpi-2)-0.5*log(pConvMap[iOrient].ForLogProbfromConv[iConv])+0.5*log(PI)+(1-Ntotpi*0.5)*(log(2*PI)+1);
@@ -70,7 +70,7 @@ __device__ static inline myfloat_t calc_logpro(const bioem_param_device& param,
 	return(logpro);
 }
 
-__device__ static inline void calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC, float value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
+__device__ static inline void calProb(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, float value, int disx, int disy, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
 {
 	/********************************************************/
 	/*********** Calculates the BioEM probability ***********/
@@ -80,8 +80,8 @@ __device__ static inline void calProb(int iRefMap,int iOrient, int iConv,myfloat
 
 	//update_prob<-1>(logpro, iRefMap, iOrient, iConv, disx, disy, pProb);
 	//GCC is too stupid to inline properly, so the code is copied here
-    if(pProb[iRefMap].Constoadd < logpro)
-    {
+	if(pProb[iRefMap].Constoadd < logpro)
+	{
 		pProb[iRefMap].Total = pProb[iRefMap].Total * exp(-logpro + pProb[iRefMap].Constoadd);
 		pProb[iRefMap].Constoadd = logpro;
 	}
@@ -106,33 +106,33 @@ __device__ static inline void calProb(int iRefMap,int iOrient, int iConv,myfloat
 
 __device__ static inline void doRefMapFFT(const int iRefMap, const int iOrient, const int iConv, const myfloat_t* lCC, const myfloat_t sumC, const myfloat_t sumsquareC, bioem_Probability* pProb, const bioem_param_device& param, const bioem_RefMap& RefMap)
 {
-	for (int cent_x = 0; cent_x <= param.maxDisplaceCenter; cent_x=cent_x+param.GridSpaceCenter)
+	for (int cent_x = 0; cent_x <= param.maxDisplaceCenter; cent_x = cent_x + param.GridSpaceCenter)
 	{
-		for (int cent_y = 0; cent_y <= param.maxDisplaceCenter; cent_y=cent_y+param.GridSpaceCenter)
+		for (int cent_y = 0; cent_y <= param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y, pProb, param, RefMap);
+			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y, pProb, param, RefMap);
 		}
-		for (int cent_y = param.NumberPixels-param.maxDisplaceCenter; cent_y < param.NumberPixels; cent_y=cent_y+param.GridSpaceCenter)
+		for (int cent_y = param.NumberPixels - param.maxDisplaceCenter; cent_y < param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels*param.NumberPixels), cent_x, cent_y - param.NumberPixels, pProb, param, RefMap);
+			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x, cent_y - param.NumberPixels, pProb, param, RefMap);
 		}
 	}
-	for (int cent_x = param.NumberPixels-param.maxDisplaceCenter; cent_x < param.NumberPixels; cent_x=cent_x+param.GridSpaceCenter)
+	for (int cent_x = param.NumberPixels - param.maxDisplaceCenter; cent_x < param.NumberPixels; cent_x = cent_x + param.GridSpaceCenter)
 	{
-		for (int cent_y = 0; cent_y < param.maxDisplaceCenter; cent_y=cent_y+param.GridSpaceCenter)
+		for (int cent_y = 0; cent_y < param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels*param.NumberPixels), cent_x - param.NumberPixels, cent_y, pProb, param, RefMap);
+			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x - param.NumberPixels, cent_y, pProb, param, RefMap);
 		}
-		for (int cent_y = param.NumberPixels-param.maxDisplaceCenter; cent_y <= param.NumberPixels; cent_y=cent_y+param.GridSpaceCenter)
+		for (int cent_y = param.NumberPixels - param.maxDisplaceCenter; cent_y <= param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
 		{
-			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.NumberPixels+cent_y]/ (myfloat_t) (param.NumberPixels*param.NumberPixels), cent_x - param.NumberPixels, cent_y - param.NumberPixels, pProb, param, RefMap);
+			calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] / (myfloat_t) (param.NumberPixels * param.NumberPixels), cent_x - param.NumberPixels, cent_y - param.NumberPixels, pProb, param, RefMap);
 		}
 	}
 }
 
 template <int GPUAlgo, class RefT>
 __device__ static inline void compareRefMap(const int iRefMap, const int iOrient, const int iConv, const bioem_map& Mapconv, bioem_Probability* pProb, const bioem_param_device& param, const RefT& RefMap,
-	const int cent_x, const int cent_y, const int myShift = 0, const int nShifts2 = 0, const int myRef = 0, const bool threadActive = true)
+		const int cent_x, const int cent_y, const int myShift = 0, const int nShifts2 = 0, const int myRef = 0, const bool threadActive = true)
 {
 	/**************************************************************************************/
 	/**********************  Calculating BioEM Probability ********************************/
@@ -145,9 +145,9 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
 	myfloat_t sum, sumsquare, crossproMapConv;
 	__m128 sum_v = _mm_setzero_ps(), sumsquare_v = _mm_setzero_ps(), cross_v = _mm_setzero_ps(), d1, d2;
 #else
-	myfloat_t sum=0.0;
-	myfloat_t sumsquare=0.0;
-	myfloat_t crossproMapConv=0.0;
+	myfloat_t sum = 0.0;
+	myfloat_t sumsquare = 0.0;
+	myfloat_t crossproMapConv = 0.0;
 #endif
 	/****** Loop over Pixels to calculate dot product and cross-correlations of displaced Ref Conv. Map***/
 	myfloat_t logpro;
@@ -182,7 +182,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
 			const float* ptr2 = RefMap.getp(iRefMap, i, jStart);
 			int j;
 			const int count = jEnd - jStart;
-			for (j = 0;j <= count - 4;j += 4)
+			for (j = 0; j <= count - 4; j += 4)
 			{
 				d1 = _mm_loadu_ps(ptr1);
 				d2 = _mm_loadu_ps(ptr2);
@@ -201,7 +201,7 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
 				// Crosscorrelation of calculated displaced map
 				sum += pointMap;
 				// Calculate Sum of pixels squared
-				sumsquare += pointMap*pointMap;
+				sumsquare += pointMap * pointMap;
 			}
 #endif
 		}
@@ -380,18 +380,18 @@ __device__ static inline void compareRefMap(const int iRefMap, const int iOrient
 	else
 #endif
 
-	/***** Summing & Storing total/Orientation Probabilites for each map ************/
+		/***** Summing & Storing total/Orientation Probabilites for each map ************/
 	{
-		update_prob<-1>(logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb);
+		update_prob < -1 > (logpro, iRefMap, iOrient, iConv, cent_x, cent_y, pProb);
 	}
 }
 
 template <int GPUAlgo, class RefT>
 __device__ static inline void compareRefMapShifted(const int iRefMap, const int iOrient, const int iConv, const bioem_map& Mapconv, bioem_Probability* pProb, const bioem_param_device& param, const RefT& RefMap)
 {
-	for (int cent_x = -param.maxDisplaceCenter; cent_x <= param.maxDisplaceCenter; cent_x=cent_x+param.GridSpaceCenter)
+	for (int cent_x = -param.maxDisplaceCenter; cent_x <= param.maxDisplaceCenter; cent_x = cent_x + param.GridSpaceCenter)
 	{
-		for (int cent_y = -param.maxDisplaceCenter; cent_y <= param.maxDisplaceCenter; cent_y=cent_y+param.GridSpaceCenter)
+		for (int cent_y = -param.maxDisplaceCenter; cent_y <= param.maxDisplaceCenter; cent_y = cent_y + param.GridSpaceCenter)
 		{
 			compareRefMap<GPUAlgo>(iRefMap, iOrient, iConv, Mapconv, pProb, param, RefMap, cent_x, cent_y);
 		}
diff --git a/bioem_cuda.cu b/bioem_cuda.cu
index 4e2439754bd4df795d08885728d68ab61623d3fc..e4358cc6242718bfac60810bdedb9de590bca1b1 100644
--- a/bioem_cuda.cu
+++ b/bioem_cuda.cu
@@ -39,7 +39,7 @@ __global__ void compareRefMap_kernel(const int iOrient, const int iConv, const b
 	const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
 	if (iRefMap < maxRef)
 	{
-		compareRefMap<0>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap, cent_x, cent_y);
+		compareRefMap<0>(iRefMap, iOrient, iConv, *pMap, pProb, param, *RefMap, cent_x, cent_y);
 	}
 }
 
@@ -48,7 +48,7 @@ __global__ void compareRefMapShifted_kernel(const int iOrient, const int iConv,
 	const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX;
 	if (iRefMap < maxRef)
 	{
-		compareRefMapShifted<1>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap);
+		compareRefMapShifted<1>(iRefMap, iOrient, iConv, *pMap, pProb, param, *RefMap);
 	}
 }
 
@@ -58,7 +58,7 @@ __global__ void cudaZeroMem(void* ptr, size_t size)
 	int mysize = size / sizeof(int);
 	int myid = myBlockDimX * myBlockIdxX + myThreadIdxX;
 	int mygrid = myBlockDimX * myGridDimX;
-	for (int i = myid;i < mysize;i += mygrid) myptr[i] = 0;
+	for (int i = myid; i < mysize; i += mygrid) myptr[i] = 0;
 }
 
 __global__ void compareRefMapLoopShifts_kernel(const int iOrient, const int iConv, const bioem_map* pMap, bioem_Probability* pProb, const bioem_param_device param, const bioem_RefMap* RefMap, const int blockoffset, const int nShifts, const int nShiftBits, const int maxRef)
@@ -74,7 +74,7 @@ __global__ void compareRefMapLoopShifts_kernel(const int iOrient, const int iCon
 
 	const bool threadActive = myShiftIdx < nShifts && myShiftIdy < nShifts && iRefMap < maxRef;
 
-	compareRefMap<2>(iRefMap,iOrient,iConv,*pMap, pProb, param, *RefMap, cent_x, cent_y, myShift, nShifts * nShifts, myRef, threadActive);
+	compareRefMap<2>(iRefMap, iOrient, iConv, *pMap, pProb, param, *RefMap, cent_x, cent_y, myShift, nShifts * nShifts, myRef, threadActive);
 }
 
 __global__ void multComplexMap(const mycomplex_t* convmap, const mycomplex_t* refmap, mycuComplex_t* out, const int NumberPixelsTotal, const int MapSize, const int NumberMaps, const int Offset)
@@ -82,7 +82,7 @@ __global__ void multComplexMap(const mycomplex_t* convmap, const mycomplex_t* re
 	if (myBlockIdxX >= NumberMaps) return;
 	const mycomplex_t* myin = &refmap[myBlockIdxX * MapSize + Offset];
 	mycuComplex_t* myout = &out[(myBlockIdxX * MapSize)];
-	for(int i = myThreadIdxX;i < NumberPixelsTotal;i += myBlockDimX)
+	for(int i = myThreadIdxX; i < NumberPixelsTotal; i += myBlockDimX)
 	{
 		myout[i].x = convmap[i][0] * myin[i][0] + convmap[i][1] * myin[i][1];
 		myout[i].y = convmap[i][1] * myin[i][0] - convmap[i][0] * myin[i][1];
@@ -126,7 +126,7 @@ int bioem_cuda::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& c
 	if (FFTAlgo)
 	{
 		checkCudaErrors(cudaMemcpyAsync(&pConvMapFFT[(iConv & 1) * param.RefMapSize], localmultFFT, param.RefMapSize * sizeof(mycomplex_t), cudaMemcpyHostToDevice, cudaStream));
-		for (int i = 0;i < maxRef;i += CUDA_FFTS_AT_ONCE)
+		for (int i = 0; i < maxRef; i += CUDA_FFTS_AT_ONCE)
 		{
 			const int num = min(CUDA_FFTS_AT_ONCE, maxRef - i);
 			multComplexMap<<<num, CUDA_THREAD_COUNT, 0, cudaStream>>>(&pConvMapFFT[(iConv & 1) * param.RefMapSize], pRefMapsFFT, pFFTtmp2, param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D, param.RefMapSize, num, i);
@@ -164,24 +164,24 @@ int bioem_cuda::compareRefMaps(int iProjectionOut, int iConv, const bioem_map& c
 			const int nShiftBits = ilog2(nShifts);
 			size_t totalBlocks = divup((size_t) maxRef * (size_t) nShifts * (size_t) nShifts, (size_t) CUDA_THREAD_COUNT);
 			size_t nBlocks = CUDA_BLOCK_COUNT;
-			for (size_t i = 0;i < totalBlocks;i += nBlocks)
+			for (size_t i = 0; i < totalBlocks; i += nBlocks)
 			{
-				compareRefMapLoopShifts_kernel <<<min(nBlocks, totalBlocks - i), CUDA_THREAD_COUNT, (CUDA_THREAD_COUNT * 2 + CUDA_THREAD_COUNT / (nShifts * nShifts) * 4) * sizeof(myfloat_t), cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device, i, nShifts, nShiftBits, maxRef);
+				compareRefMapLoopShifts_kernel<<<min(nBlocks, totalBlocks - i), CUDA_THREAD_COUNT, (CUDA_THREAD_COUNT * 2 + CUDA_THREAD_COUNT / (nShifts * nShifts) * 4) * sizeof(myfloat_t), cudaStream >>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device, i, nShifts, nShiftBits, maxRef);
 			}
 		}
 		else if (GPUAlgo == 1) //Split shifts in multiple kernels
 		{
-			for (int cent_x = -param.param_device.maxDisplaceCenter; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter)
+			for (int cent_x = -param.param_device.maxDisplaceCenter; cent_x <= param.param_device.maxDisplaceCenter; cent_x = cent_x + param.param_device.GridSpaceCenter)
 			{
-				for (int cent_y = -param.param_device.maxDisplaceCenter; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
+				for (int cent_y = -param.param_device.maxDisplaceCenter; cent_y <= param.param_device.maxDisplaceCenter; cent_y = cent_y + param.param_device.GridSpaceCenter)
 				{
-					compareRefMap_kernel <<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, cent_x, cent_y, maxRef);
+					compareRefMap_kernel<<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, cent_x, cent_y, maxRef);
 				}
 			}
 		}
 		else if (GPUAlgo == 0) //All shifts in one kernel
 		{
-			compareRefMapShifted_kernel <<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, maxRef);
+			compareRefMapShifted_kernel<<<divup(maxRef, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream>>> (iProjectionOut, iConv, pConvMap_device[iConv & 1], pProb_device, param.param_device, pRefMap_device_Mod, maxRef);
 		}
 		else
 		{
@@ -215,7 +215,7 @@ int bioem_cuda::deviceInit()
 	checkCudaErrors(cudaMalloc(&pRefMap_device, sizeof(bioem_RefMap)));
 	cout << "\tSize Probability\t" << sizeof(bioem_Probability) * RefMap.ntotRefMap << "\n";
 	checkCudaErrors(cudaMalloc(&pProb_device, sizeof(bioem_Probability) * RefMap.ntotRefMap));
-	for (int i = 0;i < 2;i++)
+	for (int i = 0; i < 2; i++)
 	{
 		checkCudaErrors(cudaEventCreate(&cudaEvent[i]));
 		checkCudaErrors(cudaMalloc(&pConvMap_device[i], sizeof(bioem_map)));
@@ -256,7 +256,7 @@ int bioem_cuda::deviceExit()
 	cudaStreamDestroy(cudaStream);
 	cudaFree(pRefMap_device);
 	cudaFree(pProb_device);
-	for (int i = 0;i < 2;i++)
+	for (int i = 0; i < 2; i++)
 	{
 		cudaEventDestroy(cudaEvent[i]);
 		cudaFree(pConvMap_device);
@@ -282,7 +282,7 @@ int bioem_cuda::deviceStartRun()
 
 	if (FFTAlgo)
 	{
-		for (int i = 0;i < 2;i++)
+		for (int i = 0; i < 2; i++)
 		{
 			int n[2] = {param.param_device.NumberPixels, param.param_device.NumberPixels};
 			if (cufftPlanMany(&plan[i], 2, n, NULL, 1, 0, NULL, 1, 0, CUFFT_C2R, i ? (maxRef % CUDA_FFTS_AT_ONCE) : CUDA_FFTS_AT_ONCE) != CUFFT_SUCCESS)
@@ -312,7 +312,7 @@ int bioem_cuda::deviceFinishRun()
 
 	if (FFTAlgo)
 	{
-		for (int i = 0;i < 2;i++) cufftDestroy(plan[i]);
+		for (int i = 0; i < 2; i++) cufftDestroy(plan[i]);
 	}
 
 	return(0);
diff --git a/include/bioem.h b/include/bioem.h
index 18929476a204a9475a31662d403a56f83cfdab82..b4fc9a60ff7bc608891052b239c387ea735f83e3 100644
--- a/include/bioem.h
+++ b/include/bioem.h
@@ -21,13 +21,13 @@ public:
 	int dopreCalCrossCorrelation(int iRefMap, int iRefMapLocal);
 	int run();
 	int doProjections(int iMap);
-	int createConvolutedProjectionMap(int iOreint,int iMap, mycomplex_t* lproj,bioem_map& Mapconv,mycomplex_t* localmultFFT,myfloat_t& sumC,myfloat_t& sumsquareC);
+	int createConvolutedProjectionMap(int iOreint, int iMap, mycomplex_t* lproj, bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC);
 
 	virtual int compareRefMaps(int iProjectionOut, int iConv, const bioem_map& conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap = 0);
 
 	int createProjection(int iMap, mycomplex_t* map);
-	int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare);
-	void calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,myfloat_t* lCC);
+	int calcross_cor(bioem_map& localmap, myfloat_t& sum, myfloat_t& sumsquare);
+	void calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC);
 
 	bioem_Probability* pProb;
 
diff --git a/include/map.h b/include/map.h
index 5af0d9883aad073188362b172b7cc632e8881b4e..5bb761f189deb452a4f041782dc92839c90f5922 100644
--- a/include/map.h
+++ b/include/map.h
@@ -22,13 +22,13 @@ public:
 	}
 	int readRefMaps(bioem_param& param);
 	int PreCalculateMapsFFT(bioem_param& param);
-        int  read_int(int *currlong, FILE *fin, int swap);
-        int  read_float(float *currfloat, FILE *fin, int swap);
-        int  read_float_empty (FILE *fin);
-        int  read_char_float (float *currfloat, FILE *fin) ;
-        int  test_mrc (const char *vol_file, int swap);
-	int  read_MRC(const char* filename,bioem_param& param);
-	
+	int  read_int(int *currlong, FILE *fin, int swap);
+	int  read_float(float *currfloat, FILE *fin, int swap);
+	int  read_float_empty (FILE *fin);
+	int  read_char_float (float *currfloat, FILE *fin) ;
+	int  test_mrc (const char *vol_file, int swap);
+	int  read_MRC(const char* filename, bioem_param& param);
+
 	mycomplex_t* RefMapsFFT;
 
 	const char* filemap;
@@ -38,7 +38,7 @@ public:
 	myfloat_t sumsquare_RefMap[BIOEM_MAX_MAPS];
 	myfloat_t ForLogProbfromRef[BIOEM_MAX_MAPS];
 
-	bool dumpMap, loadMap, readMRC,readMultMRC;
+	bool dumpMap, loadMap, readMRC, readMultMRC;
 
 	__host__ __device__ inline myfloat_t get(int map, int x, int y) const {return(Ref[map].points[x][y]);}
 	__host__ __device__ inline const myfloat_t* getp(int map, int x, int y) const {return(&Ref[map].points[x][y]);}
@@ -65,12 +65,12 @@ public:
 		memcpy(sum_RefMap, map.sum_RefMap, sizeof(sum_RefMap));
 		memcpy(sumsquare_RefMap, map.sumsquare_RefMap, sizeof(sumsquare_RefMap));
 		memcpy(ForLogProbfromRef, map.ForLogProbfromRef, sizeof(ForLogProbfromRef));
-#pragma omp parallel for
-		for (int i = 0;i < ntotRefMap;i++)
+		#pragma omp parallel for
+		for (int i = 0; i < ntotRefMap; i++)
 		{
-			for (int j = 0;j < BIOEM_MAP_SIZE_X;j++)
+			for (int j = 0; j < BIOEM_MAP_SIZE_X; j++)
 			{
-				for (int k = 0;k < BIOEM_MAP_SIZE_Y;k++)
+				for (int k = 0; k < BIOEM_MAP_SIZE_Y; k++)
 				{
 					Ref[j][k][i] = map.get(i, j, k);
 				}
diff --git a/include/model.h b/include/model.h
index bc08678e5d5d613e76dd9d7819c7a27efd386d48..f325e148ca5654f9bc45b75c69a5423e75861f0b 100644
--- a/include/model.h
+++ b/include/model.h
@@ -15,7 +15,7 @@ public:
 
 	myfloat_t getAminoAcidRad(char *name);
 	myfloat_t getAminoAcidDensity(char *name);
-	myfloat_t NormDen;
+	myfloat_t NormDen;
 
 	const char* filemodel;
 
diff --git a/main.cpp b/main.cpp
index 555d974e8f2c9b1550dfa45469d88899b73c8728..09877445faae16265a7150067f639f5461a35fb2 100644
--- a/main.cpp
+++ b/main.cpp
@@ -26,7 +26,7 @@ int main(int argc, char* argv[])
 	//*************************************************************************************
 
 #ifdef _MM_SET_DENORMALS_ZERO_MODE
-#pragma omp parallel
+	#pragma omp parallel
 	{
 		_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); //Flush denormals to zero in all OpenMP threads
 	}
@@ -47,7 +47,7 @@ int main(int argc, char* argv[])
 
 	//************  Configuration and Pre-calculating necessary objects *********************
 	printf("Configuring\n");
-	bio->configure(argc,argv);
+	bio->configure(argc, argv);
 
 	//*******************************  Run BioEM routine ************************************
 	printf("Running\n");
diff --git a/map.cpp b/map.cpp
index 488173b9ab3c4ec0f0b6ed16ec02c33917a8d6fd..dd13060e676312e488ebf6c78301dd8f7cecb242 100644
--- a/map.cpp
+++ b/map.cpp
@@ -13,225 +13,225 @@ using namespace std;
 
 int bioem_RefMap::readRefMaps(bioem_param& param)
 {
-  //**************************************************************************************
-  //***********************Reading reference Particle Maps************************
-  //**************************************************************************************
-  if (loadMap)
-    {
-      FILE* fp = fopen("maps.dump", "rb");
-      if (fp == NULL)
+	//**************************************************************************************
+	//***********************Reading reference Particle Maps************************
+	//**************************************************************************************
+	if (loadMap)
 	{
-	  cout << "Error opening dump file\n";
-	  exit(1);
+		FILE* fp = fopen("maps.dump", "rb");
+		if (fp == NULL)
+		{
+			cout << "Error opening dump file\n";
+			exit(1);
+		}
+		fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
+		if (ntotRefMap > BIOEM_MAX_MAPS)
+		{
+			cout << "BIOEM_MAX_MAPS too small, some maps dropped\n";
+			ntotRefMap = BIOEM_MAX_MAPS;
+		}
+		fread(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
+		fclose(fp);
+
+		cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
+		cout << "\n+++++++++++++++++++++++++++++++++++++++ \n";
 	}
-      fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
-      if (ntotRefMap > BIOEM_MAX_MAPS)
+	else
 	{
-	  cout << "BIOEM_MAX_MAPS too small, some maps dropped\n";
-	  ntotRefMap = BIOEM_MAX_MAPS;
-	}
-      fread(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
-      fclose(fp);
-
-      cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
-      cout << "\n+++++++++++++++++++++++++++++++++++++++ \n";
-    }
-  else
-    {
-
-      if(readMRC){
-	
-             ntotRefMap=0;
-	     
-	      if(readMultMRC)
-	      {
-		
-	      ifstream input(filemap);
-	      
-	      if (!input.good())
-	      {
-		cout << "Failed to open file contaning MRC names: " << filemap << "\n";
-		exit(0);
-	      }
-	      
-	      char line[512] = {0};
-	      char mapname[100];
-	      char tmpm[10] = {0};
-              const char* indifile;
-	      
-	      while (!input.eof())
-	      {
-		input.getline(line,512);
-		char tmpVals[100]  = {' '};
-
-		strncpy (tmpVals,line,100);
-		sscanf (tmpVals,"%99c",mapname);
-		
-		// Check for last line
-		strncpy (tmpm,mapname,3);		
-	        if(strcmp(tmpm,"XXX")!=0){
-		indifile = mapname;  
-		//Reading Multiple MRC
-		read_MRC(indifile,param);
-		
+
+		if(readMRC) {
+
+			ntotRefMap = 0;
+
+			if(readMultMRC)
+			{
+
+				ifstream input(filemap);
+
+				if (!input.good())
+				{
+					cout << "Failed to open file contaning MRC names: " << filemap << "\n";
+					exit(0);
+				}
+
+				char line[512] = {0};
+				char mapname[100];
+				char tmpm[10] = {0};
+				const char* indifile;
+
+				while (!input.eof())
+				{
+					input.getline(line, 512);
+					char tmpVals[100]  = {' '};
+
+					strncpy (tmpVals, line, 100);
+					sscanf (tmpVals, "%99c", mapname);
+
+					// Check for last line
+					strncpy (tmpm, mapname, 3);
+					if(strcmp(tmpm, "XXX") != 0) {
+						indifile = mapname;
+						//Reading Multiple MRC
+						read_MRC(indifile, param);
+
+					}
+					for(int i = 0; i < 3; i++)mapname[i] = 'X';
+					for(int i = 3; i < 100; i++)mapname[i] = {0};
+
+				}
+				cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
+				cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ;
+
+			} else {
+
+				read_MRC(filemap, param);
+				cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
+				cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ;
+
+			}
+
+		}
+		else {
+
+			int nummap = -1;
+			int lasti = 0;
+			int lastj = 0;
+			ifstream input(filemap);
+			if (!input.good())
+			{
+				cout << "Particle Maps Failed to open file" << endl ;
+				exit(1);
+			}
+
+			char line[512] = {' '};
+			char tmpLine[512] = {' '};
+
+			while (!input.eof())
+			{
+				input.getline(line, 512);
+
+				strncpy(tmpLine, line, strlen(line));
+				char *token = strtok(tmpLine, " ");
+
+				if (strcmp(token, "PARTICLE") == 0) // to count the number of maps
+				{
+					nummap++;
+					if (nummap % 128 == 0)
+					{
+						cout << "..." << nummap << "\n";
+					}
+					if (nummap == BIOEM_MAX_MAPS)
+					{
+						cout << "BIOEM_MAX_MAPS too small\n";
+						exit(1);
+					}
+					if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels && nummap > 0)
+					{
+						cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
+						exit(1);
+					}
+				}
+				else
+				{
+					int i, j;
+					float z;
+
+					char tmpVals[36]  = {' '};
+
+					strncpy (tmpVals, line, 8);
+					sscanf (tmpVals, "%d", &i);
+
+					strncpy (tmpVals, line + 8, 8);
+					sscanf (tmpVals, "%d", &j);
+
+					strncpy (tmpVals, line + 16, 12);
+					sscanf (tmpVals, "%f", &z);
+					//checking for Map limits
+					if(i > 0 && i - 1 < BIOEM_MAP_SIZE_X && j > 0 && j - 1 < BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
+					{
+						Ref[nummap].points[i - 1][j - 1] = (myfloat_t) z;
+						lasti = i;
+						lastj = j;
+					}
+					else
+					{
+						cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
+						exit(1);
+					}
+				}
+			}
+			cout << ".";
+			ntotRefMap = nummap + 1;
+			cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
+			cout << "Particle Maps read from Standard File\n" ;
+		}
+		cout << "Total Number of particles: " << ntotRefMap ;
+		cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
+
+		if (dumpMap)
+		{
+			FILE* fp = fopen("maps.dump", "w+b");
+			if (fp == NULL)
+			{
+				cout << "Error opening dump file\n";
+				exit(1);
+			}
+			fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
+			fwrite(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
+			fclose(fp);
 		}
-		for(int i=0;i<3;i++)mapname[i] = 'X';	
-		for(int i=3;i<100;i++)mapname[i] = {0};	
-		
-	      }
-	      cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
-	      cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ;
-	      
-	      } else {
-	      
-		read_MRC(filemap,param);
-		cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
-		cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ;
-		
-	      }
-      
-      }
-      else {
-
-	int nummap=-1;
-	int lasti=0;
-	int lastj=0;
-	ifstream input(filemap);
-	if (!input.good())
-	  {
-	    cout << "Particle Maps Failed to open file" << endl ;
-	    exit(1);
-	  }
-
-	char line[512] = {' '};
-	char tmpLine[512] = {' '};
-
-	while (!input.eof())
-	  {
-	    input.getline(line,512);
-
-	    strncpy(tmpLine,line,strlen(line));
-	    char *token = strtok(tmpLine," ");
-
-	    if (strcmp(token,"PARTICLE")==0) // to count the number of maps
-	      {
-		nummap++;
-		if (nummap % 128 == 0)
-		  {
-		    cout << "..." << nummap << "\n";
-		  }
-		if (nummap == BIOEM_MAX_MAPS)
-		  {
-		    cout << "BIOEM_MAX_MAPS too small\n";
-		    exit(1);
-		  }
-		if(lasti!=param.param_device.NumberPixels && lastj!=param.param_device.NumberPixels && nummap>0)
-		  {
-		    cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
-		    exit(1);
-		  }
-	      }
-	    else
-	      {
-		int i,j;
-		float z;
-
-		char tmpVals[36]  = {' '};
-
-		strncpy (tmpVals,line,8);
-		sscanf (tmpVals,"%d",&i);
-
-		strncpy (tmpVals,line+8,8);
-		sscanf (tmpVals,"%d",&j);
-
-		strncpy (tmpVals,line+16,12);
-		sscanf (tmpVals,"%f",&z);
-		//checking for Map limits
-		if(i>0 && i-1 <BIOEM_MAP_SIZE_X&& j>0 && j-1<BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
-		  {
-		    Ref[nummap].points[i-1][j-1] = (myfloat_t) z;
-		    lasti=i;
-		    lastj=j;
-		  }
-		else
-		  {
-		    cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
-		    exit(1);
-		  }
-	      }
-	  }
-	cout << ".";
-	ntotRefMap=nummap+1;
-	cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
-	cout << "Particle Maps read from Standard File\n" ;	
-      }
-        cout << "Total Number of particles: " << ntotRefMap ;
-	cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
-	
-      if (dumpMap)
-	{
-	  FILE* fp = fopen("maps.dump", "w+b");
-	  if (fp == NULL)
-	    {
-	      cout << "Error opening dump file\n";
-	      exit(1);
-	    }
-	  fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
-	  fwrite(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
-	  fclose(fp);
 	}
-    }
 
-  return(0);
+	return(0);
 }
 
 int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
 {
-  //**************************************************************************************
-  //********** Routine that pre-calculates Kernels for Convolution ***********************
-  //**************************************************************************************
+	//**************************************************************************************
+	//********** Routine that pre-calculates Kernels for Convolution ***********************
+	//**************************************************************************************
 
-  myfloat_t* localMap;
+	myfloat_t* localMap;
 
-  localMap= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
+	localMap = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
 
-  RefMapsFFT = new mycomplex_t[ntotRefMap * param.RefMapSize];
+	RefMapsFFT = new mycomplex_t[ntotRefMap * param.RefMapSize];
 
-  mycomplex_t* localout;
-  localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
+	mycomplex_t* localout;
+	localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
 
-  for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
-    {
-      //Assigning localMap values to padded Map
-      for(int i=0; i< param.param_device.NumberPixels; i++)
+	for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
 	{
-	  for(int j=0; j< param.param_device.NumberPixels; j++)
-	    {
-	      localMap[i*param.param_device.NumberPixels+j]=Ref[iRefMap].points[i][j];
-	    }
-	}
+		//Assigning localMap values to padded Map
+		for(int i = 0; i < param.param_device.NumberPixels; i++)
+		{
+			for(int j = 0; j < param.param_device.NumberPixels; j++)
+			{
+				localMap[i * param.param_device.NumberPixels + j] = Ref[iRefMap].points[i][j];
+			}
+		}
 
-      //Calling FFT_Forward
-      myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localMap,localout);
+		//Calling FFT_Forward
+		myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localMap, localout);
 
-      // Normalizing and saving the Reference CTFs
-      mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.RefMapSize];
-      for(int i=0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D ; i++ )
-	{
-	  RefMap[i][0]=localout[i][0];
-	  RefMap[i][1]=localout[i][1];
+		// Normalizing and saving the Reference CTFs
+		mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.RefMapSize];
+		for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D ; i++ )
+		{
+			RefMap[i][0] = localout[i][0];
+			RefMap[i][1] = localout[i][1];
+		}
 	}
-    }
 
-  myfftw_free(localMap);
-  myfftw_free(localout);
+	myfftw_free(localMap);
+	myfftw_free(localout);
 
-  return(0);
+	return(0);
 }
 
-int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
-  {
-        myfloat_t st,st2; 
+int bioem_RefMap::read_MRC(const char* filename, bioem_param& param)
+{
+	myfloat_t st, st2;
 	unsigned long count;
 	FILE *fin;
 	float currfloat;
@@ -239,63 +239,63 @@ int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
 	float xlen, ylen, zlen;
 	int mode, ncstart, nrstart, nsstart, ispg, nsymbt, lskflg;
 	float a_tmp, b_tmp, g_tmp;
-	int  mx, my, mz,mapc, mapr, maps;
+	int  mx, my, mz, mapc, mapr, maps;
 	float dmin, dmax, dmean;
-        int n_range_viol0, n_range_viol1;
+	int n_range_viol0, n_range_viol1;
 
 	fin = fopen(filename, "rb");
 	if( fin == NULL ) {
-	  cout << "ERROR opening MRC: " << filename;
-	  exit(1);
+		cout << "ERROR opening MRC: " << filename;
+		exit(1);
 	}
-	n_range_viol0 = test_mrc(filename,0);
-	 n_range_viol1 = test_mrc(filename,1);
-
-
-  if (n_range_viol0 < n_range_viol1) { //* guess endianism 
-    swap = 0;
-    if (n_range_viol0 > 0) {
-      printf(" Warning: %i header field range violations detected in file %s \n", n_range_viol0,filename);
-    }
-  } else {
-    swap = 1;
-    if (n_range_viol1 > 0) {
-      printf("Warning: %i header field range violations detected in file %s \n", n_range_viol1,filename);
-    }
-  }
-       printf("\n+++++++++++++++++++++++++++++++++++++++++++\n");
+	n_range_viol0 = test_mrc(filename, 0);
+	n_range_viol1 = test_mrc(filename, 1);
+
+
+	if (n_range_viol0 < n_range_viol1) { //* guess endianism
+		swap = 0;
+		if (n_range_viol0 > 0) {
+			printf(" Warning: %i header field range violations detected in file %s \n", n_range_viol0, filename);
+		}
+	} else {
+		swap = 1;
+		if (n_range_viol1 > 0) {
+			printf("Warning: %i header field range violations detected in file %s \n", n_range_viol1, filename);
+		}
+	}
+	printf("\n+++++++++++++++++++++++++++++++++++++++++++\n");
 	printf("Reading Information from MRC: %s \n", filename);
-	header_ok *= read_int(&nc,fin,swap);
-	header_ok *= read_int(&nr,fin,swap);
-	header_ok *= read_int(&ns,fin,swap);
-	header_ok *= read_int(&mode,fin,swap);
-	header_ok *= read_int(&ncstart,fin,swap);
-	header_ok *= read_int(&nrstart,fin,swap);
-	header_ok *= read_int(&nsstart,fin,swap);
-	header_ok *= read_int(&mx,fin,swap);
-	header_ok *= read_int(&my,fin,swap);
-	header_ok *= read_int(&mz,fin,swap);
-	header_ok *= read_float(&xlen,fin,swap);
-	header_ok *= read_float(&ylen,fin,swap);
-	header_ok *= read_float(&zlen,fin,swap);
-	header_ok *= read_float(&a_tmp,fin,swap);
-	header_ok *= read_float(&b_tmp,fin,swap);
-	header_ok *= read_float(&g_tmp,fin,swap);
-	header_ok *= read_int(&mapc,fin,swap);
-	header_ok *= read_int(&mapr,fin,swap);
-	header_ok *= read_int(&maps,fin,swap);
-	header_ok *= read_float(&dmin,fin,swap);
-	header_ok *= read_float(&dmax,fin,swap);
-	header_ok *= read_float(&dmean,fin,swap);
-	header_ok *= read_int(&ispg,fin,swap);
-	header_ok *= read_int(&nsymbt,fin,swap);
-	header_ok *= read_int(&lskflg,fin,swap);
-
-	printf("Number Columns  = %8d \n",nc);
-	printf("Number Rows     = %8d \n",nr);
-	printf("Number Sections = %8d \n",ns);
-	printf("MODE = %4d (only data type mode 2: 32-bit)\n",mode);
-       printf("NSYMBT = %4d (# bytes symmetry operators)\n",nsymbt);
+	header_ok *= read_int(&nc, fin, swap);
+	header_ok *= read_int(&nr, fin, swap);
+	header_ok *= read_int(&ns, fin, swap);
+	header_ok *= read_int(&mode, fin, swap);
+	header_ok *= read_int(&ncstart, fin, swap);
+	header_ok *= read_int(&nrstart, fin, swap);
+	header_ok *= read_int(&nsstart, fin, swap);
+	header_ok *= read_int(&mx, fin, swap);
+	header_ok *= read_int(&my, fin, swap);
+	header_ok *= read_int(&mz, fin, swap);
+	header_ok *= read_float(&xlen, fin, swap);
+	header_ok *= read_float(&ylen, fin, swap);
+	header_ok *= read_float(&zlen, fin, swap);
+	header_ok *= read_float(&a_tmp, fin, swap);
+	header_ok *= read_float(&b_tmp, fin, swap);
+	header_ok *= read_float(&g_tmp, fin, swap);
+	header_ok *= read_int(&mapc, fin, swap);
+	header_ok *= read_int(&mapr, fin, swap);
+	header_ok *= read_int(&maps, fin, swap);
+	header_ok *= read_float(&dmin, fin, swap);
+	header_ok *= read_float(&dmax, fin, swap);
+	header_ok *= read_float(&dmean, fin, swap);
+	header_ok *= read_int(&ispg, fin, swap);
+	header_ok *= read_int(&nsymbt, fin, swap);
+	header_ok *= read_int(&lskflg, fin, swap);
+
+	printf("Number Columns  = %8d \n", nc);
+	printf("Number Rows     = %8d \n", nr);
+	printf("Number Sections = %8d \n", ns);
+	printf("MODE = %4d (only data type mode 2: 32-bit)\n", mode);
+	printf("NSYMBT = %4d (# bytes symmetry operators)\n", nsymbt);
 
 	/* printf("  NCSTART = %8d  (index of first column, counting from 0)\n",ncstart);
 	   printf(">  NRSTART = %8d  (index of first row, counting from 0)\n",nrstart);
@@ -320,179 +320,179 @@ int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
 	   printf("   LSKFLG = %8d  (skew matrix flag: 0:none, 1:follows)\n",lskflg);*/
 
 	if (header_ok == 0) {
-	  cout << "ERROR reading MRC header: " << filename;
-	  exit(1);
+		cout << "ERROR reading MRC header: " << filename;
+		exit(1);
 	}
 
-	if(nr!=param.param_device.NumberPixels || nc!=param.param_device.NumberPixels )
-	  {
-	    cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << nc << ", j " << nr << ")" << "\n";
-	    exit(1);
-	  }
-
+	if(nr != param.param_device.NumberPixels || nc != param.param_device.NumberPixels )
+	{
+		cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << nc << ", j " << nr << ")" << "\n";
+		exit(1);
+	}
 
-	if(mode!=2)
-	  {
-	    cout << "ERROR with MRC mode " << mode << "\n";
-	    cout << "Currently mode 2 is the only one allowed" << "\n";
-	    exit(1);
 
-	  } else
-	  {  
-               rewind (fin);
-	    for (count=0; count<256; ++count) if (read_float_empty(fin)==0) {
-		cout << "ERROR Converting Data: " <<  filename;
+	if(mode != 2)
+	{
+		cout << "ERROR with MRC mode " << mode << "\n";
+		cout << "Currently mode 2 is the only one allowed" << "\n";
 		exit(1);
-	      }
-
-	    for (count=0; count<(unsigned long)nsymbt; ++count) if (read_char_float(&currfloat,fin)==0) {
-		cout << "ERROR Converting Data: " <<  filename;
-		exit(1);    
-	      }
-      
-	    for ( int nmap = 0 ; nmap < ns ; nmap ++ )
-              {
-        
-               st=0.0;
-               st2=0.0;
-	      for ( int j = 0 ; j < nr ; j ++ )
-		for ( int i = 0 ; i < nc ; i ++ )
-		  {
-		    if (read_float(&currfloat,fin,swap)==0) {
-		      cout << "ERROR Converting Data: " <<  filename;
-		      exit(1);
-		    } else {
-                      Ref[nmap+ntotRefMap].points[i][j] = currfloat;
-                      st += currfloat;
-                      st2 += currfloat*currfloat;
-	
-		    }
-		  }	
-
-            //Normaling maps to zero mean and unit standard deviation
-                    st /= float(nr*nc);
-                    st2 = sqrt(st2/float(nr*nc)-st*st); 
-                    for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){
-                          Ref[nmap+ntotRefMap].points[i][j] =  Ref[nmap+ntotRefMap].points[i][j]/st2 - st/st2;
-  	  //    if(nmap+ntotRefMap==300)cout << i << " " << j << " " << nmap+ntotRefMap << " " <<  Ref[nmap+ntotRefMap].points[i][j]  << "\n";                          
-                          }
-              }
-           ntotRefMap += ns ;
-	 //  cout << ntotRefMap << "\n";
-	  }
+
+	} else
+	{
+		rewind (fin);
+		for (count = 0; count < 256; ++count) if (read_float_empty(fin) == 0) {
+				cout << "ERROR Converting Data: " <<  filename;
+				exit(1);
+			}
+
+		for (count = 0; count < (unsigned long)nsymbt; ++count) if (read_char_float(&currfloat, fin) == 0) {
+				cout << "ERROR Converting Data: " <<  filename;
+				exit(1);
+			}
+
+		for ( int nmap = 0 ; nmap < ns ; nmap ++ )
+		{
+
+			st = 0.0;
+			st2 = 0.0;
+			for ( int j = 0 ; j < nr ; j ++ )
+				for ( int i = 0 ; i < nc ; i ++ )
+				{
+					if (read_float(&currfloat, fin, swap) == 0) {
+						cout << "ERROR Converting Data: " <<  filename;
+						exit(1);
+					} else {
+						Ref[nmap + ntotRefMap].points[i][j] = currfloat;
+						st += currfloat;
+						st2 += currfloat * currfloat;
+
+					}
+				}
+
+			//Normaling maps to zero mean and unit standard deviation
+			st /= float(nr * nc);
+			st2 = sqrt(st2 / float(nr * nc) - st * st);
+			for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ) {
+					Ref[nmap + ntotRefMap].points[i][j] =  Ref[nmap + ntotRefMap].points[i][j] / st2 - st / st2;
+					//    if(nmap+ntotRefMap==300)cout << i << " " << j << " " << nmap+ntotRefMap << " " <<  Ref[nmap+ntotRefMap].points[i][j]  << "\n";
+				}
+		}
+		ntotRefMap += ns ;
+		//  cout << ntotRefMap << "\n";
+	}
 	fclose (fin);
-	
+
 	return(0);
-  }
-  
+}
+
 int bioem_RefMap::read_float(float *currfloat, FILE *fin, int swap) {
-  unsigned char *cptr, tmp;
-
-  if (fread(currfloat,4,1,fin)!=1) return 0;
-  if (swap == 1) {
-    cptr = (unsigned char *)currfloat;
-    tmp = cptr[0];
-    cptr[0]=cptr[3];
-    cptr[3]=tmp;
-    tmp = cptr[1];
-    cptr[1]=cptr[2];
-    cptr[2]=tmp;
-  }
-  return 1;
+	unsigned char *cptr, tmp;
+
+	if (fread(currfloat, 4, 1, fin) != 1) return 0;
+	if (swap == 1) {
+		cptr = (unsigned char *)currfloat;
+		tmp = cptr[0];
+		cptr[0] = cptr[3];
+		cptr[3] = tmp;
+		tmp = cptr[1];
+		cptr[1] = cptr[2];
+		cptr[2] = tmp;
+	}
+	return 1;
 }
 
 int  bioem_RefMap::read_int(int *currlong, FILE *fin, int swap) {
-  unsigned char *cptr, tmp;
-
-  if (fread(currlong,4,1,fin)!=1) return 0;
-  if (swap == 1) {
-    cptr = (unsigned char *)currlong;
-    tmp = cptr[0];
-    cptr[0]=cptr[3];
-    cptr[3]=tmp;
-    tmp = cptr[1];
-    cptr[1]=cptr[2];
-    cptr[2]=tmp;
-  }
-  return 1;
+	unsigned char *cptr, tmp;
+
+	if (fread(currlong, 4, 1, fin) != 1) return 0;
+	if (swap == 1) {
+		cptr = (unsigned char *)currlong;
+		tmp = cptr[0];
+		cptr[0] = cptr[3];
+		cptr[3] = tmp;
+		tmp = cptr[1];
+		cptr[1] = cptr[2];
+		cptr[2] = tmp;
+	}
+	return 1;
 }
 int  bioem_RefMap::read_float_empty (FILE *fin) {
-  float currfloat;
+	float currfloat;
 
-  if (fread(&currfloat,4,1,fin)!=1) return 0;
-  return 1;
+	if (fread(&currfloat, 4, 1, fin) != 1) return 0;
+	return 1;
 }
 
 int  bioem_RefMap::read_char_float (float *currfloat, FILE *fin) {
-  char currchar;
+	char currchar;
 
-  if (fread(&currchar,1,1,fin)!=1) return 0;
-  *currfloat=(float)currchar;
-  return 1;
+	if (fread(&currchar, 1, 1, fin) != 1) return 0;
+	*currfloat = (float)currchar;
+	return 1;
 }
 
 int  bioem_RefMap::test_mrc (const char *vol_file, int swap) {
-  FILE *fin;
-  int nc, nr, ns, mx, my, mz;
-  int mode, ncstart, nrstart, nsstart;
-  float xlen, ylen, zlen;
-  int i, header_ok = 1, n_range_viols = 0;
-  int mapc, mapr, maps;
-  float alpha, beta, gamma;
-  float dmin, dmax, dmean, dummy, xorigin, yorigin, zorigin;
-
-  fin = fopen(vol_file, "rb");
-  if( fin == NULL ) {
-  	 cout << "ERROR opening MRC: " << vol_file;
-          exit(1);
-  }
-
- //* read header info 
-  header_ok *= read_int(&nc,fin,swap);
-  header_ok *= read_int(&nr,fin,swap);
-  header_ok *= read_int(&ns,fin,swap);
-  header_ok *= read_int(&mode,fin,swap);
-  header_ok *= read_int(&ncstart,fin,swap);
-  header_ok *= read_int(&nrstart,fin,swap);
-  header_ok *= read_int(&nsstart,fin,swap);
-  header_ok *= read_int(&mx,fin,swap);
-  header_ok *= read_int(&my,fin,swap);
-  header_ok *= read_int(&mz,fin,swap);
-  header_ok *= read_float(&xlen,fin,swap);
-  header_ok *= read_float(&ylen,fin,swap);
-  header_ok *= read_float(&zlen,fin,swap);
-  header_ok *= read_float(&alpha,fin,swap);
-  header_ok *= read_float(&beta,fin,swap);
-  header_ok *= read_float(&gamma,fin,swap);
-  header_ok *= read_int(&mapc,fin,swap);
-  header_ok *= read_int(&mapr,fin,swap);
-  header_ok *= read_int(&maps,fin,swap);
-  header_ok *= read_float(&dmin,fin,swap);
-  header_ok *= read_float(&dmax,fin,swap);
-  header_ok *= read_float(&dmean,fin,swap);
-  for (i=23; i<50; ++i) header_ok *= read_float(&dummy,fin,swap);
-  header_ok *= read_float(&xorigin,fin,swap);
-  header_ok *= read_float(&yorigin,fin,swap);
-  header_ok *= read_float(&zorigin,fin,swap);
-  fclose (fin);
-  if (header_ok == 0) {
- 	cout << "ERROR reading MRC header: " << vol_file;
-          exit(1);
-  }
-
-  n_range_viols += (nc>5000); n_range_viols += (nc<0);
-  n_range_viols += (nr>5000); n_range_viols += (nr<0);
-  n_range_viols += (ns>5000); n_range_viols += (ns<0);
-  n_range_viols += (ncstart>5000); n_range_viols += (ncstart<-5000);
-  n_range_viols += (nrstart>5000); n_range_viols += (nrstart<-5000);
-  n_range_viols += (nsstart>5000); n_range_viols += (nsstart<-5000);
-  n_range_viols += (mx>5000); n_range_viols += (mx<0);
-  n_range_viols += (my>5000); n_range_viols += (my<0);
-  n_range_viols += (mz>5000); n_range_viols += (mz<0);
-  n_range_viols += (alpha>360.0f); n_range_viols += (alpha<-360.0f);
-  n_range_viols += (beta>360.0f); n_range_viols += (beta<-360.0f);
-  n_range_viols += (gamma>360.0f); n_range_viols += (gamma<-360.0f);
-
-  return n_range_viols;
+	FILE *fin;
+	int nc, nr, ns, mx, my, mz;
+	int mode, ncstart, nrstart, nsstart;
+	float xlen, ylen, zlen;
+	int i, header_ok = 1, n_range_viols = 0;
+	int mapc, mapr, maps;
+	float alpha, beta, gamma;
+	float dmin, dmax, dmean, dummy, xorigin, yorigin, zorigin;
+
+	fin = fopen(vol_file, "rb");
+	if( fin == NULL ) {
+		cout << "ERROR opening MRC: " << vol_file;
+		exit(1);
+	}
+
+//* read header info
+	header_ok *= read_int(&nc, fin, swap);
+	header_ok *= read_int(&nr, fin, swap);
+	header_ok *= read_int(&ns, fin, swap);
+	header_ok *= read_int(&mode, fin, swap);
+	header_ok *= read_int(&ncstart, fin, swap);
+	header_ok *= read_int(&nrstart, fin, swap);
+	header_ok *= read_int(&nsstart, fin, swap);
+	header_ok *= read_int(&mx, fin, swap);
+	header_ok *= read_int(&my, fin, swap);
+	header_ok *= read_int(&mz, fin, swap);
+	header_ok *= read_float(&xlen, fin, swap);
+	header_ok *= read_float(&ylen, fin, swap);
+	header_ok *= read_float(&zlen, fin, swap);
+	header_ok *= read_float(&alpha, fin, swap);
+	header_ok *= read_float(&beta, fin, swap);
+	header_ok *= read_float(&gamma, fin, swap);
+	header_ok *= read_int(&mapc, fin, swap);
+	header_ok *= read_int(&mapr, fin, swap);
+	header_ok *= read_int(&maps, fin, swap);
+	header_ok *= read_float(&dmin, fin, swap);
+	header_ok *= read_float(&dmax, fin, swap);
+	header_ok *= read_float(&dmean, fin, swap);
+	for (i = 23; i < 50; ++i) header_ok *= read_float(&dummy, fin, swap);
+	header_ok *= read_float(&xorigin, fin, swap);
+	header_ok *= read_float(&yorigin, fin, swap);
+	header_ok *= read_float(&zorigin, fin, swap);
+	fclose (fin);
+	if (header_ok == 0) {
+		cout << "ERROR reading MRC header: " << vol_file;
+		exit(1);
+	}
+
+	n_range_viols += (nc > 5000); n_range_viols += (nc < 0);
+	n_range_viols += (nr > 5000); n_range_viols += (nr < 0);
+	n_range_viols += (ns > 5000); n_range_viols += (ns < 0);
+	n_range_viols += (ncstart > 5000); n_range_viols += (ncstart < -5000);
+	n_range_viols += (nrstart > 5000); n_range_viols += (nrstart < -5000);
+	n_range_viols += (nsstart > 5000); n_range_viols += (nsstart < -5000);
+	n_range_viols += (mx > 5000); n_range_viols += (mx < 0);
+	n_range_viols += (my > 5000); n_range_viols += (my < 0);
+	n_range_viols += (mz > 5000); n_range_viols += (mz < 0);
+	n_range_viols += (alpha > 360.0f); n_range_viols += (alpha < -360.0f);
+	n_range_viols += (beta > 360.0f); n_range_viols += (beta < -360.0f);
+	n_range_viols += (gamma > 360.0f); n_range_viols += (gamma < -360.0f);
+
+	return n_range_viols;
 }
 
diff --git a/model.cpp b/model.cpp
index e25d8f89aa6b73793284ba4ebc76bb90b5759af9..6e46591645997dd90f8420fbce23d822ccf650ba 100644
--- a/model.cpp
+++ b/model.cpp
@@ -30,19 +30,19 @@ int bioem_model::readModel()
 
 		char line[512] = {' '};
 		char tmpLine[512] = {' '};
-		int numres=0;
-		NormDen=0.0;
+		int numres = 0;
+		NormDen = 0.0;
 
 		//  cout << " HERE	" << filemodel ;
 		// for eachline in the file
 		while (!input.eof())
 		{
-			input.getline(line,512);
+			input.getline(line, 512);
 
-			strncpy(tmpLine,line,strlen(line));
-			char *token = strtok(tmpLine," ");
+			strncpy(tmpLine, line, strlen(line));
+			char *token = strtok(tmpLine, " ");
 
-			if (strcmp(token,"ATOM")==0) // Optional,Mandatory if standard residues exist
+			if (strcmp(token, "ATOM") == 0) // Optional,Mandatory if standard residues exist
 			{
 				/*
 				1-6 			"ATOM  "
@@ -66,26 +66,26 @@ int bioem_model::readModel()
 				char tmp[6] = {' '};
 
 				// parse name
-				strncpy(tmp,line+12,4);
-				sscanf (tmp,"%s",name);
+				strncpy(tmp, line + 12, 4);
+				sscanf (tmp, "%s", name);
 
 				// parse resName
-				strncpy(tmp,line+17,3);
-				sscanf (tmp,"%s",resName);
+				strncpy(tmp, line + 17, 3);
+				sscanf (tmp, "%s", resName);
 
 				// parse x, y, z
 				char tmpVals[36]  = {' '};
 
-				strncpy (tmpVals,line+30,8);
-				sscanf (tmpVals,"%f",&x);
+				strncpy (tmpVals, line + 30, 8);
+				sscanf (tmpVals, "%f", &x);
 
-				strncpy (tmpVals,line+38,8);
-				sscanf (tmpVals,"%f",&y);
+				strncpy (tmpVals, line + 38, 8);
+				sscanf (tmpVals, "%f", &y);
 
-				strncpy (tmpVals,line+46,8);
-				sscanf (tmpVals,"%f",&z);
+				strncpy (tmpVals, line + 46, 8);
+				sscanf (tmpVals, "%f", &z);
 
-				if (strcmp(name,"CA") == 0)
+				if (strcmp(name, "CA") == 0)
 				{
 					if (numres >= BIOEM_MODEL_SIZE)
 					{
@@ -93,30 +93,30 @@ int bioem_model::readModel()
 						exit(1);
 					}
 					//Getting residue Radius and electron density
-					radiusPointsModel[numres]=getAminoAcidRad(resName);
-					densityPointsModel[numres]=getAminoAcidDensity(resName);
-					NormDen+=densityPointsModel[numres];
+					radiusPointsModel[numres] = getAminoAcidRad(resName);
+					densityPointsModel[numres] = getAminoAcidDensity(resName);
+					NormDen += densityPointsModel[numres];
 
 					//Getting the coordinates
 					PointsModel[numres].pos[0] = (myfloat_t) x;
 					PointsModel[numres].pos[1] = (myfloat_t) y;
 					PointsModel[numres].pos[2] = (myfloat_t) z;
-					exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres]  << "\n";
+					exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres]  << "\n";
 					numres++;
 				}
 			}
 
 
 		}
-		nPointsModel=numres;
+		nPointsModel = numres;
 		cout << "Protein structure read from PDB\nTotal Number of Residues " << nPointsModel;
-		
+
 	}
 	else //Reading model from FILE FORMAT x,y,z,rad,density
 	{
 		char line[128];
-		int numres=0;
-		NormDen=0.0;
+		int numres = 0;
+		NormDen = 0.0;
 		FILE *file = fopen ( filemodel , "r" );
 		if ( file != NULL )
 		{
@@ -135,13 +135,13 @@ int bioem_model::readModel()
 				radiusPointsModel[numres] = (myfloat_t) tmpval[3];
 				densityPointsModel[numres] = (myfloat_t) tmpval[4];
 
-				exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " <<densityPointsModel[numres]<< "\n";
-				NormDen+=densityPointsModel[numres];
+				exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres] << "\n";
+				NormDen += densityPointsModel[numres];
 				numres++;
 			}
-			nPointsModel=numres;
+			nPointsModel = numres;
 			cout << "Protein structure read from Standard File \nTotal Number of Voxels " << nPointsModel ;
-			
+
 		}
 	}
 	exampleReadCoor.close();
@@ -149,22 +149,22 @@ int bioem_model::readModel()
 	//Moving to Model to its center of mass:
 	myfloat3_t r_cm;
 
-	for(int n=0; n < 3; n++)r_cm.pos[n] = 0.0;
+	for(int n = 0; n < 3; n++)r_cm.pos[n] = 0.0;
 
-	for(int n=0; n < nPointsModel; n++)
+	for(int n = 0; n < nPointsModel; n++)
 	{
 		r_cm.pos[0] += PointsModel[n].pos[0];
 		r_cm.pos[1] += PointsModel[n].pos[1];
 		r_cm.pos[2] += PointsModel[n].pos[2];
 	}
-	r_cm.pos[0]=r_cm.pos[0] / (myfloat_t) nPointsModel;
-	r_cm.pos[1]=r_cm.pos[1] / (myfloat_t) nPointsModel;
-	r_cm.pos[2]=r_cm.pos[2] / (myfloat_t) nPointsModel;
-	for(int n=0; n < nPointsModel; n++)
+	r_cm.pos[0] = r_cm.pos[0] / (myfloat_t) nPointsModel;
+	r_cm.pos[1] = r_cm.pos[1] / (myfloat_t) nPointsModel;
+	r_cm.pos[2] = r_cm.pos[2] / (myfloat_t) nPointsModel;
+	for(int n = 0; n < nPointsModel; n++)
 	{
-		PointsModel[n].pos[0]-= r_cm.pos[0];
-		PointsModel[n].pos[1]-= r_cm.pos[1];
-		PointsModel[n].pos[2]-= r_cm.pos[2];
+		PointsModel[n].pos[0] -= r_cm.pos[0];
+		PointsModel[n].pos[1] -= r_cm.pos[1];
+		PointsModel[n].pos[2] -= r_cm.pos[2];
 
 	}
 
@@ -174,28 +174,28 @@ int bioem_model::readModel()
 myfloat_t bioem_model::getAminoAcidRad(char *name)
 {
 	//************** Function that gets the radius for each amino acid ****************
-	myfloat_t iaa=0;
-
-	if(std::strcmp(name,"CYS")==0)iaa=2.75;
-	if(std::strcmp(name,"PHE")==0)iaa=3.2;
-	if(std::strcmp(name,"LEU")==0)iaa=3.1;
-	if(std::strcmp(name,"TRP")==0)iaa=3.4;
-	if(std::strcmp(name,"VAL")==0)iaa=2.95;
-	if(std::strcmp(name,"ILE")==0)iaa=3.1;
-	if(std::strcmp(name,"MET")==0)iaa=3.1;
-	if(std::strcmp(name,"HIS")==0)iaa=3.05;
-	if(std::strcmp(name,"TYR")==0)iaa=3.25;
-	if(std::strcmp(name,"ALA")==0)iaa=2.5;
-	if(std::strcmp(name,"GLY")==0)iaa=2.25;
-	if(std::strcmp(name,"PRO")==0)iaa=2.8;
-	if(std::strcmp(name,"ASN")==0)iaa=2.85;
-	if(std::strcmp(name,"THR")==0)iaa=2.8;
-	if(std::strcmp(name,"SER")==0)iaa=2.6;
-	if(std::strcmp(name,"ARG")==0)iaa=3.3;
-	if(std::strcmp(name,"GLN")==0)iaa=3.0;
-	if(std::strcmp(name,"ASP")==0)iaa=2.8;
-	if(std::strcmp(name,"LYS")==0)iaa=3.2;
-	if(std::strcmp(name,"GLU")==0)iaa=2.95;
+	myfloat_t iaa = 0;
+
+	if(std::strcmp(name, "CYS") == 0)iaa = 2.75;
+	if(std::strcmp(name, "PHE") == 0)iaa = 3.2;
+	if(std::strcmp(name, "LEU") == 0)iaa = 3.1;
+	if(std::strcmp(name, "TRP") == 0)iaa = 3.4;
+	if(std::strcmp(name, "VAL") == 0)iaa = 2.95;
+	if(std::strcmp(name, "ILE") == 0)iaa = 3.1;
+	if(std::strcmp(name, "MET") == 0)iaa = 3.1;
+	if(std::strcmp(name, "HIS") == 0)iaa = 3.05;
+	if(std::strcmp(name, "TYR") == 0)iaa = 3.25;
+	if(std::strcmp(name, "ALA") == 0)iaa = 2.5;
+	if(std::strcmp(name, "GLY") == 0)iaa = 2.25;
+	if(std::strcmp(name, "PRO") == 0)iaa = 2.8;
+	if(std::strcmp(name, "ASN") == 0)iaa = 2.85;
+	if(std::strcmp(name, "THR") == 0)iaa = 2.8;
+	if(std::strcmp(name, "SER") == 0)iaa = 2.6;
+	if(std::strcmp(name, "ARG") == 0)iaa = 3.3;
+	if(std::strcmp(name, "GLN") == 0)iaa = 3.0;
+	if(std::strcmp(name, "ASP") == 0)iaa = 2.8;
+	if(std::strcmp(name, "LYS") == 0)iaa = 3.2;
+	if(std::strcmp(name, "GLU") == 0)iaa = 2.95;
 
 	if(iaa == 0)
 	{
@@ -209,28 +209,28 @@ myfloat_t bioem_model::getAminoAcidRad(char *name)
 myfloat_t bioem_model::getAminoAcidDensity(char *name)
 {
 	//************** Function that gets the number of electrons for each amino acid ****************
-	myfloat_t iaa=0.0;
-
-	if(std::strcmp(name,"CYS")==0)iaa=64.0;
-	if(std::strcmp(name,"PHE")==0)iaa=88.0;
-	if(std::strcmp(name,"LEU")==0)iaa=72.0;
-	if(std::strcmp(name,"TRP")==0)iaa=108.0;
-	if(std::strcmp(name,"VAL")==0)iaa=64.0;
-	if(std::strcmp(name,"ILE")==0)iaa=72.0;
-	if(std::strcmp(name,"MET")==0)iaa=80.0;
-	if(std::strcmp(name,"HIS")==0)iaa=82.0;
-	if(std::strcmp(name,"TYR")==0)iaa=96.0;
-	if(std::strcmp(name,"ALA")==0)iaa=48.0;
-	if(std::strcmp(name,"GLY")==0)iaa=40.0;
-	if(std::strcmp(name,"PRO")==0)iaa=62.0;
-	if(std::strcmp(name,"ASN")==0)iaa=66.0;
-	if(std::strcmp(name,"THR")==0)iaa=64.0;
-	if(std::strcmp(name,"SER")==0)iaa=56.0;
-	if(std::strcmp(name,"ARG")==0)iaa=93.0;
-	if(std::strcmp(name,"GLN")==0)iaa=78.0;
-	if(std::strcmp(name,"ASP")==0)iaa=59.0;
-	if(std::strcmp(name,"LYS")==0)iaa=79.0;
-	if(std::strcmp(name,"GLU")==0)iaa=53.0;
+	myfloat_t iaa = 0.0;
+
+	if(std::strcmp(name, "CYS") == 0)iaa = 64.0;
+	if(std::strcmp(name, "PHE") == 0)iaa = 88.0;
+	if(std::strcmp(name, "LEU") == 0)iaa = 72.0;
+	if(std::strcmp(name, "TRP") == 0)iaa = 108.0;
+	if(std::strcmp(name, "VAL") == 0)iaa = 64.0;
+	if(std::strcmp(name, "ILE") == 0)iaa = 72.0;
+	if(std::strcmp(name, "MET") == 0)iaa = 80.0;
+	if(std::strcmp(name, "HIS") == 0)iaa = 82.0;
+	if(std::strcmp(name, "TYR") == 0)iaa = 96.0;
+	if(std::strcmp(name, "ALA") == 0)iaa = 48.0;
+	if(std::strcmp(name, "GLY") == 0)iaa = 40.0;
+	if(std::strcmp(name, "PRO") == 0)iaa = 62.0;
+	if(std::strcmp(name, "ASN") == 0)iaa = 66.0;
+	if(std::strcmp(name, "THR") == 0)iaa = 64.0;
+	if(std::strcmp(name, "SER") == 0)iaa = 56.0;
+	if(std::strcmp(name, "ARG") == 0)iaa = 93.0;
+	if(std::strcmp(name, "GLN") == 0)iaa = 78.0;
+	if(std::strcmp(name, "ASP") == 0)iaa = 59.0;
+	if(std::strcmp(name, "LYS") == 0)iaa = 79.0;
+	if(std::strcmp(name, "GLU") == 0)iaa = 53.0;
 
 	if(iaa == 0.0)
 	{
diff --git a/param.cpp b/param.cpp
index e103f09858b854fa6996924dea9fa72b95ff4929..077e378f27546da0038745f0d36c59f570107692 100644
--- a/param.cpp
+++ b/param.cpp
@@ -14,7 +14,7 @@ using namespace std;
 bioem_param::bioem_param()
 {
 	//Number of Pixels
-	param_device.NumberPixels=0;
+	param_device.NumberPixels = 0;
 	param_device.NumberFFTPixels1D = 0;
 	// Euler angle grid spacing
 	angleGridPointsAlpha = 0;
@@ -41,25 +41,25 @@ int bioem_param::readParameters()
 	//***************************** Reading Input Parameters ******************************
 	//**************************************************************************************
 
-        // Control for Parameters
-        bool yesPixSi = false;
-        bool yesNumPix = false;
-        bool yesGPal = false;
+	// Control for Parameters
+	bool yesPixSi = false;
+	bool yesNumPix = false;
+	bool yesGPal = false;
 	bool yesGPbe = false;
-        bool yesGPEnv = false;
-        bool yesGPamp = false;
-        bool yesGPpha = false;
-        bool yesSTEnv = false;
+	bool yesGPEnv = false;
+	bool yesGPamp = false;
+	bool yesGPpha = false;
+	bool yesSTEnv = false;
 	bool yesSTamp = false;
 	bool yesSTpha = false;
 	bool yesGSPamp = false ;
 	bool yesGSPEnv = false ;
 	bool yesGSPpha = false ;
 	bool yesMDC = false ;
-        bool yesGCen = false ;	
-	
-	
-	
+	bool yesGCen = false ;
+
+
+
 	ifstream input(fileinput);
 	if (!input.good())
 	{
@@ -76,186 +76,258 @@ int bioem_param::readParameters()
 	cout << "+++++++++++++++++++++++++++++++++++++++++++ \n";
 	while (!input.eof())
 	{
-		input.getline(line,512);
-		strcpy(saveline,line);
-		char *token = strtok(line," ");
+		input.getline(line, 512);
+		strcpy(saveline, line);
+		char *token = strtok(line, " ");
 
-		if (token==NULL || line[0] == '#' || strlen(token)==0)
+		if (token == NULL || line[0] == '#' || strlen(token) == 0)
 		{
 			// comment or blank line
 		}
-		else if (strcmp(token,"PIXEL_SIZE")==0)
+		else if (strcmp(token, "PIXEL_SIZE") == 0)
 		{
-			token = strtok(NULL," ");
-			pixelSize=atof(token);
-			if (pixelSize < 0 ) { cout << "*** Error: Negative pixelSize "; exit(1);}
+			token = strtok(NULL, " ");
+			pixelSize = atof(token);
+			if (pixelSize < 0 ) {
+				cout << "*** Error: Negative pixelSize ";
+				exit(1);
+			}
 			cout << "Pixel Sixe " << pixelSize << "\n";
-                        yesPixSi= true;
+			yesPixSi = true;
 
 		}
-		else if (strcmp(token,"NUMBER_PIXELS")==0)
+		else if (strcmp(token, "NUMBER_PIXELS") == 0)
 		{
-			token = strtok(NULL," ");
-			param_device.NumberPixels=int(atoi(token));
-			if (param_device.NumberPixels < 0 ) { cout << "*** Error: Negative Number of Pixels "; exit(1);}
-			cout << "Number of Pixels " << param_device.NumberPixels << "\n";			
-                        yesNumPix= true ;
+			token = strtok(NULL, " ");
+			param_device.NumberPixels = int(atoi(token));
+			if (param_device.NumberPixels < 0 ) {
+				cout << "*** Error: Negative Number of Pixels ";
+				exit(1);
+			}
+			cout << "Number of Pixels " << param_device.NumberPixels << "\n";
+			yesNumPix = true ;
 		}
-		else if (strcmp(token,"GRIDPOINTS_ALPHA")==0)
+		else if (strcmp(token, "GRIDPOINTS_ALPHA") == 0)
 		{
-			token = strtok(NULL," ");
-			angleGridPointsAlpha=int(atoi(token));
-			if (angleGridPointsAlpha < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ALPHA "; exit(1);}
+			token = strtok(NULL, " ");
+			angleGridPointsAlpha = int(atoi(token));
+			if (angleGridPointsAlpha < 0 ) {
+				cout << "*** Error: Negative GRIDPOINTS_ALPHA ";
+				exit(1);
+			}
 			cout << "Grid points alpha " << angleGridPointsAlpha << "\n";
-                        yesGPal= true;
+			yesGPal = true;
 		}
-		else if (strcmp(token,"GRIDPOINTS_BETA")==0)
+		else if (strcmp(token, "GRIDPOINTS_BETA") == 0)
 		{
-			token = strtok(NULL," ");
-			angleGridPointsBeta=int(atoi(token));			
-			if (angleGridPointsBeta < 0 ) { cout << "*** Error: Negative GRIDPOINTS_BETA "; exit(1);}
+			token = strtok(NULL, " ");
+			angleGridPointsBeta = int(atoi(token));
+			if (angleGridPointsBeta < 0 ) {
+				cout << "*** Error: Negative GRIDPOINTS_BETA ";
+				exit(1);
+			}
 			cout << "Grid points beta " << angleGridPointsBeta << "\n";
-                        yesGPbe= true;
+			yesGPbe = true;
 		}
 
-		else if (strcmp(token,"GRIDPOINTS_ENVELOPE")==0)
+		else if (strcmp(token, "GRIDPOINTS_ENVELOPE") == 0)
 		{
-			token = strtok(NULL," ");
-			numberGridPointsEnvelop=int(atoi(token));
-			if (numberGridPointsDisplaceCenter < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ENVELOPE "; exit(1);}
+			token = strtok(NULL, " ");
+			numberGridPointsEnvelop = int(atoi(token));
+			if (numberGridPointsDisplaceCenter < 0 ) {
+				cout << "*** Error: Negative GRIDPOINTS_ENVELOPE ";
+				exit(1);
+			}
 			cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";
 			yesGPEnv = true;
 		}
-		else if (strcmp(token,"START_ENVELOPE")==0)
+		else if (strcmp(token, "START_ENVELOPE") == 0)
 		{
-			token = strtok(NULL," ");
-			startGridEnvelop=atof(token);
-			if (startGridEnvelop < 0 ) { cout << "*** Error: Negative START_ENVELOPE "; exit(1);}
+			token = strtok(NULL, " ");
+			startGridEnvelop = atof(token);
+			if (startGridEnvelop < 0 ) {
+				cout << "*** Error: Negative START_ENVELOPE ";
+				exit(1);
+			}
 			cout << "Start Envelope " << startGridEnvelop << "\n";
 			yesSTEnv = true ;
 
 		}
-		else if (strcmp(token,"GRIDSPACE_ENVELOPE")==0)
+		else if (strcmp(token, "GRIDSPACE_ENVELOPE") == 0)
 		{
-			token = strtok(NULL," ");
-			gridEnvelop=atof(token);
-			if (gridEnvelop < 0 ) { cout << "*** Error: Negative GRIDSPACE_ENVELOPE "; exit(1);}
+			token = strtok(NULL, " ");
+			gridEnvelop = atof(token);
+			if (gridEnvelop < 0 ) {
+				cout << "*** Error: Negative GRIDSPACE_ENVELOPE ";
+				exit(1);
+			}
 			cout << "Grid spacing Envelope " << gridEnvelop << "\n";
 			yesGSPEnv = true ;
 
 		}
-		else if (strcmp(token,"GRIDPOINTS_PSF_PHASE")==0)
+		else if (strcmp(token, "GRIDPOINTS_PSF_PHASE") == 0)
 		{
-			token = strtok(NULL," ");
-			numberGridPointsCTF_phase=int(atoi(token));
+			token = strtok(NULL, " ");
+			numberGridPointsCTF_phase = int(atoi(token));
 			cout << "Grid points PSF " << numberGridPointsCTF_phase << "\n";
 			yesGPpha = true;
 
 		}
-		else if (strcmp(token,"START_PSF_PHASE")==0)
+		else if (strcmp(token, "START_PSF_PHASE") == 0)
 		{
-			token = strtok(NULL," ");
-			startGridCTF_phase=atof(token);
+			token = strtok(NULL, " ");
+			startGridCTF_phase = atof(token);
 			cout << "Start PSF " << startGridCTF_phase << "\n";
 			yesSTpha = true ;
 
 		}
-		else if (strcmp(token,"GRIDSPACE_PSF_PHASE")==0)
+		else if (strcmp(token, "GRIDSPACE_PSF_PHASE") == 0)
 		{
-			token = strtok(NULL," ");
-			gridCTF_phase=atof(token);
+			token = strtok(NULL, " ");
+			gridCTF_phase = atof(token);
 			cout << "Grid Space PSF " << gridCTF_phase << "\n";
 			yesGSPpha = true ;
 
-		} else if (strcmp(token,"GRIDPOINTS_PSF_AMP")==0)
+		} else if (strcmp(token, "GRIDPOINTS_PSF_AMP") == 0)
 		{
-			token = strtok(NULL," ");
-			numberGridPointsCTF_amp=int(atoi(token));
+			token = strtok(NULL, " ");
+			numberGridPointsCTF_amp = int(atoi(token));
 			cout << "Grid points PSF " << numberGridPointsCTF_amp << "\n";
 			yesGPamp = true ;
 
 		}
-		else if (strcmp(token,"START_PSF_AMP")==0)
+		else if (strcmp(token, "START_PSF_AMP") == 0)
 		{
-			token = strtok(NULL," ");
-			startGridCTF_amp=atof(token);
+			token = strtok(NULL, " ");
+			startGridCTF_amp = atof(token);
 			cout << "Start PSF " << startGridCTF_amp << "\n";
 			yesSTamp = true ;
 
 		}
-		else if (strcmp(token,"GRIDSPACE_PSF_AMP")==0)
+		else if (strcmp(token, "GRIDSPACE_PSF_AMP") == 0)
 		{
-			token = strtok(NULL," ");
-			gridCTF_amp=atof(token);
+			token = strtok(NULL, " ");
+			gridCTF_amp = atof(token);
 			cout << "Grid Space PSF " << gridCTF_amp << "\n";
 			yesGSPamp = true ;
 
 		}
-		else if (strcmp(token,"MAX_D_CENTER")==0)
+		else if (strcmp(token, "MAX_D_CENTER") == 0)
 		{
-			token = strtok(NULL," ");
-			param_device.maxDisplaceCenter=int(atoi(token));
-			if (param_device.maxDisplaceCenter < 0 ) { cout << "*** Error: Negative MAX_D_CENTER "; exit(1);}
+			token = strtok(NULL, " ");
+			param_device.maxDisplaceCenter = int(atoi(token));
+			if (param_device.maxDisplaceCenter < 0 ) {
+				cout << "*** Error: Negative MAX_D_CENTER ";
+				exit(1);
+			}
 			cout << "Maximum displacement Center " <<  param_device.maxDisplaceCenter << "\n";
 			yesMDC = true;
 
 		}
-		else if (strcmp(token,"PIXEL_GRID_CENTER")==0)
+		else if (strcmp(token, "PIXEL_GRID_CENTER") == 0)
 		{
-			token = strtok(NULL," ");
-			param_device.GridSpaceCenter=int(atoi(token));
-			if (param_device.GridSpaceCenter < 0 ) { cout << "*** Error: Negative PIXEL_GRID_CENTER "; exit(1);}
+			token = strtok(NULL, " ");
+			param_device.GridSpaceCenter = int(atoi(token));
+			if (param_device.GridSpaceCenter < 0 ) {
+				cout << "*** Error: Negative PIXEL_GRID_CENTER ";
+				exit(1);
+			}
 			cout << "Grid space displacement center " <<   param_device.GridSpaceCenter << "\n";
 			yesGCen = true;
 
 		}
-		else if (strcmp(token,"WRITE_PROB_ANGLES")==0)
+		else if (strcmp(token, "WRITE_PROB_ANGLES") == 0)
 		{
-			writeAngles=true;
+			writeAngles = true;
 			cout << "Writing Probabilies of each angle \n";
 
 		}
 
 	}
 	input.close();
-	
-	// Checks for ALL INPUT 
-	if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);};
-        if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);};
-        if( not ( yesGPal ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);};
-	if( not ( yesGPbe ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ; exit (1);};
-        if( not (  yesGPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ; exit (1);};
-        if( not (  yesGPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ; exit (1);};
-        if( not (  yesGPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ; exit (1);};
-        if( not (  yesSTEnv  ) ) { cout << "**** INPUT MISSING: Please provide START_ENVELOPE \n" ; exit (1);};
-	if( not (  yesSTamp  ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_AMP \n" ; exit (1);};
-	if( not (  yesSTpha  ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_PHASE \n" ; exit (1);};
-	if( not (  yesGSPamp  ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_AMP \n" ; exit (1);};
-	if( not ( yesGSPEnv  ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_ENVELOPE \n" ; exit (1);};
-	if( not ( yesGSPpha  ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_PHASE \n" ; exit (1);};
-	if( not (  yesMDC  ) ) { cout << "**** INPUT MISSING: Please provide MAX_D_CENTER \n" ; exit (1);};
-        if( not (  yesGCen  ) ) { cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ; exit (1);};
-	
+
+	// Checks for ALL INPUT
+	if( not ( yesPixSi ) ) {
+		cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ;
+		exit (1);
+	};
+	if( not ( yesNumPix ) ) {
+		cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ;
+		exit (1);
+	};
+	if( not ( yesGPal ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ;
+		exit (1);
+	};
+	if( not ( yesGPbe ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ;
+		exit (1);
+	};
+	if( not (  yesGPEnv ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ;
+		exit (1);
+	};
+	if( not (  yesGPamp ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ;
+		exit (1);
+	};
+	if( not (  yesGPpha ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ;
+		exit (1);
+	};
+	if( not (  yesSTEnv  ) ) {
+		cout << "**** INPUT MISSING: Please provide START_ENVELOPE \n" ;
+		exit (1);
+	};
+	if( not (  yesSTamp  ) ) {
+		cout << "**** INPUT MISSING: Please provide START_PSF_AMP \n" ;
+		exit (1);
+	};
+	if( not (  yesSTpha  ) ) {
+		cout << "**** INPUT MISSING: Please provide START_PSF_PHASE \n" ;
+		exit (1);
+	};
+	if( not (  yesGSPamp  ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_AMP \n" ;
+		exit (1);
+	};
+	if( not ( yesGSPEnv  ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDSPACE_ENVELOPE \n" ;
+		exit (1);
+	};
+	if( not ( yesGSPpha  ) ) {
+		cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_PHASE \n" ;
+		exit (1);
+	};
+	if( not (  yesMDC  ) ) {
+		cout << "**** INPUT MISSING: Please provide MAX_D_CENTER \n" ;
+		exit (1);
+	};
+	if( not (  yesGCen  ) ) {
+		cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ;
+		exit (1);
+	};
+
 	//More checks with input parameters
-	
-	// Envelope should not have a standard deviation greater than Npix/2 
-	if(sqrt(1./( (myfloat_t) numberGridPointsDisplaceCenter  * gridEnvelop + startGridEnvelop))> float(param_device.NumberPixels)/2.0) {
-	  cout << "MAX Standar deviation of envelope is larger than Allowed KERNEL Length " ;
-	  exit(1);
+
+	// Envelope should not have a standard deviation greater than Npix/2
+	if(sqrt(1. / ( (myfloat_t) numberGridPointsDisplaceCenter  * gridEnvelop + startGridEnvelop)) > float(param_device.NumberPixels) / 2.0) {
+		cout << "MAX Standar deviation of envelope is larger than Allowed KERNEL Length " ;
+		exit(1);
 	}
 	// Envelop param should be positive
-	if( startGridCTF_amp < 0 || startGridCTF_amp > 1){
-	   cout << "PSF Amplitud should be between 0 and 1  \n" ;	   
-	   exit(1);
+	if( startGridCTF_amp < 0 || startGridCTF_amp > 1) {
+		cout << "PSF Amplitud should be between 0 and 1  \n" ;
+		exit(1);
 	}
-	
-	if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp > 1){
-	   cout << "PSF Amplitud should be between 0 and 1  \n" ;	   
-	   exit(1);
+
+	if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp > 1) {
+		cout << "PSF Amplitud should be between 0 and 1  \n" ;
+		exit(1);
 	}
-	
-	
-	  
+
+
+
 	param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
 	RefMapSize = param_device.NumberPixels * param_device.NumberFFTPixels1D;
 	cout << "+++++++++++++++++++++++++++++++++++++++++++ \n";
@@ -300,14 +372,14 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
 	//************************************************************************************
 	//**************** Routine that pre-calculates Euler angle grids *********************
 	//************************************************************************************
-	myfloat_t grid_alpha,cos_grid_beta;
-	int n=0;
+	myfloat_t grid_alpha, cos_grid_beta;
+	int n = 0;
 
 	//alpha and gamma are uniform in -PI,PI
-	grid_alpha=2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
+	grid_alpha = 2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
 
 	//cosine beta is uniform in -1,1
-	cos_grid_beta=2.f / (myfloat_t) angleGridPointsBeta;
+	cos_grid_beta = 2.f / (myfloat_t) angleGridPointsBeta;
 
 	// Euler Angle Array
 	for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
@@ -316,25 +388,25 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
 		{
 			for (int igamma = 0; igamma < angleGridPointsAlpha; igamma ++)
 			{
-				angles[n].pos[0]= (myfloat_t) ialpha * grid_alpha - M_PI + grid_alpha * 0.5f; //ALPHA centered in the middle
-				angles[n].pos[1]= acos((myfloat_t) ibeta * cos_grid_beta - 1 + cos_grid_beta * 0.5f); //BETA centered in the middle
-				angles[n].pos[2]= (myfloat_t) igamma * grid_alpha - M_PI + grid_alpha * 0.5f; //GAMMA centered in the middle
+				angles[n].pos[0] = (myfloat_t) ialpha * grid_alpha - M_PI + grid_alpha * 0.5f; //ALPHA centered in the middle
+				angles[n].pos[1] = acos((myfloat_t) ibeta * cos_grid_beta - 1 + cos_grid_beta * 0.5f); //BETA centered in the middle
+				angles[n].pos[2] = (myfloat_t) igamma * grid_alpha - M_PI + grid_alpha * 0.5f; //GAMMA centered in the middle
 				n++;
 			}
 		}
 	}
-	nTotGridAngles=n;
+	nTotGridAngles = n;
 
 	//********** Calculating normalized volumen element *********
 
 	param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
-					  * gridCTF_phase * gridCTF_amp * gridEnvelop / (2.f * M_PI) / (2.f * M_PI) / 2.f / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp)
-					  / ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase + startGridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop + startGridEnvelop);
+						* gridCTF_phase * gridCTF_amp * gridEnvelop / (2.f * M_PI) / (2.f * M_PI) / 2.f / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp)
+						/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase + startGridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop + startGridEnvelop);
 
 	//*** Number of total pixels***
 
-	param_device.Ntotpi= (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
-	param_device.NtotDist=(2 * (int) (param_device.maxDisplaceCenter/param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter/param_device.GridSpaceCenter) + 1);
+	param_device.Ntotpi = (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
+	param_device.NtotDist = (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter / param_device.GridSpaceCenter) + 1);
 
 	return(0);
 
@@ -346,12 +418,12 @@ int bioem_param::CalculateRefCTF()
 	//********** Routine that pre-calculates Kernels for Convolution *********************
 	//************************************************************************************
 
-	myfloat_t amp,env,phase,ctf,radsq;
+	myfloat_t amp, env, phase, ctf, radsq;
 	myfloat_t* localCTF;
-	int nctfmax= param_device.NumberPixels / 2;
-	int n=0;
+	int nctfmax = param_device.NumberPixels / 2;
+	int n = 0;
 
-	localCTF= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels*param_device.NumberPixels);
+	localCTF = (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels * param_device.NumberPixels);
 
 	nTotCTFs = numberGridPointsCTF_amp * numberGridPointsCTF_phase * numberGridPointsEnvelop;
 	delete[] refCTF;
@@ -362,47 +434,47 @@ int bioem_param::CalculateRefCTF()
 	for (int iamp = 0; iamp <  numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
 	{
 		amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
-		
+
 		for (int iphase = 0; iphase <  numberGridPointsCTF_phase ; iphase++)//Loop over phase
 		{
 			phase = (myfloat_t) iphase * gridCTF_phase + startGridCTF_phase;
 
 			for (int ienv = 0; ienv <  numberGridPointsEnvelop ; ienv++)//Loop over envelope
 			{
-				env= (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
+				env = (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
 
-				memset(localCTF,0,param_device.NumberPixels*param_device.NumberPixels*sizeof(myfloat_t));
+				memset(localCTF, 0, param_device.NumberPixels * param_device.NumberPixels * sizeof(myfloat_t));
 
 				//Assigning CTF values
-				for(int i=0; i< nctfmax; i++)
+				for(int i = 0; i < nctfmax; i++)
 				{
-					for(int j=0; j< nctfmax; j++)
+					for(int j = 0; j < nctfmax; j++)
 					{
-						radsq=(myfloat_t) (i*i+j*j) * pixelSize * pixelSize;
-						ctf=exp(-radsq*env/2.0)*(amp*cos(radsq*phase/2.0)-sqrt((1-amp*amp))*sin(radsq*phase/2.0));
-						
-						localCTF[i*param_device.NumberPixels+j]=(myfloat_t) ctf;
-						localCTF[i*param_device.NumberPixels+param_device.NumberPixels-j-1]=(myfloat_t) ctf;
-						localCTF[(param_device.NumberPixels-i-1)*param_device.NumberPixels+j]=(myfloat_t) ctf;
-						localCTF[(param_device.NumberPixels-i-1)*param_device.NumberPixels+param_device.NumberPixels-j-1]=(myfloat_t) ctf;
+						radsq = (myfloat_t) (i * i + j * j) * pixelSize * pixelSize;
+						ctf = exp(-radsq * env / 2.0) * (amp * cos(radsq * phase / 2.0) - sqrt((1 - amp * amp)) * sin(radsq * phase / 2.0));
+
+						localCTF[i * param_device.NumberPixels + j] = (myfloat_t) ctf;
+						localCTF[i * param_device.NumberPixels + param_device.NumberPixels - j - 1] = (myfloat_t) ctf;
+						localCTF[(param_device.NumberPixels - i - 1)*param_device.NumberPixels + j] = (myfloat_t) ctf;
+						localCTF[(param_device.NumberPixels - i - 1)*param_device.NumberPixels + param_device.NumberPixels - j - 1] = (myfloat_t) ctf;
 					}
 				}
 				//Creatting FFT_Forward of Kernel to store
 				mycomplex_t* localout;
-				localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param_device.NumberPixels*param_device.NumberFFTPixels1D);
+				localout = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberFFTPixels1D);
 				//Calling FFT_Forward
-				myfftw_execute_dft_r2c(fft_plan_r2c_forward,localCTF,localout);
+				myfftw_execute_dft_r2c(fft_plan_r2c_forward, localCTF, localout);
 
 				// Normalizing and saving the Reference CTFs
 				mycomplex_t* curRef = &refCTF[n * RefMapSize];
-				for(int i=0; i < param_device.NumberPixels * param_device.NumberFFTPixels1D; i++ )
+				for(int i = 0; i < param_device.NumberPixels * param_device.NumberFFTPixels1D; i++ )
 				{
 					curRef[i][0] = localout[i][0];
 					curRef[i][1] = localout[i][1];
 				}
-				CtfParam[n].pos[0]=amp;
-				CtfParam[n].pos[1]=phase;
-				CtfParam[n].pos[2]=env;
+				CtfParam[n].pos[0] = amp;
+				CtfParam[n].pos[1] = phase;
+				CtfParam[n].pos[2] = env;
 				n++;
 				myfftw_free(localout);
 			}
@@ -423,7 +495,7 @@ int bioem_param::CalculateRefCTF()
 bioem_param::~bioem_param()
 {
 	releaseFFTPlans();
-	param_device.NumberPixels=0;
+	param_device.NumberPixels = 0;
 	angleGridPointsAlpha = 0;
 	angleGridPointsBeta = 0;
 	numberGridPointsEnvelop = 0;