function [zeta,gradzeta] = myzeta(h,r)

q=abs(r)/h;

if(q>=2.0)
    zeta = 0.0;
    gradzeta = 0.0;
else
    fact = 1/h; % normalization factor
    dqdr = r/(h*abs(r));

    if(q>1.0)
        zeta = 1./6.*fact*(2.0-q)*(2.0-q)*(2.0-q);
        gradzeta =  -dqdr*fact*0.5*(2.0-q)*(2.0-q);
    else
        zeta = 2./3. * fact * (1.0 - 3./2.*q*q + 3./4.*q*q*q);
        gradzeta = dqdr*2./3. * fact*(-3.0*q + 9./4.*q*q);
    end
end