Calculating Normals across a sphere with a wave-li

I've been trying to get the correct normals for a sphere I'm messing with using a vertex shader. The algorithm can be boiled down simply to

vert.xyz += max(0, sin(time + 0.004*vert.x))*10*normal.xyz

This causes a wave to roll across the sphere. In order to make my normals correct, I need to transform them as well. I can take the tangent vector at a given x,y,z, get a perpendicular vector (0, -vert.z, vert.y), and then cross the tangent with the perp vector.

I've been having some issue with the math though, and it's become a personal vendetta at this point. I've solved for the derivative hundreds of times but I keep getting it incorrect. How can I get the tangent? Breaking down the above line, I can make a math function

f(x,y,z) = max(0, sin(time + 0.004*x))*10*Norm(x,y,z) + (x,y,z)

where Norm(..) is Normalize((x,y,z) - CenterOfSphere)

After applying f(x,y,z), unchanged normals

What is the correct f '(x,y,z)?

I've accounted for the weirdness caused by the max in f(...), so that's not the issue.

Edit: The most successful algorithm I have right now is as follows:

Tangent vector.x = 0.004*10*cos(0.004*vert.x + time)*norm.x + 10*sin(0.004*vert.x + time) + 1

Tangent vector.y = 10*sin(0.004*vert.x + time) + 1

Tangent vector.z = 10*sin(0.004*vert.x + time) + 1

2nd Tangent vector.x = 0

2nd Tangent vector.y = -norm.z

2nd Tangent vector.z = norm.y

Normalize both, and perform Cross(Tangent2, Tangent1). Normalize again, and done (it should be Cross(Tangent1, Tangent2), but this seems to have better results... more hints of an issue in my math!).

This yields this

Get tangent/normal by derivate of function can sometimes fail if your surface points are nonlinearly distributed and or some math singularity is present or if you make a math mistake (which is the case in 99.99%). Anyway you can always use the geometric approach:

1. you can get the tangents easy by

U(x,y,z)=f(x+d,y,z)-f(x,y,z);
V(x,y,z)=f(x,y+d,z)-f(x,y,z);
where d is some small enough step
and f(x,y,z) is you current surface point computation
not sure why you use 3 input variables I would use just 2
but therefore if the shifted point is the same as unshifted
use this instead =f(x,y,z+d)-f(x,y,z);
at the end do not forget to normalize U,V size to unit vector

2. next step

if bullet 1 leads to correct normals
then you can simply solve the U,V algebraically
so rewrite U(x,y,z)=f(x+d,y,z)-f(x,y,z); to full equation
by substituting f(x,y,z) with the surface point equation
and simplify

[notes]

sometimes well selected d can simplify normalization to multipliyng by a constant
you should add normals visualization for example like this:
to actually see what is really happening (for debug purposses)