def visco_P():
    global P; global p; global pos; global N; global m
    ret = np.zeros([3, N])
    for i in range(N):
        for j in range(N):
            v_pis=v_pi(i,j)
            b=m*(2*np.sqrt(P[j]*P[i])/(p[j]*p[i])+v_pis)*1/2*(gradker(pos[:, i], pos[:, j], i) + gradker(pos[:, i], pos[:, j], j))
            ret[:, i] += b
            ret[:, j] -= b
    return ret
#speed sound
def v_pis(i, j):
    global p; global pos; global vel; global N; global c
    global h;
    alpha=1
    beta=2
    nu=0.01
    pij=(p[j]+p[i])/2
    mu=0
    h=(h[i]+h[j])/2
    r=pos[:, i]-pos[:, j]
    v=vel[:, i]-vel[:, j]
    if(np.dot(r, v)<0):
        mu=np.dot(r, v)/(h*((r/h)**2+nu))
    return (-alpha*mu*c+beta*(mu**2))/pij

def grad_v():
    global P; global p; global pos; global vel; global N; global m;
    grad = np.zeros([1, N])
    for i in range(N):
        for j in range(i):
            vij=vel[:,i]-vel[:, j]
            b1=-1/p[i]*m*1/2*(np.dot(gradker(pos[:, i], pos[:, j], i), vel) + np.dot(gradker(pos[:, i], pos[:, j], j), vel))
            b2=-1/p[i]*m*1/2*(np.dot(gradker(pos[:, j], pos[:, i], i), -vel) + np.dot(gradker(pos[:, j], pos[:, i], j)), -vel)
            grad[i] += b1
            grad[j] += b2
    return grad



def v_pis_al(i,j):
    global P; global p; global gradv; global pos; global vel; global N; global m;
    global c; global h;
    qi=0
    qj=0
    alpha=1
    beta=2
    if gradv[i]<0:
        qi=alpha*h[i]*p[i]*c*abs(gradv[i])+ beta*(h[i]**2)*p[i]*(gradv[i]**2)
        qi=qi/(p[i]**2)
    if gradv[j]<0:
        qj=alpha*h[j]*p[j]*c*abs(gradv[j])+ beta*(h[j]**2)*p[j]*(gradv[j]**2)
        qj=qj/(p[j]**2)
    return qi+qj