Reproducing plane waves with Huygens principle and phase Question:
In my physics class we were asked to derive Snell's law from Huygens principle in the context of this hypothetical situation:
There are plane waves moving through the ocean in a direction parallel to our x-axis. If I draw a line $x=c$, then I can see the state of the water (its oscillation) at any point and at any time along that line, then approximately re-create the original plane wave using $n$ point sources stretched between two points on the line $x=c$. If the original waves are parallel to $x=c$ then the waves created by each point source will be in phase with each other. But if the waves are moving at an angle $\theta$ relative to the $x$-axis, then the phase $\phi$ of the waves at each point source is not the same. I needed to find $\phi$.
After working with classmates I got the answer $\phi = k n \sin(\theta)$ and was able to use this to derive Snell's law. But I don't understand why. I was able to "work it backwards" to sort of understand why that solution works, but I want to know how to get there, and that is something that my classmates could not explain to me. Can someone please help me understand how to find $\phi$ as a function of $\theta$?
My backwards working of the solution:
Let $y'$ be the distance between crests along a cross-section parallel to the $y$-axis of the plane wave and $λ \le y' < \infty$.
$$k = 2\pi / \lambda$$
$$y' = \lambda / \sin(\theta)$$
$$\phi = k n \sin(\theta) = (2\pi / \lambda) n \sin(\theta) = (2\pi / y') n$$
As $y'$ increases towards infinity, the plane waves become more and more parallel to $x=c$, and so $\phi$ gets closer and closer to zero. As $y'$ decreases towards $\lambda$, $\phi$ increases towards its maximum value (though I don't know how to express/find its maximum value). So I can see that $\phi$ is inversely proportional to $y'$, and so also to $\sin(\theta)/k$. I understand that the wave is periodic, and multiplying by $2\pi$ is used often to convert to angular units, so that makes sense. And $n$ is the number of point sources used, so that it would be involved makes sense. But beyond that I am lost. How can I get to $\phi = k n \sin(\theta)$ working it "the right way around"?
I already tried asking classmates, reading our text books, and googling/searching on this site, but everything I find is either too advanced, too simple, or not related.
Edit:
Here is a photo of some revised work based on @mikuszefski's answer. It makes more sense now, but I'm still not sure why or how the $k$ comes in.
 A: This is not complete, but I would start like this:

If all is in phase the wavefront is parallel to $x$. If we now want an angle $\theta$ we have $d=\xi \sin\theta$. The wave propagates with $\cos(\omega t-k r)$, and we may set $t=0$. So if we want the point at the end of $d$ being in phase the source at $\xi$ must be at phase $\phi=k d$. (According to OP's comment/question, let me elaborate this point). The left point in the sketch I set as phase $\phi=0$. I am free to do so and I decide to measure all with respect to this point. This point, hence, also defines the wave front.
Now the point at $\xi$ gives me a wave $A \cos(k r)$, were $r =\sqrt{ (x-\xi)^2 +y^2}$, but wait: I want it to have a phase so it is $A \cos(k r-\phi(x))$. Now I want this wave to have a phase zero at the wave front (sort of by definition). In other words the argument of the cosine must be zero for $r=d$, i.e. $k d(\xi)-\phi(\xi)=0$. Only then this is part of the wavefront. So $k d(\xi)=\phi(\xi)$ and plugging $d=\xi \sin\theta$ gives, hence, 
$\phi=k \xi \sin\theta$.
Here a small python test for your point approximation:

First image is just a single point source. Second shows 100 sources in phase. They are distributed between $-20 \le x\le 20$. Last image introduces in $x$ dependent phase.Phases are calculated for $15^\circ$. (Note, here the color of the sources is just showing a change in phase, you could color it properly and use $\phi(x)\mod 2\pi$)
(code comes here, should be self-explaining)
import matplotlib
matplotlib.use('Qt4Agg')
import matplotlib.pyplot as plt
import numpy as np

def f(x,y,phi=0): return np.cos(2*np.pi*np.sqrt(x**2+y**2)-phi)

def wavephase(x,y,theta=0,pnts=100):
    out=0
    for pt in np.linspace(-20,20,pnts):
        out+=f(x-pt,y,2*np.pi*pt*np.sin(theta))
    return out

degree=np.pi/180.

n = 120
x = np.linspace(-10,10,2*n)
y = np.linspace(-10,0,n)
X,Y = np.meshgrid(x,y)

fff=15
fig=plt.figure()
ax=fig.add_subplot(3,1,1)
bx=fig.add_subplot(3,1,2)
cx=fig.add_subplot(3,1,3)
ax.contourf(X, Y, f(X,Y), 22, alpha=1, cmap="Blues")
bx.contourf(X, Y, wavephase(X,Y, theta=0*degree), 22, alpha=1, cmap="Blues")
bx.scatter(np.linspace(-20,20,100)[25:-25],np.linspace(-20,20,100)[25:-25]*0,c='b',s=45)
cx.contourf(X, Y, wavephase(X,Y, theta=fff*degree), 22, alpha=1, cmap="Blues")
cx.plot([0,10*np.cos(fff*degree)],[0,-10*np.sin(fff*degree)],color='#ffaa00')
cx.scatter(np.linspace(-20,20,100)[25:-25],np.linspace(-20,20,100)[25:-25]*0,c=np.linspace(-20,20,100)[25:-25],s=45)
for ma in [ax,bx,cx]:
    ma.set_xlim([-10,10])
    ma.set_ylim([-10,0])
plt.show()

A: I am not sure if this the correct diagram but what seems to be missing from your explanation is the separation of the point sources?

The phase, $\phi$ ,  between each adjacent point source is found from 
$\dfrac{\phi}{2 \pi}=\dfrac{d\sin \theta}{\lambda}$
and so $\phi_{\rm n}= \dfrac{2 \pi d \sin \theta}{\lambda} n$
