/* first puncture */
xp = xx - par_b;
yp = yy;
zp = zz;
rp = sqrt (xp*xp + yp*yp + zp*zp);
rp = pow (pow (rp, 4) + pow (TP_epsilon, 4), 0.25);
if (rp < TP_Tiny)
rp = TP_Tiny;
ir = 1.0/rp;
if (rp < TP_Extend_Radius) {
ir = EXTEND(1., rp);
}
s1 = 0.5* *mp *ir;
s3 = -s1*ir*ir;
s5 = -3.0*s3*ir*ir;
p += s1;
px += xp*s3;
py += yp*s3;
pz += zp*s3;
pxx += xp*xp*s5 + s3;
pxy += xp*yp*s5;
pxz += xp*zp*s5;
pyy += yp*yp*s5 + s3;
pyz += yp*zp*s5;
pzz += zp*zp*s5 + s3;
/* second puncture */
xp = xx + par_b;
yp = yy;
zp = zz;
rp = sqrt (xp*xp + yp*yp + zp*zp);
rp = pow (pow (rp, 4) + pow (TP_epsilon, 4), 0.25);
if (rp < TP_Tiny)
rp = TP_Tiny;
ir = 1.0/rp;
if (rp < TP_Extend_Radius) {
ir = EXTEND(1., rp);
}
s1 = 0.5* *mm *ir;
s3 = -s1*ir*ir;
s5 = -3.0*s3*ir*ir;
p += s1;
px += xp*s3;
py += yp*s3;
pz += zp*s3;
pxx += xp*xp*s5 + s3;
pxy += xp*yp*s5;
pxz += xp*zp*s5;
pyy += yp*yp*s5 + s3;
pyz += yp*zp*s5;
pzz += zp*zp*s5 + s3;
if (*conformal_state >= 1) {
static_psi = p;
psi[ind] = static_psi;
}
if (*conformal_state >= 2) {
psix[ind] = px / static_psi;
psiy[ind] = py / static_psi;
psiz[ind] = pz / static_psi;
}
if (*conformal_state >= 3) {
psixx[ind] = pxx / static_psi;
psixy[ind] = pxy / static_psi;
psixz[ind] = pxz / static_psi;
psiyy[ind] = pyy / static_psi;
psiyz[ind] = pyz / static_psi;
psizz[ind] = pzz / static_psi;
}
if (pmn_lapse)
alp[ind] = pow(p, initial_lapse_psi_exponent);
if (brownsville_lapse)
alp[ind] = 2.0/(1.0+pow(p, initial_lapse_psi_exponent));
} /* if conformal-state > 0 */