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
154 /* #include <rpc/rpc.h>
155 #include <rpc/xdr.h> */
161 int ftocstr(char *, int, char *, int);
162 int ctofstr(char *, int, char *);
165 static FILE *xdrfiles[MAXID];
166 static XDR *xdridptr[MAXID];
167 static char xdrmodes[MAXID];
168 static unsigned int cnt;
170 typedef void (* FUNCTION(xdrfproc)) (int *, void *, int *);
173 FUNCTION(xdrfbool) ARGS(`xdrid, pb, ret')
177 *ret = xdr_bool(xdridptr[*xdrid], (bool_t *) pb);
182 FUNCTION(xdrfchar) ARGS(`xdrid, cp, ret')
186 *ret = xdr_char(xdridptr[*xdrid], cp);
191 FUNCTION(xdrfdouble) ARGS(`xdrid, dp, ret')
195 *ret = xdr_double(xdridptr[*xdrid], dp);
196 cnt += sizeof(double);
200 FUNCTION(xdrffloat) ARGS(`xdrid, fp, ret')
204 *ret = xdr_float(xdridptr[*xdrid], fp);
205 cnt += sizeof(float);
209 FUNCTION(xdrfint) ARGS(`xdrid, ip, ret')
213 *ret = xdr_int(xdridptr[*xdrid], ip);
218 FUNCTION(xdrflong) ARGS(`xdrid, lp, ret')
222 *ret = xdr_long(xdridptr[*xdrid], lp);
227 FUNCTION(xdrfshort) ARGS(`xdrid, sp, ret')
231 *ret = xdr_short(xdridptr[*xdrid], sp);
236 FUNCTION(xdrfuchar) ARGS(`xdrid, ucp, ret')
240 *ret = xdr_u_char(xdridptr[*xdrid], ucp);
245 FUNCTION(xdrfulong) ARGS(`xdrid, ulp, ret')
249 *ret = xdr_u_long(xdridptr[*xdrid], ulp);
250 cnt += sizeof(unsigned long);
254 FUNCTION(xdrfushort) ARGS(`xdrid, usp, ret')
258 *ret = xdr_u_short(xdridptr[*xdrid], usp);
259 cnt += sizeof(unsigned short);
263 FUNCTION(xdrf3dfcoord) ARGS(`xdrid, fp, size, precision, ret')
269 *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
273 FUNCTION(xdrfstring) ARGS(`xdrid, STRING_ARG(sp), maxsize, ret')
280 tsp = (char*) malloc(((STRING_LEN(sp)) + 1) * sizeof(char));
285 if (ftocstr(tsp, *maxsize+1, STRING_PTR(sp), STRING_LEN(sp))) {
290 *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int) *maxsize);
291 ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
297 FUNCTION(xdrfwrapstring) ARGS(`xdrid, STRING_ARG(sp), ret')
303 maxsize = (STRING_LEN(sp)) + 1;
304 tsp = (char*) malloc(maxsize * sizeof(char));
309 if (ftocstr(tsp, maxsize, STRING_PTR(sp), STRING_LEN(sp))) {
314 *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
315 ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
321 FUNCTION(xdrfopaque) ARGS(`xdrid, cp, ccnt, ret')
326 *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
331 FUNCTION(xdrfsetpos) ARGS(`xdrid, pos, ret')
335 *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
339 FUNCTION(xdrf) ARGS(`xdrid, pos')
342 *pos = xdr_getpos(xdridptr[*xdrid]);
346 FUNCTION(xdrfvector) ARGS(`xdrid, cp, size, elproc, ret')
350 FUNCTION(xdrfproc) elproc;
354 for (lcnt = 0; lcnt < *size; lcnt++) {
355 elproc(xdrid, (cp+cnt) , ret);
361 FUNCTION(xdrfclose) ARGS(`xdrid, ret')
365 *ret = xdrclose(xdridptr[*xdrid]);
370 FUNCTION(xdrfopen) ARGS(`xdrid, STRING_ARG(fp), STRING_ARG(mode), ret')
373 STRING_ARG_DECL(mode);
379 if (ftocstr(fname, sizeof(fname), STRING_PTR(fp), STRING_LEN(fp))) {
382 if (ftocstr(fmode, sizeof(fmode), STRING_PTR(mode),
387 *xdrid = xdropen(NULL, fname, fmode);
394 /*___________________________________________________________________________
396 | what follows are the C routines for opening, closing xdr streams
397 | and the routine to read/write compressed coordinates together
398 | with some routines to assist in this task (those are marked
399 | static and cannot be called from user programs)
401 #define MAXABS INT_MAX-2
404 #define MIN(x,y) ((x) < (y) ? (x):(y))
407 #define MAX(x,y) ((x) > (y) ? (x):(y))
410 #define SQR(x) ((x)*(x))
412 static int magicints[] = {
413 0, 0, 0, 0, 0, 0, 0, 0, 0,
414 8, 10, 12, 16, 20, 25, 32, 40, 50, 64,
415 80, 101, 128, 161, 203, 256, 322, 406, 512, 645,
416 812, 1024, 1290, 1625, 2048, 2580, 3250, 4096, 5060, 6501,
417 8192, 10321, 13003, 16384, 20642, 26007, 32768, 41285, 52015, 65536,
418 82570, 104031, 131072, 165140, 208063, 262144, 330280, 416127, 524287, 660561,
419 832255, 1048576, 1321122, 1664510, 2097152, 2642245, 3329021, 4194304, 5284491, 6658042,
420 8388607, 10568983, 13316085, 16777216 };
423 /* note that magicints[FIRSTIDX-1] == 0 */
424 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
427 /*__________________________________________________________________________
429 | xdropen - open xdr file
431 | This versions differs from xdrstdio_create, because I need to know
432 | the state of the file (read or write) so I can use xdr3dfcoord
433 | in eigther read or write mode, and the file descriptor
434 | so I can close the file (something xdr_destroy doesn't do).
438 int xdropen(XDR *xdrs, const char *filename, const char *type) {
439 static int init_done = 0;
444 if (init_done == 0) {
445 for (xdrid = 1; xdrid < MAXID; xdrid++) {
446 xdridptr[xdrid] = NULL;
451 while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
454 if (xdrid == MAXID) {
457 if (*type == 'w' || *type == 'W') {
461 } else if (*type == 'a' || *type == 'A') {
470 xdrfiles[xdrid] = fopen(filename, type1);
471 if (xdrfiles[xdrid] == NULL) {
475 xdrmodes[xdrid] = *type;
476 /* next test isn't usefull in the case of C language
477 * but is used for the Fortran interface
478 * (C users are expected to pass the address of an already allocated
482 xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
483 xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
485 xdridptr[xdrid] = xdrs;
486 xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
491 /*_________________________________________________________________________
493 | xdrclose - close a xdr file
495 | This will flush the xdr buffers, and destroy the xdr stream.
496 | It also closes the associated file descriptor (this is *not*
497 | done by xdr_destroy).
501 int xdrclose(XDR *xdrs) {
505 fprintf(stderr, "xdrclose: passed a NULL pointer\n");
508 for (xdrid = 1; xdrid < MAXID; xdrid++) {
509 if (xdridptr[xdrid] == xdrs) {
512 fclose(xdrfiles[xdrid]);
513 xdridptr[xdrid] = NULL;
517 fprintf(stderr, "xdrclose: no such open xdr file\n");
522 /*____________________________________________________________________________
524 | sendbits - encode num into buf using the specified number of bits
526 | This routines appends the value of num to the bits already present in
527 | the array buf. You need to give it the number of bits to use and you
528 | better make sure that this number of bits is enough to hold the value
529 | Also num must be positive.
533 static void sendbits(int buf[], int num_of_bits, int num) {
535 unsigned int cnt, lastbyte;
537 unsigned char * cbuf;
539 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
540 cnt = (unsigned int) buf[0];
542 lastbyte =(unsigned int) buf[2];
543 while (num_of_bits >= 8) {
544 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
545 cbuf[cnt++] = lastbyte >> lastbits;
548 if (num_of_bits > 0) {
549 lastbyte = (lastbyte << num_of_bits) | num;
550 lastbits += num_of_bits;
553 cbuf[cnt++] = lastbyte >> lastbits;
560 cbuf[cnt] = lastbyte << (8 - lastbits);
564 /*_________________________________________________________________________
566 | sizeofint - calculate bitsize of an integer
568 | return the number of bits needed to store an integer with given max size
572 static int sizeofint(const int size) {
573 unsigned int num = 1;
576 while (size >= num && num_of_bits < 32) {
583 /*___________________________________________________________________________
585 | sizeofints - calculate 'bitsize' of compressed ints
587 | given the number of small unsigned integers and the maximum value
588 | return the number of bits needed to read or write them with the
589 | routines receiveints and sendints. You need this parameter when
590 | calling these routines. Note that for many calls I can use
591 | the variable 'smallidx' which is exactly the number of bits, and
592 | So I don't need to call 'sizeofints for those calls.
595 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
597 unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
601 for (i=0; i < num_of_ints; i++) {
603 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
604 tmp = bytes[bytecnt] * sizes[i] + tmp;
605 bytes[bytecnt] = tmp & 0xff;
609 bytes[bytecnt++] = tmp & 0xff;
612 num_of_bytes = bytecnt;
616 while (bytes[num_of_bytes] >= num) {
620 return num_of_bits + num_of_bytes * 8;
624 /*____________________________________________________________________________
626 | sendints - send a small set of small integers in compressed format
628 | this routine is used internally by xdr3dfcoord, to send a set of
629 | small integers to the buffer.
630 | Multiplication with fixed (specified maximum ) sizes is used to get
631 | to one big, multibyte integer. Allthough the routine could be
632 | modified to handle sizes bigger than 16777216, or more than just
633 | a few integers, this is not done, because the gain in compression
634 | isn't worth the effort. Note that overflowing the multiplication
635 | or the byte buffer (32 bytes) is unchecked and causes bad results.
639 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
640 unsigned int sizes[], unsigned int nums[]) {
643 unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
648 bytes[num_of_bytes++] = tmp & 0xff;
652 for (i = 1; i < num_of_ints; i++) {
653 if (nums[i] >= sizes[i]) {
654 fprintf(stderr,"major breakdown in sendints num %d doesn't "
655 "match size %d\n", nums[i], sizes[i]);
658 /* use one step multiply */
660 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
661 tmp = bytes[bytecnt] * sizes[i] + tmp;
662 bytes[bytecnt] = tmp & 0xff;
666 bytes[bytecnt++] = tmp & 0xff;
669 num_of_bytes = bytecnt;
671 if (num_of_bits >= num_of_bytes * 8) {
672 for (i = 0; i < num_of_bytes; i++) {
673 sendbits(buf, 8, bytes[i]);
675 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
677 for (i = 0; i < num_of_bytes-1; i++) {
678 sendbits(buf, 8, bytes[i]);
680 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
685 /*___________________________________________________________________________
687 | receivebits - decode number from buf using specified number of bits
689 | extract the number of bits from the array buf and construct an integer
690 | from it. Return that value.
694 static int receivebits(int buf[], int num_of_bits) {
697 unsigned int lastbits, lastbyte;
698 unsigned char * cbuf;
699 int mask = (1 << num_of_bits) -1;
701 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
703 lastbits = (unsigned int) buf[1];
704 lastbyte = (unsigned int) buf[2];
707 while (num_of_bits >= 8) {
708 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
709 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
712 if (num_of_bits > 0) {
713 if (lastbits < num_of_bits) {
715 lastbyte = (lastbyte << 8) | cbuf[cnt++];
717 lastbits -= num_of_bits;
718 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
727 /*____________________________________________________________________________
729 | receiveints - decode 'small' integers from the buf array
731 | this routine is the inverse from sendints() and decodes the small integers
732 | written to buf by calculating the remainder and doing divisions with
733 | the given sizes[]. You need to specify the total number of bits to be
734 | used from buf in num_of_bits.
738 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
739 unsigned int sizes[], int nums[]) {
741 int i, j, num_of_bytes, p, num;
743 bytes[1] = bytes[2] = bytes[3] = 0;
745 while (num_of_bits > 8) {
746 bytes[num_of_bytes++] = receivebits(buf, 8);
749 if (num_of_bits > 0) {
750 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
752 for (i = num_of_ints-1; i > 0; i--) {
754 for (j = num_of_bytes-1; j >=0; j--) {
755 num = (num << 8) | bytes[j];
758 num = num - p * sizes[i];
762 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
765 /*____________________________________________________________________________
767 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
769 | this routine reads or writes (depending on how you opened the file with
770 | xdropen() ) a large number of 3d coordinates (stored in *fp).
771 | The number of coordinates triplets to write is given by *size. On
772 | read this number may be zero, in which case it reads as many as were written
773 | or it may specify the number if triplets to read (which should match the
775 | Compression is achieved by first converting all floating numbers to integer
776 | using multiplication by *precision and rounding to the nearest integer.
777 | Then the minimum and maximum value are calculated to determine the range.
778 | The limited range of integers so found, is used to compress the coordinates.
779 | In addition the differences between succesive coordinates is calculated.
780 | If the difference happens to be 'small' then only the difference is saved,
781 | compressing the data even more. The notion of 'small' is changed dynamically
782 | and is enlarged or reduced whenever needed or possible.
783 | Extra compression is achieved in the case of GROMOS and coordinates of
784 | water molecules. GROMOS first writes out the Oxygen position, followed by
785 | the two hydrogens. In order to make the differences smaller (and thereby
786 | compression the data better) the order is changed into first one hydrogen
787 | then the oxygen, followed by the other hydrogen. This is rather special, but
788 | it shouldn't harm in the general case.
792 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
795 static int *ip = NULL;
799 int minint[3], maxint[3], mindiff, *lip, diff;
800 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
802 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
804 int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
806 int tmp, *thiscoord, prevcoord[3];
807 unsigned int tmpcoord[30];
809 int bufsize, xdrid, lsize;
810 unsigned int bitsize;
814 /* find out if xdrs is opened for reading or for writing */
816 while (xdridptr[xdrid] != xdrs) {
818 if (xdrid >= MAXID) {
819 fprintf(stderr, "xdr error. no open xdr stream\n");
823 if (xdrmodes[xdrid] == 'w') {
825 /* xdrs is open for writing */
827 if (xdr_int(xdrs, size) == 0)
830 /* when the number of coordinates is small, don't try to compress; just
831 * write them as floats using xdr_vector
834 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
835 (xdrproc_t)xdr_float));
838 xdr_float(xdrs, precision);
840 ip = (int *)malloc(size3 * sizeof(*ip));
842 fprintf(stderr,"malloc failed\n");
845 bufsize = size3 * 1.2;
846 buf = (int *)malloc(bufsize * sizeof(*buf));
848 fprintf(stderr,"malloc failed\n");
852 } else if (*size > oldsize) {
853 ip = (int *)realloc(ip, size3 * sizeof(*ip));
855 fprintf(stderr,"malloc failed\n");
858 bufsize = size3 * 1.2;
859 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
861 fprintf(stderr,"malloc failed\n");
866 /* buf[0-2] are special and do not contain actual data */
867 buf[0] = buf[1] = buf[2] = 0;
868 minint[0] = minint[1] = minint[2] = INT_MAX;
869 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
874 oldlint1 = oldlint2 = oldlint3 = 0;
875 while(lfp < fp + size3 ) {
876 /* find nearest integer */
878 lf = *lfp * *precision + 0.5;
880 lf = *lfp * *precision - 0.5;
881 if (fabs(lf) > MAXABS) {
882 /* scaling would cause overflow */
886 if (lint1 < minint[0]) minint[0] = lint1;
887 if (lint1 > maxint[0]) maxint[0] = lint1;
891 lf = *lfp * *precision + 0.5;
893 lf = *lfp * *precision - 0.5;
894 if (fabs(lf) > MAXABS) {
895 /* scaling would cause overflow */
899 if (lint2 < minint[1]) minint[1] = lint2;
900 if (lint2 > maxint[1]) maxint[1] = lint2;
904 lf = *lfp * *precision + 0.5;
906 lf = *lfp * *precision - 0.5;
907 if (fabs(lf) > MAXABS) {
908 /* scaling would cause overflow */
912 if (lint3 < minint[2]) minint[2] = lint3;
913 if (lint3 > maxint[2]) maxint[2] = lint3;
916 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
917 if (diff < mindiff && lfp > fp + 3)
923 xdr_int(xdrs, &(minint[0]));
924 xdr_int(xdrs, &(minint[1]));
925 xdr_int(xdrs, &(minint[2]));
927 xdr_int(xdrs, &(maxint[0]));
928 xdr_int(xdrs, &(maxint[1]));
929 xdr_int(xdrs, &(maxint[2]));
931 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
932 (float)maxint[1] - (float)minint[1] >= MAXABS ||
933 (float)maxint[2] - (float)minint[2] >= MAXABS) {
934 /* turning value in unsigned by subtracting minint
935 * would cause overflow
939 sizeint[0] = maxint[0] - minint[0]+1;
940 sizeint[1] = maxint[1] - minint[1]+1;
941 sizeint[2] = maxint[2] - minint[2]+1;
943 /* check if one of the sizes is to big to be multiplied */
944 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
945 bitsizeint[0] = sizeofint(sizeint[0]);
946 bitsizeint[1] = sizeofint(sizeint[1]);
947 bitsizeint[2] = sizeofint(sizeint[2]);
948 bitsize = 0; /* flag the use of large sizes */
950 bitsize = sizeofints(3, sizeint);
953 luip = (unsigned int *) ip;
955 while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
958 xdr_int(xdrs, &smallidx);
959 maxidx = MIN(LASTIDX, smallidx + 8) ;
960 minidx = maxidx - 8; /* often this equal smallidx */
961 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
962 small = magicints[smallidx] / 2;
963 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
964 larger = magicints[maxidx] / 2;
968 thiscoord = (int *)(luip) + i * 3;
969 if (smallidx < maxidx && i >= 1 &&
970 abs(thiscoord[0] - prevcoord[0]) < larger &&
971 abs(thiscoord[1] - prevcoord[1]) < larger &&
972 abs(thiscoord[2] - prevcoord[2]) < larger) {
974 } else if (smallidx > minidx) {
980 if (abs(thiscoord[0] - thiscoord[3]) < small &&
981 abs(thiscoord[1] - thiscoord[4]) < small &&
982 abs(thiscoord[2] - thiscoord[5]) < small) {
983 /* interchange first with second atom for better
984 * compression of water molecules
986 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
988 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
990 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
996 tmpcoord[0] = thiscoord[0] - minint[0];
997 tmpcoord[1] = thiscoord[1] - minint[1];
998 tmpcoord[2] = thiscoord[2] - minint[2];
1000 sendbits(buf, bitsizeint[0], tmpcoord[0]);
1001 sendbits(buf, bitsizeint[1], tmpcoord[1]);
1002 sendbits(buf, bitsizeint[2], tmpcoord[2]);
1004 sendints(buf, 3, bitsize, sizeint, tmpcoord);
1006 prevcoord[0] = thiscoord[0];
1007 prevcoord[1] = thiscoord[1];
1008 prevcoord[2] = thiscoord[2];
1009 thiscoord = thiscoord + 3;
1013 if (is_small == 0 && is_smaller == -1)
1015 while (is_small && run < 8*3) {
1016 if (is_smaller == -1 && (
1017 SQR(thiscoord[0] - prevcoord[0]) +
1018 SQR(thiscoord[1] - prevcoord[1]) +
1019 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
1023 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1024 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1025 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1027 prevcoord[0] = thiscoord[0];
1028 prevcoord[1] = thiscoord[1];
1029 prevcoord[2] = thiscoord[2];
1032 thiscoord = thiscoord + 3;
1035 abs(thiscoord[0] - prevcoord[0]) < small &&
1036 abs(thiscoord[1] - prevcoord[1]) < small &&
1037 abs(thiscoord[2] - prevcoord[2]) < small) {
1041 if (run != prevrun || is_smaller != 0) {
1043 sendbits(buf, 1, 1); /* flag the change in run-length */
1044 sendbits(buf, 5, run+is_smaller+1);
1046 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1048 for (k=0; k < run; k+=3) {
1049 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
1051 if (is_smaller != 0) {
1052 smallidx += is_smaller;
1053 if (is_smaller < 0) {
1055 smaller = magicints[smallidx-1] / 2;
1058 small = magicints[smallidx] / 2;
1060 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1063 if (buf[1] != 0) buf[0]++;;
1064 xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1065 return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1068 /* xdrs is open for reading */
1070 if (xdr_int(xdrs, &lsize) == 0)
1072 if (*size != 0 && lsize != *size) {
1073 fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1074 "%d arg vs %d in file", *size, lsize);
1079 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
1080 (xdrproc_t)xdr_float));
1082 xdr_float(xdrs, precision);
1084 ip = (int *)malloc(size3 * sizeof(*ip));
1086 fprintf(stderr,"malloc failed\n");
1089 bufsize = size3 * 1.2;
1090 buf = (int *)malloc(bufsize * sizeof(*buf));
1092 fprintf(stderr,"malloc failed\n");
1096 } else if (*size > oldsize) {
1097 ip = (int *)realloc(ip, size3 * sizeof(*ip));
1099 fprintf(stderr,"malloc failed\n");
1102 bufsize = size3 * 1.2;
1103 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1105 fprintf(stderr,"malloc failed\n");
1110 buf[0] = buf[1] = buf[2] = 0;
1112 xdr_int(xdrs, &(minint[0]));
1113 xdr_int(xdrs, &(minint[1]));
1114 xdr_int(xdrs, &(minint[2]));
1116 xdr_int(xdrs, &(maxint[0]));
1117 xdr_int(xdrs, &(maxint[1]));
1118 xdr_int(xdrs, &(maxint[2]));
1120 sizeint[0] = maxint[0] - minint[0]+1;
1121 sizeint[1] = maxint[1] - minint[1]+1;
1122 sizeint[2] = maxint[2] - minint[2]+1;
1124 /* check if one of the sizes is to big to be multiplied */
1125 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1126 bitsizeint[0] = sizeofint(sizeint[0]);
1127 bitsizeint[1] = sizeofint(sizeint[1]);
1128 bitsizeint[2] = sizeofint(sizeint[2]);
1129 bitsize = 0; /* flag the use of large sizes */
1131 bitsize = sizeofints(3, sizeint);
1134 xdr_int(xdrs, &smallidx);
1135 maxidx = MIN(LASTIDX, smallidx + 8) ;
1136 minidx = maxidx - 8; /* often this equal smallidx */
1137 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1138 small = magicints[smallidx] / 2;
1139 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1140 larger = magicints[maxidx];
1142 /* buf[0] holds the length in bytes */
1144 if (xdr_int(xdrs, &(buf[0])) == 0)
1146 if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1148 buf[0] = buf[1] = buf[2] = 0;
1151 inv_precision = 1.0 / * precision;
1155 while ( i < lsize ) {
1156 thiscoord = (int *)(lip) + i * 3;
1159 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1160 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1161 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1163 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1167 thiscoord[0] += minint[0];
1168 thiscoord[1] += minint[1];
1169 thiscoord[2] += minint[2];
1171 prevcoord[0] = thiscoord[0];
1172 prevcoord[1] = thiscoord[1];
1173 prevcoord[2] = thiscoord[2];
1176 flag = receivebits(buf, 1);
1179 run = receivebits(buf, 5);
1180 is_smaller = run % 3;
1186 for (k = 0; k < run; k+=3) {
1187 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1189 thiscoord[0] += prevcoord[0] - small;
1190 thiscoord[1] += prevcoord[1] - small;
1191 thiscoord[2] += prevcoord[2] - small;
1193 /* interchange first with second atom for better
1194 * compression of water molecules
1196 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1198 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1200 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1202 *lfp++ = prevcoord[0] * inv_precision;
1203 *lfp++ = prevcoord[1] * inv_precision;
1204 *lfp++ = prevcoord[2] * inv_precision;
1206 prevcoord[0] = thiscoord[0];
1207 prevcoord[1] = thiscoord[1];
1208 prevcoord[2] = thiscoord[2];
1210 *lfp++ = thiscoord[0] * inv_precision;
1211 *lfp++ = thiscoord[1] * inv_precision;
1212 *lfp++ = thiscoord[2] * inv_precision;
1215 *lfp++ = thiscoord[0] * inv_precision;
1216 *lfp++ = thiscoord[1] * inv_precision;
1217 *lfp++ = thiscoord[2] * inv_precision;
1219 smallidx += is_smaller;
1220 if (is_smaller < 0) {
1222 if (smallidx > FIRSTIDX) {
1223 smaller = magicints[smallidx - 1] /2;
1227 } else if (is_smaller > 0) {
1229 small = magicints[smallidx] / 2;
1231 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;