 	vec2 N_v = compute_N_min_max(v);
 	vec2 N_uv[2] = vec2[2](N_u, N_v);
 	
+
+	vec2 modUV = getRotation(u,v,-phi);
 	
 	// only specular contribution within epsilon range: i.e. fixed number of lambdas
 	if(abs(u) < eps && abs(v) < eps){
-
-
-		
+		for(int iter = 0; iter < 16; iter++){
+			float lambda_iter = fixed_lambdas[iter]*pow(10.0,-9.0);
+			k = 2.0*PI / lambda_iter;
+			vec2 coords = vec2((k*modUV.x/Omega) + bias, (k*modUV.y/Omega) + bias); //2d
+
+			if(coords.x < 0.0 || coords.x > 1.0 || coords.y < 0.0 || coords.y > 1.0) continue;
+			
+			float w_u = k*modUV.x;
+			float w_v = k*modUV.y;
+			
+			P = taylorApproximation(coords, k, w);
+			float pq_scale = compute_pq_scale_factor(w_u,w_v);
+			P *= pq_scale;
+			
+			float abs_P_Sq = P.x*P.x + P.y*P.y;
+			
+			float diffractionCoeff = getFactor(k, F, G, w);
+			vec3 waveColor = avgWeighted_XYZ_weight(lambda_iter);
+			brdf += vec4(diffractionCoeff * abs_P_Sq * waveColor, 1.0);
+			maxBRDF += vec4(diffractionCoeff * brdfMax * waveColor, 1.0);		
+		}
 	}else{
-		float bias = 50.0/99.0;
-		vec2 modUV = getRotation(u,v,-phi);
+
 		
 		// iterate twice: once for N_u and once for N_v lower,upper
 		for(int variant = 0; variant < 2; variant++){


// error constants
const float eps_pq = 0.0001; 
const float eps = 1.0*pow(10.0, -3.0);
const float tolerance = 0.999999999; 

// period constants

const float t_0 = dx / N_1;
const float T_1 = t_0 * N_1;
const float T_2 = t_0 * N_1;
const float periods = 26.0-1.0; // 26 // number of patch periods along surface
const 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);
const float bias = (N_2/2.0)/(N_2-1.0); // old: 50.0/99.0;

// transformation constant
const mat3 M_Adobe_XR = mat3(
