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') {
466 xdrfiles[xdrid] = fopen(filename, type1);
467 if (xdrfiles[xdrid] == NULL) {
471 xdrmodes[xdrid] = *type;
472 /* next test isn't usefull in the case of C language
473 * but is used for the Fortran interface
474 * (C users are expected to pass the address of an already allocated
478 xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
479 xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
481 xdridptr[xdrid] = xdrs;
482 xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
487 /*_________________________________________________________________________
489 | xdrclose - close a xdr file
491 | This will flush the xdr buffers, and destroy the xdr stream.
492 | It also closes the associated file descriptor (this is *not*
493 | done by xdr_destroy).
497 int xdrclose(XDR *xdrs) {
501 fprintf(stderr, "xdrclose: passed a NULL pointer\n");
504 for (xdrid = 1; xdrid < MAXID; xdrid++) {
505 if (xdridptr[xdrid] == xdrs) {
508 fclose(xdrfiles[xdrid]);
509 xdridptr[xdrid] = NULL;
513 fprintf(stderr, "xdrclose: no such open xdr file\n");
518 /*____________________________________________________________________________
520 | sendbits - encode num into buf using the specified number of bits
522 | This routines appends the value of num to the bits already present in
523 | the array buf. You need to give it the number of bits to use and you
524 | better make sure that this number of bits is enough to hold the value
525 | Also num must be positive.
529 static void sendbits(int buf[], int num_of_bits, int num) {
531 unsigned int cnt, lastbyte;
533 unsigned char * cbuf;
535 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
536 cnt = (unsigned int) buf[0];
538 lastbyte =(unsigned int) buf[2];
539 while (num_of_bits >= 8) {
540 lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
541 cbuf[cnt++] = lastbyte >> lastbits;
544 if (num_of_bits > 0) {
545 lastbyte = (lastbyte << num_of_bits) | num;
546 lastbits += num_of_bits;
549 cbuf[cnt++] = lastbyte >> lastbits;
556 cbuf[cnt] = lastbyte << (8 - lastbits);
560 /*_________________________________________________________________________
562 | sizeofint - calculate bitsize of an integer
564 | return the number of bits needed to store an integer with given max size
568 static int sizeofint(const int size) {
569 unsigned int num = 1;
572 while (size >= num && num_of_bits < 32) {
579 /*___________________________________________________________________________
581 | sizeofints - calculate 'bitsize' of compressed ints
583 | given the number of small unsigned integers and the maximum value
584 | return the number of bits needed to read or write them with the
585 | routines receiveints and sendints. You need this parameter when
586 | calling these routines. Note that for many calls I can use
587 | the variable 'smallidx' which is exactly the number of bits, and
588 | So I don't need to call 'sizeofints for those calls.
591 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
593 unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
597 for (i=0; i < num_of_ints; i++) {
599 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
600 tmp = bytes[bytecnt] * sizes[i] + tmp;
601 bytes[bytecnt] = tmp & 0xff;
605 bytes[bytecnt++] = tmp & 0xff;
608 num_of_bytes = bytecnt;
612 while (bytes[num_of_bytes] >= num) {
616 return num_of_bits + num_of_bytes * 8;
620 /*____________________________________________________________________________
622 | sendints - send a small set of small integers in compressed format
624 | this routine is used internally by xdr3dfcoord, to send a set of
625 | small integers to the buffer.
626 | Multiplication with fixed (specified maximum ) sizes is used to get
627 | to one big, multibyte integer. Allthough the routine could be
628 | modified to handle sizes bigger than 16777216, or more than just
629 | a few integers, this is not done, because the gain in compression
630 | isn't worth the effort. Note that overflowing the multiplication
631 | or the byte buffer (32 bytes) is unchecked and causes bad results.
635 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
636 unsigned int sizes[], unsigned int nums[]) {
639 unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
644 bytes[num_of_bytes++] = tmp & 0xff;
648 for (i = 1; i < num_of_ints; i++) {
649 if (nums[i] >= sizes[i]) {
650 fprintf(stderr,"major breakdown in sendints num %d doesn't "
651 "match size %d\n", nums[i], sizes[i]);
654 /* use one step multiply */
656 for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
657 tmp = bytes[bytecnt] * sizes[i] + tmp;
658 bytes[bytecnt] = tmp & 0xff;
662 bytes[bytecnt++] = tmp & 0xff;
665 num_of_bytes = bytecnt;
667 if (num_of_bits >= num_of_bytes * 8) {
668 for (i = 0; i < num_of_bytes; i++) {
669 sendbits(buf, 8, bytes[i]);
671 sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
673 for (i = 0; i < num_of_bytes-1; i++) {
674 sendbits(buf, 8, bytes[i]);
676 sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
681 /*___________________________________________________________________________
683 | receivebits - decode number from buf using specified number of bits
685 | extract the number of bits from the array buf and construct an integer
686 | from it. Return that value.
690 static int receivebits(int buf[], int num_of_bits) {
693 unsigned int lastbits, lastbyte;
694 unsigned char * cbuf;
695 int mask = (1 << num_of_bits) -1;
697 cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
699 lastbits = (unsigned int) buf[1];
700 lastbyte = (unsigned int) buf[2];
703 while (num_of_bits >= 8) {
704 lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
705 num |= (lastbyte >> lastbits) << (num_of_bits - 8);
708 if (num_of_bits > 0) {
709 if (lastbits < num_of_bits) {
711 lastbyte = (lastbyte << 8) | cbuf[cnt++];
713 lastbits -= num_of_bits;
714 num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
723 /*____________________________________________________________________________
725 | receiveints - decode 'small' integers from the buf array
727 | this routine is the inverse from sendints() and decodes the small integers
728 | written to buf by calculating the remainder and doing divisions with
729 | the given sizes[]. You need to specify the total number of bits to be
730 | used from buf in num_of_bits.
734 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
735 unsigned int sizes[], int nums[]) {
737 int i, j, num_of_bytes, p, num;
739 bytes[1] = bytes[2] = bytes[3] = 0;
741 while (num_of_bits > 8) {
742 bytes[num_of_bytes++] = receivebits(buf, 8);
745 if (num_of_bits > 0) {
746 bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
748 for (i = num_of_ints-1; i > 0; i--) {
750 for (j = num_of_bytes-1; j >=0; j--) {
751 num = (num << 8) | bytes[j];
754 num = num - p * sizes[i];
758 nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
761 /*____________________________________________________________________________
763 | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
765 | this routine reads or writes (depending on how you opened the file with
766 | xdropen() ) a large number of 3d coordinates (stored in *fp).
767 | The number of coordinates triplets to write is given by *size. On
768 | read this number may be zero, in which case it reads as many as were written
769 | or it may specify the number if triplets to read (which should match the
771 | Compression is achieved by first converting all floating numbers to integer
772 | using multiplication by *precision and rounding to the nearest integer.
773 | Then the minimum and maximum value are calculated to determine the range.
774 | The limited range of integers so found, is used to compress the coordinates.
775 | In addition the differences between succesive coordinates is calculated.
776 | If the difference happens to be 'small' then only the difference is saved,
777 | compressing the data even more. The notion of 'small' is changed dynamically
778 | and is enlarged or reduced whenever needed or possible.
779 | Extra compression is achieved in the case of GROMOS and coordinates of
780 | water molecules. GROMOS first writes out the Oxygen position, followed by
781 | the two hydrogens. In order to make the differences smaller (and thereby
782 | compression the data better) the order is changed into first one hydrogen
783 | then the oxygen, followed by the other hydrogen. This is rather special, but
784 | it shouldn't harm in the general case.
788 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
791 static int *ip = NULL;
795 int minint[3], maxint[3], mindiff, *lip, diff;
796 int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
798 unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
800 int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
802 int tmp, *thiscoord, prevcoord[3];
803 unsigned int tmpcoord[30];
805 int bufsize, xdrid, lsize;
806 unsigned int bitsize;
810 /* find out if xdrs is opened for reading or for writing */
812 while (xdridptr[xdrid] != xdrs) {
814 if (xdrid >= MAXID) {
815 fprintf(stderr, "xdr error. no open xdr stream\n");
819 if (xdrmodes[xdrid] == 'w') {
821 /* xdrs is open for writing */
823 if (xdr_int(xdrs, size) == 0)
826 /* when the number of coordinates is small, don't try to compress; just
827 * write them as floats using xdr_vector
830 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
831 (xdrproc_t)xdr_float));
834 xdr_float(xdrs, precision);
836 ip = (int *)malloc(size3 * sizeof(*ip));
838 fprintf(stderr,"malloc failed\n");
841 bufsize = size3 * 1.2;
842 buf = (int *)malloc(bufsize * sizeof(*buf));
844 fprintf(stderr,"malloc failed\n");
848 } else if (*size > oldsize) {
849 ip = (int *)realloc(ip, size3 * sizeof(*ip));
851 fprintf(stderr,"malloc failed\n");
854 bufsize = size3 * 1.2;
855 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
857 fprintf(stderr,"malloc failed\n");
862 /* buf[0-2] are special and do not contain actual data */
863 buf[0] = buf[1] = buf[2] = 0;
864 minint[0] = minint[1] = minint[2] = INT_MAX;
865 maxint[0] = maxint[1] = maxint[2] = INT_MIN;
870 oldlint1 = oldlint2 = oldlint3 = 0;
871 while(lfp < fp + size3 ) {
872 /* find nearest integer */
874 lf = *lfp * *precision + 0.5;
876 lf = *lfp * *precision - 0.5;
877 if (fabs(lf) > MAXABS) {
878 /* scaling would cause overflow */
882 if (lint1 < minint[0]) minint[0] = lint1;
883 if (lint1 > maxint[0]) maxint[0] = lint1;
887 lf = *lfp * *precision + 0.5;
889 lf = *lfp * *precision - 0.5;
890 if (fabs(lf) > MAXABS) {
891 /* scaling would cause overflow */
895 if (lint2 < minint[1]) minint[1] = lint2;
896 if (lint2 > maxint[1]) maxint[1] = lint2;
900 lf = *lfp * *precision + 0.5;
902 lf = *lfp * *precision - 0.5;
903 if (fabs(lf) > MAXABS) {
904 /* scaling would cause overflow */
908 if (lint3 < minint[2]) minint[2] = lint3;
909 if (lint3 > maxint[2]) maxint[2] = lint3;
912 diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
913 if (diff < mindiff && lfp > fp + 3)
919 xdr_int(xdrs, &(minint[0]));
920 xdr_int(xdrs, &(minint[1]));
921 xdr_int(xdrs, &(minint[2]));
923 xdr_int(xdrs, &(maxint[0]));
924 xdr_int(xdrs, &(maxint[1]));
925 xdr_int(xdrs, &(maxint[2]));
927 if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
928 (float)maxint[1] - (float)minint[1] >= MAXABS ||
929 (float)maxint[2] - (float)minint[2] >= MAXABS) {
930 /* turning value in unsigned by subtracting minint
931 * would cause overflow
935 sizeint[0] = maxint[0] - minint[0]+1;
936 sizeint[1] = maxint[1] - minint[1]+1;
937 sizeint[2] = maxint[2] - minint[2]+1;
939 /* check if one of the sizes is to big to be multiplied */
940 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
941 bitsizeint[0] = sizeofint(sizeint[0]);
942 bitsizeint[1] = sizeofint(sizeint[1]);
943 bitsizeint[2] = sizeofint(sizeint[2]);
944 bitsize = 0; /* flag the use of large sizes */
946 bitsize = sizeofints(3, sizeint);
949 luip = (unsigned int *) ip;
951 while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
954 xdr_int(xdrs, &smallidx);
955 maxidx = MIN(LASTIDX, smallidx + 8) ;
956 minidx = maxidx - 8; /* often this equal smallidx */
957 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
958 small = magicints[smallidx] / 2;
959 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
960 larger = magicints[maxidx] / 2;
964 thiscoord = (int *)(luip) + i * 3;
965 if (smallidx < maxidx && i >= 1 &&
966 abs(thiscoord[0] - prevcoord[0]) < larger &&
967 abs(thiscoord[1] - prevcoord[1]) < larger &&
968 abs(thiscoord[2] - prevcoord[2]) < larger) {
970 } else if (smallidx > minidx) {
976 if (abs(thiscoord[0] - thiscoord[3]) < small &&
977 abs(thiscoord[1] - thiscoord[4]) < small &&
978 abs(thiscoord[2] - thiscoord[5]) < small) {
979 /* interchange first with second atom for better
980 * compression of water molecules
982 tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
984 tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
986 tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
992 tmpcoord[0] = thiscoord[0] - minint[0];
993 tmpcoord[1] = thiscoord[1] - minint[1];
994 tmpcoord[2] = thiscoord[2] - minint[2];
996 sendbits(buf, bitsizeint[0], tmpcoord[0]);
997 sendbits(buf, bitsizeint[1], tmpcoord[1]);
998 sendbits(buf, bitsizeint[2], tmpcoord[2]);
1000 sendints(buf, 3, bitsize, sizeint, tmpcoord);
1002 prevcoord[0] = thiscoord[0];
1003 prevcoord[1] = thiscoord[1];
1004 prevcoord[2] = thiscoord[2];
1005 thiscoord = thiscoord + 3;
1009 if (is_small == 0 && is_smaller == -1)
1011 while (is_small && run < 8*3) {
1012 if (is_smaller == -1 && (
1013 SQR(thiscoord[0] - prevcoord[0]) +
1014 SQR(thiscoord[1] - prevcoord[1]) +
1015 SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
1019 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1020 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1021 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1023 prevcoord[0] = thiscoord[0];
1024 prevcoord[1] = thiscoord[1];
1025 prevcoord[2] = thiscoord[2];
1028 thiscoord = thiscoord + 3;
1031 abs(thiscoord[0] - prevcoord[0]) < small &&
1032 abs(thiscoord[1] - prevcoord[1]) < small &&
1033 abs(thiscoord[2] - prevcoord[2]) < small) {
1037 if (run != prevrun || is_smaller != 0) {
1039 sendbits(buf, 1, 1); /* flag the change in run-length */
1040 sendbits(buf, 5, run+is_smaller+1);
1042 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1044 for (k=0; k < run; k+=3) {
1045 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);
1047 if (is_smaller != 0) {
1048 smallidx += is_smaller;
1049 if (is_smaller < 0) {
1051 smaller = magicints[smallidx-1] / 2;
1054 small = magicints[smallidx] / 2;
1056 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1059 if (buf[1] != 0) buf[0]++;;
1060 xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1061 return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1064 /* xdrs is open for reading */
1066 if (xdr_int(xdrs, &lsize) == 0)
1068 if (*size != 0 && lsize != *size) {
1069 fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1070 "%d arg vs %d in file", *size, lsize);
1075 return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
1076 (xdrproc_t)xdr_float));
1078 xdr_float(xdrs, precision);
1080 ip = (int *)malloc(size3 * sizeof(*ip));
1082 fprintf(stderr,"malloc failed\n");
1085 bufsize = size3 * 1.2;
1086 buf = (int *)malloc(bufsize * sizeof(*buf));
1088 fprintf(stderr,"malloc failed\n");
1092 } else if (*size > oldsize) {
1093 ip = (int *)realloc(ip, size3 * sizeof(*ip));
1095 fprintf(stderr,"malloc failed\n");
1098 bufsize = size3 * 1.2;
1099 buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1101 fprintf(stderr,"malloc failed\n");
1106 buf[0] = buf[1] = buf[2] = 0;
1108 xdr_int(xdrs, &(minint[0]));
1109 xdr_int(xdrs, &(minint[1]));
1110 xdr_int(xdrs, &(minint[2]));
1112 xdr_int(xdrs, &(maxint[0]));
1113 xdr_int(xdrs, &(maxint[1]));
1114 xdr_int(xdrs, &(maxint[2]));
1116 sizeint[0] = maxint[0] - minint[0]+1;
1117 sizeint[1] = maxint[1] - minint[1]+1;
1118 sizeint[2] = maxint[2] - minint[2]+1;
1120 /* check if one of the sizes is to big to be multiplied */
1121 if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1122 bitsizeint[0] = sizeofint(sizeint[0]);
1123 bitsizeint[1] = sizeofint(sizeint[1]);
1124 bitsizeint[2] = sizeofint(sizeint[2]);
1125 bitsize = 0; /* flag the use of large sizes */
1127 bitsize = sizeofints(3, sizeint);
1130 xdr_int(xdrs, &smallidx);
1131 maxidx = MIN(LASTIDX, smallidx + 8) ;
1132 minidx = maxidx - 8; /* often this equal smallidx */
1133 smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1134 small = magicints[smallidx] / 2;
1135 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1136 larger = magicints[maxidx];
1138 /* buf[0] holds the length in bytes */
1140 if (xdr_int(xdrs, &(buf[0])) == 0)
1142 if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1144 buf[0] = buf[1] = buf[2] = 0;
1147 inv_precision = 1.0 / * precision;
1151 while ( i < lsize ) {
1152 thiscoord = (int *)(lip) + i * 3;
1155 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1156 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1157 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1159 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1163 thiscoord[0] += minint[0];
1164 thiscoord[1] += minint[1];
1165 thiscoord[2] += minint[2];
1167 prevcoord[0] = thiscoord[0];
1168 prevcoord[1] = thiscoord[1];
1169 prevcoord[2] = thiscoord[2];
1172 flag = receivebits(buf, 1);
1175 run = receivebits(buf, 5);
1176 is_smaller = run % 3;
1182 for (k = 0; k < run; k+=3) {
1183 receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1185 thiscoord[0] += prevcoord[0] - small;
1186 thiscoord[1] += prevcoord[1] - small;
1187 thiscoord[2] += prevcoord[2] - small;
1189 /* interchange first with second atom for better
1190 * compression of water molecules
1192 tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1194 tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1196 tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1198 *lfp++ = prevcoord[0] * inv_precision;
1199 *lfp++ = prevcoord[1] * inv_precision;
1200 *lfp++ = prevcoord[2] * inv_precision;
1202 prevcoord[0] = thiscoord[0];
1203 prevcoord[1] = thiscoord[1];
1204 prevcoord[2] = thiscoord[2];
1206 *lfp++ = thiscoord[0] * inv_precision;
1207 *lfp++ = thiscoord[1] * inv_precision;
1208 *lfp++ = thiscoord[2] * inv_precision;
1211 *lfp++ = thiscoord[0] * inv_precision;
1212 *lfp++ = thiscoord[1] * inv_precision;
1213 *lfp++ = thiscoord[2] * inv_precision;
1215 smallidx += is_smaller;
1216 if (is_smaller < 0) {
1218 if (smallidx > FIRSTIDX) {
1219 smaller = magicints[smallidx - 1] /2;
1223 } else if (is_smaller > 0) {
1225 small = magicints[smallidx] / 2;
1227 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;