00001 typedef unsigned long size_t;
00002 typedef long fpos_t;
00003 #pragma pack( 8 )
00004 struct _Kee1 {
00005 int _cnt;
00006 unsigned char *_ptr;
00007 unsigned char *_base;
00008 int _bufsiz;
00009 short _flag;
00010 short _file;
00011 char *__newbase;
00012 void *_lock;
00013 unsigned char *_bufendp;
00014 } ;
00015 typedef struct _Kee1 FILE;
00016 extern struct _Kee1 _iob[];
00017 extern int fread( );
00018 extern int fwrite( );
00019 extern int _flsbuf( );
00020 extern int _filbuf( );
00021 extern int ferror( );
00022 extern int feof( );
00023 extern void clearerr( );
00024 extern int putchar( );
00025 extern int getchar( );
00026 extern int putc( );
00027 extern int getc( );
00028 extern int remove( );
00029 extern int rename( );
00030 extern struct _Kee1 *tmpfile( );
00031 extern char *tmpnam( );
00032 extern int fclose( );
00033 extern int fflush( );
00034 extern struct _Kee1 *fopen( );
00035 extern struct _Kee1 *freopen( );
00036 extern void setbuf( );
00037 extern int setvbuf( );
00038 extern int fprintf( );
00039 extern int fscanf( );
00040 extern int printf( );
00041 extern int scanf( );
00042 extern int sprintf( );
00043 extern int sscanf( );
00044 struct _Kee2 {
00045 char * *_a0;
00046 int _offset;
00047 } ;
00048 typedef struct _Kee2 va_list;
00049 extern int vfprintf( );
00050 extern int vprintf( );
00051 extern int vsprintf( );
00052 extern int fgetc( );
00053 extern char *fgets( );
00054 extern int fputc( );
00055 extern int fputs( );
00056 extern char *gets( );
00057 extern int puts( );
00058 extern int ungetc( );
00059 extern int fgetpos( );
00060 extern int fseek( );
00061 extern int fsetpos( );
00062 extern long ftell( );
00063 extern void rewind( );
00064 extern void perror( );
00065 typedef signed long ptrdiff_t;
00066 typedef unsigned int wchar_t;
00067 typedef unsigned int wctype_t;
00068 typedef int time_t;
00069 typedef int clock_t;
00070 typedef long ssize_t;
00071 typedef unsigned char uchar_t;
00072 typedef unsigned short ushort_t;
00073 typedef unsigned int uint_t;
00074 typedef unsigned long ulong_t;
00075 typedef volatile unsigned char vuchar_t;
00076 typedef volatile unsigned short vushort_t;
00077 typedef volatile unsigned int vuint_t;
00078 typedef volatile unsigned long vulong_t;
00079 struct _Kee3 {
00080 long r[1];
00081 } ;
00082 typedef struct _Kee3 *physadr_t;
00083 struct label_t {
00084 long val[10];
00085 } ;
00086 typedef struct label_t label_t;
00087 typedef int level_t;
00088 typedef int daddr_t;
00089 typedef char *caddr_t;
00090 typedef long *qaddr_t;
00091 typedef char *addr_t;
00092 typedef unsigned int ino_t;
00093 typedef short cnt_t;
00094 typedef int dev_t;
00095 typedef int chan_t;
00096 typedef long off_t;
00097 typedef unsigned long rlim_t;
00098 typedef int paddr_t;
00099 typedef unsigned short nlink_t;
00100 typedef int key_t;
00101 typedef unsigned int mode_t;
00102 typedef unsigned int uid_t;
00103 typedef unsigned int gid_t;
00104 typedef void *mid_t;
00105 typedef int pid_t;
00106 typedef char slab_t[12];
00107 typedef unsigned long shmatt_t;
00108 typedef unsigned long msgqnum_t;
00109 typedef unsigned long msglen_t;
00110 typedef unsigned int wint_t;
00111 typedef unsigned long sigset_t;
00112 typedef long timer_t;
00113 typedef void (*sig_t)( );
00114 typedef int id_t;
00115 typedef unsigned int major_t;
00116 typedef unsigned int minor_t;
00117 typedef unsigned int devs_t;
00118 typedef unsigned int unit_t;
00119 typedef unsigned long vm_offset_t;
00120 typedef unsigned long vm_size_t;
00121 typedef unsigned char uchar;
00122 typedef unsigned short ushort;
00123 typedef unsigned int uint;
00124 typedef unsigned long ulong;
00125 typedef struct _Kee3 *physadr;
00126 typedef unsigned char u_char;
00127 typedef unsigned short u_short;
00128 typedef unsigned int u_int;
00129 typedef unsigned long u_long;
00130 typedef volatile unsigned char vu_char;
00131 typedef volatile unsigned short vu_short;
00132 typedef volatile unsigned int vu_int;
00133 typedef volatile unsigned long vu_long;
00134 struct _quad {
00135 int val[2];
00136 } ;
00137 typedef struct _quad quad;
00138 typedef long swblk_t;
00139 typedef unsigned long fixpt_t;
00140 typedef int fd_mask;
00141 struct fd_set {
00142 int fds_bits[128];
00143 } ;
00144 typedef struct fd_set fd_set;
00145 extern void bzero( );
00146 #pragma pack( 0 )
00147 struct timeval ;
00148 int select( );
00149 extern int fileno( );
00150 extern struct _Kee1 *fdopen( );
00151 extern char *cuserid( );
00152 extern int getopt( );
00153 extern char *optarg;
00154 extern int optind;
00155 extern int optopt;
00156 extern int opterr;
00157 extern char *ctermid( );
00158 extern int getw( );
00159 extern int pclose( );
00160 extern int putw( );
00161 extern struct _Kee1 *popen( );
00162 extern char *tempnam( );
00163 extern void setbuffer( );
00164 extern void setlinebuf( );
00165 unsigned int h3wait_( );
00166 int npipes = 3;
00167 double xjmem[1000][3];
00168 double vjmem[1000][3];
00169 double xjpred[1000][3];
00170 double vjpred[1000][3];
00171 double ajmem[1000][3];
00172 double jjmem[1000][3];
00173 double mjmem[1000];
00174 double tjmem[1000];
00175 double xjbuf[50][26][3];
00176 double vjbuf[50][26][3];
00177 double ajbuf[50][26][3];
00178 double jjbuf[50][26][3];
00179 double mjbuf[50][26];
00180 double tjbuf[50][26];
00181 int indexbuf[50][26];
00182 int njbuf[50];
00183 double timem;
00184 double ximem[3][3];
00185 double vimem[3][3];
00186 double hmem[3];
00187 double epsmem[3];
00188 int nnb[3];
00189 int nb[3][1024];
00190 static int debug_level = 0;
00191 void h3open_( )
00192
00193 {
00194
00195 }
00196 void h3close_( )
00197
00198 {
00199
00200 }
00201 int h3npipe_( )
00202
00203 {
00204
00205 return npipes;
00206 }
00207 void h3setnboards_( nb )
00208 int *nb;
00209
00210 {
00211
00212 }
00213 int h3getnboards_( )
00214
00215 {
00216
00217 return 1;
00218 }
00219 unsigned int h3wait_( )
00220
00221 {
00222
00223 return 0;
00224 }
00225 void h3setti_( ti )
00226 double *ti;
00227
00228 {
00229
00230 timem = ti[0];
00231 }
00232 int h3jpmax_( )
00233
00234 {
00235
00236 return 26;
00237 }
00238 void h3setmode_( )
00239
00240 {
00241
00242 }
00243 void h3jpdma_indirect_( int *, int *, double [][3], double [][3], double [][3], double [][3], double *, double *, int * );
00244 void h3mjpdma_indirect_( int *, int *, double [][3], double [][3], double [][3], double [][3], double *, double *, int *,
00245 int * );
00246 void h3mjpdma_start_( int * );
00247 void h3mjpdma_flush_( );
00248 void h3jpdma_indirect_( nj, hostindex, xj, vj, aj, jj, mj, tj, mode )
00249 int *nj;
00250 int *hostindex;
00251 double (*xj)[3];
00252 double (*vj)[3];
00253 double (*aj)[3];
00254 double (*jj)[3];
00255 double *mj;
00256 double *tj;
00257 int *mode;
00258
00259 {
00260 int zero;
00261
00262 zero = 0;
00263 h3mjpdma_indirect_( nj, hostindex, xj, vj, aj, jj, mj, tj, mode, &zero );
00264 h3mjpdma_start_( &zero );
00265 }
00266 static int buff_used = 0;
00267 void h3mjpdma_indirect_( nj, hostindex, xj, vj, aj, jj, mj, tj, mode, buff_id )
00268 int *nj;
00269 int *hostindex;
00270 double (*xj)[3];
00271 double (*vj)[3];
00272 double (*aj)[3];
00273 double (*jj)[3];
00274 double *mj;
00275 double *tj;
00276 int *mode;
00277 int *buff_id;
00278
00279 {
00280 register int i;
00281 int bufno;
00282 int k;
00283 int hostid;
00284 long _Kii1;
00285
00286 if (debug_level) {
00287 fprintf( (_iob + 2), "(h3mjpdma) called %d %d\n", (nj[0]), (mode[0]) );
00288 }
00289 if (nj[0] > 26) {
00290 fprintf( (_iob + 2), "(h3jpdma) Too large nj (%d > %d\n", nj, 26 );
00291 exit( ((-1)) );
00292 }
00293 if (buff_id[0] >= 0) {
00294 bufno = buff_id[0];
00295 } else {
00296 bufno = buff_used;
00297 buff_used++;
00298 }
00299 njbuf[bufno] = nj[0];
00300 i = 0;
00301 _Kii1 = ((void *)(xjbuf[bufno][i] + 2) < (void *)(vjbuf[bufno][i]) || (void *)(xjbuf[bufno][i]) > (void *)(vjbuf[bufno][i]
00302 + 2)) && ((void *)(xjbuf[bufno][i] + 2) < (void *)(ajbuf[bufno][i]) || (void *)(xjbuf[bufno][i]) > (void *)(ajbuf[
00303 bufno][i] + 2)) && ((void *)(xjbuf[bufno][i] + 2) < (void *)(jjbuf[bufno][i]) || (void *)(xjbuf[bufno][i]) > (void *)(
00304 jjbuf[bufno][i] + 2)) && ((void *)(xjbuf[bufno][i] + 2) < (void *)(xj[hostid]) || (void *)(xjbuf[bufno][i]) > (void *)
00305 (xj[hostid] + 2)) && ((void *)(xjbuf[bufno][i] + 2) < (void *)(vj[hostid]) || (void *)(xjbuf[bufno][i]) > (void *)(vj[
00306 hostid] + 2)) && ((void *)(xjbuf[bufno][i] + 2) < (void *)(aj[hostid]) || (void *)(xjbuf[bufno][i]) > (void *)(aj[
00307 hostid] + 2)) && ((void *)(xjbuf[bufno][i] + 2) < (void *)(jj[hostid]) || (void *)(xjbuf[bufno][i]) > (void *)(jj[
00308 hostid] + 2)) && ((void *)(vjbuf[bufno][i] + 2) < (void *)(ajbuf[bufno][i]) || (void *)(vjbuf[bufno][i]) > (void *)(
00309 ajbuf[bufno][i] + 2)) && ((void *)(vjbuf[bufno][i] + 2) < (void *)(jjbuf[bufno][i]) || (void *)(vjbuf[bufno][i]) > (
00310 void *)(jjbuf[bufno][i] + 2)) && ((void *)(vjbuf[bufno][i] + 2) < (void *)(xj[hostid]) || (void *)(vjbuf[bufno][i]) >
00311 (void *)(xj[hostid] + 2)) && ((void *)(vjbuf[bufno][i] + 2) < (void *)(vj[hostid]) || (void *)(vjbuf[bufno][i]) > (
00312 void *)(vj[hostid] + 2)) && ((void *)(vjbuf[bufno][i] + 2) < (void *)(aj[hostid]) || (void *)(vjbuf[bufno][i]) > (
00313 void *)(aj[hostid] + 2)) && ((void *)(vjbuf[bufno][i] + 2) < (void *)(jj[hostid]) || (void *)(vjbuf[bufno][i]) > (
00314 void *)(jj[hostid] + 2)) && ((void *)(ajbuf[bufno][i] + 2) < (void *)(jjbuf[bufno][i]) || (void *)(ajbuf[bufno][i]) >
00315 (void *)(jjbuf[bufno][i] + 2)) && ((void *)(ajbuf[bufno][i] + 2) < (void *)(xj[hostid]) || (void *)(ajbuf[bufno][i])
00316 > (void *)(xj[hostid] + 2)) && ((void *)(ajbuf[bufno][i] + 2) < (void *)(vj[hostid]) || (void *)(ajbuf[bufno][i]) > (
00317 void *)(vj[hostid] + 2)) && ((void *)(ajbuf[bufno][i] + 2) < (void *)(aj[hostid]) || (void *)(ajbuf[bufno][i]) > (
00318 void *)(aj[hostid] + 2)) && ((void *)(ajbuf[bufno][i] + 2) < (void *)(jj[hostid]) || (void *)(ajbuf[bufno][i]) > (
00319 void *)(jj[hostid] + 2)) && ((void *)(jjbuf[bufno][i] + 2) < (void *)(xj[hostid]) || (void *)(jjbuf[bufno][i]) > (
00320 void *)(xj[hostid] + 2)) && ((void *)(jjbuf[bufno][i] + 2) < (void *)(vj[hostid]) || (void *)(jjbuf[bufno][i]) > (
00321 void *)(vj[hostid] + 2)) && ((void *)(jjbuf[bufno][i] + 2) < (void *)(aj[hostid]) || (void *)(jjbuf[bufno][i]) > (
00322 void *)(aj[hostid] + 2)) && ((void *)(jjbuf[bufno][i] + 2) < (void *)(jj[hostid]) || (void *)(jjbuf[bufno][i]) > (
00323 void *)(jj[hostid] + 2));
00324 if (!_Kii1) {
00325 while ( i < nj[0] ) {
00326 hostid = hostindex[i] - 1;
00327 xjbuf[bufno][i][0] = xj[hostid][0];
00328 vjbuf[bufno][i][0] = vj[hostid][0];
00329 ajbuf[bufno][i][0] = aj[hostid][0];
00330 jjbuf[bufno][i][0] = jj[hostid][0];
00331 xjbuf[bufno][i][1] = xj[hostid][1];
00332 vjbuf[bufno][i][1] = vj[hostid][1];
00333 ajbuf[bufno][i][1] = aj[hostid][1];
00334 jjbuf[bufno][i][1] = jj[hostid][1];
00335 xjbuf[bufno][i][2] = xj[hostid][2];
00336 vjbuf[bufno][i][2] = vj[hostid][2];
00337 ajbuf[bufno][i][2] = aj[hostid][2];
00338 jjbuf[bufno][i][2] = jj[hostid][2];
00339 tjbuf[bufno][i] = tj[hostid];
00340 mjbuf[bufno][i] = mj[hostid];
00341 indexbuf[bufno][i] = hostid;
00342 i++;
00343 }
00344 } else {
00345 while ( i < nj[0] ) {
00346 hostid = hostindex[i] - 1;
00347 xjbuf[bufno][i][0] = xj[hostid][0];
00348 xjbuf[bufno][i][1] = xj[hostid][1];
00349 xjbuf[bufno][i][2] = xj[hostid][2];
00350 vjbuf[bufno][i][0] = vj[hostid][0];
00351 vjbuf[bufno][i][1] = vj[hostid][1];
00352 vjbuf[bufno][i][2] = vj[hostid][2];
00353 ajbuf[bufno][i][0] = aj[hostid][0];
00354 ajbuf[bufno][i][1] = aj[hostid][1];
00355 ajbuf[bufno][i][2] = aj[hostid][2];
00356 jjbuf[bufno][i][0] = jj[hostid][0];
00357 jjbuf[bufno][i][1] = jj[hostid][1];
00358 jjbuf[bufno][i][2] = jj[hostid][2];
00359 tjbuf[bufno][i] = tj[hostid];
00360 mjbuf[bufno][i] = mj[hostid];
00361 indexbuf[bufno][i] = hostid;
00362 i++;
00363 }
00364 }
00365 }
00366 void h3mjpdma_start_( buff_id )
00367 int *buff_id;
00368
00369 {
00370 register int i;
00371 register int k;
00372 register int j;
00373 register int buff;
00374 int _Kii2;
00375 int _Kii1;
00376
00377 if (debug_level > 0) {
00378 fprintf( (_iob + 2), "(h3mjpdma_start) called %d \n", (buff_id[0]) );
00379 }
00380 buff = buff_id[0];
00381 for ( i = 0; i<njbuf[buff]; i++ ) {
00382 _Kii1 = indexbuf[buff][i];
00383 if (debug_level > 2) {
00384 printf( "jpdma_start %d %d %d %f\n", i, _Kii1, buff, xjbuf[buff][i][0] );
00385 }
00386 xjmem[_Kii1][0] = xjbuf[buff][i][0];
00387 xjmem[_Kii1][1] = xjbuf[buff][i][1];
00388 xjmem[_Kii1][2] = xjbuf[buff][i][2];
00389 vjmem[_Kii1][0] = vjbuf[buff][i][0];
00390 vjmem[_Kii1][1] = vjbuf[buff][i][1];
00391 vjmem[_Kii1][2] = vjbuf[buff][i][2];
00392 ajmem[_Kii1][0] = ajbuf[buff][i][0];
00393 ajmem[_Kii1][1] = ajbuf[buff][i][1];
00394 ajmem[_Kii1][2] = ajbuf[buff][i][2];
00395 jjmem[_Kii1][0] = jjbuf[buff][i][0];
00396 jjmem[_Kii1][1] = jjbuf[buff][i][1];
00397 jjmem[_Kii1][2] = jjbuf[buff][i][2];
00398 tjmem[_Kii1] = tjbuf[buff][i];
00399 mjmem[_Kii1] = mjbuf[buff][i];
00400 }
00401 }
00402 void h3jpdma_flush_( )
00403
00404 {
00405 int i;
00406
00407 i = 0;
00408 while ( i < buff_used ) {
00409 h3mjpdma_start_( &i );
00410 i++;
00411 }
00412 buff_used = 0;
00413 }
00414 void h3jpdma_reset_( )
00415
00416 {
00417
00418 buff_used = 0;
00419 }
00420 static int njset;
00421 int h3calc_firsthalf_( nj, ni, xi, vi, eps2, h2 )
00422 int *nj;
00423 int *ni;
00424 double (*xi)[3];
00425 double (*vi)[3];
00426 double *eps2;
00427 double *h2;
00428
00429 {
00430 int i;
00431 int j;
00432 int k;
00433 double s;
00434 double s1;
00435 double s2;
00436 long _Kii1;
00437 long _Kii2;
00438
00439 njset = nj[0];
00440 if (ni[0] > 3) {
00441 fprintf( (_iob + 2), "(h3ipdma) Too large ni (%d > %d\n", ni, 3 );
00442 exit( ((-1)) );
00443 }
00444 i = 0;
00445 _Kii1 = (void *)(ximem[i] + 2) < (void *)(xi[i]) || (void *)(ximem[i]) > (void *)(xi[i] + 2);
00446 _Kii2 = (void *)(vimem[i] + 2) < (void *)(vi[i]) || (void *)(vimem[i]) > (void *)(vi[i] + 2);
00447 if (!_Kii1) {
00448 if (!_Kii2) {
00449 while ( i < ni[0] ) {
00450 ximem[i][0] = xi[i][0];
00451 ximem[i][1] = xi[i][1];
00452 ximem[i][2] = xi[i][2];
00453 vimem[i][0] = vi[i][0];
00454 vimem[i][1] = vi[i][1];
00455 vimem[i][2] = vi[i][2];
00456 hmem[i] = h2[i] * 2;
00457 epsmem[i] = eps2[i] * 2;
00458 i++;
00459 }
00460 } else {
00461 while ( i < ni[0] ) {
00462 ximem[i][0] = xi[i][0];
00463 ximem[i][1] = xi[i][1];
00464 ximem[i][2] = xi[i][2];
00465 vimem[i][0] = vi[i][0];
00466 vimem[i][1] = vi[i][1];
00467 vimem[i][2] = vi[i][2];
00468 hmem[i] = h2[i] * 2;
00469 epsmem[i] = eps2[i] * 2;
00470 i++;
00471 }
00472 }
00473 } else {
00474 if (!_Kii2) {
00475 while ( i < ni[0] ) {
00476 ximem[i][0] = xi[i][0];
00477 ximem[i][1] = xi[i][1];
00478 ximem[i][2] = xi[i][2];
00479 vimem[i][0] = vi[i][0];
00480 vimem[i][1] = vi[i][1];
00481 vimem[i][2] = vi[i][2];
00482 hmem[i] = h2[i] * 2;
00483 epsmem[i] = eps2[i] * 2;
00484 i++;
00485 }
00486 } else {
00487 while ( i < ni[0] ) {
00488 ximem[i][0] = xi[i][0];
00489 ximem[i][1] = xi[i][1];
00490 ximem[i][2] = xi[i][2];
00491 vimem[i][0] = vi[i][0];
00492 vimem[i][1] = vi[i][1];
00493 vimem[i][2] = vi[i][2];
00494 hmem[i] = h2[i] * 2;
00495 epsmem[i] = eps2[i] * 2;
00496 i++;
00497 }
00498 }
00499 }
00500 if (debug_level > 0) {
00501 i = 0;
00502 while ( i < ni[0] ) {
00503 if (hmem[i] != hmem[i+1] != 0) {
00504 fprintf( (_iob + 2), "h3calc_firsthalf, eps/h2 error for %d\n", i );
00505 }
00506 i += 2;
00507 }
00508 }
00509 j = 0;
00510 while ( j < nj[0] ) {
00511 s = timem - tjmem[j];
00512 s1 = 1.5 * s;
00513 s2 = s * 2.e0;
00514 for ( k = 0; k<=2; k++ ) {
00515 xjpred[j][k] = xjmem[j][k] + (vjmem[j][k] + (ajmem[j][k] + jjmem[j][k] * s) * s) * s;
00516 vjpred[j][k] = vjmem[j][k] + (ajmem[j][k] + jjmem[j][k] * s1) * s2;
00517 if (debug_level > 2) {
00518 printf( "predict %d %d %22.14g %22.14g\n", j, k, xjpred[j][k], vjpred[j][k] );
00519 }
00520 }
00521 j++;
00522 }
00523 }
00524 int calculate_force_on_one_particle( i, acc, jerk, pot )
00525 int i;
00526 double (*acc)[3];
00527 double (*jerk)[3];
00528 double *pot;
00529
00530 {
00531 register int j;
00532 register int k;
00533 double sqrt( );
00534 register double xi;
00535 register double yi;
00536 register double zi;
00537 register double vxi;
00538 register double vyi;
00539 register double vzi;
00540 register double dx;
00541 register double dy;
00542 register double dz;
00543 register double dvx;
00544 register double dvy;
00545 register double dvz;
00546 register double dr2;
00547 register double drdv;
00548 register double dr2i;
00549 register double dr3i;
00550 register double ax;
00551 register double ay;
00552 register double az;
00553 register double jx;
00554 register double jy;
00555 register double jz;
00556 register double epsi;
00557
00558 jz = 0.e0;
00559 jy = 0.e0;
00560 jx = 0.e0;
00561 az = 0.e0;
00562 ay = 0.e0;
00563 ax = 0.e0;
00564 nnb[i] = 0;
00565 pot[i] = 0.e0;
00566 xi = ximem[i][0];
00567 yi = ximem[i][1];
00568 zi = ximem[i][2];
00569 vxi = vimem[i][0];
00570 vyi = vimem[i][1];
00571 vzi = vimem[i][2];
00572 epsi = epsmem[i];
00573 j = 0;
00574 while ( j < njset ) {
00575 dx = xjpred[j][0] - xi;
00576 dy = xjpred[j][1] - yi;
00577 dz = xjpred[j][2] - zi;
00578 dvx = vjpred[j][0] - vxi;
00579 dvy = vjpred[j][1] - vyi;
00580 dvz = vjpred[j][2] - vzi;
00581 dr2 = epsi + dx * dx + dy * dy + dz * dz;
00582 drdv = dx * dvx + dy * dvy + dz * dvz;
00583 if (dr2 < hmem[i]) {
00584 nb[i][nnb[i]] = j + 1;
00585 nnb[i]++;
00586 }
00587 if (dr2 > 0.e0) {
00588 dr2i = 1.e0 / dr2;
00589 dr3i = mjmem[j] * dr2i * sqrt( dr2i );
00590 drdv *= 3.e0 * dr2i;
00591 ax += dx * dr3i;
00592 ay += dy * dr3i;
00593 az += dz * dr3i;
00594 jx += (dvx - dx * drdv) * dr3i;
00595 jy += (dvy - dy * drdv) * dr3i;
00596 jz += (dvz - dz * drdv) * dr3i;
00597 }
00598 j++;
00599 }
00600 acc[i][0] = ax;
00601 acc[i][1] = ay;
00602 acc[i][2] = az;
00603 jerk[i][0] = jx;
00604 jerk[i][1] = jy;
00605 jerk[i][2] = jz;
00606 }
00607 int h3calc_lasthalf_( ni, acc, jerk, pot )
00608 int *ni;
00609 double (*acc)[3];
00610 double (*jerk)[3];
00611 double *pot;
00612
00613 {
00614 register int i;
00615 register int j;
00616 register int k;
00617 int _Kii1;
00618
00619 if (ni[0] > 3) {
00620 fprintf( (_iob + 2), "(h3ipdma) Too large ni (%d > %d\n", ni, 3 );
00621 exit( ((-1)) );
00622 }
00623 if (debug_level) {
00624 fprintf( (_iob + 2), "lasthalf ni=%d\n", (ni[0]) );
00625 }
00626 for ( _Kii1 = 0; _Kii1<ni[0]; _Kii1++ ) {
00627 calculate_force_on_one_particle( _Kii1, acc, jerk, pot );
00628 }
00629 }
00630 int h3calc_lasthalf_old( ni, acc, jerk, pot )
00631 int *ni;
00632 double (*acc)[3];
00633 double (*jerk)[3];
00634 double *pot;
00635
00636 {
00637 register int i;
00638 register int j;
00639 register int k;
00640 double sqrt( );
00641 register double dx[3];
00642 register double dv[3];
00643 register double dr2;
00644 register double drdv;
00645 register double dr2i;
00646 register double dr3i;
00647 long _Kii1;
00648 int _Kii2;
00649
00650 if (ni[0] > 3) {
00651 fprintf( (_iob + 2), "(h3ipdma) Too large ni (%d > %d\n", ni, 3 );
00652 exit( ((-1)) );
00653 }
00654 if (debug_level) {
00655 fprintf( (_iob + 2), "lasthalf ni=%d\n", (ni[0]) );
00656 }
00657 i = 0;
00658 _Kii1 = (void *)(acc[i] + 2) < (void *)(jerk[i]) || (void *)(acc[i]) > (void *)(jerk[i] + 2);
00659 _Kii2 = !_Kii1;
00660 while ( i < ni[0] ) {
00661 if (_Kii2) {
00662 acc[i][0] = 0.e0;
00663 acc[i][1] = 0.e0;
00664 acc[i][2] = 0.e0;
00665 jerk[i][0] = 0.e0;
00666 jerk[i][1] = 0.e0;
00667 jerk[i][2] = 0.e0;
00668 } else {
00669 acc[i][0] = 0.e0;
00670 acc[i][1] = 0.e0;
00671 acc[i][2] = 0.e0;
00672 jerk[i][0] = 0.e0;
00673 jerk[i][1] = 0.e0;
00674 jerk[i][2] = 0.e0;
00675 }
00676 nnb[i] = 0;
00677 pot[i] = 0.e0;
00678 j = 0;
00679 while ( j < njset ) {
00680 dr2 = epsmem[i];
00681 drdv = 0.e0;
00682 for ( k = 0; k<=2; k++ ) {
00683 dx[k] = xjpred[j][k] - ximem[i][k];
00684 dv[k] = vjpred[j][k] - vimem[i][k];
00685 dr2 += dx[k] * dx[k];
00686 drdv += dx[k] * dv[k];
00687 if (debug_level > 2) {
00688 fprintf( (_iob + 2), "rijk %d %d %d %g %21.14g %21.14g\n", i, j, k, dx[k], xjpred[j][k], ximem[i][k] );
00689 }
00690 }
00691 if (dr2 < hmem[i]) {
00692 nb[i][nnb[i]] = j + 1;
00693 nnb[i]++;
00694 }
00695 if (debug_level > 1) {
00696 fprintf( (_iob + 2), "i, j, nnb[i], dr2, h = %d %d %d %g %g\n", i, j, nnb[i], dr2, hmem[i] );
00697 }
00698 if (dr2 > 0.e0) {
00699 dr2i = 1.e0 / dr2;
00700 dr3i = mjmem[j] * dr2i * sqrt( dr2i );
00701 drdv *= 3.e0 * dr2i;
00702 for ( k = 0; k<=2; k++ ) {
00703 acc[i][k] += dx[k] * dr3i;
00704 jerk[i][k] += (dv[k] - dx[k] * drdv) * dr3i;
00705 if (debug_level > 1) {
00706 fprintf( (_iob + 2), "force i,j,k %d %d %d %20.10g \n", i, j, k, acc[i][k] );
00707 }
00708 }
00709 } else {
00710 if (debug_level > 1) {
00711 fprintf( (_iob + 2), "force skip i,j %d %d \n", i, j );
00712 }
00713 }
00714 j++;
00715 }
00716 i++;
00717 }
00718 }
00719 int h3calc_( nj, ni, xi, vi, eps2, h2, acc, jerk, pot )
00720 int *nj;
00721 int *ni;
00722 double (*xi)[3];
00723 double (*vi)[3];
00724 double *eps2;
00725 double *h2;
00726 double (*acc)[3];
00727 double (*jerk)[3];
00728 double *pot;
00729
00730 {
00731
00732 h3calc_firsthalf_( nj, ni, xi, vi, eps2, h2 );
00733 h3calc_lasthalf_( ni, acc, jerk, pot );
00734 }
00735 void h3nbread_( nboards )
00736 int *nboards;
00737
00738 {
00739
00740 }
00741 void set_debug_level_( int *level )
00742 {
00743
00744 debug_level = level[0];
00745 }
00746 void h3setdebuglevel_( int *level )
00747 {
00748
00749 set_debug_level_( level );
00750 }
00751 int compare( i, j )
00752 int *i;
00753 int *j;
00754
00755 {
00756
00757 return i[0] - j[0];
00758 }
00759 int h3nblist_( board, chip, nblist )
00760 int *board;
00761 int *chip;
00762 int *nblist;
00763
00764 {
00765 int i;
00766
00767 i = 0;
00768 while ( i < nnb[chip[0]] ) {
00769 nblist[i] = nb[chip[0]][i];
00770 i++;
00771 }
00772 return nnb[chip[0]];
00773 }
