 				float uu_N_base = uu_N_n - neighborRadius;
 				float uv_N_base = uv_N_n - neighborRadius;
 				
-				// xyz value of color for current wavelength.
+				// xyz value of color for current wavelength (regarding current wavenumber k).
 				vec3 waveColor = avgWeighted_XYZ_weight(lambda_iter);
-				
 				
 				// iterte over neighbor winodw of size (2*neighborRadius+1)x(2*neighborRadius+1)
 				for(float ind1 = uu_N_base; ind1 <= uu_N_base + 2.0*neighborRadius; ind1 = ind1 + 1.0){
 					for(float ind2 = uv_N_base; ind2 <= uv_N_base + 2.0*neighborRadius; ind2 = ind2 + 1.0){
+							
+						// get L1 distance to neighbor from current pixel
+						float dist2 = pow(ind1-uu_N_base, 2.0) + pow(ind2-uv_N_n_hat, 2.0);
 						
-						
-						// get L1 distance to neighbor
-						float dist2 = pow(ind1-uu_N_base, 2.0) + pow(ind2-uv_N_n_hat, 2.0);
-						 
-						// case distinction: depending on value of t, i.e. value of variance
-						if(variant == 0.0){
-							coords = vec2( (ind1/(dimN-1)) + bias, (ind2/(dimN-1)) + bias); //2d case
-						}else{
-							coords = vec2( (ind2/(dimN-1)) + bias, (ind1/(dimN-1)) + bias); //2d case
-						}
-
+						// get lookup coordinates for current pixel
+						coords = getLookupCoordinates(variant, ind1, ind2);
 						
 						// clip if coordiantes are not within bound [0,1]x[0,1]
 						if(coords.x < 0.0 || coords.x > 1.0 || coords.y < 0.0 || coords.y > 1.0) continue;
-						
 						
 						// complex valued frequency contribution of current pixel.
 						P = taylorApproximation(coords, k, w);
 						
 						// compute gaussion weight of current pixel
-						float exponent = -dist2/(sigma_f_pix);
-						float w_ij = exp(exponent);
-							
-						if(abs(exponent) < 1.0*pow(10.0, -18.0) ){
-							w_ij = 1.0;
-						}
-								
-//						if(abs(norm_fact) < 1.0*pow(10.0, -10.0)) norm_fact = 1.0;
-						w_ij /= norm_fact;
+						float w_ij = getGaussianWeight(dist2, sigma_f_pix);
 						
 						// phaser of current pixel contribution.
 						float pq_scale = compute_pq_scale_factor(w_u,w_v);
-						if(isinf(pq_scale)) pq_scale = 1.0;
-						if(isnan(pq_scale)) pq_scale = 1.0;
+						P *= pq_scale;
 						
-						P *= pq_scale;
-							
 						// amplitute of current pixel contribution
 						float abs_P_Sq = P.x*P.x + P.y*P.y;
+						
+						// weight with window
 						abs_P_Sq *= w_ij;
-						
-						// brdf coefficient after squaring them. 
-						float diffractionCoeff = getFactor(k, F, G, w);
-						diffractionCoeff = 1.0;
 
 						// summation of brdf over all wavelengths under consideration.
-						brdf += vec4(diffractionCoeff * abs_P_Sq * waveColor, 0.0);	
-						maxBRDF += vec4(waveColor, 0.0);					
+						brdf += vec4(abs_P_Sq * waveColor, 0.0);					
 					}
 				}
 			}
 		}
 	}
-	
-
-
-
-//	brdf.xyz = fac2*fac2*fac2*fac2*brdf.xyz;
-		
-	
+
 	float ambient = 0.0;	
 	// remove negative values
 	if(brdf.x < 0.0 ) brdf.x = 0.0;

float periods = periodCount-1.0; // 26 // number of patch periods along surface
float Omega = ((N_1/N_2)*2.0*PI)/t_0; // (N_1/N_2)*2*PI/t_0, before 8.0*PI*pow(10.0,7.0);
float bias = (N_2/2.0)/(N_2-1.0); // old: 50.0/99.0;
float neighborRadius = (neigh_rad < 5 && neigh_rad > -1) ? float(neigh_rad) : 0.0;

//transformation constant

	float phi = -PI/2.0;
	phi = 0.0;
	
	vec4 o_col = vec4(0.0, 0.0, 0.0, 0.0);
	vec4 brdf = vec4(0.0, 0.0, 0.0, 1.0);
	vec4 maxBRDF = vec4(0.0, 0.0, 0.0, 1.0);
	vec2 P = vec2(0.0, 0.0);
	vec2 coords = vec2(0.0);
	
	
	// okay from here on
	thetaR = asin(sqrt(o_org_pos.x * o_org_pos.x + o_org_pos.y * o_org_pos.y ));
	float phiR = atan(o_org_pos.y, o_org_pos.x);

	vec3 _k1 = vec3(0.0f);
	vec3 _k2 = vec3(0.0f);
	
	
	_k1.x = - sin(thetaI)*cos(phiI);
	_k1.y = - sin(thetaI)*sin(phiI);
	_k1.z = - cos(thetaI);
	
	_k2.x = sin(thetaR)*cos(phiR);
	_k2.y = sin(thetaR)*sin(phiR);
	_k2.z = cos(thetaR);

	float uu0 = _k1.x - _k2.x;
	float vv0 = _k1.y - _k2.y;
	float ww = _k1.z - _k2.z;	
	
	float fF = getFresnelFactorAbsolute(_k1, _k2);
	float nSkin = 1.5;
	float R0 = pow( (nSkin - 1.0) / (nSkin + 1.0) , 2.0); 
	
	fFByR0 = fF/R0;
	float shadowF = getShadowMaskFactor(_k1, _k2);

	vec3 V = vec3(uu0, vv0, ww);
	float u = V.x; float v = V.y; float w = V.z;
	
	// get iteration bounds for given (u,v)
	vec2 N_u = compute_N_min_max(u);
	vec2 N_v = compute_N_min_max(v);
	vec2 N_uv[2] = vec2[2](N_u, N_v);
//	vec2 modUV = getRotation(u,v,phi);
	
	float uv_sqr = pow(u*u+v*v, 0.5);

	float iterMax = 100.0;
	float lambdaStep = (lambda_max - lambda_min)/(iterMax-1.0);
	float F2 = fFByR0*fFByR0;
	
	
	float stepSize = 1.0;
