I have read some threads here asking a similar question. But the examples are different, or more easy than this. I have tested first the most simple wave and I got the normal working properly.
For me this is the case:
float theta =dot(K, wp) -w;
float2 XY=(K/k)* A * sin(theta );
float Z= A* cos (theta );
return float3(XY.x.x,XY.y, Z);
As explained in tessendorf, K is the wind vector and L is lambda.
So the question here would be how to derivate this and finally how to get the normal.