1 /*____________________________________________________________________________
3 | libxdrf - portable fortran interface to xdr. some xdr routines
4 | are C routines for compressed coordinates
8 | This collection of routines is intended to write and read
9 | data in a portable way to a file, so data written on one type
10 | of machine can be read back on a different type.
12 | all fortran routines use an integer 'xdrid', which is an id to the
13 | current xdr file, and is set by xdrfopen.
14 | most routines have in integer 'ret' which is the return value.
15 | The value of 'ret' is zero on failure, and most of the time one
18 | There are three routines useful for C users:
19 | xdropen(), xdrclose(), xdr3dfcoord().
20 | The first two replace xdrstdio_create and xdr_destroy, and *must* be
21 | used when you plan to use xdr3dfcoord(). (they are also a bit
22 | easier to interface). For writing data other than compressed coordinates
23 | you should use the standard C xdr routines (see xdr man page)
25 | xdrfopen(xdrid, filename, mode, ret)
26 | character *(*) filename
29 | this will open the file with the given filename (string)
30 | and the given mode, it returns an id in xdrid, which is
31 | to be used in all other calls to xdrf routines.
32 | mode is 'w' to create, or update an file, for all other
33 | values of mode the file is opened for reading
35 | you need to call xdrfclose to flush the output and close
37 | Note that you should not use xdrstdio_create, which comes with the
38 | standard xdr library
40 | xdrfclose(xdrid, ret)
41 | flush the data to the file, and closes the file;
42 | You should not use xdr_destroy (which comes standard with
45 | xdrfbool(xdrid, bp, ret)
48 | This filter produces values of either 1 or 0
50 | xdrfchar(xdrid, cp, ret)
53 | filter that translate between characters and their xdr representation
54 | Note that the characters in not compressed and occupies 4 bytes.
56 | xdrfdouble(xdrid, dp, ret)
59 | read/write a double.
61 | xdrffloat(xdrid, fp, ret)
66 | xdrfint(xdrid, ip, ret)
71 | xdrflong(xdrid, lp, ret)
74 | this routine has a possible portablility problem due to 64 bits longs.
76 | xdrfshort(xdrid, sp, ret)
79 | xdrfstring(xdrid, sp, maxsize, ret)
83 | read/write a string, with maximum length given by maxsize
85 | xdrfwrapstring(xdris, sp, ret)
88 | read/write a string (it is the same as xdrfstring accept that it finds
89 | the stringlength itself.
91 | xdrfvector(xdrid, cp, size, xdrfproc, ret)
96 | read/write an array pointed to by cp, with number of elements
97 | defined by 'size'. the routine 'xdrfproc' is the name
98 | of one of the above routines to read/write data (like xdrfdouble)
99 | In contrast with the c-version you don't need to specify the
100 | byte size of an element.
101 | xdrfstring is not allowed here (it is in the c version)
103 | xdrf3dfcoord(xdrid, fp, size, precision, ret)
108 | this is *NOT* a standard xdr routine. I named it this way, because
109 | it invites people to use the other xdr routines.
110 | It is introduced to store specifically 3d coordinates of molecules
111 | (as found in molecular dynamics) and it writes it in a compressed way.
112 | It starts by multiplying all numbers by precision and
113 | rounding the result to integer. effectively converting
114 | all floating point numbers to fixed point.
115 | it uses an algorithm for compression that is optimized for
116 | molecular data, but could be used for other 3d coordinates
117 | as well. There is subtantial overhead involved, so call this
118 | routine only if you have a large number of coordinates to read/write
120 | ________________________________________________________________________
122 | Below are the routines to be used by C programmers. Use the 'normal'
123 | xdr routines to write integers, floats, etc (see man xdr)
125 | int xdropen(XDR *xdrs, const char *filename, const char *type)
126 | This will open the file with the given filename and the
127 | given mode. You should pass it an allocated XDR struct
128 | in xdrs, to be used in all other calls to xdr routines.
129 | Mode is 'w' to create, or update an file, and for all
130 | other values of mode the file is opened for reading.
131 | You need to call xdrclose to flush the output and close
134 | Note that you should not use xdrstdio_create, which
135 | comes with the standard xdr library.
137 | int xdrclose(XDR *xdrs)
138 | Flush the data to the file, and close the file;
139 | You should not use xdr_destroy (which comes standard
140 | with the xdr libraries).
142 | int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision)
143 | This is \fInot\fR a standard xdr routine. I named it this
144 | way, because it invites people to use the other xdr
147 | (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl
160 int ftocstr(char *, int, char *, int);
161 int ctofstr(char *, int, char *);
164 static FILE *xdrfiles[MAXID];
165 static XDR *xdridptr[MAXID];
166 static char xdrmodes[MAXID];
167 static unsigned int cnt;
169 typedef void (* FUNCTION(xdrfproc)) (int *, void *, int *);
172 FUNCTION(xdrfbool) ARGS(`xdrid, pb, ret')
176 *ret = xdr_bool(xdridptr[*xdrid], (bool_t *) pb);
181 FUNCTION(xdrfchar) ARGS(`xdrid, cp, ret')
185 *ret = xdr_char(xdridptr[*xdrid], cp);
190 FUNCTION(xdrfdouble) ARGS(`xdrid, dp, ret')
194 *ret = xdr_double(xdridptr[*xdrid], dp);
195 cnt += sizeof(double);
199 FUNCTION(xdrffloat) ARGS(`xdrid, fp, ret')
203 *ret = xdr_float(xdridptr[*xdrid], fp);
204 cnt += sizeof(float);
208 FUNCTION(xdrfint) ARGS(`xdrid, ip, ret')
212 *ret = xdr_int(xdridptr[*xdrid], ip);
217 FUNCTION(xdrflong) ARGS(`xdrid, lp, ret')
221 *ret = xdr_long(xdridptr[*xdrid], lp);
226 FUNCTION(xdrfshort) ARGS(`xdrid, sp, ret')
230 *ret = xdr_short(xdridptr[*xdrid], sp);
235 FUNCTION(xdrfuchar) ARGS(`xdrid, ucp, ret')
239 *ret = xdr_u_char(xdridptr[*xdrid], ucp);
244 FUNCTION(xdrfulong) ARGS(`xdrid, ulp, ret')
248 *ret = xdr_u_long(xdridptr[*xdrid], ulp);
249 cnt += sizeof(unsigned long);
253 FUNCTION(xdrfushort) ARGS(`xdrid, usp, ret')
257 *ret = xdr_u_short(xdridptr[*xdrid], usp);
258 cnt += sizeof(unsigned short);
262 FUNCTION(xdrf3dfcoord) ARGS(`xdrid, fp, size, precision, ret')
268 *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
272 FUNCTION(xdrfstring) ARGS(`xdrid, STRING_ARG(sp), maxsize, ret')
279 tsp = (char*) malloc(((STRING_LEN(sp)) + 1) * sizeof(char));
284 if (ftocstr(tsp, *maxsize+1, STRING_PTR(sp), STRING_LEN(sp))) {
289 *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int) *maxsize);
290 ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
296 FUNCTION(xdrfwrapstring) ARGS(`xdrid, STRING_ARG(sp), ret')
302 maxsize = (STRING_LEN(sp)) + 1;
303 tsp = (char*) malloc(maxsize * sizeof(char));
308 if (ftocstr(tsp, maxsize, STRING_PTR(sp), STRING_LEN(sp))) {
313 *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
314 ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
320 FUNCTION(xdrfopaque) ARGS(`xdrid, cp, ccnt, ret')
325 *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
330 FUNCTION(xdrfsetpos) ARGS(`xdrid, pos, ret')
334 *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
338 FUNCTION(xdrf) ARGS(`xdrid, pos')
341 *pos = xdr_getpos(xdridptr[*xdrid]);
345 FUNCTION(xdrfvector) ARGS(`xdrid, cp, size, elproc, ret')
349 FUNCTION(xdrfproc) elproc;
353 for (lcnt = 0; lcnt < *size; lcnt++) {
354 elproc(xdrid, (cp+cnt) , ret);
360 FUNCTION(xdrfclose) ARGS(`xdrid, ret')
364 *ret = xdrclose(xdridptr[*xdrid]);
369 FUNCTION(xdrfopen) ARGS(`xdrid, STRING_ARG(fp), STRING_ARG(mode), ret')
372 STRING_ARG_DECL(mode);
378 if (ftocstr(fname, sizeof(fname), STRING_PTR(fp), STRING_LEN(fp))) {
381 if (ftocstr(fmode, sizeof(fmode), STRING_PTR(mode),
386 *xdrid = xdropen(NULL, fname, fmode);
393 /*___________________________________________________________________________
395 | what follows are the C routines for opening, closing xdr streams
396 | and the routine to read/write compressed coordinates together
397 | with some routines to assist in this task (those are marked
398 | static and cannot be called from user programs)
400 #define MAXABS INT_MAX-2
403 #define MIN(x,y) ((x) < (y) ? (x):(y))
406 #define MAX(x,y) ((x) > (y) ? (x):(y))
409 #define SQR(x) ((x)*(x))
411 static int magicints[] = {
412 0, 0, 0, 0, 0, 0, 0, 0, 0,
413 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
414 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
415 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
416 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
417 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
418 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
419 8388607, 10568983, 13316085, 16777216 };
422 /* note that magicints[FIRSTIDX-1] == 0 */
423 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
426 /*__________________________________________________________________________
428 | xdropen - open xdr file
430 | This versions differs from xdrstdio_create, because I need to know
431 | the state of the file (read or write) so I can use xdr3dfcoord
432 | in eigther read or write mode, and the file descriptor
433 | so I can close the file (something xdr_destroy doesn't do).
437 int xdropen(XDR *xdrs, const char *filename, const char *type) {
438 static int init_done = 0;
443 if (init_done == 0) {
444 for (xdrid = 1; xdrid < MAXID; xdrid++) {
445 xdridptr[xdrid] = NULL;
450 while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
453 if (xdrid == MAXID) {
456 if (*type == 'w' || *type == 'W') {
465 xdrfiles[xdrid] = fopen(filename, type1);
466 if (xdrfiles[xdrid] == NULL) {
470 xdrmodes[xdrid] = *type;
471 /* next test isn't usefull in the case of C language
472 * but is used for the Fortran interface
473 * (C users are expected to pass the address of an already allocated
477 xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
478 xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
480 xdridptr[xdrid] = xdrs;
481 xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
486 /*_________________________________________________________________________
488 | xdrclose - close a xdr file
490 | This will flush the xdr buffers, and destroy the xdr stream.
491 | It also closes the associated file descriptor (this is *not*
492 | done by xdr_destroy).
496 int xdrclose(XDR *xdrs) {
500 fprintf(stderr, "xdrclose: passed a NULL pointer\n");
503 for (xdrid = 1; xdrid < MAXID; xdrid++) {
504 if (xdridptr[xdrid] == xdrs) {
507 fclose(xdrfiles[xdrid]);
508 xdridptr[xdrid] = NULL;
512 fprintf(stderr, "xdrclose: no such open xdr file\n");
517 /*____________________________________________________________________________
519 | sendbits - encode num into buf using the specified number of bits
521 | This routines appends the value of num to the bits already present in
522 | the array buf. You need to give it the number of bits to use and you
523 | better make sure that this number of bits is enough to hold the value
524 | Also num must be positive.
528 static void sendbits(int buf[], int num_of_bits, int num) {
530 unsigned int cnt, lastbyte;
532 unsigned char * cbuf;
534 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
535 cnt = (unsigned int) buf[0];
537 lastbyte =(unsigned int) buf[2];
538 while (num_of_bits >= 8) {
539 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
540 cbuf[cnt++] = lastbyte >> lastbits;
543 if (num_of_bits > 0) {
544 lastbyte = (lastbyte << num_of_bits) | num;
545 lastbits += num_of_bits;
548 cbuf[cnt++] = lastbyte >> lastbits;
555 cbuf[cnt] = lastbyte << (8 - lastbits);
559 /*_________________________________________________________________________
561 | sizeofint - calculate bitsize of an integer
563 | return the number of bits needed to store an integer with given max size
567 static int sizeofint(const int size) {
568 unsigned int num = 1;
571 while (size >= num && num_of_bits < 32) {
578 /*___________________________________________________________________________
580 | sizeofints - calculate 'bitsize' of compressed ints
582 | given the number of small unsigned integers and the maximum value
583 | return the number of bits needed to read or write them with the
584 | routines receiveints and sendints. You need this parameter when
585 | calling these routines. Note that for many calls I can use
586 | the variable 'smallidx' which is exactly the number of bits, and
587 | So I don't need to call 'sizeofints for those calls.
590 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
592 unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
596 for (i=0; i < num_of_ints; i++) {
598 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
599 tmp = bytes[bytecnt] * sizes[i] + tmp;
600 bytes[bytecnt] = tmp & 0xff;
604 bytes[bytecnt++] = tmp & 0xff;
607 num_of_bytes = bytecnt;
611 while (bytes[num_of_bytes] >= num) {
615 return num_of_bits + num_of_bytes * 8;
619 /*____________________________________________________________________________
621 | sendints - send a small set of small integers in compressed format
623 | this routine is used internally by xdr3dfcoord, to send a set of
624 | small integers to the buffer.
625 | Multiplication with fixed (specified maximum ) sizes is used to get
626 | to one big, multibyte integer. Allthough the routine could be
627 | modified to handle sizes bigger than 16777216, or more than just
628 | a few integers, this is not done, because the gain in compression
629 | isn't worth the effort. Note that overflowing the multiplication
630 | or the byte buffer (32 bytes) is unchecked and causes bad results.
634 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
635 unsigned int sizes[], unsigned int nums[]) {
638 unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
643 bytes[num_of_bytes++] = tmp & 0xff;
647 for (i = 1; i < num_of_ints; i++) {
648 if (nums[i] >= sizes[i]) {
649 fprintf(stderr,"major breakdown in sendints num %d doesn't "
650 "match size %d\n", nums[i], sizes[i]);
653 /* use one step multiply */
655 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
656 tmp = bytes[bytecnt] * sizes[i] + tmp;
657 bytes[bytecnt] = tmp & 0xff;
661 bytes[bytecnt++] = tmp & 0xff;
664 num_of_bytes = bytecnt;
666 if (num_of_bits >= num_of_bytes * 8) {
667 for (i = 0; i < num_of_bytes; i++) {
668 sendbits(buf, 8, bytes[i]);
670 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
672 for (i = 0; i < num_of_bytes-1; i++) {
673 sendbits(buf, 8, bytes[i]);
675 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
680 /*___________________________________________________________________________
682 | receivebits - decode number from buf using specified number of bits
684 | extract the number of bits from the array buf and construct an integer
685 | from it. Return that value.
689 static int receivebits(int buf[], int num_of_bits) {
692 unsigned int lastbits, lastbyte;
693 unsigned char * cbuf;
694 int mask = (1 << num_of_bits) -1;
696 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
698 lastbits = (unsigned int) buf[1];
699 lastbyte = (unsigned int) buf[2];
702 while (num_of_bits >= 8) {
703 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
704 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
707 if (num_of_bits > 0) {
708 if (lastbits < num_of_bits) {
710 lastbyte = (lastbyte << 8) | cbuf[cnt++];
712 lastbits -= num_of_bits;
713 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
722 /*____________________________________________________________________________
724 | receiveints - decode 'small' integers from the buf array
726 | this routine is the inverse from sendints() and decodes the small integers
727 | written to buf by calculating the remainder and doing divisions with
728 | the given sizes[]. You need to specify the total number of bits to be
729 | used from buf in num_of_bits.
733 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
734 unsigned int sizes[], int nums[]) {
736 int i, j, num_of_bytes, p, num;
738 bytes[1] = bytes[2] = bytes[3] = 0;
740 while (num_of_bits > 8) {
741 bytes[num_of_bytes++] = receivebits(buf, 8);
744 if (num_of_bits > 0) {
745 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
747 for (i = num_of_ints-1; i > 0; i--) {
749 for (j = num_of_bytes-1; j >=0; j--) {
750 num = (num << 8) | bytes[j];
753 num = num - p * sizes[i];
757 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
760 /*____________________________________________________________________________
762 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
764 | this routine reads or writes (depending on how you opened the file with
765 | xdropen() ) a large number of 3d coordinates (stored in *fp).
766 | The number of coordinates triplets to write is given by *size. On
767 | read this number may be zero, in which case it reads as many as were written
768 | or it may specify the number if triplets to read (which should match the
770 | Compression is achieved by first converting all floating numbers to integer
771 | using multiplication by *precision and rounding to the nearest integer.
772 | Then the minimum and maximum value are calculated to determine the range.
773 | The limited range of integers so found, is used to compress the coordinates.
774 | In addition the differences between succesive coordinates is calculated.
775 | If the difference happens to be 'small' then only the difference is saved,
776 | compressing the data even more. The notion of 'small' is changed dynamically
777 | and is enlarged or reduced whenever needed or possible.
778 | Extra compression is achieved in the case of GROMOS and coordinates of
779 | water molecules. GROMOS first writes out the Oxygen position, followed by
780 | the two hydrogens. In order to make the differences smaller (and thereby
781 | compression the data better) the order is changed into first one hydrogen
782 | then the oxygen, followed by the other hydrogen. This is rather special, but
783 | it shouldn't harm in the general case.
787 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
790 static int *ip = NULL;
794 int minint[3], maxint[3], mindiff, *lip, diff;
795 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
797 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
799 int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
801 int tmp, *thiscoord, prevcoord[3];
802 unsigned int tmpcoord[30];
804 int bufsize, xdrid, lsize;
805 unsigned int bitsize;
809 /* find out if xdrs is opened for reading or for writing */
811 while (xdridptr[xdrid] != xdrs) {
813 if (xdrid >= MAXID) {
814 fprintf(stderr, "xdr error. no open xdr stream\n");
818 if (xdrmodes[xdrid] == 'w') {
820 /* xdrs is open for writing */
822 if (xdr_int(xdrs, size) == 0)
825 /* when the number of coordinates is small, don't try to compress; just
826 * write them as floats using xdr_vector
829 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
830 (xdrproc_t)xdr_float));
833 xdr_float(xdrs, precision);
835 ip = (int *)malloc(size3 * sizeof(*ip));
837 fprintf(stderr,"malloc failed\n");
840 bufsize = size3 * 1.2;
841 buf = (int *)malloc(bufsize * sizeof(*buf));
843 fprintf(stderr,"malloc failed\n");
847 } else if (*size > oldsize) {
848 ip = (int *)realloc(ip, size3 * sizeof(*ip));
850 fprintf(stderr,"malloc failed\n");
853 bufsize = size3 * 1.2;
854 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
856 fprintf(stderr,"malloc failed\n");
861 /* buf[0-2] are special and do not contain actual data */
862 buf[0] = buf[1] = buf[2] = 0;
863 minint[0] = minint[1] = minint[2] = INT_MAX;
864 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
869 oldlint1 = oldlint2 = oldlint3 = 0;
870 while(lfp < fp + size3 ) {
871 /* find nearest integer */
873 lf = *lfp * *precision + 0.5;
875 lf = *lfp * *precision - 0.5;
876 if (fabs(lf) > MAXABS) {
877 /* scaling would cause overflow */
881 if (lint1 < minint[0]) minint[0] = lint1;
882 if (lint1 > maxint[0]) maxint[0] = lint1;
886 lf = *lfp * *precision + 0.5;
888 lf = *lfp * *precision - 0.5;
889 if (fabs(lf) > MAXABS) {
890 /* scaling would cause overflow */
894 if (lint2 < minint[1]) minint[1] = lint2;
895 if (lint2 > maxint[1]) maxint[1] = lint2;
899 lf = *lfp * *precision + 0.5;
901 lf = *lfp * *precision - 0.5;
902 if (fabs(lf) > MAXABS) {
903 /* scaling would cause overflow */
907 if (lint3 < minint[2]) minint[2] = lint3;
908 if (lint3 > maxint[2]) maxint[2] = lint3;
911 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
912 if (diff < mindiff && lfp > fp + 3)
918 xdr_int(xdrs, &(minint[0]));
919 xdr_int(xdrs, &(minint[1]));
920 xdr_int(xdrs, &(minint[2]));
922 xdr_int(xdrs, &(maxint[0]));
923 xdr_int(xdrs, &(maxint[1]));
924 xdr_int(xdrs, &(maxint[2]));
926 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
927 (float)maxint[1] - (float)minint[1] >= MAXABS ||
928 (float)maxint[2] - (float)minint[2] >= MAXABS) {
929 /* turning value in unsigned by subtracting minint
930 * would cause overflow
934 sizeint[0] = maxint[0] - minint[0]+1;
935 sizeint[1] = maxint[1] - minint[1]+1;
936 sizeint[2] = maxint[2] - minint[2]+1;
938 /* check if one of the sizes is to big to be multiplied */
939 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
940 bitsizeint[0] = sizeofint(sizeint[0]);
941 bitsizeint[1] = sizeofint(sizeint[1]);
942 bitsizeint[2] = sizeofint(sizeint[2]);
943 bitsize = 0; /* flag the use of large sizes */
945 bitsize = sizeofints(3, sizeint);
948 luip = (unsigned int *) ip;
950 while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
953 xdr_int(xdrs, &smallidx);
954 maxidx = MIN(LASTIDX, smallidx + 8) ;
955 minidx = maxidx - 8; /* often this equal smallidx */
956 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
957 small = magicints[smallidx] / 2;
958 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
959 larger = magicints[maxidx] / 2;
963 thiscoord = (int *)(luip) + i * 3;
964 if (smallidx < maxidx && i >= 1 &&
965 abs(thiscoord[0] - prevcoord[0]) < larger &&
966 abs(thiscoord[1] - prevcoord[1]) < larger &&
967 abs(thiscoord[2] - prevcoord[2]) < larger) {
969 } else if (smallidx > minidx) {
975 if (abs(thiscoord[0] - thiscoord[3]) < small &&
976 abs(thiscoord[1] - thiscoord[4]) < small &&
977 abs(thiscoord[2] - thiscoord[5]) < small) {
978 /* interchange first with second atom for better
979 * compression of water molecules
981 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
983 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
985 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
991 tmpcoord[0] = thiscoord[0] - minint[0];
992 tmpcoord[1] = thiscoord[1] - minint[1];
993 tmpcoord[2] = thiscoord[2] - minint[2];
995 sendbits(buf, bitsizeint[0], tmpcoord[0]);
996 sendbits(buf, bitsizeint[1], tmpcoord[1]);
997 sendbits(buf, bitsizeint[2], tmpcoord[2]);
999 sendints(buf, 3, bitsize, sizeint, tmpcoord);
1001 prevcoord[0] = thiscoord[0];
1002 prevcoord[1] = thiscoord[1];
1003 prevcoord[2] = thiscoord[2];
1004 thiscoord = thiscoord + 3;
1008 if (is_small == 0 && is_smaller == -1)
1010 while (is_small && run < 8*3) {
1011 if (is_smaller == -1 && (
1012 SQR(thiscoord[0] - prevcoord[0]) +
1013 SQR(thiscoord[1] - prevcoord[1]) +
1014 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
1018 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1019 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1020 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1022 prevcoord[0] = thiscoord[0];
1023 prevcoord[1] = thiscoord[1];
1024 prevcoord[2] = thiscoord[2];
1027 thiscoord = thiscoord + 3;
1030 abs(thiscoord[0] - prevcoord[0]) < small &&
1031 abs(thiscoord[1] - prevcoord[1]) < small &&
1032 abs(thiscoord[2] - prevcoord[2]) < small) {
1036 if (run != prevrun || is_smaller != 0) {
1038 sendbits(buf, 1, 1); /* flag the change in run-length */
1039 sendbits(buf, 5, run+is_smaller+1);
1041 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1043 for (k=0; k < run; k+=3) {
1044 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
1046 if (is_smaller != 0) {
1047 smallidx += is_smaller;
1048 if (is_smaller < 0) {
1050 smaller = magicints[smallidx-1] / 2;
1053 small = magicints[smallidx] / 2;
1055 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1058 if (buf[1] != 0) buf[0]++;;
1059 xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1060 return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1063 /* xdrs is open for reading */
1065 if (xdr_int(xdrs, &lsize) == 0)
1067 if (*size != 0 && lsize != *size) {
1068 fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1069 "%d arg vs %d in file", *size, lsize);
1074 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
1075 (xdrproc_t)xdr_float));
1077 xdr_float(xdrs, precision);
1079 ip = (int *)malloc(size3 * sizeof(*ip));
1081 fprintf(stderr,"malloc failed\n");
1084 bufsize = size3 * 1.2;
1085 buf = (int *)malloc(bufsize * sizeof(*buf));
1087 fprintf(stderr,"malloc failed\n");
1091 } else if (*size > oldsize) {
1092 ip = (int *)realloc(ip, size3 * sizeof(*ip));
1094 fprintf(stderr,"malloc failed\n");
1097 bufsize = size3 * 1.2;
1098 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1100 fprintf(stderr,"malloc failed\n");
1105 buf[0] = buf[1] = buf[2] = 0;
1107 xdr_int(xdrs, &(minint[0]));
1108 xdr_int(xdrs, &(minint[1]));
1109 xdr_int(xdrs, &(minint[2]));
1111 xdr_int(xdrs, &(maxint[0]));
1112 xdr_int(xdrs, &(maxint[1]));
1113 xdr_int(xdrs, &(maxint[2]));
1115 sizeint[0] = maxint[0] - minint[0]+1;
1116 sizeint[1] = maxint[1] - minint[1]+1;
1117 sizeint[2] = maxint[2] - minint[2]+1;
1119 /* check if one of the sizes is to big to be multiplied */
1120 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1121 bitsizeint[0] = sizeofint(sizeint[0]);
1122 bitsizeint[1] = sizeofint(sizeint[1]);
1123 bitsizeint[2] = sizeofint(sizeint[2]);
1124 bitsize = 0; /* flag the use of large sizes */
1126 bitsize = sizeofints(3, sizeint);
1129 xdr_int(xdrs, &smallidx);
1130 maxidx = MIN(LASTIDX, smallidx + 8) ;
1131 minidx = maxidx - 8; /* often this equal smallidx */
1132 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1133 small = magicints[smallidx] / 2;
1134 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1135 larger = magicints[maxidx];
1137 /* buf[0] holds the length in bytes */
1139 if (xdr_int(xdrs, &(buf[0])) == 0)
1141 if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1143 buf[0] = buf[1] = buf[2] = 0;
1146 inv_precision = 1.0 / * precision;
1150 while ( i < lsize ) {
1151 thiscoord = (int *)(lip) + i * 3;
1154 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1155 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1156 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1158 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1162 thiscoord[0] += minint[0];
1163 thiscoord[1] += minint[1];
1164 thiscoord[2] += minint[2];
1166 prevcoord[0] = thiscoord[0];
1167 prevcoord[1] = thiscoord[1];
1168 prevcoord[2] = thiscoord[2];
1171 flag = receivebits(buf, 1);
1174 run = receivebits(buf, 5);
1175 is_smaller = run % 3;
1181 for (k = 0; k < run; k+=3) {
1182 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1184 thiscoord[0] += prevcoord[0] - small;
1185 thiscoord[1] += prevcoord[1] - small;
1186 thiscoord[2] += prevcoord[2] - small;
1188 /* interchange first with second atom for better
1189 * compression of water molecules
1191 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1193 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1195 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1197 *lfp++ = prevcoord[0] * inv_precision;
1198 *lfp++ = prevcoord[1] * inv_precision;
1199 *lfp++ = prevcoord[2] * inv_precision;
1201 prevcoord[0] = thiscoord[0];
1202 prevcoord[1] = thiscoord[1];
1203 prevcoord[2] = thiscoord[2];
1205 *lfp++ = thiscoord[0] * inv_precision;
1206 *lfp++ = thiscoord[1] * inv_precision;
1207 *lfp++ = thiscoord[2] * inv_precision;
1210 *lfp++ = thiscoord[0] * inv_precision;
1211 *lfp++ = thiscoord[1] * inv_precision;
1212 *lfp++ = thiscoord[2] * inv_precision;
1214 smallidx += is_smaller;
1215 if (is_smaller < 0) {
1217 if (smallidx > FIRSTIDX) {
1218 smaller = magicints[smallidx - 1] /2;
1222 } else if (is_smaller > 0) {
1224 small = magicints[smallidx] / 2;
1226 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;