Corrected WHAM and UNRES reading in SS dyn...
[unres.git] / source / lib / xdrf / libxdrf.m4
1 /*____________________________________________________________________________
2  |
3  | libxdrf - portable fortran interface to xdr. some xdr routines
4  |           are C routines for compressed coordinates
5  |
6  | version 1.1
7  |
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.
11  |
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
16  | on succes.
17  |
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)
24  |
25  | xdrfopen(xdrid, filename, mode, ret)
26  |      character *(*) filename
27  |      character *(*) mode
28  |
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
34  |
35  |      you need to call xdrfclose to flush the output and close
36  |      the file.
37  |      Note that you should not use xdrstdio_create, which comes with the
38  |      standard xdr library
39  |
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
43  |      the xdr libraries.
44  |
45  | xdrfbool(xdrid, bp, ret)
46  |      integer pb
47  |
48  |      This filter produces values of either 1 or 0    
49  |
50  | xdrfchar(xdrid, cp, ret)
51  |      character cp
52  |
53  |      filter that translate between characters and their xdr representation
54  |      Note that the characters in not compressed and occupies 4 bytes.
55  |
56  | xdrfdouble(xdrid, dp, ret)
57  |      double dp
58  |
59  |      read/write a double.
60  |
61  | xdrffloat(xdrid, fp, ret)
62  |      float fp
63  |
64  |      read/write a float.
65  |
66  | xdrfint(xdrid, ip, ret)
67  |      integer ip
68  |
69  |      read/write integer.
70  |
71  | xdrflong(xdrid, lp, ret)
72  |      integer lp
73  |
74  |      this routine has a possible portablility problem due to 64 bits longs.
75  |
76  | xdrfshort(xdrid, sp, ret)
77  |      integer *2 sp
78  |
79  | xdrfstring(xdrid, sp, maxsize, ret)
80  |      character *(*)
81  |      integer maxsize
82  |
83  |      read/write a string, with maximum length given by maxsize
84  |
85  | xdrfwrapstring(xdris, sp, ret)
86  |      character *(*)
87  |
88  |      read/write a string (it is the same as xdrfstring accept that it finds
89  |      the stringlength itself.
90  |
91  | xdrfvector(xdrid, cp, size, xdrfproc, ret)
92  |      character *(*)
93  |      integer size
94  |      external xdrfproc
95  |
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)
102  |      
103  | xdrf3dfcoord(xdrid, fp, size, precision, ret)
104  |      real (*) fp
105  |      real precision
106  |      integer size
107  |
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
119  |
120  | ________________________________________________________________________
121  |
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)    
124  |
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
132  |      the file.
133  |
134  |      Note that you should not use xdrstdio_create, which
135  |      comes with the standard xdr library.
136  |
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).
141  |       
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 
145  |      routines.
146  |
147  |      (c) 1995 Frans van Hoesel, hoesel@chem.rug.nl
148 */      
149
150
151 #include <limits.h>
152 #include <malloc.h>
153 #include <math.h>
154 #include <rpc/rpc.h>
155 #include <rpc/xdr.h> 
156 #include <stdio.h>
157 #include <stdlib.h>
158 #include "xdrf.h"
159
160 int ftocstr(char *, int, char *, int);
161 int ctofstr(char *, int, char *);
162
163 #define MAXID 20
164 static FILE *xdrfiles[MAXID];
165 static XDR *xdridptr[MAXID];
166 static char xdrmodes[MAXID];
167 static unsigned int cnt;
168
169 typedef void (* FUNCTION(xdrfproc)) (int *, void *, int *);
170
171 void
172 FUNCTION(xdrfbool) ARGS(`xdrid, pb, ret')
173 int *xdrid, *ret;
174 int *pb;
175 {
176         *ret = xdr_bool(xdridptr[*xdrid], (bool_t *) pb);
177         cnt += sizeof(int);
178 }
179
180 void
181 FUNCTION(xdrfchar) ARGS(`xdrid, cp, ret')
182 int *xdrid, *ret;
183 char *cp;
184 {
185         *ret = xdr_char(xdridptr[*xdrid], cp);
186         cnt += sizeof(char);
187 }
188
189 void
190 FUNCTION(xdrfdouble) ARGS(`xdrid, dp, ret')
191 int *xdrid, *ret;
192 double *dp;
193 {
194         *ret = xdr_double(xdridptr[*xdrid], dp);
195         cnt += sizeof(double);
196 }
197
198 void
199 FUNCTION(xdrffloat) ARGS(`xdrid, fp, ret')
200 int *xdrid, *ret;
201 float *fp;
202 {
203         *ret = xdr_float(xdridptr[*xdrid], fp);
204         cnt += sizeof(float);
205 }
206
207 void
208 FUNCTION(xdrfint) ARGS(`xdrid, ip, ret')
209 int *xdrid, *ret;
210 int *ip;
211 {
212         *ret = xdr_int(xdridptr[*xdrid], ip);
213         cnt += sizeof(int);
214 }
215
216 void
217 FUNCTION(xdrflong) ARGS(`xdrid, lp, ret')
218 int *xdrid, *ret;
219 long *lp;
220 {
221         *ret = xdr_long(xdridptr[*xdrid], lp);
222         cnt += sizeof(long);
223 }
224
225 void
226 FUNCTION(xdrfshort) ARGS(`xdrid, sp, ret')
227 int *xdrid, *ret;
228 short *sp;
229 {
230         *ret = xdr_short(xdridptr[*xdrid], sp);
231         cnt += sizeof(sp);
232 }
233
234 void
235 FUNCTION(xdrfuchar) ARGS(`xdrid, ucp, ret')
236 int *xdrid, *ret;
237 char *ucp;
238 {
239         *ret = xdr_u_char(xdridptr[*xdrid], ucp);
240         cnt += sizeof(char);
241 }
242
243 void
244 FUNCTION(xdrfulong) ARGS(`xdrid, ulp, ret')
245 int *xdrid, *ret;
246 unsigned long *ulp;
247 {
248         *ret = xdr_u_long(xdridptr[*xdrid], ulp);
249         cnt += sizeof(unsigned long);
250 }
251
252 void
253 FUNCTION(xdrfushort) ARGS(`xdrid, usp, ret')
254 int *xdrid, *ret;
255 unsigned short *usp;
256 {
257         *ret = xdr_u_short(xdridptr[*xdrid], usp);
258         cnt += sizeof(unsigned short);
259 }
260
261 void 
262 FUNCTION(xdrf3dfcoord) ARGS(`xdrid, fp, size, precision, ret')
263 int *xdrid, *ret;
264 float *fp;
265 int *size;
266 float *precision;
267 {
268         *ret = xdr3dfcoord(xdridptr[*xdrid], fp, size, precision);
269 }
270
271 void
272 FUNCTION(xdrfstring) ARGS(`xdrid, STRING_ARG(sp), maxsize, ret')
273 int *xdrid, *ret;
274 STRING_ARG_DECL(sp);
275 int *maxsize;
276 {
277         char *tsp;
278
279         tsp = (char*) malloc(((STRING_LEN(sp)) + 1) * sizeof(char));
280         if (tsp == NULL) {
281             *ret = -1;
282             return;
283         }
284         if (ftocstr(tsp, *maxsize+1, STRING_PTR(sp), STRING_LEN(sp))) {
285             *ret = -1;
286             free(tsp);
287             return;
288         }
289         *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int) *maxsize);
290         ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
291         cnt += *maxsize;
292         free(tsp);
293 }
294
295 void
296 FUNCTION(xdrfwrapstring) ARGS(`xdrid,  STRING_ARG(sp), ret')
297 int *xdrid, *ret;
298 STRING_ARG_DECL(sp);
299 {
300         char *tsp;
301         int maxsize;
302         maxsize = (STRING_LEN(sp)) + 1;
303         tsp = (char*) malloc(maxsize * sizeof(char));
304         if (tsp == NULL) {
305             *ret = -1;
306             return;
307         }
308         if (ftocstr(tsp, maxsize, STRING_PTR(sp), STRING_LEN(sp))) {
309             *ret = -1;
310             free(tsp);
311             return;
312         }
313         *ret = xdr_string(xdridptr[*xdrid], (char **) &tsp, (u_int)maxsize);
314         ctofstr( STRING_PTR(sp), STRING_LEN(sp), tsp);
315         cnt += maxsize;
316         free(tsp);
317 }
318
319 void
320 FUNCTION(xdrfopaque) ARGS(`xdrid, cp, ccnt, ret')
321 int *xdrid, *ret;
322 caddr_t *cp;
323 int *ccnt;
324 {
325         *ret = xdr_opaque(xdridptr[*xdrid], (caddr_t)*cp, (u_int)*ccnt);
326         cnt += *ccnt;
327 }
328
329 void
330 FUNCTION(xdrfsetpos) ARGS(`xdrid, pos, ret')
331 int *xdrid, *ret;
332 int *pos;
333 {
334         *ret = xdr_setpos(xdridptr[*xdrid], (u_int) *pos);
335 }
336
337 void
338 FUNCTION(xdrf) ARGS(`xdrid, pos')
339 int *xdrid, *pos;
340 {
341         *pos = xdr_getpos(xdridptr[*xdrid]);
342 }
343
344 void
345 FUNCTION(xdrfvector) ARGS(`xdrid, cp, size, elproc, ret')
346 int *xdrid, *ret;
347 char *cp;
348 int *size;
349 FUNCTION(xdrfproc) elproc;
350 {
351         int lcnt;
352         cnt = 0;
353         for (lcnt = 0; lcnt < *size; lcnt++) {
354                 elproc(xdrid, (cp+cnt) , ret);
355         }
356 }
357
358
359 void
360 FUNCTION(xdrfclose) ARGS(`xdrid, ret')
361 int *xdrid;
362 int *ret;
363 {
364         *ret = xdrclose(xdridptr[*xdrid]);
365         cnt = 0;
366 }
367
368 void
369 FUNCTION(xdrfopen) ARGS(`xdrid,  STRING_ARG(fp), STRING_ARG(mode), ret')
370 int *xdrid;
371 STRING_ARG_DECL(fp);
372 STRING_ARG_DECL(mode);
373 int *ret;
374 {
375         char fname[512];
376         char fmode[3];
377
378         if (ftocstr(fname, sizeof(fname), STRING_PTR(fp), STRING_LEN(fp))) {
379                 *ret = 0;
380         }
381         if (ftocstr(fmode, sizeof(fmode), STRING_PTR(mode),
382                         STRING_LEN(mode))) {
383                 *ret = 0;
384         }
385
386         *xdrid = xdropen(NULL, fname, fmode);
387         if (*xdrid == 0)
388                 *ret = 0;
389         else 
390                 *ret = 1;       
391 }
392
393 /*___________________________________________________________________________
394  |
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)
399 */
400 #define MAXABS INT_MAX-2
401
402 #ifndef MIN
403 #define MIN(x,y) ((x) < (y) ? (x):(y))
404 #endif
405 #ifndef MAX
406 #define MAX(x,y) ((x) > (y) ? (x):(y))
407 #endif
408 #ifndef SQR
409 #define SQR(x) ((x)*(x))
410 #endif
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 };
420
421 #define FIRSTIDX 9
422 /* note that magicints[FIRSTIDX-1] == 0 */
423 #define LASTIDX (sizeof(magicints) / sizeof(*magicints))
424
425
426 /*__________________________________________________________________________
427  |
428  | xdropen - open xdr file
429  |
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).
434  |
435 */
436
437 int xdropen(XDR *xdrs, const char *filename, const char *type) {
438     static int init_done = 0;
439     enum xdr_op lmode;
440     const char *type1;
441     int xdrid;
442     
443     if (init_done == 0) {
444         for (xdrid = 1; xdrid < MAXID; xdrid++) {
445             xdridptr[xdrid] = NULL;
446         }
447         init_done = 1;
448     }
449     xdrid = 1;
450     while (xdrid < MAXID && xdridptr[xdrid] != NULL) {
451         xdrid++;
452     }
453     if (xdrid == MAXID) {
454         return 0;
455     }
456     if (*type == 'w' || *type == 'W') {
457             type = "w+";
458             type1 = "w+";
459             lmode = XDR_ENCODE;
460     } else if (*type == 'a' || *type == 'A') {
461             type = "w+";
462             type1 = "a+";
463             lmode = XDR_ENCODE;
464     } else {
465             type = "r";
466             type1 = "r";
467             lmode = XDR_DECODE;
468     }
469     xdrfiles[xdrid] = fopen(filename, type1);
470     if (xdrfiles[xdrid] == NULL) {
471         xdrs = NULL;
472         return 0;
473     }
474     xdrmodes[xdrid] = *type;
475     /* next test isn't usefull in the case of C language
476      * but is used for the Fortran interface
477      * (C users are expected to pass the address of an already allocated
478      * XDR staructure)
479      */
480     if (xdrs == NULL) {
481         xdridptr[xdrid] = (XDR *) malloc(sizeof(XDR));
482         xdrstdio_create(xdridptr[xdrid], xdrfiles[xdrid], lmode);
483     } else {
484         xdridptr[xdrid] = xdrs;
485         xdrstdio_create(xdrs, xdrfiles[xdrid], lmode);
486     }
487     return xdrid;
488 }
489
490 /*_________________________________________________________________________
491  |
492  | xdrclose - close a xdr file
493  |
494  | This will flush the xdr buffers, and destroy the xdr stream.
495  | It also closes the associated file descriptor (this is *not*
496  | done by xdr_destroy).
497  |
498 */
499  
500 int xdrclose(XDR *xdrs) {
501     int xdrid;
502     
503     if (xdrs == NULL) {
504         fprintf(stderr, "xdrclose: passed a NULL pointer\n");
505         exit(1);
506     }
507     for (xdrid = 1; xdrid < MAXID; xdrid++) {
508         if (xdridptr[xdrid] == xdrs) {
509             
510             xdr_destroy(xdrs);
511             fclose(xdrfiles[xdrid]);
512             xdridptr[xdrid] = NULL;
513             return 1;
514         }
515     } 
516     fprintf(stderr, "xdrclose: no such open xdr file\n");
517     exit(1);
518     
519 }
520
521 /*____________________________________________________________________________
522  |
523  | sendbits - encode num into buf using the specified number of bits
524  |
525  | This routines appends the value of num to the bits already present in
526  | the array buf. You need to give it the number of bits to use and you
527  | better make sure that this number of bits is enough to hold the value
528  | Also num must be positive.
529  |
530 */
531
532 static void sendbits(int buf[], int num_of_bits, int num) {
533     
534     unsigned int cnt, lastbyte;
535     int lastbits;
536     unsigned char * cbuf;
537     
538     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
539     cnt = (unsigned int) buf[0];
540     lastbits = buf[1];
541     lastbyte =(unsigned int) buf[2];
542     while (num_of_bits >= 8) {
543         lastbyte = (lastbyte << 8) | ((num >> (num_of_bits -8)) /* & 0xff*/);
544         cbuf[cnt++] = lastbyte >> lastbits;
545         num_of_bits -= 8;
546     }
547     if (num_of_bits > 0) {
548         lastbyte = (lastbyte << num_of_bits) | num;
549         lastbits += num_of_bits;
550         if (lastbits >= 8) {
551             lastbits -= 8;
552             cbuf[cnt++] = lastbyte >> lastbits;
553         }
554     }
555     buf[0] = cnt;
556     buf[1] = lastbits;
557     buf[2] = lastbyte;
558     if (lastbits>0) {
559         cbuf[cnt] = lastbyte << (8 - lastbits);
560     }
561 }
562
563 /*_________________________________________________________________________
564  |
565  | sizeofint - calculate bitsize of an integer
566  |
567  | return the number of bits needed to store an integer with given max size
568  |
569 */
570
571 static int sizeofint(const int size) {
572     unsigned int num = 1;
573     int num_of_bits = 0;
574     
575     while (size >= num && num_of_bits < 32) {
576         num_of_bits++;
577         num <<= 1;
578     }
579     return num_of_bits;
580 }
581
582 /*___________________________________________________________________________
583  |
584  | sizeofints - calculate 'bitsize' of compressed ints
585  |
586  | given the number of small unsigned integers and the maximum value
587  | return the number of bits needed to read or write them with the
588  | routines receiveints and sendints. You need this parameter when
589  | calling these routines. Note that for many calls I can use
590  | the variable 'smallidx' which is exactly the number of bits, and
591  | So I don't need to call 'sizeofints for those calls.
592 */
593
594 static int sizeofints( const int num_of_ints, unsigned int sizes[]) {
595     int i, num;
596     unsigned int num_of_bytes, num_of_bits, bytes[32], bytecnt, tmp;
597     num_of_bytes = 1;
598     bytes[0] = 1;
599     num_of_bits = 0;
600     for (i=0; i < num_of_ints; i++) {   
601         tmp = 0;
602         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
603             tmp = bytes[bytecnt] * sizes[i] + tmp;
604             bytes[bytecnt] = tmp & 0xff;
605             tmp >>= 8;
606         }
607         while (tmp != 0) {
608             bytes[bytecnt++] = tmp & 0xff;
609             tmp >>= 8;
610         }
611         num_of_bytes = bytecnt;
612     }
613     num = 1;
614     num_of_bytes--;
615     while (bytes[num_of_bytes] >= num) {
616         num_of_bits++;
617         num *= 2;
618     }
619     return num_of_bits + num_of_bytes * 8;
620
621 }
622     
623 /*____________________________________________________________________________
624  |
625  | sendints - send a small set of small integers in compressed format
626  |
627  | this routine is used internally by xdr3dfcoord, to send a set of
628  | small integers to the buffer. 
629  | Multiplication with fixed (specified maximum ) sizes is used to get
630  | to one big, multibyte integer. Allthough the routine could be
631  | modified to handle sizes bigger than 16777216, or more than just
632  | a few integers, this is not done, because the gain in compression
633  | isn't worth the effort. Note that overflowing the multiplication
634  | or the byte buffer (32 bytes) is unchecked and causes bad results.
635  |
636  */
637  
638 static void sendints(int buf[], const int num_of_ints, const int num_of_bits,
639         unsigned int sizes[], unsigned int nums[]) {
640
641     int i;
642     unsigned int bytes[32], num_of_bytes, bytecnt, tmp;
643
644     tmp = nums[0];
645     num_of_bytes = 0;
646     do {
647         bytes[num_of_bytes++] = tmp & 0xff;
648         tmp >>= 8;
649     } while (tmp != 0);
650
651     for (i = 1; i < num_of_ints; i++) {
652         if (nums[i] >= sizes[i]) {
653             fprintf(stderr,"major breakdown in sendints num %d doesn't "
654                     "match size %d\n", nums[i], sizes[i]);
655             exit(1);
656         }
657         /* use one step multiply */    
658         tmp = nums[i];
659         for (bytecnt = 0; bytecnt < num_of_bytes; bytecnt++) {
660             tmp = bytes[bytecnt] * sizes[i] + tmp;
661             bytes[bytecnt] = tmp & 0xff;
662             tmp >>= 8;
663         }
664         while (tmp != 0) {
665             bytes[bytecnt++] = tmp & 0xff;
666             tmp >>= 8;
667         }
668         num_of_bytes = bytecnt;
669     }
670     if (num_of_bits >= num_of_bytes * 8) {
671         for (i = 0; i < num_of_bytes; i++) {
672             sendbits(buf, 8, bytes[i]);
673         }
674         sendbits(buf, num_of_bits - num_of_bytes * 8, 0);
675     } else {
676         for (i = 0; i < num_of_bytes-1; i++) {
677             sendbits(buf, 8, bytes[i]);
678         }
679         sendbits(buf, num_of_bits- (num_of_bytes -1) * 8, bytes[i]);
680     }
681 }
682
683
684 /*___________________________________________________________________________
685  |
686  | receivebits - decode number from buf using specified number of bits
687  | 
688  | extract the number of bits from the array buf and construct an integer
689  | from it. Return that value.
690  |
691 */
692
693 static int receivebits(int buf[], int num_of_bits) {
694
695     int cnt, num; 
696     unsigned int lastbits, lastbyte;
697     unsigned char * cbuf;
698     int mask = (1 << num_of_bits) -1;
699
700     cbuf = ((unsigned char *)buf) + 3 * sizeof(*buf);
701     cnt = buf[0];
702     lastbits = (unsigned int) buf[1];
703     lastbyte = (unsigned int) buf[2];
704     
705     num = 0;
706     while (num_of_bits >= 8) {
707         lastbyte = ( lastbyte << 8 ) | cbuf[cnt++];
708         num |=  (lastbyte >> lastbits) << (num_of_bits - 8);
709         num_of_bits -=8;
710     }
711     if (num_of_bits > 0) {
712         if (lastbits < num_of_bits) {
713             lastbits += 8;
714             lastbyte = (lastbyte << 8) | cbuf[cnt++];
715         }
716         lastbits -= num_of_bits;
717         num |= (lastbyte >> lastbits) & ((1 << num_of_bits) -1);
718     }
719     num &= mask;
720     buf[0] = cnt;
721     buf[1] = lastbits;
722     buf[2] = lastbyte;
723     return num; 
724 }
725
726 /*____________________________________________________________________________
727  |
728  | receiveints - decode 'small' integers from the buf array
729  |
730  | this routine is the inverse from sendints() and decodes the small integers
731  | written to buf by calculating the remainder and doing divisions with
732  | the given sizes[]. You need to specify the total number of bits to be
733  | used from buf in num_of_bits.
734  |
735 */
736
737 static void receiveints(int buf[], const int num_of_ints, int num_of_bits,
738         unsigned int sizes[], int nums[]) {
739     int bytes[32];
740     int i, j, num_of_bytes, p, num;
741     
742     bytes[1] = bytes[2] = bytes[3] = 0;
743     num_of_bytes = 0;
744     while (num_of_bits > 8) {
745         bytes[num_of_bytes++] = receivebits(buf, 8);
746         num_of_bits -= 8;
747     }
748     if (num_of_bits > 0) {
749         bytes[num_of_bytes++] = receivebits(buf, num_of_bits);
750     }
751     for (i = num_of_ints-1; i > 0; i--) {
752         num = 0;
753         for (j = num_of_bytes-1; j >=0; j--) {
754             num = (num << 8) | bytes[j];
755             p = num / sizes[i];
756             bytes[j] = p;
757             num = num - p * sizes[i];
758         }
759         nums[i] = num;
760     }
761     nums[0] = bytes[0] | (bytes[1] << 8) | (bytes[2] << 16) | (bytes[3] << 24);
762 }
763     
764 /*____________________________________________________________________________
765  |
766  | xdr3dfcoord - read or write compressed 3d coordinates to xdr file.
767  |
768  | this routine reads or writes (depending on how you opened the file with
769  | xdropen() ) a large number of 3d coordinates (stored in *fp).
770  | The number of coordinates triplets to write is given by *size. On
771  | read this number may be zero, in which case it reads as many as were written
772  | or it may specify the number if triplets to read (which should match the
773  | number written).
774  | Compression is achieved by first converting all floating numbers to integer
775  | using multiplication by *precision and rounding to the nearest integer.
776  | Then the minimum and maximum value are calculated to determine the range.
777  | The limited range of integers so found, is used to compress the coordinates.
778  | In addition the differences between succesive coordinates is calculated.
779  | If the difference happens to be 'small' then only the difference is saved,
780  | compressing the data even more. The notion of 'small' is changed dynamically
781  | and is enlarged or reduced whenever needed or possible.
782  | Extra compression is achieved in the case of GROMOS and coordinates of
783  | water molecules. GROMOS first writes out the Oxygen position, followed by
784  | the two hydrogens. In order to make the differences smaller (and thereby
785  | compression the data better) the order is changed into first one hydrogen
786  | then the oxygen, followed by the other hydrogen. This is rather special, but
787  | it shouldn't harm in the general case.
788  |
789  */
790  
791 int xdr3dfcoord(XDR *xdrs, float *fp, int *size, float *precision) {
792     
793
794     static int *ip = NULL;
795     static int oldsize;
796     static int *buf;
797
798     int minint[3], maxint[3], mindiff, *lip, diff;
799     int lint1, lint2, lint3, oldlint1, oldlint2, oldlint3, smallidx;
800     int minidx, maxidx;
801     unsigned sizeint[3], sizesmall[3], bitsizeint[3], size3, *luip;
802     int flag, k;
803     int small, smaller, larger, i, is_small, is_smaller, run, prevrun;
804     float *lfp, lf;
805     int tmp, *thiscoord,  prevcoord[3];
806     unsigned int tmpcoord[30];
807
808     int bufsize, xdrid, lsize;
809     unsigned int bitsize;
810     float inv_precision;
811     int errval = 1;
812
813     /* find out if xdrs is opened for reading or for writing */
814     xdrid = 0;
815     while (xdridptr[xdrid] != xdrs) {
816         xdrid++;
817         if (xdrid >= MAXID) {
818             fprintf(stderr, "xdr error. no open xdr stream\n");
819             exit (1);
820         }
821     }
822     if (xdrmodes[xdrid] == 'w') {
823
824         /* xdrs is open for writing */
825
826         if (xdr_int(xdrs, size) == 0)
827             return 0;
828         size3 = *size * 3;
829         /* when the number of coordinates is small, don't try to compress; just
830          * write them as floats using xdr_vector
831          */
832         if (*size <= 9 ) {
833             return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
834                 (xdrproc_t)xdr_float));
835         }
836         
837         xdr_float(xdrs, precision);
838         if (ip == NULL) {
839             ip = (int *)malloc(size3 * sizeof(*ip));
840             if (ip == NULL) {
841                 fprintf(stderr,"malloc failed\n");
842                 exit(1);
843             }
844             bufsize = size3 * 1.2;
845             buf = (int *)malloc(bufsize * sizeof(*buf));
846             if (buf == NULL) {
847                 fprintf(stderr,"malloc failed\n");
848                 exit(1);
849             }
850             oldsize = *size;
851         } else if (*size > oldsize) {
852             ip = (int *)realloc(ip, size3 * sizeof(*ip));
853             if (ip == NULL) {
854                 fprintf(stderr,"malloc failed\n");
855                 exit(1);
856             }
857             bufsize = size3 * 1.2;
858             buf = (int *)realloc(buf, bufsize * sizeof(*buf));
859             if (buf == NULL) {
860                 fprintf(stderr,"malloc failed\n");
861                 exit(1);
862             }
863             oldsize = *size;
864         }
865         /* buf[0-2] are special and do not contain actual data */
866         buf[0] = buf[1] = buf[2] = 0;
867         minint[0] = minint[1] = minint[2] = INT_MAX;
868         maxint[0] = maxint[1] = maxint[2] = INT_MIN;
869         prevrun = -1;
870         lfp = fp;
871         lip = ip;
872         mindiff = INT_MAX;
873         oldlint1 = oldlint2 = oldlint3 = 0;
874         while(lfp < fp + size3 ) {
875             /* find nearest integer */
876             if (*lfp >= 0.0)
877                 lf = *lfp * *precision + 0.5;
878             else
879                 lf = *lfp * *precision - 0.5;
880             if (fabs(lf) > MAXABS) {
881                 /* scaling would cause overflow */
882                 errval = 0;
883             }
884             lint1 = lf;
885             if (lint1 < minint[0]) minint[0] = lint1;
886             if (lint1 > maxint[0]) maxint[0] = lint1;
887             *lip++ = lint1;
888             lfp++;
889             if (*lfp >= 0.0)
890                 lf = *lfp * *precision + 0.5;
891             else
892                 lf = *lfp * *precision - 0.5;
893             if (fabs(lf) > MAXABS) {
894                 /* scaling would cause overflow */
895                 errval = 0;
896             }
897             lint2 = lf;
898             if (lint2 < minint[1]) minint[1] = lint2;
899             if (lint2 > maxint[1]) maxint[1] = lint2;
900             *lip++ = lint2;
901             lfp++;
902             if (*lfp >= 0.0)
903                 lf = *lfp * *precision + 0.5;
904             else
905                 lf = *lfp * *precision - 0.5;
906             if (fabs(lf) > MAXABS) {
907                 /* scaling would cause overflow */
908                 errval = 0;
909             }
910             lint3 = lf;
911             if (lint3 < minint[2]) minint[2] = lint3;
912             if (lint3 > maxint[2]) maxint[2] = lint3;
913             *lip++ = lint3;
914             lfp++;
915             diff = abs(oldlint1-lint1)+abs(oldlint2-lint2)+abs(oldlint3-lint3);
916             if (diff < mindiff && lfp > fp + 3)
917                 mindiff = diff;
918             oldlint1 = lint1;
919             oldlint2 = lint2;
920             oldlint3 = lint3;
921         }
922         xdr_int(xdrs, &(minint[0]));
923         xdr_int(xdrs, &(minint[1]));
924         xdr_int(xdrs, &(minint[2]));
925         
926         xdr_int(xdrs, &(maxint[0]));
927         xdr_int(xdrs, &(maxint[1]));
928         xdr_int(xdrs, &(maxint[2]));
929         
930         if ((float)maxint[0] - (float)minint[0] >= MAXABS ||
931                 (float)maxint[1] - (float)minint[1] >= MAXABS ||
932                 (float)maxint[2] - (float)minint[2] >= MAXABS) {
933             /* turning value in unsigned by subtracting minint
934              * would cause overflow
935              */
936             errval = 0;
937         }
938         sizeint[0] = maxint[0] - minint[0]+1;
939         sizeint[1] = maxint[1] - minint[1]+1;
940         sizeint[2] = maxint[2] - minint[2]+1;
941         
942         /* check if one of the sizes is to big to be multiplied */
943         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
944             bitsizeint[0] = sizeofint(sizeint[0]);
945             bitsizeint[1] = sizeofint(sizeint[1]);
946             bitsizeint[2] = sizeofint(sizeint[2]);
947             bitsize = 0; /* flag the use of large sizes */
948         } else {
949             bitsize = sizeofints(3, sizeint);
950         }
951         lip = ip;
952         luip = (unsigned int *) ip;
953         smallidx = FIRSTIDX;
954         while (smallidx < LASTIDX && magicints[smallidx] < mindiff) {
955             smallidx++;
956         }
957         xdr_int(xdrs, &smallidx);
958         maxidx = MIN(LASTIDX, smallidx + 8) ;
959         minidx = maxidx - 8; /* often this equal smallidx */
960         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
961         small = magicints[smallidx] / 2;
962         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
963         larger = magicints[maxidx] / 2;
964         i = 0;
965         while (i < *size) {
966             is_small = 0;
967             thiscoord = (int *)(luip) + i * 3;
968             if (smallidx < maxidx && i >= 1 &&
969                     abs(thiscoord[0] - prevcoord[0]) < larger &&
970                     abs(thiscoord[1] - prevcoord[1]) < larger &&
971                     abs(thiscoord[2] - prevcoord[2]) < larger) {
972                 is_smaller = 1;
973             } else if (smallidx > minidx) {
974                 is_smaller = -1;
975             } else {
976                 is_smaller = 0;
977             }
978             if (i + 1 < *size) {
979                 if (abs(thiscoord[0] - thiscoord[3]) < small &&
980                         abs(thiscoord[1] - thiscoord[4]) < small &&
981                         abs(thiscoord[2] - thiscoord[5]) < small) {
982                     /* interchange first with second atom for better
983                      * compression of water molecules
984                      */
985                     tmp = thiscoord[0]; thiscoord[0] = thiscoord[3];
986                         thiscoord[3] = tmp;
987                     tmp = thiscoord[1]; thiscoord[1] = thiscoord[4];
988                         thiscoord[4] = tmp;
989                     tmp = thiscoord[2]; thiscoord[2] = thiscoord[5];
990                         thiscoord[5] = tmp;
991                     is_small = 1;
992                 }
993     
994             }
995             tmpcoord[0] = thiscoord[0] - minint[0];
996             tmpcoord[1] = thiscoord[1] - minint[1];
997             tmpcoord[2] = thiscoord[2] - minint[2];
998             if (bitsize == 0) {
999                 sendbits(buf, bitsizeint[0], tmpcoord[0]);
1000                 sendbits(buf, bitsizeint[1], tmpcoord[1]);
1001                 sendbits(buf, bitsizeint[2], tmpcoord[2]);
1002             } else {
1003                 sendints(buf, 3, bitsize, sizeint, tmpcoord);
1004             }
1005             prevcoord[0] = thiscoord[0];
1006             prevcoord[1] = thiscoord[1];
1007             prevcoord[2] = thiscoord[2];
1008             thiscoord = thiscoord + 3;
1009             i++;
1010             
1011             run = 0;
1012             if (is_small == 0 && is_smaller == -1)
1013                 is_smaller = 0;
1014             while (is_small && run < 8*3) {
1015                 if (is_smaller == -1 && (
1016                         SQR(thiscoord[0] - prevcoord[0]) +
1017                         SQR(thiscoord[1] - prevcoord[1]) +
1018                         SQR(thiscoord[2] - prevcoord[2]) >= smaller * smaller)) {
1019                     is_smaller = 0;
1020                 }
1021
1022                 tmpcoord[run++] = thiscoord[0] - prevcoord[0] + small;
1023                 tmpcoord[run++] = thiscoord[1] - prevcoord[1] + small;
1024                 tmpcoord[run++] = thiscoord[2] - prevcoord[2] + small;
1025                 
1026                 prevcoord[0] = thiscoord[0];
1027                 prevcoord[1] = thiscoord[1];
1028                 prevcoord[2] = thiscoord[2];
1029
1030                 i++;
1031                 thiscoord = thiscoord + 3;
1032                 is_small = 0;
1033                 if (i < *size &&
1034                         abs(thiscoord[0] - prevcoord[0]) < small &&
1035                         abs(thiscoord[1] - prevcoord[1]) < small &&
1036                         abs(thiscoord[2] - prevcoord[2]) < small) {
1037                     is_small = 1;
1038                 }
1039             }
1040             if (run != prevrun || is_smaller != 0) {
1041                 prevrun = run;
1042                 sendbits(buf, 1, 1); /* flag the change in run-length */
1043                 sendbits(buf, 5, run+is_smaller+1);
1044             } else {
1045                 sendbits(buf, 1, 0); /* flag the fact that runlength did not change */
1046             }
1047             for (k=0; k < run; k+=3) {
1048                 sendints(buf, 3, smallidx, sizesmall, &tmpcoord[k]);    
1049             }
1050             if (is_smaller != 0) {
1051                 smallidx += is_smaller;
1052                 if (is_smaller < 0) {
1053                     small = smaller;
1054                     smaller = magicints[smallidx-1] / 2;
1055                 } else {
1056                     smaller = small;
1057                     small = magicints[smallidx] / 2;
1058                 }
1059                 sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx];
1060             }
1061         }
1062         if (buf[1] != 0) buf[0]++;;
1063         xdr_int(xdrs, &(buf[0])); /* buf[0] holds the length in bytes */
1064         return errval * (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]));
1065     } else {
1066         
1067         /* xdrs is open for reading */
1068         
1069         if (xdr_int(xdrs, &lsize) == 0) 
1070             return 0;
1071         if (*size != 0 && lsize != *size) {
1072             fprintf(stderr, "wrong number of coordinates in xdr3dfcoor; "
1073                     "%d arg vs %d in file", *size, lsize);
1074         }
1075         *size = lsize;
1076         size3 = *size * 3;
1077         if (*size <= 9) {
1078             return (xdr_vector(xdrs, (char *) fp, size3, sizeof(*fp),
1079                 (xdrproc_t)xdr_float));
1080         }
1081         xdr_float(xdrs, precision);
1082         if (ip == NULL) {
1083             ip = (int *)malloc(size3 * sizeof(*ip));
1084             if (ip == NULL) {
1085                 fprintf(stderr,"malloc failed\n");
1086                 exit(1);
1087             }
1088             bufsize = size3 * 1.2;
1089             buf = (int *)malloc(bufsize * sizeof(*buf));
1090             if (buf == NULL) {
1091                 fprintf(stderr,"malloc failed\n");
1092                 exit(1);
1093             }
1094             oldsize = *size;
1095         } else if (*size > oldsize) {
1096             ip = (int *)realloc(ip, size3 * sizeof(*ip));
1097             if (ip == NULL) {
1098                 fprintf(stderr,"malloc failed\n");
1099                 exit(1);
1100             }
1101             bufsize = size3 * 1.2;
1102             buf = (int *)realloc(buf, bufsize * sizeof(*buf));
1103             if (buf == NULL) {
1104                 fprintf(stderr,"malloc failed\n");
1105                 exit(1);
1106             }
1107             oldsize = *size;
1108         }
1109         buf[0] = buf[1] = buf[2] = 0;
1110         
1111         xdr_int(xdrs, &(minint[0]));
1112         xdr_int(xdrs, &(minint[1]));
1113         xdr_int(xdrs, &(minint[2]));
1114
1115         xdr_int(xdrs, &(maxint[0]));
1116         xdr_int(xdrs, &(maxint[1]));
1117         xdr_int(xdrs, &(maxint[2]));
1118                 
1119         sizeint[0] = maxint[0] - minint[0]+1;
1120         sizeint[1] = maxint[1] - minint[1]+1;
1121         sizeint[2] = maxint[2] - minint[2]+1;
1122         
1123         /* check if one of the sizes is to big to be multiplied */
1124         if ((sizeint[0] | sizeint[1] | sizeint[2] ) > 0xffffff) {
1125             bitsizeint[0] = sizeofint(sizeint[0]);
1126             bitsizeint[1] = sizeofint(sizeint[1]);
1127             bitsizeint[2] = sizeofint(sizeint[2]);
1128             bitsize = 0; /* flag the use of large sizes */
1129         } else {
1130             bitsize = sizeofints(3, sizeint);
1131         }
1132         
1133         xdr_int(xdrs, &smallidx);
1134         maxidx = MIN(LASTIDX, smallidx + 8) ;
1135         minidx = maxidx - 8; /* often this equal smallidx */
1136         smaller = magicints[MAX(FIRSTIDX, smallidx-1)] / 2;
1137         small = magicints[smallidx] / 2;
1138         sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1139         larger = magicints[maxidx];
1140
1141         /* buf[0] holds the length in bytes */
1142
1143         if (xdr_int(xdrs, &(buf[0])) == 0)
1144             return 0;
1145         if (xdr_opaque(xdrs, (caddr_t)&(buf[3]), (u_int)buf[0]) == 0)
1146             return 0;
1147         buf[0] = buf[1] = buf[2] = 0;
1148         
1149         lfp = fp;
1150         inv_precision = 1.0 / * precision;
1151         run = 0;
1152         i = 0;
1153         lip = ip;
1154         while ( i < lsize ) {
1155             thiscoord = (int *)(lip) + i * 3;
1156
1157             if (bitsize == 0) {
1158                 thiscoord[0] = receivebits(buf, bitsizeint[0]);
1159                 thiscoord[1] = receivebits(buf, bitsizeint[1]);
1160                 thiscoord[2] = receivebits(buf, bitsizeint[2]);
1161             } else {
1162                 receiveints(buf, 3, bitsize, sizeint, thiscoord);
1163             }
1164             
1165             i++;
1166             thiscoord[0] += minint[0];
1167             thiscoord[1] += minint[1];
1168             thiscoord[2] += minint[2];
1169             
1170             prevcoord[0] = thiscoord[0];
1171             prevcoord[1] = thiscoord[1];
1172             prevcoord[2] = thiscoord[2];
1173             
1174            
1175             flag = receivebits(buf, 1);
1176             is_smaller = 0;
1177             if (flag == 1) {
1178                 run = receivebits(buf, 5);
1179                 is_smaller = run % 3;
1180                 run -= is_smaller;
1181                 is_smaller--;
1182             }
1183             if (run > 0) {
1184                 thiscoord += 3;
1185                 for (k = 0; k < run; k+=3) {
1186                     receiveints(buf, 3, smallidx, sizesmall, thiscoord);
1187                     i++;
1188                     thiscoord[0] += prevcoord[0] - small;
1189                     thiscoord[1] += prevcoord[1] - small;
1190                     thiscoord[2] += prevcoord[2] - small;
1191                     if (k == 0) {
1192                         /* interchange first with second atom for better
1193                          * compression of water molecules
1194                          */
1195                         tmp = thiscoord[0]; thiscoord[0] = prevcoord[0];
1196                                 prevcoord[0] = tmp;
1197                         tmp = thiscoord[1]; thiscoord[1] = prevcoord[1];
1198                                 prevcoord[1] = tmp;
1199                         tmp = thiscoord[2]; thiscoord[2] = prevcoord[2];
1200                                 prevcoord[2] = tmp;
1201                         *lfp++ = prevcoord[0] * inv_precision;
1202                         *lfp++ = prevcoord[1] * inv_precision;
1203                         *lfp++ = prevcoord[2] * inv_precision;
1204                     } else {
1205                         prevcoord[0] = thiscoord[0];
1206                         prevcoord[1] = thiscoord[1];
1207                         prevcoord[2] = thiscoord[2];
1208                     }
1209                     *lfp++ = thiscoord[0] * inv_precision;
1210                     *lfp++ = thiscoord[1] * inv_precision;
1211                     *lfp++ = thiscoord[2] * inv_precision;
1212                 }
1213             } else {
1214                 *lfp++ = thiscoord[0] * inv_precision;
1215                 *lfp++ = thiscoord[1] * inv_precision;
1216                 *lfp++ = thiscoord[2] * inv_precision;          
1217             }
1218             smallidx += is_smaller;
1219             if (is_smaller < 0) {
1220                 small = smaller;
1221                 if (smallidx > FIRSTIDX) {
1222                     smaller = magicints[smallidx - 1] /2;
1223                 } else {
1224                     smaller = 0;
1225                 }
1226             } else if (is_smaller > 0) {
1227                 smaller = small;
1228                 small = magicints[smallidx] / 2;
1229             }
1230             sizesmall[0] = sizesmall[1] = sizesmall[2] = magicints[smallidx] ;
1231         }
1232     }
1233     return 1;
1234 }
1235
1236
1237