input

sph1.q

/VARI  xi, yi, zi, vxi, vyi, vzi, hi;
/VARJ  xj, yj, zj, vxj, vyj, vzj, hj, mj;
/VARF  rho, rox, roy, roz, dd;

dx = xi - xj;
dy = yi - yj;
dz = zi - zj;
r2 = dx*dx + dy*dy + dz*dz;
r1  = sqrt(r2);

hh  = 0.5*(hi+hj);
hhi = powm11(hh);
q   = r1*hhi;
qi = powm11(q);

######################## Kernel
## 1 < q <= 2
d1 = 2.0 - q;
d12 = d1*d1;
k1 = 0.25*d12*d1;

## 0 <= q <= 1
d2 = q*q;
k2 = 1.0 - 1.5*d2 + 0.75*d2*q;

kk1 = mux(q, 2.0, 0.0, k1);
ker = mux(q, 1.0, kk1, k2);

######################## dKernel
## 1 < q <= 2
dk1 = -0.75*d12*qi;

## 2/3 < q <= 1
dk2 = -3.0 + 2.25*q;

## 0 < q <= 2/3
dk3 = -1.0*qi;

dkk1 = mux(q, 2.0, 0.0,  dk1);
dkk2 = mux(q, 1.0, dkk1, dk2);
dker = mux(q, 0.66666666666666667, dkk2, dk3);

hhi2  = hhi*hhi;
mhhi2 = mj*hhi2;
mhhi3 = mhhi2*hhi;
mhhi5 = mhhi3*hhi2;

## density
rho = ker*mhhi3;

#############
dvx = vxi - vxj;
dvy = vyi - vyj;
dvz = vzi - vzj;

dum1 = dker*mhhi5;
## rotation
rox = (dvy*dz - dvz*dy)*dum1;
roy = (dvz*dx - dvx*dz)*dum1;
roz = (dvx*dy - dvy*dx)*dum1;

## time derivative of density
xv = dx*dvx + dy*dvy + dz*dvz;
dd = xv*dum1;