| /*************************************************************************/ |
| /* */ |
| /* Copyright (c) 1994 Stanford University */ |
| /* */ |
| /* All rights reserved. */ |
| /* */ |
| /* Permission is given to use, copy, and modify this software for any */ |
| /* non-commercial purpose as long as this copyright notice is not */ |
| /* removed. All other uses, including redistribution in whole or in */ |
| /* part, are forbidden without prior written permission. */ |
| /* */ |
| /* This software is provided with absolutely no warranty and no */ |
| /* support. */ |
| /* */ |
| /*************************************************************************/ |
| |
| /* |
| * GRAV.C: |
| */ |
| |
| EXTERN_ENV |
| #define global extern |
| |
| #include "code.h" |
| |
| /* |
| * HACKGRAV: evaluate grav field at a given particle. |
| */ |
| |
| hackgrav(p,ProcessId) |
| bodyptr p; |
| unsigned ProcessId; |
| |
| { |
| extern gravsub(); |
| |
| Local[ProcessId].pskip = p; |
| SETV(Local[ProcessId].pos0, Pos(p)); |
| Local[ProcessId].phi0 = 0.0; |
| CLRV(Local[ProcessId].acc0); |
| Local[ProcessId].myn2bterm = 0; |
| Local[ProcessId].mynbcterm = 0; |
| Local[ProcessId].skipself = FALSE; |
| hackwalk(gravsub, ProcessId); |
| Phi(p) = Local[ProcessId].phi0; |
| SETV(Acc(p), Local[ProcessId].acc0); |
| #ifdef QUADPOLE |
| Cost(p) = Local[ProcessId].myn2bterm + NDIM * Local[ProcessId].mynbcterm; |
| #else |
| Cost(p) = Local[ProcessId].myn2bterm + Local[ProcessId].mynbcterm; |
| #endif |
| } |
| |
| |
| |
| /* |
| * GRAVSUB: compute a single body-body or body-cell interaction. |
| */ |
| |
| gravsub(p, ProcessId, level) |
| register nodeptr p; /* body or cell to interact with */ |
| unsigned ProcessId; |
| int level; |
| { |
| double sqrt(); |
| real drabs, phii, mor3; |
| vector ai, quaddr; |
| real dr5inv, phiquad, drquaddr; |
| |
| if (p != Local[ProcessId].pmem) { |
| SUBV(Local[ProcessId].dr, Pos(p), Local[ProcessId].pos0); |
| DOTVP(Local[ProcessId].drsq, Local[ProcessId].dr, Local[ProcessId].dr); |
| } |
| |
| Local[ProcessId].drsq += epssq; |
| drabs = sqrt((double) Local[ProcessId].drsq); |
| phii = Mass(p) / drabs; |
| Local[ProcessId].phi0 -= phii; |
| mor3 = phii / Local[ProcessId].drsq; |
| MULVS(ai, Local[ProcessId].dr, mor3); |
| ADDV(Local[ProcessId].acc0, Local[ProcessId].acc0, ai); |
| if(Type(p) != BODY) { /* a body-cell/leaf interaction? */ |
| Local[ProcessId].mynbcterm++; |
| #ifdef QUADPOLE |
| dr5inv = 1.0/(Local[ProcessId].drsq * Local[ProcessId].drsq * drabs); |
| MULMV(quaddr, Quad(p), Local[ProcessId].dr); |
| DOTVP(drquaddr, Local[ProcessId].dr, quaddr); |
| phiquad = -0.5 * dr5inv * drquaddr; |
| Local[ProcessId].phi0 += phiquad; |
| phiquad = 5.0 * phiquad / Local[ProcessId].drsq; |
| MULVS(ai, Local[ProcessId].dr, phiquad); |
| SUBV(Local[ProcessId].acc0, Local[ProcessId].acc0, ai); |
| MULVS(quaddr, quaddr, dr5inv); |
| SUBV(Local[ProcessId].acc0, Local[ProcessId].acc0, quaddr); |
| #endif |
| } |
| else { /* a body-body interaction */ |
| Local[ProcessId].myn2bterm++; |
| } |
| } |
| |
| /* |
| * HACKWALK: walk the tree opening cells too close to a given point. |
| */ |
| |
| local proced hacksub; |
| |
| hackwalk(sub, ProcessId) |
| proced sub; /* routine to do calculation */ |
| unsigned ProcessId; |
| { |
| walksub(Global->G_root, Global->rsize * Global->rsize, ProcessId); |
| } |
| |
| /* |
| * WALKSUB: recursive routine to do hackwalk operation. |
| */ |
| |
| walksub(n, dsq, ProcessId) |
| nodeptr n; /* pointer into body-tree */ |
| real dsq; /* size of box squared */ |
| unsigned ProcessId; |
| { |
| bool subdivp(); |
| nodeptr* nn; |
| leafptr l; |
| bodyptr p; |
| int i; |
| |
| if (subdivp(n, dsq, ProcessId)) { |
| if (Type(n) == CELL) { |
| for (nn = Subp(n); nn < Subp(n) + NSUB; nn++) { |
| if (*nn != NULL) { |
| walksub(*nn, dsq / 4.0, ProcessId); |
| } |
| } |
| } |
| else { |
| l = (leafptr) n; |
| for (i = 0; i < l->num_bodies; i++) { |
| p = Bodyp(l)[i]; |
| if (p != Local[ProcessId].pskip) { |
| gravsub(p, ProcessId); |
| } |
| else { |
| Local[ProcessId].skipself = TRUE; |
| } |
| } |
| } |
| } |
| else { |
| gravsub(n, ProcessId); |
| } |
| } |
| |
| /* |
| * SUBDIVP: decide if a node should be opened. |
| * Side effects: sets pmem,dr, and drsq. |
| */ |
| |
| bool subdivp(p, dsq, ProcessId) |
| register nodeptr p; /* body/cell to be tested */ |
| real dsq; /* size of cell squared */ |
| unsigned ProcessId; |
| { |
| SUBV(Local[ProcessId].dr, Pos(p), Local[ProcessId].pos0); |
| DOTVP(Local[ProcessId].drsq, Local[ProcessId].dr, Local[ProcessId].dr); |
| Local[ProcessId].pmem = p; |
| return (tolsq * Local[ProcessId].drsq < dsq); |
| } |