 			(*simulations)[i].bulks[ii].flc=ftemp[3];
 			(*simulations)[i].bulks[ii].g=ftemp[7];
 			if(ftemp[2]==0.0f)(*simulations)[i].bulks[ii].mutr=FLT_MAX; //Glass
-			else(*simulations)[i].bulks[ii].mutr=1.0f/(ftemp[1]+ftemp[2]);
-			(*simulations)[i].bulks[ii].albedof =  ftemp[2]/(ftemp[3]*ftemp[4] + ftemp[2]); //mus/(flc*AbsCExi + mus)
-			printf("Bulk desc. %i\n n=%f\n mua=%f\n muaf=%f\n mutr=%f\n eY=%f\n flc=%e\n g=%f\n albedof=%f\n AbsCExi=%f\n AbsCEmi=%f\n\n", ii,
+				else(*simulations)[i].bulks[ii].mutr=1.0f/(ftemp[1]+ftemp[2]);
+			if ((ftemp[3] >0 && ftemp[4]>0) || ftemp[2]>0) (*simulations)[i].bulks[ii].albedof =  ftemp[2]/(ftemp[3]*ftemp[4] + ftemp[2]); //mus/(flc*AbsCExi + mus)
+				else (*simulations)[i].bulks[ii].albedof=0;
+		}//end ii<n_bulks
+
+		for(ii=1;ii<=n_bulks;ii++)	{
+			printf("Bulk desc. %i\n n=%f\n mua=%f\n muaf=%f\n mutr=%f\n eY=%f\n flc=%e\n g=%f\n albedof=%f\n \n", ii,
 				(*simulations)[i].bulks[ii].n,
 				(*simulations)[i].bulks[ii].mua,
 				(*simulations)[i].bulks[ii].muaf,

		(*simulations)[i].bulks[n_bulks+1].g=1.;
		(*simulations)[i].bulks[n_bulks+1].mutr=FLT_MAX; //Glass
		(*simulations)[i].bulks[n_bulks+1].albedof = 0.; //mus/(flc*AbsCExi + mus)

		//Read matrix filename
		ii=0;
		while(ii<=0){

		(*simulations)[i].bulks[0].mutr=FLT_MAX; //Glass
		(*simulations)[i].bulks[0].albedof = 0.; //mus/(flc*AbsCExi + mus)

		for(ii=1;ii<=n_bulks+1;ii++)	{
			// Read Layer data (9x float)
			if(!readfloats(8, ftemp, pFile)){perror ("Error reading bulk descriptors data");return 0;}
			(*simulations)[i].bulks[ii].n=ftemp[0];
