https://github.com/N-BodyShop/changa
Tip revision: 681e5936cd310d517ef6e80f4c247efae2da03db authored by Tim Haines on 10 August 2016, 22:04:28 UTC
Make single-precision for gravity calculations configurable
Make single-precision for gravity calculations configurable
Tip revision: 681e593
moments.c
/*
* Mulitpole moment routines.
* Main author: Joachim Stadel as part of PKGRAV
*/
#include "moments.h"
#include <assert.h>
#include <math.h>
#include <stdio.h>
void momClearLocr(LOCR *l) {
l->m = 0;
l->x = 0;
l->y = 0;
l->z = 0;
l->xx = 0;
l->xy = 0;
l->yy = 0;
l->xz = 0;
l->yz = 0;
l->xxx = 0;
l->xxy = 0;
l->xyy = 0;
l->yyy = 0;
l->xxz = 0;
l->xyz = 0;
l->yyz = 0;
l->xxxx = 0;
l->xxxy = 0;
l->xxyy = 0;
l->xyyy = 0;
l->yyyy = 0;
l->xxxz = 0;
l->xxyz = 0;
l->xyyz = 0;
l->yyyz = 0;
l->xxxxx = 0;
l->xxxxy = 0;
l->xxxyy = 0;
l->xxyyy = 0;
l->xyyyy = 0;
l->yyyyy = 0;
l->xxxxz = 0;
l->xxxyz = 0;
l->xxyyz = 0;
l->xyyyz = 0;
l->yyyyz = 0;
}
void momClearMomr(MOMR *mr) {
mr->m = 0.0;
mr->xx = 0.0;
mr->yy = 0.0;
mr->xy = 0.0;
mr->xz = 0.0;
mr->yz = 0.0;
mr->xxx = 0.0;
mr->xyy = 0.0;
mr->xxy = 0.0;
mr->yyy = 0.0;
mr->xxz = 0.0;
mr->yyz = 0.0;
mr->xyz = 0.0;
mr->xxxx = 0.0;
mr->xyyy = 0.0;
mr->xxxy = 0.0;
mr->yyyy = 0.0;
mr->xxxz = 0.0;
mr->yyyz = 0.0;
mr->xxyy = 0.0;
mr->xxyz = 0.0;
mr->xyyz = 0.0;
}
void momClearFmomr(FMOMR *l) {
l->m = 0;
l->xx = 0;
l->yy = 0;
l->xy = 0;
l->xz = 0;
l->yz = 0;
l->xxx = 0;
l->xyy = 0;
l->xxy = 0;
l->yyy = 0;
l->xxz = 0;
l->yyz = 0;
l->xyz = 0;
l->xxxx = 0;
l->xyyy = 0;
l->xxxy = 0;
l->yyyy = 0;
l->xxxz = 0;
l->yyyz = 0;
l->xxyy = 0;
l->xxyz = 0;
l->xyyz = 0;
}
/*
** This function adds the complete moment ma to the complete moment mc
*/
void momAddMomc(MOMC *mc, MOMC *ma) {
mc->m += ma->m;
mc->xx += ma->xx;
mc->yy += ma->yy;
mc->xy += ma->xy;
mc->xz += ma->xz;
mc->yz += ma->yz;
mc->xxx += ma->xxx;
mc->xyy += ma->xyy;
mc->xxy += ma->xxy;
mc->yyy += ma->yyy;
mc->xxz += ma->xxz;
mc->yyz += ma->yyz;
mc->xyz += ma->xyz;
mc->xxxx += ma->xxxx;
mc->xyyy += ma->xyyy;
mc->xxxy += ma->xxxy;
mc->yyyy += ma->yyyy;
mc->xxxz += ma->xxxz;
mc->yyyz += ma->yyyz;
mc->xxyy += ma->xxyy;
mc->xxyz += ma->xxyz;
mc->xyyz += ma->xyyz;
mc->zz += ma->zz;
mc->xzz += ma->xzz;
mc->yzz += ma->yzz;
mc->zzz += ma->zzz;
mc->xxzz += ma->xxzz;
mc->xyzz += ma->xyzz;
mc->xzzz += ma->xzzz;
mc->yyzz += ma->yyzz;
mc->yzzz += ma->yzzz;
mc->zzzz += ma->zzzz;
}
/*
** This function adds the reduced moment ma to the reduced moment mc
*/
void momAddMomr(MOMR *mr, MOMR *ma) {
mr->m += ma->m;
mr->xx += ma->xx;
mr->yy += ma->yy;
mr->xy += ma->xy;
mr->xz += ma->xz;
mr->yz += ma->yz;
mr->xxx += ma->xxx;
mr->xyy += ma->xyy;
mr->xxy += ma->xxy;
mr->yyy += ma->yyy;
mr->xxz += ma->xxz;
mr->yyz += ma->yyz;
mr->xyz += ma->xyz;
mr->xxxx += ma->xxxx;
mr->xyyy += ma->xyyy;
mr->xxxy += ma->xxxy;
mr->yyyy += ma->yyyy;
mr->xxxz += ma->xxxz;
mr->yyyz += ma->yyyz;
mr->xxyy += ma->xxyy;
mr->xxyz += ma->xxyz;
mr->xyyz += ma->xyyz;
}
/*
** This function adds the reduced moment ma to the reduced moment mr.
*/
void momAddFmomr(FMOMR *mr, FMOMR *ma) {
mr->m += ma->m;
mr->xx += ma->xx;
mr->yy += ma->yy;
mr->xy += ma->xy;
mr->xz += ma->xz;
mr->yz += ma->yz;
mr->xxx += ma->xxx;
mr->xyy += ma->xyy;
mr->xxy += ma->xxy;
mr->yyy += ma->yyy;
mr->xxz += ma->xxz;
mr->yyz += ma->yyz;
mr->xyz += ma->xyz;
mr->xxxx += ma->xxxx;
mr->xyyy += ma->xyyy;
mr->xxxy += ma->xxxy;
mr->yyyy += ma->yyyy;
mr->xxxz += ma->xxxz;
mr->yyyz += ma->yyyz;
mr->xxyy += ma->xxyy;
mr->xxyz += ma->xxyz;
mr->xyyz += ma->xyyz;
}
/*
** This function adds the reduced scaled moment ma to the reduced scaled moment
*mr.
** It needs to correctly rescale the moments of ma to be compatible with the
*scaling of mr.
*/
void momScaledAddFmomr(FMOMR *mr, cosmoType ur, FMOMR *ma, cosmoType ua) {
cosmoType f, s;
assert(ur > 0.0 && ua > 0);
f = ua / ur;
s = f;
mr->m += ma->m;
s *= f;
mr->xx += s * ma->xx;
mr->yy += s * ma->yy;
mr->xy += s * ma->xy;
mr->xz += s * ma->xz;
mr->yz += s * ma->yz;
s *= f;
mr->xxx += s * ma->xxx;
mr->xyy += s * ma->xyy;
mr->xxy += s * ma->xxy;
mr->yyy += s * ma->yyy;
mr->xxz += s * ma->xxz;
mr->yyz += s * ma->yyz;
mr->xyz += s * ma->xyz;
s *= f;
mr->xxxx += s * ma->xxxx;
mr->xyyy += s * ma->xyyy;
mr->xxxy += s * ma->xxxy;
mr->yyyy += s * ma->yyyy;
mr->xxxz += s * ma->xxxz;
mr->yyyz += s * ma->yyyz;
mr->xxyy += s * ma->xxyy;
mr->xxyz += s * ma->xxyz;
mr->xyyz += s * ma->xyyz;
}
/*
** This function rescales reduced scaled moment mr.
*/
void momRescaleFmomr(FMOMR *mr, cosmoType unew, cosmoType uold) {
cosmoType f, s;
f = uold / unew;
s = f;
s *= f;
mr->xx *= s;
mr->yy *= s;
mr->xy *= s;
mr->xz *= s;
mr->yz *= s;
s *= f;
mr->xxx *= s;
mr->xyy *= s;
mr->xxy *= s;
mr->yyy *= s;
mr->xxz *= s;
mr->yyz *= s;
mr->xyz *= s;
s *= f;
mr->xxxx *= s;
mr->xyyy *= s;
mr->xxxy *= s;
mr->yyyy *= s;
mr->xxxz *= s;
mr->yyyz *= s;
mr->xxyy *= s;
mr->xxyz *= s;
mr->xyyz *= s;
}
/*
** This function multiply-adds the complete moment ma
*/
void momMulAddMomc(MOMC *mc, double m, MOMC *ma) {
mc->m += m * ma->m;
mc->xx += m * ma->xx;
mc->yy += m * ma->yy;
mc->xy += m * ma->xy;
mc->xz += m * ma->xz;
mc->yz += m * ma->yz;
mc->xxx += m * ma->xxx;
mc->xyy += m * ma->xyy;
mc->xxy += m * ma->xxy;
mc->yyy += m * ma->yyy;
mc->xxz += m * ma->xxz;
mc->yyz += m * ma->yyz;
mc->xyz += m * ma->xyz;
mc->xxxx += m * ma->xxxx;
mc->xyyy += m * ma->xyyy;
mc->xxxy += m * ma->xxxy;
mc->yyyy += m * ma->yyyy;
mc->xxxz += m * ma->xxxz;
mc->yyyz += m * ma->yyyz;
mc->xxyy += m * ma->xxyy;
mc->xxyz += m * ma->xxyz;
mc->xyyz += m * ma->xyyz;
mc->zz += m * ma->zz;
mc->xzz += m * ma->xzz;
mc->yzz += m * ma->yzz;
mc->zzz += m * ma->zzz;
mc->xxzz += m * ma->xxzz;
mc->xyzz += m * ma->xyzz;
mc->xzzz += m * ma->xzzz;
mc->yyzz += m * ma->yyzz;
mc->yzzz += m * ma->yzzz;
mc->zzzz += m * ma->zzzz;
}
/*
** This function multiply-adds the reduced moment ma
*/
void momMulAddMomr(MOMR *mr, double m, MOMR *ma) {
mr->m += m * ma->m;
mr->xx += m * ma->xx;
mr->yy += m * ma->yy;
mr->xy += m * ma->xy;
mr->xz += m * ma->xz;
mr->yz += m * ma->yz;
mr->xxx += m * ma->xxx;
mr->xyy += m * ma->xyy;
mr->xxy += m * ma->xxy;
mr->yyy += m * ma->yyy;
mr->xxz += m * ma->xxz;
mr->yyz += m * ma->yyz;
mr->xyz += m * ma->xyz;
mr->xxxx += m * ma->xxxx;
mr->xyyy += m * ma->xyyy;
mr->xxxy += m * ma->xxxy;
mr->yyyy += m * ma->yyyy;
mr->xxxz += m * ma->xxxz;
mr->yyyz += m * ma->yyyz;
mr->xxyy += m * ma->xxyy;
mr->xxyz += m * ma->xxyz;
mr->xyyz += m * ma->xyyz;
}
/*
** This function multiply-adds the reduced scaled moment ma
*/
void momMulAddFmomr(FMOMR *mr, cosmoType ur, cosmoType m, FMOMR *ma,
cosmoType ua) {
cosmoType f;
assert(ua > 0.0 && ur > 0.0);
f = ua / ur;
mr->m += m * ma->m;
m *= f;
m *= f;
mr->xx += m * ma->xx;
mr->yy += m * ma->yy;
mr->xy += m * ma->xy;
mr->xz += m * ma->xz;
mr->yz += m * ma->yz;
m *= f;
mr->xxx += m * ma->xxx;
mr->xyy += m * ma->xyy;
mr->xxy += m * ma->xxy;
mr->yyy += m * ma->yyy;
mr->xxz += m * ma->xxz;
mr->yyz += m * ma->yyz;
mr->xyz += m * ma->xyz;
m *= f;
mr->xxxx += m * ma->xxxx;
mr->xyyy += m * ma->xyyy;
mr->xxxy += m * ma->xxxy;
mr->yyyy += m * ma->yyyy;
mr->xxxz += m * ma->xxxz;
mr->yyyz += m * ma->yyyz;
mr->xxyy += m * ma->xxyy;
mr->xxyz += m * ma->xxyz;
mr->xyyz += m * ma->xyyz;
}
/*
** This function adds the reduced local expansion la to the reduced local
*explansion lr.
*/
void momAddFlocr(FLOCR *lr, FLOCR *la) {
lr->m += la->m;
lr->x += la->x;
lr->y += la->y;
lr->z += la->z;
lr->xx += la->xx;
lr->yy += la->yy;
lr->xy += la->xy;
lr->xz += la->xz;
lr->yz += la->yz;
lr->xxx += la->xxx;
lr->xyy += la->xyy;
lr->xxy += la->xxy;
lr->yyy += la->yyy;
lr->xxz += la->xxz;
lr->yyz += la->yyz;
lr->xyz += la->xyz;
lr->xxxx += la->xxxx;
lr->xyyy += la->xyyy;
lr->xxxy += la->xxxy;
lr->yyyy += la->yyyy;
lr->xxxz += la->xxxz;
lr->yyyz += la->yyyz;
lr->xxyy += la->xxyy;
lr->xxyz += la->xxyz;
lr->xyyz += la->xyyz;
lr->xxxxx += la->xxxxx;
lr->xyyyy += la->xyyyy;
lr->xxxxy += la->xxxxy;
lr->yyyyy += la->yyyyy;
lr->xxxxz += la->xxxxz;
lr->yyyyz += la->yyyyz;
lr->xxxyy += la->xxxyy;
lr->xxyyy += la->xxyyy;
lr->xxxyz += la->xxxyz;
lr->xyyyz += la->xyyyz;
lr->xxyyz += la->xxyyz;
}
/*
** This function adds the reduced scaled local expansion la to the reduced
*scaled local expansion lr.
** It needs to correctly rescale the elements of la to be compatible with the
*scaling of lr.
*/
void momScaledAddFlocr(FLOCR *lr, cosmoType vr, FLOCR *la, cosmoType va) {
cosmoType f, s;
assert(vr > 0.0 && va > 0);
f = va / vr;
s = f;
lr->m += la->m;
lr->x += s * la->x;
lr->y += s * la->y;
lr->z += s * la->z;
s *= f;
lr->xx += s * la->xx;
lr->yy += s * la->yy;
lr->xy += s * la->xy;
lr->xz += s * la->xz;
lr->yz += s * la->yz;
s *= f;
lr->xxx += s * la->xxx;
lr->xyy += s * la->xyy;
lr->xxy += s * la->xxy;
lr->yyy += s * la->yyy;
lr->xxz += s * la->xxz;
lr->yyz += s * la->yyz;
lr->xyz += s * la->xyz;
s *= f;
lr->xxxx += s * la->xxxx;
lr->xyyy += s * la->xyyy;
lr->xxxy += s * la->xxxy;
lr->yyyy += s * la->yyyy;
lr->xxxz += s * la->xxxz;
lr->yyyz += s * la->yyyz;
lr->xxyy += s * la->xxyy;
lr->xxyz += s * la->xxyz;
lr->xyyz += s * la->xyyz;
s *= f;
lr->xxxxx += s * la->xxxxx;
lr->xyyyy += s * la->xyyyy;
lr->xxxxy += s * la->xxxxy;
lr->yyyyy += s * la->yyyyy;
lr->xxxxz += s * la->xxxxz;
lr->yyyyz += s * la->yyyyz;
lr->xxxyy += s * la->xxxyy;
lr->xxyyy += s * la->xxyyy;
lr->xxxyz += s * la->xxxyz;
lr->xyyyz += s * la->xyyyz;
lr->xxyyz += s * la->xxyyz;
}
/*
** This function rescales the reduced scaled local expansion lr.
*/
void momRescaleFlocr(FLOCR *lr, cosmoType vnew, cosmoType vold) {
cosmoType f, s;
assert(vnew > 0.0 && vold > 0);
f = vnew / vold;
s = f;
lr->x *= s;
lr->y *= s;
lr->z *= s;
s *= f;
lr->xx *= s;
lr->yy *= s;
lr->xy *= s;
lr->xz *= s;
lr->yz *= s;
s *= f;
lr->xxx *= s;
lr->xyy *= s;
lr->xxy *= s;
lr->yyy *= s;
lr->xxz *= s;
lr->yyz *= s;
lr->xyz *= s;
s *= f;
lr->xxxx *= s;
lr->xyyy *= s;
lr->xxxy *= s;
lr->yyyy *= s;
lr->xxxz *= s;
lr->yyyz *= s;
lr->xxyy *= s;
lr->xxyz *= s;
lr->xyyz *= s;
s *= f;
lr->xxxxx *= s;
lr->xyyyy *= s;
lr->xxxxy *= s;
lr->yyyyy *= s;
lr->xxxxz *= s;
lr->yyyyz *= s;
lr->xxxyy *= s;
lr->xxyyy *= s;
lr->xxxyz *= s;
lr->xyyyz *= s;
lr->xxyyz *= s;
}
/*
** This function subtracts the complete moment ma from the complete moment mc
*/
void momSubMomc(MOMC *mc, MOMC *ma) {
mc->m -= ma->m;
mc->xx -= ma->xx;
mc->yy -= ma->yy;
mc->xy -= ma->xy;
mc->xz -= ma->xz;
mc->yz -= ma->yz;
mc->xxx -= ma->xxx;
mc->xyy -= ma->xyy;
mc->xxy -= ma->xxy;
mc->yyy -= ma->yyy;
mc->xxz -= ma->xxz;
mc->yyz -= ma->yyz;
mc->xyz -= ma->xyz;
mc->xxxx -= ma->xxxx;
mc->xyyy -= ma->xyyy;
mc->xxxy -= ma->xxxy;
mc->yyyy -= ma->yyyy;
mc->xxxz -= ma->xxxz;
mc->yyyz -= ma->yyyz;
mc->xxyy -= ma->xxyy;
mc->xxyz -= ma->xxyz;
mc->xyyz -= ma->xyyz;
mc->zz -= ma->zz;
mc->xzz -= ma->xzz;
mc->yzz -= ma->yzz;
mc->zzz -= ma->zzz;
mc->xxzz -= ma->xxzz;
mc->xyzz -= ma->xyzz;
mc->xzzz -= ma->xzzz;
mc->yyzz -= ma->yyzz;
mc->yzzz -= ma->yzzz;
mc->zzzz -= ma->zzzz;
}
/*
** This function subtracts the reduced moment ma from the reduced moment mr
*/
void momSubMomr(MOMR *mr, MOMR *ma) {
mr->m -= ma->m;
mr->xx -= ma->xx;
mr->yy -= ma->yy;
mr->xy -= ma->xy;
mr->xz -= ma->xz;
mr->yz -= ma->yz;
mr->xxx -= ma->xxx;
mr->xyy -= ma->xyy;
mr->xxy -= ma->xxy;
mr->yyy -= ma->yyy;
mr->xxz -= ma->xxz;
mr->yyz -= ma->yyz;
mr->xyz -= ma->xyz;
mr->xxxx -= ma->xxxx;
mr->xyyy -= ma->xyyy;
mr->xxxy -= ma->xxxy;
mr->yyyy -= ma->yyyy;
mr->xxxz -= ma->xxxz;
mr->yyyz -= ma->yyyz;
mr->xxyy -= ma->xxyy;
mr->xxyz -= ma->xxyz;
mr->xyyz -= ma->xyyz;
}
/*
** This function subtracts the reduced scaled moment ma to the
** reduced scaled moment mr.
** It needs to correctly rescale the moments of ma to be compatible
** with the scaling of mr.
*/
void momScaledSubFmomr(FMOMR *mr, cosmoType ur, FMOMR *ma, cosmoType ua) {
cosmoType f, s;
assert(ur > 0.0 && ua > 0);
f = ua / ur;
s = f;
mr->m -= ma->m;
s *= f;
mr->xx -= s * ma->xx;
mr->yy -= s * ma->yy;
mr->xy -= s * ma->xy;
mr->xz -= s * ma->xz;
mr->yz -= s * ma->yz;
s *= f;
mr->xxx -= s * ma->xxx;
mr->xyy -= s * ma->xyy;
mr->xxy -= s * ma->xxy;
mr->yyy -= s * ma->yyy;
mr->xxz -= s * ma->xxz;
mr->yyz -= s * ma->yyz;
mr->xyz -= s * ma->xyz;
s *= f;
mr->xxxx -= s * ma->xxxx;
mr->xyyy -= s * ma->xyyy;
mr->xxxy -= s * ma->xxxy;
mr->yyyy -= s * ma->yyyy;
mr->xxxz -= s * ma->xxxz;
mr->yyyz -= s * ma->yyyz;
mr->xxyy -= s * ma->xxyy;
mr->xxyz -= s * ma->xxyz;
mr->xyyz -= s * ma->xyyz;
}
/*
** This function calculates a complete multipole from a single
** particle at position <x,y,z> from the center of mass.
** The strange order of evaluation reduces the number of
** multiplications to a minimum.
** <x,y,z> := d := r(particle) - rcm.
**
** OpCount (*,+) = (34,0)
**
*/
void momMakeMomc(MOMC *mc, double m, double x, double y, double z) {
double tx, ty, tz;
mc->m = m;
/*
** Calculate the Quadrupole Moment.
*/
tx = m * x;
ty = m * y;
mc->xy = tx * y;
mc->xz = tx * z;
mc->yz = ty * z;
tx *= x;
ty *= y;
tz = m * z * z;
mc->xx = tx;
mc->yy = ty;
mc->zz = tz;
/*
** Calculate the Octopole Moment.
*/
mc->xxy = tx * y;
mc->xxz = tx * z;
mc->yyz = ty * z;
mc->xyy = ty * x;
mc->xzz = tz * x;
mc->yzz = tz * y;
mc->xyz = mc->xy * z;
tx *= x;
ty *= y;
tz *= z;
mc->xxx = tx;
mc->yyy = ty;
mc->zzz = tz;
/*
** Calculate the Hexadecapole Moment.
*/
mc->xxxx = tx * x;
mc->xxxy = tx * y;
mc->xxxz = tx * z;
mc->xyyy = ty * x;
mc->yyyy = ty * y;
mc->yyyz = ty * z;
mc->xzzz = tz * x;
mc->yzzz = tz * y;
mc->zzzz = tz * z;
mc->xxyy = mc->xxy * y;
mc->xxyz = mc->xxy * z;
mc->xyyz = mc->yyz * x;
mc->yyzz = mc->yyz * z;
mc->xxzz = mc->xzz * x;
mc->xyzz = mc->xzz * y;
}
/*
** This function calculates a reduced multipole from a single
** particle at position <x,y,z> from the center of mass.
** The strange order of evaluation reduces the number of
** multiplications to a minimum.
** <x,y,z> := d := r(particle) - rcm.
** returns: d^2
**
** OpCount (*,+) = (43,18) = ~60
*/
double momMakeMomr(MOMR *mr, double m, double x, double y, double z) {
double tx, ty, t, dx, dy;
double x2 = x * x;
double y2 = y * y;
double d2 = x2 + y2 + z * z;
mr->m = m;
/*
** Calculate the Quadrupole Moment.
*/
tx = m * x;
ty = m * y;
mr->xy = tx * y;
mr->xz = tx * z;
mr->yz = ty * z;
tx *= x;
ty *= y;
m *= d2;
t = m / 3;
mr->xx = tx - t;
mr->yy = ty - t;
/*
** Calculate the Octopole Moment.
*/
t = 0.2 * m;
dx = tx - t;
dy = ty - t;
mr->xxy = dx * y;
mr->xxz = dx * z;
mr->yyz = dy * z;
mr->xyy = dy * x;
mr->xyz = mr->xy * z;
t *= 3;
mr->xxx = (tx - t) * x;
mr->yyy = (ty - t) * y;
/*
** Calculate the Hexadecapole Moment.
*/
t = m / 7;
mr->xxyz = (tx - t) * y * z;
mr->xyyz = (ty - t) * x * z;
dx = (tx - 3 * t) * x;
dy = (ty - 3 * t) * y;
mr->xxxy = dx * y;
mr->xxxz = dx * z;
mr->xyyy = dy * x;
mr->yyyz = dy * z;
dx = t * (x2 - 0.1 * d2);
dy = t * (y2 - 0.1 * d2);
mr->xxxx = tx * x2 - 6 * dx;
mr->yyyy = ty * y2 - 6 * dy;
mr->xxyy = tx * y2 - dx - dy;
return (d2);
}
/*
** This function calculates a reduced scaled multipole from a single
** particle at position <x,y,z> from the any type of "center". A scaling
** factor 'u' for the positions must also be specified, which should
** typically be the value of b_max.
**
** The strange order of evaluation reduces the number of
** multiplications to a minimum.
** <x,y,z> := d := r(particle) - rcm.
** returns: d^2 scaled by u^2.
**
** OpCount (*,+) = (43,18) = ~60
*/
cosmoType momMakeFmomr(FMOMR *mr, cosmoType m, cosmoType u, cosmoType x,
cosmoType y, cosmoType z) {
cosmoType tx, ty, t, dx, dy;
cosmoType x2;
cosmoType y2;
cosmoType d2, iu;
assert(u > 0.0);
iu = 1.0f / u;
x *= iu;
y *= iu;
z *= iu;
x2 = x * x;
y2 = y * y;
d2 = x2 + y2 + z * z;
mr->m = m;
tx = m * x;
ty = m * y;
/*
** Calculate the Quadrupole Moment.
*/
mr->xy = tx * y;
mr->xz = tx * z;
mr->yz = ty * z;
tx *= x;
ty *= y;
m *= d2;
t = (1.0f / 3.0f) * m;
mr->xx = tx - t;
mr->yy = ty - t;
/*
** Calculate the Octopole Moment.
*/
t = 0.2f * m;
dx = tx - t;
dy = ty - t;
mr->xxy = dx * y;
mr->xxz = dx * z;
mr->yyz = dy * z;
mr->xyy = dy * x;
mr->xyz = mr->xy * z;
t *= 3.0f;
mr->xxx = (tx - t) * x;
mr->yyy = (ty - t) * y;
/*
** Calculate the Hexadecapole Moment.
*/
t = (1.0f / 7.0f) * m;
mr->xxyz = (tx - t) * y * z;
mr->xyyz = (ty - t) * x * z;
dx = (tx - 3.0f * t) * x;
dy = (ty - 3.0f * t) * y;
mr->xxxy = dx * y;
mr->xxxz = dx * z;
mr->xyyy = dy * x;
mr->yyyz = dy * z;
dx = t * (x2 - 0.1f * d2);
dy = t * (y2 - 0.1f * d2);
mr->xxxx = tx * x2 - 6.0f * dx;
mr->yyyy = ty * y2 - 6.0f * dy;
mr->xxyy = tx * y2 - dx - dy;
return (d2);
}
/*
** This function calculates a reduced multipole from a single
** particle at position <x,y,z> from the center of mass.
** This is the "straight forward" implementation which we
** used in the original version of PKDGRAV. It remains a good
** test of more peculiar looking code.
**
** <x,y,z> := d := r(particle) - rcm.
**
** OpCount (*,+) = (115,20) = ~135
*/
void momOldMakeMomr(MOMR *mr, double m, double x, double y, double z) {
double d2 = x * x + y * y + z * z;
mr->xxxx = m * (x * x * x * x - 6.0 / 7.0 * d2 * (x * x - 0.1 * d2));
mr->xyyy = m * (x * y * y * y - 3.0 / 7.0 * d2 * x * y);
mr->xxxy = m * (x * x * x * y - 3.0 / 7.0 * d2 * x * y);
mr->yyyy = m * (y * y * y * y - 6.0 / 7.0 * d2 * (y * y - 0.1 * d2));
mr->xxxz = m * (x * x * x * z - 3.0 / 7.0 * d2 * x * z);
mr->yyyz = m * (y * y * y * z - 3.0 / 7.0 * d2 * y * z);
mr->xxyy =
m * (x * x * y * y - 1.0 / 7.0 * d2 * (x * x + y * y - 0.2 * d2));
mr->xxyz = m * (x * x * y * z - 1.0 / 7.0 * d2 * y * z);
mr->xyyz = m * (x * y * y * z - 1.0 / 7.0 * d2 * x * z);
/*
** Calculate reduced octopole moment...
*/
mr->xxx = m * (x * x * x - 0.6 * d2 * x);
mr->xyy = m * (x * y * y - 0.2 * d2 * x);
mr->xxy = m * (x * x * y - 0.2 * d2 * y);
mr->yyy = m * (y * y * y - 0.6 * d2 * y);
mr->xxz = m * (x * x * z - 0.2 * d2 * z);
mr->yyz = m * (y * y * z - 0.2 * d2 * z);
mr->xyz = m * x * y * z;
/*
** Calculate quadrupole moment...
*/
mr->xx = m * (x * x - 1.0 / 3.0 * d2);
mr->yy = m * (y * y - 1.0 / 3.0 * d2);
mr->xy = m * x * y;
mr->xz = m * x * z;
mr->yz = m * y * z;
mr->m = m;
}
/*
** This function shifts a complete multipole (MOMC) to a new center of mass.
** <x,y,z> := d := rcm(old) - rcm(new).
**
** OpCount ShiftMomc (*,+) = (111,84)
** MakeMomc (*,+) = (34,0)
** MulAddMomc (*,+) = (32,32)
** Total (*,+) = (177,116) = 293
*/
void momShiftMomc(MOMC *m, double x, double y, double z) {
MOMC f;
momMakeMomc(&f, 1, x, y, z);
/*
** Shift the Hexadecapole.
*/
m->xxxx += 4 * m->xxx * x + 6 * m->xx * f.xx;
m->yyyy += 4 * m->yyy * y + 6 * m->yy * f.yy;
m->zzzz += 4 * m->zzz * z + 6 * m->zz * f.zz;
m->xyyy += m->yyy * x + 3 * (m->xyy * y + m->yy * f.xy + m->xy * f.yy);
m->xxxy += m->xxx * y + 3 * (m->xxy * x + m->xx * f.xy + m->xy * f.xx);
m->xxxz += m->xxx * z + 3 * (m->xxz * x + m->xx * f.xz + m->xz * f.xx);
m->yyyz += m->yyy * z + 3 * (m->yyz * y + m->yy * f.yz + m->yz * f.yy);
m->xzzz += m->zzz * x + 3 * (m->xzz * z + m->zz * f.xz + m->xz * f.zz);
m->yzzz += m->zzz * y + 3 * (m->yzz * z + m->zz * f.yz + m->yz * f.zz);
m->xxyy += 2 * (m->xxy * y + m->xyy * x) + m->xx * f.yy + m->yy * f.xx +
4 * m->xy * f.xy;
m->xxzz += 2 * (m->xxz * z + m->xzz * x) + m->xx * f.zz + m->zz * f.xx +
4 * m->xz * f.xz;
m->yyzz += 2 * (m->yyz * z + m->yzz * y) + m->yy * f.zz + m->zz * f.yy +
4 * m->yz * f.yz;
m->xxyz += m->xxy * z + m->xxz * y + m->xx * f.yz + m->yz * f.xx +
2 * (m->xyz * x + m->xy * f.xz + m->xz * f.xy);
m->xyyz += m->xyy * z + m->yyz * x + m->yy * f.xz + m->xz * f.yy +
2 * (m->xyz * y + m->xy * f.yz + m->yz * f.xy);
m->xyzz += m->yzz * x + m->xzz * y + m->zz * f.xy + m->xy * f.zz +
2 * (m->xyz * z + m->yz * f.xz + m->xz * f.yz);
/*
** Now shift the Octopole.
*/
m->xxx += 3 * m->xx * x;
m->yyy += 3 * m->yy * y;
m->zzz += 3 * m->zz * z;
m->xyy += 2 * m->xy * y + m->yy * x;
m->xzz += 2 * m->xz * z + m->zz * x;
m->yzz += 2 * m->yz * z + m->zz * y;
m->xxy += 2 * m->xy * x + m->xx * y;
m->xxz += 2 * m->xz * x + m->xx * z;
m->yyz += 2 * m->yz * y + m->yy * z;
m->xyz += m->xy * z + m->xz * y + m->yz * x;
/*
** Now deal with the monopole terms.
*/
f.m = 0;
momMulAddMomc(m, m->m, &f);
}
/*
** This function shifts a reduced multipole (MOMR) to a new center of mass.
** <x,y,z> := d := rcm(old) - rcm(new).
**
** OpCount ShiftMomr (*,+) = (128,111)
** MakeMomr (*,+) = (43,18)
** MulAddMomr (*,+) = (22,22)
** Total (*,+) = (193,151) = 344
*/
void momShiftMomr(MOMR *m, double x, double y, double z) {
MOMR f;
double t, tx, ty, tz, txx, tyy, txy, tyz, txz;
const double twosevenths = 2.0 / 7.0;
momMakeMomr(&f, 1, x, y, z);
/*
** Calculate the correction terms.
*/
tx = 0.4 * (m->xx * x + m->xy * y + m->xz * z);
ty = 0.4 * (m->xy * x + m->yy * y + m->yz * z);
tz = 0.4 * (m->xz * x + m->yz * y - (m->xx + m->yy) * z);
t = tx * x + ty * y + tz * z;
txx = twosevenths *
(m->xxx * x + m->xxy * y + m->xxz * z +
2 * (m->xx * f.xx + m->xy * f.xy + m->xz * f.xz) - 0.5 * t);
tyy = twosevenths *
(m->xyy * x + m->yyy * y + m->yyz * z +
2 * (m->xy * f.xy + m->yy * f.yy + m->yz * f.yz) - 0.5 * t);
txy = twosevenths *
(m->xxy * x + m->xyy * y + m->xyz * z + m->xy * (f.xx + f.yy) +
(m->xx + m->yy) * f.xy + m->yz * f.xz + m->xz * f.yz);
tyz = twosevenths *
(m->xyz * x + m->yyz * y - (m->xxy + m->yyy) * z - m->yz * f.xx -
m->xx * f.yz + m->xz * f.xy + m->xy * f.xz);
txz = twosevenths *
(m->xxz * x + m->xyz * y - (m->xxx + m->xyy) * z - m->xz * f.yy -
m->yy * f.xz + m->yz * f.xy + m->xy * f.yz);
/*
** Shift the Hexadecapole.
*/
m->xxxx += 4 * m->xxx * x + 6 * (m->xx * f.xx - txx);
m->yyyy += 4 * m->yyy * y + 6 * (m->yy * f.yy - tyy);
m->xyyy +=
m->yyy * x + 3 * (m->xyy * y + m->yy * f.xy + m->xy * f.yy - txy);
m->xxxy +=
m->xxx * y + 3 * (m->xxy * x + m->xx * f.xy + m->xy * f.xx - txy);
m->xxxz +=
m->xxx * z + 3 * (m->xxz * x + m->xx * f.xz + m->xz * f.xx - txz);
m->yyyz +=
m->yyy * z + 3 * (m->yyz * y + m->yy * f.yz + m->yz * f.yy - tyz);
m->xxyy += 2 * (m->xxy * y + m->xyy * x) + m->xx * f.yy + m->yy * f.xx +
4 * m->xy * f.xy - txx - tyy;
m->xxyz += m->xxy * z + m->xxz * y + m->xx * f.yz + m->yz * f.xx +
2 * (m->xyz * x + m->xy * f.xz + m->xz * f.xy) - tyz;
m->xyyz += m->xyy * z + m->yyz * x + m->yy * f.xz + m->xz * f.yy +
2 * (m->xyz * y + m->xy * f.yz + m->yz * f.xy) - txz;
/*
** Now shift the Octopole.
*/
m->xxx += 3 * (m->xx * x - tx);
m->xyy += 2 * m->xy * y + m->yy * x - tx;
m->yyy += 3 * (m->yy * y - ty);
m->xxy += 2 * m->xy * x + m->xx * y - ty;
m->xxz += 2 * m->xz * x + m->xx * z - tz;
m->yyz += 2 * m->yz * y + m->yy * z - tz;
m->xyz += m->xy * z + m->xz * y + m->yz * x;
/*
** Now deal with the monopole terms.
*/
f.m = 0;
momMulAddMomr(m, m->m, &f);
}
/*
** This function shifts a reduced scaled multipole (FMOMR) to a new center of
*mass.
** <x,y,z> := d := rcm(old) - rcm(new).
**
** OpCount ShiftMomr (*,+) = (128,111)
** MakeMomr (*,+) = (43,18)
** MulAddMomr (*,+) = (22,22)
** Total (*,+) = (193,151) = 344
*/
void momShiftFmomr(FMOMR *m, cosmoType u, cosmoType x, cosmoType y,
cosmoType z) {
FMOMR f;
cosmoType t, tx, ty, tz, txx, tyy, txy, tyz, txz, iu;
const cosmoType twosevenths = 2.0f / 7.0f;
momMakeFmomr(&f, 1.0f, u, x, y, z);
iu = 1.0f / u;
x *= iu;
y *= iu;
z *= iu;
/*
** Calculate the correction terms.
*/
tx = 0.4f * (m->xx * x + m->xy * y + m->xz * z);
ty = 0.4f * (m->xy * x + m->yy * y + m->yz * z);
tz = 0.4f * (m->xz * x + m->yz * y - (m->xx + m->yy) * z);
t = tx * x + ty * y + tz * z;
txx = twosevenths *
(m->xxx * x + m->xxy * y + m->xxz * z +
2.0f * (m->xx * f.xx + m->xy * f.xy + m->xz * f.xz) - 0.5f * t);
tyy = twosevenths *
(m->xyy * x + m->yyy * y + m->yyz * z +
2.0f * (m->xy * f.xy + m->yy * f.yy + m->yz * f.yz) - 0.5f * t);
txy = twosevenths *
(m->xxy * x + m->xyy * y + m->xyz * z + m->xy * (f.xx + f.yy) +
(m->xx + m->yy) * f.xy + m->yz * f.xz + m->xz * f.yz);
tyz = twosevenths *
(m->xyz * x + m->yyz * y - (m->xxy + m->yyy) * z - m->yz * f.xx -
m->xx * f.yz + m->xz * f.xy + m->xy * f.xz);
txz = twosevenths *
(m->xxz * x + m->xyz * y - (m->xxx + m->xyy) * z - m->xz * f.yy -
m->yy * f.xz + m->yz * f.xy + m->xy * f.yz);
/*
** Shift the Hexadecapole.
*/
m->xxxx += 4.0f * m->xxx * x + 6.0f * (m->xx * f.xx - txx);
m->yyyy += 4.0f * m->yyy * y + 6.0f * (m->yy * f.yy - tyy);
m->xyyy +=
m->yyy * x + 3.0f * (m->xyy * y + m->yy * f.xy + m->xy * f.yy - txy);
m->xxxy +=
m->xxx * y + 3.0f * (m->xxy * x + m->xx * f.xy + m->xy * f.xx - txy);
m->xxxz +=
m->xxx * z + 3.0f * (m->xxz * x + m->xx * f.xz + m->xz * f.xx - txz);
m->yyyz +=
m->yyy * z + 3.0f * (m->yyz * y + m->yy * f.yz + m->yz * f.yy - tyz);
m->xxyy += 2.0f * (m->xxy * y + m->xyy * x) + m->xx * f.yy + m->yy * f.xx +
4.0f * m->xy * f.xy - txx - tyy;
m->xxyz += m->xxy * z + m->xxz * y + m->xx * f.yz + m->yz * f.xx +
2.0f * (m->xyz * x + m->xy * f.xz + m->xz * f.xy) - tyz;
m->xyyz += m->xyy * z + m->yyz * x + m->yy * f.xz + m->xz * f.yy +
2.0f * (m->xyz * y + m->xy * f.yz + m->yz * f.xy) - txz;
/*
** Now shift the Octopole.
*/
m->xxx += 3.0f * (m->xx * x - tx);
m->xyy += 2.0f * m->xy * y + m->yy * x - tx;
m->yyy += 3.0f * (m->yy * y - ty);
m->xxy += 2.0f * m->xy * x + m->xx * y - ty;
m->xxz += 2.0f * m->xz * x + m->xx * z - tz;
m->yyz += 2.0f * m->yz * y + m->yy * z - tz;
m->xyz += m->xy * z + m->xz * y + m->yz * x;
/*
** Now deal with the monopole terms.
*/
f.m = 0;
momMulAddFmomr(m, 1.0, m->m, &f, 1.0);
}
/*
** This function shifts a reduced local expansion (LOCR) to a new center of
*expansion.
** <x,y,z> := d := rexp(new) - rexp(old).
**
** Op Count (*,+) = (159,173) = 332
*/
/* IBM brain damage */
#undef hz
double momShiftLocr(LOCR *l, double x, double y, double z) {
const double onethird = 1.0 / 3.0;
double hx, hy, hz, tx, ty, tz;
double L, Lx, Ly, Lz, Lxx, Lxy, Lxz, Lyy, Lyz, Lxxx, Lxxy, Lxxz, Lxyy, Lxyz,
Lyyy, Lyyz;
double Lxxxx, Lxxxy, Lxxyy, Lxyyy, Lyyyy, Lxxxz, Lxxyz, Lxyyz, Lyyyz;
hx = 0.5 * x;
hy = 0.5 * y;
hz = 0.5 * z;
tx = onethird * x;
ty = onethird * y;
tz = onethird * z;
L = l->x * x + l->y * y + l->z * z;
l->m += L;
Lx = l->xx * x + l->xy * y + l->xz * z;
Ly = l->xy * x + l->yy * y + l->yz * z;
Lz = l->xz * x + l->yz * y - (l->xx + l->yy) * z;
L = Lx * hx + Ly * hy + Lz * hz;
l->x += Lx;
l->y += Ly;
l->z += Lz;
l->m += L;
Lxx = l->xxx * x + l->xxy * y + l->xxz * z;
Lxy = l->xxy * x + l->xyy * y + l->xyz * z;
Lxz = l->xxz * x + l->xyz * y - (l->xxx + l->xyy) * z;
Lyy = l->xyy * x + l->yyy * y + l->yyz * z;
Lyz = l->xyz * x + l->yyz * y - (l->xxy + l->yyy) * z;
Lx = Lxx * hx + Lxy * hy + Lxz * hz;
Ly = Lxy * hx + Lyy * hy + Lyz * hz;
Lz = Lxz * hx + Lyz * hy - (Lxx + Lyy) * hz;
L = Lx * tx + Ly * ty + Lz * tz;
l->xx += Lxx;
l->xy += Lxy;
l->xz += Lxz;
l->yy += Lyy;
l->yz += Lyz;
l->x += Lx;
l->y += Ly;
l->z += Lz;
l->m += L;
Lxxx = l->xxxx * x + l->xxxy * y + l->xxxz * z;
Lxxy = l->xxxy * x + l->xxyy * y + l->xxyz * z;
Lxxz = l->xxxz * x + l->xxyz * y - (l->xxxx + l->xxyy) * z;
Lxyy = l->xxyy * x + l->xyyy * y + l->xyyz * z;
Lxyz = l->xxyz * x + l->xyyz * y - (l->xxxy + l->xyyy) * z;
Lyyy = l->xyyy * x + l->yyyy * y + l->yyyz * z;
Lyyz = l->xyyz * x + l->yyyz * y - (l->xxyy + l->yyyy) * z;
Lxx = Lxxx * hx + Lxxy * hy + Lxxz * hz;
Lxy = Lxxy * hx + Lxyy * hy + Lxyz * hz;
Lxz = Lxxz * hx + Lxyz * hy - (Lxxx + Lxyy) * hz;
Lyy = Lxyy * hx + Lyyy * hy + Lyyz * hz;
Lyz = Lxyz * hx + Lyyz * hy - (Lxxy + Lyyy) * hz;
Lx = Lxx * tx + Lxy * ty + Lxz * tz;
Ly = Lxy * tx + Lyy * ty + Lyz * tz;
Lz = Lxz * tx + Lyz * ty - (Lxx + Lyy) * tz;
L = Lx * hx + Ly * hy + Lz * hz;
l->xxx += Lxxx;
l->xxy += Lxxy;
l->xxz += Lxxz;
l->xyy += Lxyy;
l->xyz += Lxyz;
l->yyy += Lyyy;
l->yyz += Lyyz;
l->xx += Lxx;
l->xy += Lxy;
l->xz += Lxz;
l->yy += Lyy;
l->yz += Lyz;
l->x += Lx;
l->y += Ly;
l->z += Lz;
l->m += 0.5 * L;
Lxxxx = l->xxxxx * x + l->xxxxy * y + l->xxxxz * z;
Lxxxy = l->xxxxy * x + l->xxxyy * y + l->xxxyz * z;
Lxxyy = l->xxxyy * x + l->xxyyy * y + l->xxyyz * z;
Lxyyy = l->xxyyy * x + l->xyyyy * y + l->xyyyz * z;
Lyyyy = l->xyyyy * x + l->yyyyy * y + l->yyyyz * z;
Lxxxz = l->xxxxz * x + l->xxxyz * y - (l->xxxxx + l->xxxyy) * z;
Lxxyz = l->xxxyz * x + l->xxyyz * y - (l->xxxxy + l->xxyyy) * z;
Lxyyz = l->xxyyz * x + l->xyyyz * y - (l->xxxyy + l->xyyyy) * z;
Lyyyz = l->xyyyz * x + l->yyyyz * y - (l->xxyyy + l->yyyyy) * z;
Lxxx = Lxxxx * hx + Lxxxy * hy + Lxxxz * hz;
Lxxy = Lxxxy * hx + Lxxyy * hy + Lxxyz * hz;
Lxyy = Lxxyy * hx + Lxyyy * hy + Lxyyz * hz;
Lyyy = Lxyyy * hx + Lyyyy * hy + Lyyyz * hz;
Lxxz = Lxxxz * hx + Lxxyz * hy - (Lxxxx + Lxxyy) * hz;
Lxyz = Lxxyz * hx + Lxyyz * hy - (Lxxxy + Lxyyy) * hz;
Lyyz = Lxyyz * hx + Lyyyz * hy - (Lxxyy + Lyyyy) * hz;
Lxx = Lxxx * tx + Lxxy * ty + Lxxz * tz;
Lxy = Lxxy * tx + Lxyy * ty + Lxyz * tz;
Lyy = Lxyy * tx + Lyyy * ty + Lyyz * tz;
Lxz = Lxxz * tx + Lxyz * ty - (Lxxx + Lxyy) * tz;
Lyz = Lxyz * tx + Lyyz * ty - (Lxxy + Lyyy) * tz;
Lx = Lxx * hx + Lxy * hy + Lxz * hz;
Ly = Lxy * hx + Lyy * hy + Lyz * hz;
Lz = Lxz * hx + Lyz * hy - (Lxx + Lyy) * hz;
L = Lx * hx + Ly * hy + Lz * hz;
l->xxxx += Lxxxx;
l->xxxy += Lxxxy;
l->xxyy += Lxxyy;
l->xyyy += Lxyyy;
l->yyyy += Lyyyy;
l->xxxz += Lxxxz;
l->xxyz += Lxxyz;
l->xyyz += Lxyyz;
l->yyyz += Lyyyz;
l->xxx += Lxxx;
l->xxy += Lxxy;
l->xxz += Lxxz;
l->xyy += Lxyy;
l->xyz += Lxyz;
l->yyy += Lyyy;
l->yyz += Lyyz;
l->xx += Lxx;
l->xy += Lxy;
l->xz += Lxz;
l->yy += Lyy;
l->yz += Lyz;
l->x += 0.5 * Lx;
l->y += 0.5 * Ly;
l->z += 0.5 * Lz;
l->m += 0.2 * L;
return 332.0;
}
/*
** This function shifts a reduced scaled local expansion (FLOCR) to a new
*center of expansion.
** <x,y,z> := d := rexp(new) - rexp(old).
**
** Op Count (/,*,+) = (1,162,173) = 336
*/
double momShiftFlocr(FLOCR *l, cosmoType v, cosmoType x, cosmoType y,
cosmoType z) {
const cosmoType onethird = 1.0f / 3.0f;
cosmoType iv, hx, hy, hz, tx, ty, tz;
cosmoType L, Lx, Ly, Lz, Lxx, Lxy, Lxz, Lyy, Lyz, Lxxx, Lxxy, Lxxz, Lxyy,
Lxyz, Lyyy, Lyyz;
cosmoType Lxxxx, Lxxxy, Lxxyy, Lxyyy, Lyyyy, Lxxxz, Lxxyz, Lxyyz, Lyyyz;
iv = 1.0f / v;
x *= iv;
y *= iv;
z *= iv;
hx = 0.5 * x;
hy = 0.5 * y;
hz = 0.5 * z;
tx = onethird * x;
ty = onethird * y;
tz = onethird * z;
L = l->x * x + l->y * y + l->z * z;
l->m += L;
Lx = l->xx * x + l->xy * y + l->xz * z;
Ly = l->xy * x + l->yy * y + l->yz * z;
Lz = l->xz * x + l->yz * y - (l->xx + l->yy) * z;
L = Lx * hx + Ly * hy + Lz * hz;
l->x += Lx;
l->y += Ly;
l->z += Lz;
l->m += L;
Lxx = l->xxx * x + l->xxy * y + l->xxz * z;
Lxy = l->xxy * x + l->xyy * y + l->xyz * z;
Lxz = l->xxz * x + l->xyz * y - (l->xxx + l->xyy) * z;
Lyy = l->xyy * x + l->yyy * y + l->yyz * z;
Lyz = l->xyz * x + l->yyz * y - (l->xxy + l->yyy) * z;
Lx = Lxx * hx + Lxy * hy + Lxz * hz;
Ly = Lxy * hx + Lyy * hy + Lyz * hz;
Lz = Lxz * hx + Lyz * hy - (Lxx + Lyy) * hz;
L = Lx * tx + Ly * ty + Lz * tz;
l->xx += Lxx;
l->xy += Lxy;
l->xz += Lxz;
l->yy += Lyy;
l->yz += Lyz;
l->x += Lx;
l->y += Ly;
l->z += Lz;
l->m += L;
Lxxx = l->xxxx * x + l->xxxy * y + l->xxxz * z;
Lxxy = l->xxxy * x + l->xxyy * y + l->xxyz * z;
Lxxz = l->xxxz * x + l->xxyz * y - (l->xxxx + l->xxyy) * z;
Lxyy = l->xxyy * x + l->xyyy * y + l->xyyz * z;
Lxyz = l->xxyz * x + l->xyyz * y - (l->xxxy + l->xyyy) * z;
Lyyy = l->xyyy * x + l->yyyy * y + l->yyyz * z;
Lyyz = l->xyyz * x + l->yyyz * y - (l->xxyy + l->yyyy) * z;
Lxx = Lxxx * hx + Lxxy * hy + Lxxz * hz;
Lxy = Lxxy * hx + Lxyy * hy + Lxyz * hz;
Lxz = Lxxz * hx + Lxyz * hy - (Lxxx + Lxyy) * hz;
Lyy = Lxyy * hx + Lyyy * hy + Lyyz * hz;
Lyz = Lxyz * hx + Lyyz * hy - (Lxxy + Lyyy) * hz;
Lx = Lxx * tx + Lxy * ty + Lxz * tz;
Ly = Lxy * tx + Lyy * ty + Lyz * tz;
Lz = Lxz * tx + Lyz * ty - (Lxx + Lyy) * tz;
L = Lx * hx + Ly * hy + Lz * hz;
l->xxx += Lxxx;
l->xxy += Lxxy;
l->xxz += Lxxz;
l->xyy += Lxyy;
l->xyz += Lxyz;
l->yyy += Lyyy;
l->yyz += Lyyz;
l->xx += Lxx;
l->xy += Lxy;
l->xz += Lxz;
l->yy += Lyy;
l->yz += Lyz;
l->x += Lx;
l->y += Ly;
l->z += Lz;
l->m += 0.5 * L;
Lxxxx = l->xxxxx * x + l->xxxxy * y + l->xxxxz * z;
Lxxxy = l->xxxxy * x + l->xxxyy * y + l->xxxyz * z;
Lxxyy = l->xxxyy * x + l->xxyyy * y + l->xxyyz * z;
Lxyyy = l->xxyyy * x + l->xyyyy * y + l->xyyyz * z;
Lyyyy = l->xyyyy * x + l->yyyyy * y + l->yyyyz * z;
Lxxxz = l->xxxxz * x + l->xxxyz * y - (l->xxxxx + l->xxxyy) * z;
Lxxyz = l->xxxyz * x + l->xxyyz * y - (l->xxxxy + l->xxyyy) * z;
Lxyyz = l->xxyyz * x + l->xyyyz * y - (l->xxxyy + l->xyyyy) * z;
Lyyyz = l->xyyyz * x + l->yyyyz * y - (l->xxyyy + l->yyyyy) * z;
Lxxx = Lxxxx * hx + Lxxxy * hy + Lxxxz * hz;
Lxxy = Lxxxy * hx + Lxxyy * hy + Lxxyz * hz;
Lxyy = Lxxyy * hx + Lxyyy * hy + Lxyyz * hz;
Lyyy = Lxyyy * hx + Lyyyy * hy + Lyyyz * hz;
Lxxz = Lxxxz * hx + Lxxyz * hy - (Lxxxx + Lxxyy) * hz;
Lxyz = Lxxyz * hx + Lxyyz * hy - (Lxxxy + Lxyyy) * hz;
Lyyz = Lxyyz * hx + Lyyyz * hy - (Lxxyy + Lyyyy) * hz;
Lxx = Lxxx * tx + Lxxy * ty + Lxxz * tz;
Lxy = Lxxy * tx + Lxyy * ty + Lxyz * tz;
Lyy = Lxyy * tx + Lyyy * ty + Lyyz * tz;
Lxz = Lxxz * tx + Lxyz * ty - (Lxxx + Lxyy) * tz;
Lyz = Lxyz * tx + Lyyz * ty - (Lxxy + Lyyy) * tz;
Lx = Lxx * hx + Lxy * hy + Lxz * hz;
Ly = Lxy * hx + Lyy * hy + Lyz * hz;
Lz = Lxz * hx + Lyz * hy - (Lxx + Lyy) * hz;
L = Lx * hx + Ly * hy + Lz * hz;
l->xxxx += Lxxxx;
l->xxxy += Lxxxy;
l->xxyy += Lxxyy;
l->xyyy += Lxyyy;
l->yyyy += Lyyyy;
l->xxxz += Lxxxz;
l->xxyz += Lxxyz;
l->xyyz += Lxyyz;
l->yyyz += Lyyyz;
l->xxx += Lxxx;
l->xxy += Lxxy;
l->xxz += Lxxz;
l->xyy += Lxyy;
l->xyz += Lxyz;
l->yyy += Lyyy;
l->yyz += Lyyz;
l->xx += Lxx;
l->xy += Lxy;
l->xz += Lxz;
l->yy += Lyy;
l->yz += Lyz;
l->x += 0.5 * Lx;
l->y += 0.5 * Ly;
l->z += 0.5 * Lz;
l->m += 0.2 * L;
return 332.0;
}
/*
** This function converts a complete multipole (MOMC) to a reduced one (MOMR).
*/
void momReduceMomc(MOMC *mc, MOMR *mr) {
double t, tx, ty, tz, txx, txy, txz, tyy, tyz, tzz;
/*
** First reduce Hexadecapole.
*/
txx = (mc->xxxx + mc->xxyy + mc->xxzz) / 7;
txy = (mc->xxxy + mc->xyyy + mc->xyzz) / 7;
txz = (mc->xxxz + mc->xyyz + mc->xzzz) / 7;
tyy = (mc->xxyy + mc->yyyy + mc->yyzz) / 7;
tyz = (mc->xxyz + mc->yyyz + mc->yzzz) / 7;
tzz = (mc->xxzz + mc->yyzz + mc->zzzz) / 7;
t = 0.1 * (txx + tyy + tzz);
mr->xxxx = mc->xxxx - 6 * (txx - t);
mr->xyyy = mc->xyyy - 3 * txy;
mr->xxxy = mc->xxxy - 3 * txy;
mr->yyyy = mc->yyyy - 6 * (tyy - t);
mr->xxxz = mc->xxxz - 3 * txz;
mr->yyyz = mc->yyyz - 3 * tyz;
mr->xxyy = mc->xxyy - (txx + tyy - 2 * t);
mr->xxyz = mc->xxyz - tyz;
mr->xyyz = mc->xyyz - txz;
/*
** Now reduce the Octopole.
*/
tx = (mc->xxx + mc->xyy + mc->xzz) / 5;
ty = (mc->xxy + mc->yyy + mc->yzz) / 5;
tz = (mc->xxz + mc->yyz + mc->zzz) / 5;
mr->xxx = mc->xxx - 3 * tx;
mr->xyy = mc->xyy - tx;
mr->xxy = mc->xxy - ty;
mr->yyy = mc->yyy - 3 * ty;
mr->xxz = mc->xxz - tz;
mr->yyz = mc->yyz - tz;
mr->xyz = mc->xyz;
/*
** Now reduce the Quadrupole.
*/
t = (mc->xx + mc->yy + mc->zz) / 3;
mr->xx = mc->xx - t;
mr->yy = mc->yy - t;
mr->xy = mc->xy;
mr->xz = mc->xz;
mr->yz = mc->yz;
/*
** Finally the mass remains the same.
*/
mr->m = mc->m;
}
/*
** This is a new fast version of QEVAL which evaluates
** the interaction due to the reduced moment 'm'.
** This version is nearly two times as fast as a naive
** implementation.
**
** OpCount = (*,+) = (103,72) = 175 - 8 = 167
*/
void momEvalMomr(MOMR *m, double dir0, double x, double y, double z,
double *fPot, double *ax, double *ay, double *az) {
const double onethird = 1.0 / 3.0;
double xx, xy, xz, yy, yz, zz;
double xxx, xxy, xxz, xyy, yyy, yyz, xyz;
double tx, ty, tz, dir2, g2, g3, g4;
double dir = -dir0;
dir2 = dir * dir;
g2 = 3 * dir * dir2 * dir2;
g3 = -5 * g2 * dir2;
g4 = -7 * g3 * dir2;
/*
** Calculate the funky distance terms.
*/
xx = 0.5 * x * x;
xy = x * y;
xz = x * z;
yy = 0.5 * y * y;
yz = y * z;
zz = 0.5 * z * z;
xxx = x * (onethird * xx - zz);
xxz = z * (xx - onethird * zz);
yyy = y * (onethird * yy - zz);
yyz = z * (yy - onethird * zz);
xx -= zz;
yy -= zz;
xxy = y * xx;
xyy = x * yy;
xyz = xy * z;
/*
** Now calculate the interaction up to Hexadecapole order.
*/
tx = g4 * (m->xxxx * xxx + m->xyyy * yyy + m->xxxy * xxy + m->xxxz * xxz +
m->xxyy * xyy + m->xxyz * xyz + m->xyyz * yyz);
ty = g4 * (m->xyyy * xyy + m->xxxy * xxx + m->yyyy * yyy + m->yyyz * yyz +
m->xxyy * xxy + m->xxyz * xxz + m->xyyz * xyz);
tz = g4 * (-m->xxxx * xxz - (m->xyyy + m->xxxy) * xyz - m->yyyy * yyz +
m->xxxz * xxx + m->yyyz * yyy - m->xxyy * (xxz + yyz) +
m->xxyz * xxy + m->xyyz * xyy);
g4 = 0.25 * (tx * x + ty * y + tz * z);
xxx = g3 *
(m->xxx * xx + m->xyy * yy + m->xxy * xy + m->xxz * xz + m->xyz * yz);
xxy = g3 *
(m->xyy * xy + m->xxy * xx + m->yyy * yy + m->yyz * yz + m->xyz * xz);
xxz = g3 * (-(m->xxx + m->xyy) * xz - (m->xxy + m->yyy) * yz + m->xxz * xx +
m->yyz * yy + m->xyz * xy);
g3 = onethird * (xxx * x + xxy * y + xxz * z);
xx = g2 * (m->xx * x + m->xy * y + m->xz * z);
xy = g2 * (m->yy * y + m->xy * x + m->yz * z);
xz = g2 * (-(m->xx + m->yy) * z + m->xz * x + m->yz * y);
g2 = 0.5 * (xx * x + xy * y + xz * z);
dir *= m->m;
dir2 *= -(dir + 5 * g2 + 7 * g3 + 9 * g4);
*fPot += dir + g2 + g3 + g4;
*ax += xx + xxx + tx + x * dir2;
*ay += xy + xxy + ty + y * dir2;
*az += xz + xxz + tz + z * dir2;
}
/*
** This is a new fast version of QEVAL which evaluates
** the interaction due to the reduced moment 'm'.
** This version is nearly two times as fast as a naive
** implementation.
**
** March 23, 2007: This function now uses unit vectors
** which reduces the required precision in the exponent
** since the highest power of r is now 5 (g4 ~ r^(-5)).
**
** OpCount = (*,+) = (106,72) = 178 - 8 = 170
**
*/
void momEvalFmomrcm(FMOMR *m, cosmoType u, cosmoType dir, cosmoType x,
cosmoType y, cosmoType z, cosmoType *fPot, cosmoType *ax,
cosmoType *ay, cosmoType *az, cosmoType *magai) {
const cosmoType onethird = 1.0f / 3.0f;
cosmoType xx, xy, xz, yy, yz, zz;
cosmoType xxx, xxy, xxz, xyy, yyy, yyz, xyz;
cosmoType tx, ty, tz, g0, g2, g3, g4;
u *= dir;
g0 = dir;
g2 = 3 * dir * u * u;
g3 = 5 * g2 * u;
g4 = 7 * g3 * u;
/*
** Calculate the trace-free distance terms.
*/
x *= dir;
y *= dir;
z *= dir;
xx = 0.5 * x * x;
xy = x * y;
xz = x * z;
yy = 0.5 * y * y;
yz = y * z;
zz = 0.5 * z * z;
xxx = x * (onethird * xx - zz);
xxz = z * (xx - onethird * zz);
yyy = y * (onethird * yy - zz);
yyz = z * (yy - onethird * zz);
xx -= zz;
yy -= zz;
xxy = y * xx;
xyy = x * yy;
xyz = xy * z;
/*
** Now calculate the interaction up to Hexadecapole order.
*/
tx = g4 * (m->xxxx * xxx + m->xyyy * yyy + m->xxxy * xxy + m->xxxz * xxz +
m->xxyy * xyy + m->xxyz * xyz + m->xyyz * yyz);
ty = g4 * (m->xyyy * xyy + m->xxxy * xxx + m->yyyy * yyy + m->yyyz * yyz +
m->xxyy * xxy + m->xxyz * xxz + m->xyyz * xyz);
tz = g4 * (-m->xxxx * xxz - (m->xyyy + m->xxxy) * xyz - m->yyyy * yyz +
m->xxxz * xxx + m->yyyz * yyy - m->xxyy * (xxz + yyz) +
m->xxyz * xxy + m->xyyz * xyy);
g4 = 0.25 * (tx * x + ty * y + tz * z);
xxx = g3 *
(m->xxx * xx + m->xyy * yy + m->xxy * xy + m->xxz * xz + m->xyz * yz);
xxy = g3 *
(m->xyy * xy + m->xxy * xx + m->yyy * yy + m->yyz * yz + m->xyz * xz);
xxz = g3 * (-(m->xxx + m->xyy) * xz - (m->xxy + m->yyy) * yz + m->xxz * xx +
m->yyz * yy + m->xyz * xy);
g3 = onethird * (xxx * x + xxy * y + xxz * z);
xx = g2 * (m->xx * x + m->xy * y + m->xz * z);
xy = g2 * (m->yy * y + m->xy * x + m->yz * z);
xz = g2 * (-(m->xx + m->yy) * z + m->xz * x + m->yz * y);
g2 = 0.5 * (xx * x + xy * y + xz * z);
g0 *= m->m;
*fPot += -(g0 + g2 + g3 + g4);
g0 += 5 * g2 + 7 * g3 + 9 * g4;
*ax += dir * (xx + xxx + tx - x * g0);
*ay += dir * (xy + xxy + ty - y * g0);
*az += dir * (xz + xxz + tz - z * g0);
*magai = g0 * dir;
}
void momMomr2Momc(MOMR *ma, MOMC *mc) {
mc->m = ma->m;
mc->xx = ma->xx;
mc->yy = ma->yy;
mc->xy = ma->xy;
mc->xz = ma->xz;
mc->yz = ma->yz;
mc->xxx = ma->xxx;
mc->xyy = ma->xyy;
mc->xxy = ma->xxy;
mc->yyy = ma->yyy;
mc->xxz = ma->xxz;
mc->yyz = ma->yyz;
mc->xyz = ma->xyz;
mc->xxxx = ma->xxxx;
mc->xyyy = ma->xyyy;
mc->xxxy = ma->xxxy;
mc->yyyy = ma->yyyy;
mc->xxxz = ma->xxxz;
mc->yyyz = ma->yyyz;
mc->xxyy = ma->xxyy;
mc->xxyz = ma->xxyz;
mc->xyyz = ma->xyyz;
mc->zz = -(ma->xx + ma->yy);
mc->xzz = -(ma->xxx + ma->xyy);
mc->yzz = -(ma->xxy + ma->yyy);
mc->zzz = -(ma->xxz + ma->yyz);
mc->xxzz = -(ma->xxxx + ma->xxyy);
mc->xyzz = -(ma->xxxy + ma->xyyy);
mc->xzzz = -(ma->xxxz + ma->xyyz);
mc->yyzz = -(ma->xxyy + ma->yyyy);
mc->yzzz = -(ma->xxyz + ma->yyyz);
mc->zzzz = -(mc->xxzz + mc->yyzz);
}
void momFmomr2Momc(FMOMR *ma, MOMC *mc) {
mc->m = ma->m;
mc->xx = ma->xx;
mc->yy = ma->yy;
mc->xy = ma->xy;
mc->xz = ma->xz;
mc->yz = ma->yz;
mc->xxx = ma->xxx;
mc->xyy = ma->xyy;
mc->xxy = ma->xxy;
mc->yyy = ma->yyy;
mc->xxz = ma->xxz;
mc->yyz = ma->yyz;
mc->xyz = ma->xyz;
mc->xxxx = ma->xxxx;
mc->xyyy = ma->xyyy;
mc->xxxy = ma->xxxy;
mc->yyyy = ma->yyyy;
mc->xxxz = ma->xxxz;
mc->yyyz = ma->yyyz;
mc->xxyy = ma->xxyy;
mc->xxyz = ma->xxyz;
mc->xyyz = ma->xyyz;
mc->zz = -(ma->xx + ma->yy);
mc->xzz = -(ma->xxx + ma->xyy);
mc->yzz = -(ma->xxy + ma->yyy);
mc->zzz = -(ma->xxz + ma->yyz);
mc->xxzz = -(ma->xxxx + ma->xxyy);
mc->xyzz = -(ma->xxxy + ma->xyyy);
mc->xzzz = -(ma->xxxz + ma->xyyz);
mc->yyzz = -(ma->xxyy + ma->yyyy);
mc->yzzz = -(ma->xxyz + ma->yyyz);
mc->zzzz = -(mc->xxzz + mc->yyzz);
}
/*
** Op Count = (*,+) = (302,207) = 509
*/
double momLocrAddMomr5(LOCR *l, MOMR *m, double dir, double x, double y,
double z, double *tax, double *tay, double *taz) {
const double onethird = 1.0 / 3.0;
double xx, xy, xz, yy, yz, zz;
double xxx, xxy, xyy, yyy, xxz, xyz, yyz;
double Ax, Ay, Az, A, Bxx, Bxy, Byy, Bxz, Byz, Bx, By, Bz, B, Cx, Cy, Cz, C;
double R1, R2, R3, T2, T3;
double g0, g1, g2, g3, g4, g5;
double g4xx, g4yy, g5xx, g5yy, fxx, fyy;
x *= dir;
y *= dir;
z *= dir;
g0 = -dir;
g1 = -g0 * dir;
g2 = -3 * g1 * dir;
g3 = -5 * g2 * dir;
g4 = -7 * g3 * dir;
g5 = -9 * g4 * dir;
/*
** Calculate the funky distance terms.
*/
xx = 0.5 * x * x;
xy = x * y;
yy = 0.5 * y * y;
xz = x * z;
yz = y * z;
zz = 0.5 * z * z;
xxx = x * (onethird * xx - zz);
xxz = z * (xx - onethird * zz);
yyy = y * (onethird * yy - zz);
yyz = z * (yy - onethird * zz);
xx -= zz;
yy -= zz;
xxy = y * xx;
xyy = x * yy;
xyz = xy * z;
Bxx = x * m->xxx + y * m->xxy + z * m->xxz;
Bxy = x * m->xxy + y * m->xyy + z * m->xyz;
Byy = x * m->xyy + y * m->yyy + z * m->yyz;
Bxz = x * m->xxz + y * m->xyz - z * (m->xxx + m->xyy);
Byz = x * m->xyz + y * m->yyz - z * (m->xxy + m->yyy);
Cx = m->xxxx * xxx + m->xyyy * yyy + m->xxxy * xxy + m->xxxz * xxz +
m->xxyy * xyy + m->xxyz * xyz + m->xyyz * yyz;
Cy = m->xyyy * xyy + m->xxxy * xxx + m->yyyy * yyy + m->yyyz * yyz +
m->xxyy * xxy + m->xxyz * xxz + m->xyyz * xyz;
Cz = -m->xxxx * xxz - (m->xyyy + m->xxxy) * xyz - m->yyyy * yyz +
m->xxxz * xxx + m->yyyz * yyy - m->xxyy * (xxz + yyz) + m->xxyz * xxy +
m->xyyz * xyy;
Ax = x * m->xx + y * m->xy + z * m->xz;
Ay = x * m->xy + y * m->yy + z * m->yz;
Az = x * m->xz + y * m->yz - z * (m->xx + m->yy);
Bx = 0.5 * (x * Bxx + y * Bxy + z * Bxz);
By = 0.5 * (x * Bxy + y * Byy + z * Byz);
Bz = 0.5 * (x * Bxz + y * Byz - z * (Bxx + Byy));
C = 0.25 * (x * Cx + y * Cy + z * Cz);
A = 0.5 * (x * Ax + y * Ay + z * Az);
B = onethird * (x * Bx + y * By + z * Bz);
xx = x * x;
yy = y * y;
l->m += g0 * m->m + g2 * A - g3 * B + g4 * C;
R1 = g1 * m->m + g3 * A - g4 * B + g5 * C;
R2 = g2 * m->m + g4 * A - g5 * B;
R3 = g3 * m->m + g5 * A;
g1 *= dir;
g2 *= dir;
g3 *= dir;
T2 = g1 * m->m + g3 * A;
T3 = g2 * m->m;
*tax = -(g2 * Ax - g3 * Bx + x * R1);
*tay = -(g2 * Ay - g3 * By + y * R1);
*taz = -(g2 * Az - g3 * Bz + z * R1);
g2 *= dir;
g4xx = g4 * xx;
g4yy = g4 * yy;
l->xxxx += (3 * g2 + (6 * g3 + g4xx) * xx) * m->m;
l->yyyy += (3 * g2 + (6 * g3 + g4yy) * yy) * m->m;
fxx = (3 * g3 + g4xx) * m->m;
l->xxxy += fxx * xy;
l->xxxz += fxx * xz;
fyy = (3 * g3 + g4yy) * m->m;
l->xyyy += fyy * xy;
l->yyyz += fyy * yz;
fxx = (g3 + g4xx);
l->xxyz += fxx * yz * m->m;
l->xxyy += (g2 + g3 * xx + fxx * yy) * m->m;
l->xyyz += (g3 + g4yy) * xz * m->m;
g4 *= dir;
T2 -= g4 * B;
T3 += g4 * A;
*tax -= g4 * Cx;
*tay -= g4 * Cy;
*taz -= g4 * Cz;
l->x -= *tax;
l->y -= *tay;
l->z -= *taz;
l->xx += g2 * m->xx + T2 + R2 * xx + 2 * g3 * Ax * x - 2 * g4 * Bx * x;
l->yy += g2 * m->yy + T2 + R2 * yy + 2 * g3 * Ay * y - 2 * g4 * By * y;
l->xy +=
g2 * m->xy + R2 * xy + g3 * (Ax * y + Ay * x) - g4 * (Bx * y + By * x);
l->xz +=
g2 * m->xz + R2 * xz + g3 * (Ax * z + Az * x) - g4 * (Bx * z + Bz * x);
l->yz +=
g2 * m->yz + R2 * yz + g3 * (Ay * z + Az * y) - g4 * (By * z + Bz * y);
g3 *= dir;
g5xx = g5 * xx;
g5yy = g5 * yy;
l->xx -= g3 * Bxx;
l->xy -= g3 * Bxy;
l->yy -= g3 * Byy;
l->xz -= g3 * Bxz;
l->yz -= g3 * Byz;
fxx = T3 + R3 * xx;
fyy = T3 + R3 * yy;
l->xxy += fxx * y + g4 * (2 * xy * Ax + xx * Ay) +
g3 * (Ay + 2 * m->xy * x + m->xx * y);
l->xxz += fxx * z + g4 * (2 * xz * Ax + xx * Az) +
g3 * (Az + 2 * m->xz * x + m->xx * z);
l->xyy += fyy * x + g4 * (2 * xy * Ay + yy * Ax) +
g3 * (Ax + 2 * m->xy * y + m->yy * x);
l->yyz += fyy * z + g4 * (2 * yz * Ay + yy * Az) +
g3 * (Az + 2 * m->yz * y + m->yy * z);
l->xyz += R3 * xy * z + g4 * (yz * Ax + xz * Ay + xy * Az) +
g3 * (m->xy * z + m->xz * y + m->yz * x);
fxx += 2 * T3;
fyy += 2 * T3;
l->xxx += fxx * x + 3 * (g4 * xx * Ax + g3 * (Ax + m->xx * x));
l->yyy += fyy * y + 3 * (g4 * yy * Ay + g3 * (Ay + m->yy * y));
fxx = 3 * g3 + (6 * g4 + g5xx) * xx;
fyy = 3 * g3 + (6 * g4 + g5yy) * yy;
x *= m->m;
y *= m->m;
z *= m->m;
l->xxxxx += (15 * g3 + (10 * g4 + g5xx) * xx) * x;
l->yyyyy += (15 * g3 + (10 * g4 + g5yy) * yy) * y;
l->xxxyy += (3 * g3 + 3 * g4 * yy + (g4 + g5yy) * xx) * x;
l->xxyyy += (3 * g3 + 3 * g4 * xx + (g4 + g5xx) * yy) * y;
l->xxyyz += (g3 + g4 * xx + (g4 + g5xx) * yy) * z;
l->xxxxy += fxx * y;
l->xxxxz += fxx * z;
l->xyyyy += fyy * x;
l->yyyyz += fyy * z;
l->xxxyz += (3 * g4 + g5xx) * xy * z;
l->xyyyz += (3 * g4 + g5yy) * xy * z;
return 509.0; /* a guess right now */
}
/*
** Op Count = (*,+,-) = (265,150,49) = 464
*/
double momFlocrAddFmomr5cm(FLOCR *l, cosmoType v, FMOMR *m, cosmoType u,
cosmoType dir, cosmoType x, cosmoType y, cosmoType z,
cosmoType *tax, cosmoType *tay, cosmoType *taz) {
const cosmoType onethird = 1.0f / 3.0f;
cosmoType u2, u3, u4;
cosmoType xx, xy, xz, yy, yz, zz;
cosmoType xxx, xxy, xyy, yyy, xxz, xyz, yyz;
cosmoType R2xx, R2xy, R2xz, R2yy, R2yz, R2x, R2y, R2z, R2, R3xx, R3xy, R3yy,
R3xz, R3yz, R3x, R3y, R3z, R3, R4x, R4y, R4z, R4;
cosmoType T0, txx, tyy, t1, t1x, t1y, t1z, t1xx, t1yy, t2x, t2y, t2z, t2xx,
t2yy, txxxx, tyyyy;
u *= dir;
x *= dir;
y *= dir;
z *= dir;
dir = -dir;
v *= dir;
u2 = 15.0f * u * u; /* becomes 15.0f*u2! */
/*
** Calculate the funky distance terms.
*/
xx = 0.5f * x * x;
xy = x * y;
yy = 0.5f * y * y;
xz = x * z;
yz = y * z;
zz = 0.5f * z * z;
xxx = x * (onethird * xx - zz);
xxz = z * (xx - onethird * zz);
yyy = y * (onethird * yy - zz);
yyz = z * (yy - onethird * zz);
xx -= zz;
yy -= zz;
xxy = y * xx;
xyy = x * yy;
xyz = xy * z;
u3 = u2 * u; /* becomes 5.0f*u3! */
R2xx = u2 * m->xx;
R2xy = u2 * m->xy;
R2xz = u2 * m->xz;
R2yy = u2 * m->yy;
R2yz = u2 * m->yz;
u4 = 7.0f * u3 * u; /* becomes 7.0f*5.0f*u4! */
R2x = x * R2xx + y * R2xy + z * R2xz;
R2y = x * R2xy + y * R2yy + z * R2yz;
R2z = x * R2xz + y * R2yz - z * (R2xx + R2yy);
R3xx = u3 * (x * m->xxx + y * m->xxy + z * m->xxz);
R3xy = u3 * (x * m->xxy + y * m->xyy + z * m->xyz);
R3yy = u3 * (x * m->xyy + y * m->yyy + z * m->yyz);
R3xz = u3 * (x * m->xxz + y * m->xyz - z * (m->xxx + m->xyy));
R3yz = u3 * (x * m->xyz + y * m->yyz - z * (m->xxy + m->yyy));
R4x = u4 * (m->xxxx * xxx + m->xyyy * yyy + m->xxxy * xxy + m->xxxz * xxz +
m->xxyy * xyy + m->xxyz * xyz + m->xyyz * yyz);
R4y = u4 * (m->xyyy * xyy + m->xxxy * xxx + m->yyyy * yyy + m->yyyz * yyz +
m->xxyy * xxy + m->xxyz * xxz + m->xyyz * xyz);
R4z = u4 * (-m->xxxx * xxz - (m->xyyy + m->xxxy) * xyz - m->yyyy * yyz +
m->xxxz * xxx + m->yyyz * yyy - m->xxyy * (xxz + yyz) +
m->xxyz * xxy + m->xyyz * xyy);
R3x = 0.5f * (x * R3xx + y * R3xy + z * R3xz);
R3y = 0.5f * (x * R3xy + y * R3yy + z * R3yz);
R3z = 0.5f * (x * R3xz + y * R3yz - z * (R3xx + R3yy));
R4 = 0.25f * (x * R4x + y * R4y + z * R4z);
R2 = 0.5f * (x * R2x + y * R2y + z * R2z);
R3 = onethird * (x * R3x + y * R3y + z * R3z);
xx = x * x;
yy = y * y;
/*
** Now we use the 'R's.
*/
l->m += dir * (m->m + 0.2f * R2 + R3 + R4);
dir *= v;
T0 = -(m->m + R2 + 7.0f * R3 + 9.0f * R4);
*tax = dir * (T0 * x + 0.2f * R2x + R3x + R4x);
*tay = dir * (T0 * y + 0.2f * R2y + R3y + R4y);
*taz = dir * (T0 * z + 0.2f * R2z + R3z + R4z);
l->x -= *tax;
l->y -= *tay;
l->z -= *taz;
dir *= v;
T0 = 3.0f * m->m + 7.0f * (R2 + 9.0f * R3);
t1 = m->m + R2 + 7.0f * R3;
t1x = R2x + 7.0f * R3x;
t1y = R2y + 7.0f * R3y;
t1z = R2z + 7.0f * R3z;
l->xx += dir * (T0 * xx - t1 - 2.0f * x * t1x + 0.2f * R2xx + R3xx);
l->yy += dir * (T0 * yy - t1 - 2.0f * y * t1y + 0.2f * R2yy + R3yy);
l->xy += dir * (T0 * xy - y * t1x - x * t1y + 0.2f * R2xy + R3xy);
l->xz += dir * (T0 * xz - z * t1x - x * t1z + 0.2f * R2xz + R3xz);
l->yz += dir * (T0 * yz - z * t1y - y * t1z + 0.2f * R2yz + R3yz);
dir *= v;
T0 = 15.0f * m->m + 63.0f * R2;
txx = T0 * xx;
tyy = T0 * yy;
t1 = 3.0f * m->m + 7.0f * R2;
t2x = -7.0f * R2x;
t2y = -7.0f * R2y;
t2z = -7.0f * R2z;
t1xx = txx - t1 + 2.0f * x * t2x + R2xx;
t1yy = tyy - t1 + 2.0f * y * t2y + R2yy;
l->xxx += dir * (x * (txx - 3.0f * (t1 - t2x * x - R2xx)) + 3.0f * R2x);
l->yyy += dir * (y * (tyy - 3.0f * (t1 - t2y * y - R2yy)) + 3.0f * R2y);
l->xxy += dir * (y * t1xx + xx * t2y + R2y + 2.0f * R2xy * x);
l->xxz += dir * (z * t1xx + xx * t2z + R2z + 2.0f * R2xz * x);
l->xyy += dir * (x * t1yy + yy * t2x + R2x + 2.0f * R2xy * y);
l->yyz += dir * (z * t1yy + yy * t2z + R2z + 2.0f * R2yz * y);
l->xyz += dir * (T0 * xyz + (yz * t2x + xz * t2y + xy * t2z) + R2xy * z +
R2yz * x + R2xz * y);
dir *= v * m->m;
txx = 105.0f * xx;
tyy = 105.0f * yy;
t2xx = txx - 90.0f;
t2yy = tyy - 90.0f;
l->xxxx += dir * (xx * t2xx + 9.0f);
l->yyyy += dir * (yy * t2yy + 9.0f);
t2xx += 45.0f;
t2yy += 45.0f;
l->xxxy += dir * xy * t2xx;
l->xxxz += dir * xz * t2xx;
l->xyyy += dir * xy * t2yy;
l->yyyz += dir * yz * t2yy;
t2xx += 30.0f;
t2yy += 30.0f;
l->xxyy += dir * (yy * t2xx - xx * 15.0f + 3.0f);
l->xxyz += dir * (yz * t2xx);
l->xyyz += dir * (xz * t2yy);
dir *= v;
x *= dir;
y *= dir;
z *= dir;
txx = 9.0f * xx - 10.0f;
tyy = 9.0f * yy - 10.0f;
xx *= 105.0f;
yy *= 105.0f;
xy *= z * 105.0f;
l->xxxxx += x * (xx * txx + 225.0f);
l->yyyyy += y * (yy * tyy + 225.0f);
txx += 4.0f;
tyy += 4.0f;
txxxx = xx * txx + 45.0f;
tyyyy = yy * tyy + 45.0f;
l->xxxxy += y * txxxx;
l->xxxxz += z * txxxx;
l->xyyyy += x * tyyyy;
l->yyyyz += z * tyyyy;
txx += 3.0f;
tyy += 3.0f;
l->xxxyz += xy * txx;
l->xyyyz += xy * tyy;
l->xxxyy += x * (yy * txx - xx + 45.0f);
l->xxyyy += y * (xx * tyy - yy + 45.0f);
tyy += 2.0f;
l->xxyyz += z * (xx * tyy - yy + 15.0f);
return (464.0);
}
/*
** Op Count = (*,+,-) = (,,) =
*/
double momFlocrAddMono5(FLOCR *l, cosmoType v, cosmoType m, cosmoType dir,
cosmoType x, cosmoType y, cosmoType z, cosmoType *tax,
cosmoType *tay, cosmoType *taz) {
cosmoType xx, xy, xz, yy, yz;
cosmoType T0, txx, tyy, t1, t1xx, t1yy, t2, txxxx, tyyyy;
x *= dir;
y *= dir;
z *= dir;
dir = -dir;
v *= dir;
/*
** Calculate the funky distance terms.
*/
xy = x * y;
xz = x * z;
yz = y * z;
xx = x * x;
yy = y * y;
l->m += dir * m;
dir *= v;
T0 = -m * dir;
*tax = T0 * x;
*tay = T0 * y;
*taz = T0 * z;
l->x -= *tax;
l->y -= *tay;
l->z -= *taz;
dir *= v;
T0 = 3.0f * m * dir;
t1 = m * dir;
l->xx += T0 * xx - t1;
l->yy += T0 * yy - t1;
l->xy += T0 * xy;
l->xz += T0 * xz;
l->yz += T0 * yz;
dir *= v;
T0 = 15.0f * m * dir;
txx = T0 * xx;
tyy = T0 * yy;
t1 = 3.0f * m * dir;
t1xx = txx - t1;
t1yy = tyy - t1;
l->xxx += x * (txx - 3.0f * t1);
l->yyy += y * (tyy - 3.0f * t1);
l->xxy += y * t1xx;
l->xxz += z * t1xx;
l->xyy += x * t1yy;
l->yyz += z * t1yy;
l->xyz += T0 * x * yz;
dir *= v;
T0 = 105.0f * m * dir;
txx = T0 * xx;
tyy = T0 * yy;
t2 = 3.0f * m * dir;
t1 = 5.0f * t2;
l->xxxx += xx * (txx - 6.0f * t1) + 3.0f * t2;
l->yyyy += yy * (tyy - 6.0f * t1) + 3.0f * t2;
t1xx = txx - 3.0f * t1;
t1yy = tyy - 3.0f * t1;
l->xxxy += xy * t1xx;
l->xxxz += xz * t1xx;
l->xyyy += xy * t1yy;
l->yyyz += yz * t1yy;
l->xxyy += yy * txx - (xx + yy) * t1 + t2;
l->xxyz += yz * (txx - t1);
l->xyyz += xz * (tyy - t1);
dir *= v * m;
x *= dir;
y *= dir;
z *= dir;
txx = 9.0f * xx - 10.0f;
tyy = 9.0f * yy - 10.0f;
xx *= 105.0f;
yy *= 105.0f;
xy *= z * 105.0f;
l->xxxxx += x * (xx * txx + 225.0f);
l->yyyyy += y * (yy * tyy + 225.0f);
txx += 4.0f;
tyy += 4.0f;
txxxx = xx * txx + 45.0f;
tyyyy = yy * tyy + 45.0f;
l->xxxxy += y * txxxx;
l->xxxxz += z * txxxx;
l->xyyyy += x * tyyyy;
l->yyyyz += z * tyyyy;
txx += 3.0f;
tyy += 3.0f;
l->xxxyz += xy * txx;
l->xyyyz += xy * tyy;
l->xxxyy += x * (yy * txx - xx + 45.0f);
l->xxyyy += y * (xx * tyy - yy + 45.0f);
tyy += 2.0f;
l->xxyyz += z * (xx * tyy - yy + 15.0f);
return (200.0);
}
void momEvalLocr(LOCR *l, double x, double y, double z, double *fPot,
double *ax, double *ay, double *az) {
const double onethird = 1.0 / 3.0;
double xx, xy, xz, yy, yz, zz, xxx, xxz, yyy, yyz, xxy, xyy, xyz, xxxx,
xxxy, xxyy, xyyy, yyyy, xxxz, xxyz, xyyz, yyyz;
double g1, A, Ax, Ay, Az, B, Bx, By, Bz, C, Cx, Cy, Cz, D, Dx, Dy, Dz;
/*
** Calculate the funky distance terms.
*/
xx = 0.5 * x * x;
xy = x * y;
xz = x * z;
yy = 0.5 * y * y;
yz = y * z;
zz = 0.5 * z * z;
xxx = x * (onethird * xx - zz);
xxz = z * (xx - onethird * zz);
yyy = y * (onethird * yy - zz);
yyz = z * (yy - onethird * zz);
xx -= zz;
yy -= zz;
xxy = y * xx;
xyy = x * yy;
xyz = xy * z;
xxxx = 0.25 * (x * xxx - z * xxz);
xxxy = y * xxx;
xxyy = xx * yy - 2 * onethird * zz * zz;
xyyy = x * yyy;
yyyy = 0.25 * (y * yyy - z * yyz);
xxxz = onethird * x * z * xx;
xxyz = y * xxz;
xyyz = x * yyz;
yyyz = y * yyz;
/*
** Now calculate the interaction.
*/
Dx = l->xxxxx * xxxx + l->xxxxy * xxxy + l->xxxyy * xxyy + l->xxyyy * xyyy +
l->xyyyy * yyyy + l->xxxxz * xxxz + l->xxxyz * xxyz + l->xxyyz * xyyz +
l->xyyyz * yyyz;
Dy = l->xxxxy * xxxx + l->xxxyy * xxxy + l->xxyyy * xxyy + l->xyyyy * xyyy +
l->yyyyy * yyyy + l->xxxyz * xxxz + l->xxyyz * xxyz + l->xyyyz * xyyz +
l->yyyyz * yyyz;
Dz = l->xxxxz * xxxx + l->xxxyz * xxxy + l->xxyyz * xxyy + l->xyyyz * xyyy +
l->yyyyz * yyyy - l->xxxxx * xxxz - l->xxxxy * xxyz -
l->xxxyy * (xxxz + xyyz) - l->xxyyy * (xxyz + yyyz) + l->xyyyy * xyyz +
l->yyyyy * yyyz;
Cx = l->xxxx * xxx + l->xyyy * yyy + l->xxxy * xxy + l->xxxz * xxz +
l->xxyy * xyy + l->xxyz * xyz + l->xyyz * yyz;
Cy = l->xyyy * xyy + l->xxxy * xxx + l->yyyy * yyy + l->yyyz * yyz +
l->xxyy * xxy + l->xxyz * xxz + l->xyyz * xyz;
Cz = -l->xxxx * xxz - (l->xyyy + l->xxxy) * xyz - l->yyyy * yyz +
l->xxxz * xxx + l->yyyz * yyy - l->xxyy * (xxz + yyz) + l->xxyz * xxy +
l->xyyz * xyy;
Bx = l->xxx * xx + l->xyy * yy + l->xxy * xy + l->xxz * xz + l->xyz * yz;
By = l->xyy * xy + l->xxy * xx + l->yyy * yy + l->yyz * yz + l->xyz * xz;
Bz = -(l->xxx + l->xyy) * xz - (l->xxy + l->yyy) * yz + l->xxz * xx +
l->yyz * yy + l->xyz * xy;
Ax = l->xx * x + l->xy * y + l->xz * z;
Ay = l->yy * y + l->xy * x + l->yz * z;
Az = -(l->xx + l->yy) * z + l->xz * x + l->yz * y;
D = 0.2 * (Dx * x + Dy * y + Dz * z);
C = 0.25 * (Cx * x + Cy * y + Cz * z);
B = onethird * (Bx * x + By * y + Bz * z);
A = 0.5 * (Ax * x + Ay * y + Az * z);
g1 = x * l->x + y * l->y + z * l->z;
*ax -= l->x + Ax + Bx + Cx + Dx;
*ay -= l->y + Ay + By + Cy + Dy;
*az -= l->z + Az + Bz + Cz + Dz;
*fPot += l->m + g1 + A + B + C + D;
}
void momEvalFlocr(FLOCR *l, cosmoType v, cosmoType x, cosmoType y, cosmoType z,
cosmoType *fPot, cosmoType *ax, cosmoType *ay,
cosmoType *az) {
const cosmoType onethird = 1.0f / 3.0f;
cosmoType xx, xy, xz, yy, yz, zz, xxx, xxz, yyy, yyz, xxy, xyy, xyz, xxxx,
xxxy, xxyy, xyyy, yyyy, xxxz, xxyz, xyyz, yyyz;
cosmoType g1, A, Ax, Ay, Az, B, Bx, By, Bz, C, Cx, Cy, Cz, D, Dx, Dy, Dz;
cosmoType iv;
/*
** Calculate the funky distance terms, but first scale x,y,z so that they
*are of order unity!
*/
iv = 1 / v;
x *= iv;
y *= iv;
z *= iv;
xx = 0.5 * x * x;
xy = x * y;
xz = x * z;
yy = 0.5 * y * y;
yz = y * z;
zz = 0.5 * z * z;
xxx = x * (onethird * xx - zz);
xxz = z * (xx - onethird * zz);
yyy = y * (onethird * yy - zz);
yyz = z * (yy - onethird * zz);
xx -= zz;
yy -= zz;
xxy = y * xx;
xyy = x * yy;
xyz = xy * z;
xxxx = 0.25 * (x * xxx - z * xxz);
xxxy = y * xxx;
xxyy = xx * yy - 2 * onethird * zz * zz;
xyyy = x * yyy;
yyyy = 0.25 * (y * yyy - z * yyz);
xxxz = onethird * x * z * xx;
xxyz = y * xxz;
xyyz = x * yyz;
yyyz = y * yyz;
/*
** Now calculate the interaction.
*/
Dx = l->xxxxx * xxxx + l->xxxxy * xxxy + l->xxxyy * xxyy + l->xxyyy * xyyy +
l->xyyyy * yyyy + l->xxxxz * xxxz + l->xxxyz * xxyz + l->xxyyz * xyyz +
l->xyyyz * yyyz;
Dy = l->xxxxy * xxxx + l->xxxyy * xxxy + l->xxyyy * xxyy + l->xyyyy * xyyy +
l->yyyyy * yyyy + l->xxxyz * xxxz + l->xxyyz * xxyz + l->xyyyz * xyyz +
l->yyyyz * yyyz;
Dz = l->xxxxz * xxxx + l->xxxyz * xxxy + l->xxyyz * xxyy + l->xyyyz * xyyy +
l->yyyyz * yyyy - l->xxxxx * xxxz - l->xxxxy * xxyz -
l->xxxyy * (xxxz + xyyz) - l->xxyyy * (xxyz + yyyz) + l->xyyyy * xyyz +
l->yyyyy * yyyz;
Cx = l->xxxx * xxx + l->xyyy * yyy + l->xxxy * xxy + l->xxxz * xxz +
l->xxyy * xyy + l->xxyz * xyz + l->xyyz * yyz;
Cy = l->xyyy * xyy + l->xxxy * xxx + l->yyyy * yyy + l->yyyz * yyz +
l->xxyy * xxy + l->xxyz * xxz + l->xyyz * xyz;
Cz = -l->xxxx * xxz - (l->xyyy + l->xxxy) * xyz - l->yyyy * yyz +
l->xxxz * xxx + l->yyyz * yyy - l->xxyy * (xxz + yyz) + l->xxyz * xxy +
l->xyyz * xyy;
Bx = l->xxx * xx + l->xyy * yy + l->xxy * xy + l->xxz * xz + l->xyz * yz;
By = l->xyy * xy + l->xxy * xx + l->yyy * yy + l->yyz * yz + l->xyz * xz;
Bz = -(l->xxx + l->xyy) * xz - (l->xxy + l->yyy) * yz + l->xxz * xx +
l->yyz * yy + l->xyz * xy;
Ax = l->xx * x + l->xy * y + l->xz * z;
Ay = l->yy * y + l->xy * x + l->yz * z;
Az = -(l->xx + l->yy) * z + l->xz * x + l->yz * y;
D = 0.2 * (Dx * x + Dy * y + Dz * z);
C = 0.25 * (Cx * x + Cy * y + Cz * z);
B = onethird * (Bx * x + By * y + Bz * z);
A = 0.5 * (Ax * x + Ay * y + Az * z);
g1 = x * l->x + y * l->y + z * l->z;
*ax -= iv * (l->x + Ax + Bx + Cx + Dx);
*ay -= iv * (l->y + Ay + By + Cy + Dy);
*az -= iv * (l->z + Az + Bz + Cz + Dz);
*fPot += l->m + g1 + A + B + C + D;
}
void momPrintMomc(MOMC *m) {
printf("MOMC:%20.15g\n", (double)m->m);
printf(" xx:%20.15g yy:%20.15g zz:%20.15g\n", (double)m->xx,
(double)m->yy, (double)m->zz);
printf(" xy:%20.15g yz:%20.15g xz:%20.15g\n", (double)m->xy,
(double)m->yz, (double)m->xz);
printf(" xxx:%20.15g xyy:%20.15g xzz:%20.15g\n", (double)m->xxx,
(double)m->xyy, (double)m->xzz);
printf(" xxy:%20.15g yyy:%20.15g yzz:%20.15g\n", (double)m->xxy,
(double)m->yyy, (double)m->yzz);
printf(" xxz:%20.15g yyz:%20.15g zzz:%20.15g\n", (double)m->xxz,
(double)m->yyz, (double)m->zzz);
printf(" xyz:%20.15g\n", (double)m->xyz);
printf("xxxx:%20.15g xxxy:%20.15g xxxz:%20.15g\n", (double)m->xxxx,
(double)m->xxxy, (double)m->xxxz);
printf("xyyy:%20.15g yyyy:%20.15g yyyz:%20.15g\n", (double)m->xyyy,
(double)m->yyyy, (double)m->yyyz);
printf("xzzz:%20.15g yzzz:%20.15g zzzz:%20.15g\n", (double)m->xzzz,
(double)m->yzzz, (double)m->zzzz);
printf("xxyy:%20.15g xxyz:%20.15g xyyz:%20.15g\n", (double)m->xxyy,
(double)m->xxyz, (double)m->xyyz);
printf("yyzz:%20.15g xxzz:%20.15g xyzz:%20.15g\n", (double)m->yyzz,
(double)m->xxzz, (double)m->xyzz);
}
void momPrintMomr(MOMR *m) {
printf("MOMR:%20.15g\n", (double)m->m);
printf(" xx:%20.15g yy:%20.15g\n", (double)m->xx, (double)m->yy);
printf(" xy:%20.15g yz:%20.15g xz:%20.15g\n", (double)m->xy,
(double)m->yz, (double)m->xz);
printf(" xxx:%20.15g xyy:%20.15g\n", (double)m->xxx, (double)m->xyy);
printf(" xxy:%20.15g yyy:%20.15g\n", (double)m->xxy, (double)m->yyy);
printf(" xxz:%20.15g yyz:%20.15g\n", (double)m->xxz, (double)m->yyz);
printf(" xyz:%20.15g\n", (double)m->xyz);
printf("xxxx:%20.15g xxxy:%20.15g xxxz:%20.15g\n", (double)m->xxxx,
(double)m->xxxy, (double)m->xxxz);
printf("xyyy:%20.15g yyyy:%20.15g yyyz:%20.15g\n", (double)m->xyyy,
(double)m->yyyy, (double)m->yyyz);
printf("xxyy:%20.15g xxyz:%20.15g xyyz:%20.15g\n", (double)m->xxyy,
(double)m->xxyz, (double)m->xyyz);
}
