ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/tcProps/readWrite.c
(Generate patch)

Comparing trunk/tcProps/readWrite.c (file contents):
Revision 1052 by mmeineke, Fri Feb 13 22:13:06 2004 UTC vs.
Revision 1085 by mmeineke, Fri Mar 5 03:10:29 2004 UTC

# Line 11 | Line 11
11   #include "readWrite.h"
12  
13  
14 + #define BUFFER_SIZE 2000
15   int isScanned;
16 + int fileOpen;
17 + int nFrames;
18 + double *frameTimes;
19   FILE* inFile;
20 + char* inName;
21  
17 –
22   struct linkedPos{
23    
24    fpos_t *myPos;
# Line 26 | Line 30 | struct staticPos{
30   struct staticPos{
31      
32    fpos_t *myPos;
29 –  double timeStamp;
33    double Hmat[3][3];
34   };
35   struct staticPos* posArray;
36  
37 + void parseDumpLine( char readBuffer[BUFFER_SIZE], int index,
38 +                    struct atomCoord* atoms );
39 + void setU( double q[4], double u[3] );
40 + void normalizeQ( double q[4] );
41  
42 + void openFile( void ){
43 +  
44 +  inFile = fopen(inName, "r");
45 +  if(inFile ==NULL){
46 +    fprintf(stderr,
47 +            "Error opening file \"%s\"\n",
48 +            inName);
49 +    exit(0);
50 +  }
51 +
52 +  fileOpen = 1;
53 + }  
54 +
55   void closeFile( void ){
56  
57    fclose( inFile );
58 +
59 +  fileOpen = 0;
60   }
61  
62 < int setFrames( char* inFile ){
62 > void setFrames( void ){
63  
42 –  int nFrames = 0;
64    int i,j,k;
65    struct linkedPos* headPos;
66    struct linkedPos* currPos;
67    fpos_t *currPT;
68 <
69 <
70 <  inFile = fopen(inName);
71 <  if(inFile ==NULL){
68 >  char readBuffer[BUFFER_SIZE];
69 >  char trash[BUFFER_SIZE];
70 >  char* foo;
71 >  int lineNum = 0;
72 >  
73 >  
74 >  if( !fileOpen ){
75      fprintf(stderr,
76 <            "Error opening file \"%s\"\n",
53 <            inName);
76 >            "Attempt to scan the file without first opening the file.\n" );
77      exit(0);
78    }
79  
80 +
81 +  nFrames = 0;
82 +  headPos = (struct linkedPos*)malloc(sizeof(struct linkedPos));
83 +  currPos = headPos;
84 +
85 +  currPT = (fpos_t *)malloc(sizeof(fpos_t));
86 +  fgetpos(inFile, currPT);
87    
88 <  headPos = (struct linkedPos*)malloc(
88 >  fgets( readBuffer, sizeof( readBuffer ), inFile );
89 >  lineNum++;
90 >  if( feof( inFile ) ){
91 >    fprintf( stderr,
92 >             "File \"%s\" ended unexpectedly at line %d\n",
93 >             inName,
94 >             lineNum );
95 >    exit(0);
96 >  }
97 >
98    while( !feof( inFile ) ){
99  
100 <    currPT
100 >    currPos->next = (struct linkedPos*)malloc(sizeof(struct linkedPos));
101 >    currPos = currPos->next;
102 >    currPos->myPos = currPT;      
103 >    nFrames++;
104 >
105 >    i = atoi(readBuffer);
106 >
107 >    fgets( readBuffer, sizeof( readBuffer ), inFile );
108 >    lineNum++;
109 >    if( feof( inFile ) ){
110 >      fprintf( stderr,
111 >               "File \"%s\" ended unexpectedly at line %d\n",
112 >               inName,
113 >               lineNum );
114 >      exit(0);
115 >    }
116 >    
117 >    
118 >    // parse the comment line
119 >    
120 >    foo = strtok(readBuffer, " ,;\t");
121 >      
122 >    if(foo == NULL){
123 >      fprintf(stderr,
124 >              "error in reading time from %s at line %d\n",
125 >              inName, lineNum );
126 >      exit(0);
127 >    }
128 >    currPos->timeStamp = atof( foo );
129 >      
130 >    // get the Hx vector
131 >      
132 >    foo = strtok(NULL, " ,;\t");
133 >    if(foo == NULL){
134 >      fprintf(stderr,
135 >              "error in reading Hx[0] from %s at line %d\n",
136 >              inName, lineNum );
137 >      exit(0);
138 >    }
139 >    currPos->Hmat[0][0] = atof( foo );
140 >      
141 >    foo = strtok(NULL, " ,;\t");
142 >    if(foo == NULL){
143 >      fprintf(stderr,
144 >              "error in reading Hx[1] from %s at line %d\n",
145 >              inName, lineNum );
146 >      exit(0);
147 >    }
148 >    currPos->Hmat[1][0] = atof( foo );
149 >      
150 >    foo = strtok(NULL, " ,;\t");
151 >    if(foo == NULL){
152 >      fprintf(stderr,
153 >              "error in reading Hx[2] from %s at line %d\n",
154 >              inName, lineNum );
155 >      exit(0);
156 >    }
157 >    currPos->Hmat[2][0] = atof( foo );    
158 >      
159 >    // get the Hy vector
160 >      
161 >    foo = strtok(NULL, " ,;\t");
162 >    if(foo == NULL){
163 >      fprintf(stderr,
164 >              "error in reading Hy[0] from %s at line %d\n",
165 >              inName, lineNum );
166 >      exit(0);
167 >    }
168 >    currPos->Hmat[0][1] = atof( foo );
169 >      
170 >    foo = strtok(NULL, " ,;\t");
171 >    if(foo == NULL){
172 >      fprintf(stderr,
173 >              "error in reading Hy[1] from %s at line %d\n",
174 >              inName, lineNum );
175 >      exit(0);
176 >    }
177 >    currPos->Hmat[1][1] = atof( foo );
178 >      
179 >    foo = strtok(NULL, " ,;\t");
180 >    if(foo == NULL){
181 >      fprintf(stderr,
182 >              "error in reading Hy[2] from %s at line %d\n",
183 >              inName, lineNum );
184 >      exit(0);
185 >    }
186 >    currPos->Hmat[2][1] = atof( foo );    
187 >      
188 >    // get the Hz vector
189 >      
190 >    foo = strtok(NULL, " ,;\t");
191 >    if(foo == NULL){
192 >      fprintf(stderr,
193 >              "error in reading Hz[0] from %s at line %d\n",
194 >              inName, lineNum );
195 >      exit(0);
196 >    }
197 >    currPos->Hmat[0][2] = atof( foo );
198 >      
199 >    foo = strtok(NULL, " ,;\t");
200 >    if(foo == NULL){
201 >      fprintf(stderr,
202 >              "error in reading Hz[1] from %s at line %d\n",
203 >              inName, lineNum );
204 >      exit(0);
205 >    }
206 >    currPos->Hmat[1][2] = atof( foo );
207 >      
208 >    foo = strtok(NULL, " ,;\t");
209 >    if(foo == NULL){
210 >      fprintf(stderr,
211 >              "error in reading Hz[2] from %s at line %d\n",
212 >              inName, lineNum );
213 >      exit(0);
214 >    }
215 >    currPos->Hmat[2][2] = atof( foo );
216 >    
217 >    // skip the atoms
218 >
219 >    for(j=0; j<i; j++){
220 >      
221 >      fgets( readBuffer, sizeof( readBuffer ), inFile );
222 >      lineNum++;    
223 >      if( feof( inFile ) ){
224 >        fprintf( stderr,
225 >                 "File \"%s\" ended unexpectedly at line %d,"
226 >                 " with atom %d\n",
227 >                 inName,
228 >                 lineNum,
229 >                 j );
230 >        exit(0);
231 >      }
232 >    }
233 >
234 >    currPT = (fpos_t *)malloc(sizeof(fpos_t));
235 >    fgetpos(inFile, currPT);
236 >    
237 >    fgets( readBuffer, sizeof( readBuffer ), inFile );
238 >    lineNum++;
239 >  }
240    
241 +  // allocate the static position array
242  
243 +  posArray = (struct staticPos*)calloc(nFrames, sizeof(struct staticPos));
244 +  frameTimes = (double *)calloc(nFrames, sizeof(double) );
245  
246 +  i=0;
247 +  currPos = headPos->next;
248 +  while( currPos != NULL ){
249 +    
250 +    if(i>=nFrames){
251 +      
252 +      fprintf( stderr,
253 +               "The number of frame pointers exceeded the number of frames.\n"
254 +               "  OOPS!\n" );
255 +      exit(0);
256 +    }
257 +    
258 +    frameTimes[i] = currPos->timeStamp;
259 +    posArray[i].myPos = currPos->myPos;
260 +    for(j=0;j<3;j++){
261 +      for(k=0;k<3;k++){
262 +        posArray[i].Hmat[j][k] = currPos->Hmat[j][k];
263 +      }
264 +    }
265 +
266 +    i++;
267 +    currPos = currPos->next;
268 +  }  
269 +
270 +  // clean up the linked list
271 +
272 +  currPos = headPos;
273 +  while( currPos != NULL ){
274 +    headPos = currPos;
275 +    currPos = headPos->next;
276 +    free(headPos);
277 +  }
278 +
279 +
280    isScanned = 1;
281 <  return nFrames;
281 > }
282 >
283 >
284 > void readFrame(int theFrame, struct atomCoord* atoms, double Hmat[3][3] ){
285 >
286 >  // list of 'a priori' constants
287 >
288 >  const int nLipAtoms = NL_ATOMS;
289 >  const int nBonds    = NBONDS;
290 >  const int nLipids   = NLIPIDS;
291 >  const int nSSD      = NSSD;
292 >  const int nAtoms    = nLipAtoms * nLipids + nSSD;
293 >
294 >  fpos_t *framePos;
295 >  char readBuffer[BUFFER_SIZE];
296 >
297 >  int i, j, k;
298 >
299 >
300 >  // quick checks
301 >
302 >  if (theFrame >= nFrames){
303 >    
304 >    fprintf(stderr,
305 >            "frame %d, is out of range\n",
306 >            theFrame );
307 >    exit(0);
308 >  }
309 >
310 >  if( !fileOpen ){
311 >    
312 >    fprintf( stderr,
313 >             "File is closed, cannot read frame %d.\n",
314 >             theFrame );
315 >    exit(0);
316 >  }
317 >
318 >  // access the frame
319 >  
320 >  framePos = posArray[theFrame].myPos;
321 >  fsetpos( inFile, framePos );
322 >
323 >  for(j=0;j<3;j++){
324 >    for(k=0;k<3;k++){
325 >      Hmat[j][k] = posArray[theFrame].Hmat[j][k];
326 >    }
327 >  }
328 >  
329 >  fgets( readBuffer, sizeof( readBuffer ), inFile );
330 >  if( feof( inFile ) ){
331 >    fprintf( stderr,
332 >             "File \"%s\" ended unexpectedly in Frame %d\n",
333 >             inName,
334 >             theFrame );
335 >    exit(0);
336 >  }
337 >  
338 >  i = atoi( readBuffer );
339 >  if( i != nAtoms ){
340 >    fprintf( stderr,
341 >             "Error in frame %d: frame atoms, %d, does not params atoms, %d\n",
342 >             theFrame,
343 >             i,
344 >             nAtoms );
345 >    exit(0);
346 >  }
347 >
348 >  // read and toss the comment line
349 >
350 >  fgets( readBuffer, sizeof( readBuffer ), inFile );
351 >  if( feof( inFile ) ){
352 >    fprintf( stderr,
353 >             "File \"%s\" ended unexpectedly in Frame %d\n",
354 >             inName,
355 >             theFrame );
356 >    exit(0);
357 >  }
358 >
359 >  // read the atoms
360 >
361 >  for( i=0; i < nAtoms; i++){
362 >    
363 >    fgets( readBuffer, sizeof( readBuffer ), inFile );
364 >    if( feof( inFile ) ){
365 >      fprintf( stderr,
366 >               "File \"%s\" ended unexpectedly in Frame %d at atom %d\n",
367 >               inName,
368 >               theFrame,
369 >               i );
370 >      exit(0);
371 >    }
372 >
373 >    
374 >    parseDumpLine( readBuffer, i, atoms );
375 >  }
376   }
377 +
378 + void parseDumpLine( char readLine[BUFFER_SIZE], int i,
379 +                    struct atomCoord* atoms ){
380 +  
381 +  const int nLipAtoms = NL_ATOMS;
382 +  const int nBonds    = NBONDS;
383 +  const int nLipids   = NLIPIDS;
384 +  const int nSSD      = NSSD;
385 +  const int nAtoms    = nLipAtoms * nLipids + nSSD;
386 +
387 +  char* foo;
388 +  double q[4];
389 +
390 +  // set the string tokenizer
391 +  
392 +  foo = strtok(readLine, " ,;\t");
393 +  
394 +  // check the atom ID
395 +
396 +  switch( atoms[i].type ){
397 +    
398 +  case HEAD:
399 +    if( strcmp(foo, "HEAD") ){
400 +      fprintf( stderr,
401 +               "Atom %s does not match master array type \"HEAD\".\n",
402 +               foo );
403 +      exit(0);
404 +    }
405 +    break;
406 +
407 +  case CH:
408 +    if( strcmp(foo, "CH") ){
409 +      fprintf( stderr,
410 +               "Atom %s does not match master array type \"CH\".\n",
411 +               foo );
412 +      exit(0);
413 +    }
414 +    break;
415 +
416 +  case CH2:
417 +    if( strcmp(foo, "CH2") ){
418 +      fprintf( stderr,
419 +               "Atom %s does not match master array type \"CH2\".\n",
420 +               foo );
421 +      exit(0);
422 +    }
423 +    break;
424 +
425 +  case CH3:
426 +    if( strcmp(foo, "CH3") ){
427 +      fprintf( stderr,
428 +               "Atom %s does not match master array type \"CH3\".\n",
429 +               foo );
430 +      exit(0);
431 +    }
432 +    break;
433 +
434 +  case SSD:
435 +    if( strcmp(foo, "SSD") ){
436 +      fprintf( stderr,
437 +               "Atom %s does not match master array type \"SSD\".\n",
438 +               foo );
439 +      exit(0);
440 +    }
441 +    break;
442 +      
443 +  default:
444 +    fprintf(stderr,
445 +            "Atom %d is an unrecognized enumeration\n",
446 +            i );
447 +    exit(0);
448 +  }
449 +  
450 +  // get the positions
451 +  
452 +  foo = strtok(NULL, " ,;\t");
453 +  if(foo == NULL){
454 +    fprintf( stderr,
455 +             "error in reading postition x from %s\n"
456 +             "natoms  = %d, i = %d\n",
457 +             inName, nAtoms, i );
458 +    exit(0);
459 +  }
460 +  atoms[i].pos[0] = atof( foo );
461 +  
462 +  foo = strtok(NULL, " ,;\t");
463 +  if(foo == NULL){
464 +    fprintf( stderr,
465 +             "error in reading postition y from %s\n"
466 +             "natoms  = %d, i = %d\n",
467 +             inName, nAtoms, i );
468 +    exit(0);
469 +  }
470 +  atoms[i].pos[1] = atof( foo );
471 +    
472 +  foo = strtok(NULL, " ,;\t");
473 +  if(foo == NULL){
474 +    fprintf( stderr,
475 +             "error in reading postition z from %s\n"
476 +             "natoms  = %d, i = %d\n",
477 +             inName, nAtoms, i );
478 +    exit(0);
479 +  }
480 +  atoms[i].pos[2] = atof( foo );    
481 +
482 +
483 +  // get the velocities
484 +
485 +  foo = strtok(NULL, " ,;\t");
486 +  if(foo == NULL){
487 +    fprintf( stderr,
488 +             "error in reading velocity x from %s\n"
489 +             "natoms  = %d, i = %d\n",
490 +             inName, nAtoms, i );
491 +    exit(0);
492 +  }
493 +  atoms[i].vel[0] = atof( foo );
494 +    
495 +  foo = strtok(NULL, " ,;\t");
496 +  if(foo == NULL){
497 +    fprintf( stderr,
498 +             "error in reading velocity y from %s\n"
499 +             "natoms  = %d, i = %d\n",
500 +             inName, nAtoms, i );
501 +    exit(0);
502 +  }
503 +  atoms[i].vel[1] = atof( foo );
504 +    
505 +  foo = strtok(NULL, " ,;\t");
506 +  if(foo == NULL){
507 +    fprintf( stderr,
508 +             "error in reading velocity z from %s\n"
509 +             "natoms  = %d, i = %d\n",
510 +             inName, nAtoms, i );
511 +    exit(0);
512 +  }
513 +  atoms[i].vel[2] = atof( foo );
514 +    
515 +    
516 +  // get the quaternions
517 +    
518 +  if( (atoms[i].type == HEAD) || (atoms[i].type == SSD) ){
519 +      
520 +    foo = strtok(NULL, " ,;\t");
521 +    if(foo == NULL){
522 +      fprintf( stderr,
523 +               "error in reading quaternion 0 from %s\n"
524 +               "natoms  = %d, i = %d\n",
525 +               inName, nAtoms, i );
526 +      exit(0);
527 +    }
528 +    q[0] = atof( foo );
529 +      
530 +    foo = strtok(NULL, " ,;\t");
531 +    if(foo == NULL){
532 +      fprintf( stderr,
533 +               "error in reading quaternion 1 from %s\n"
534 +               "natoms  = %d, i = %d\n",
535 +               inName, nAtoms, i );
536 +      exit(0);
537 +    }
538 +    q[1] = atof( foo );
539 +      
540 +    foo = strtok(NULL, " ,;\t");
541 +    if(foo == NULL){
542 +      fprintf( stderr,
543 +               "error in reading quaternion 2 from %s\n"
544 +               "natoms  = %d, i = %d\n",
545 +               inName, nAtoms, i );
546 +      exit(0);
547 +    }
548 +    q[2] = atof( foo );
549 +      
550 +    foo = strtok(NULL, " ,;\t");
551 +    if(foo == NULL){
552 +      fprintf( stderr,
553 +               "error in reading quaternion 3 from %s\n"
554 +               "natoms  = %d, i = %d\n",
555 +               inName, nAtoms, i );
556 +      exit(0);
557 +    }
558 +    q[3] = atof( foo );
559 +  
560 +    normalizeQ( q );
561 +
562 +    setU( q, atoms[i].u );
563 +  }
564 + }
565 +
566 +
567 + void setU( double the_q[4], double u[3] ){
568 +
569 +  double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
570 +  double A[3][3];
571 +  double rb[3];
572 +  int i,j,k;
573 +
574 +  // initialize the axis
575 +
576 +  rb[0] = 0.0;
577 +  rb[1] = 0.0;
578 +  rb[2] = 1.0;
579 +
580 +  // set the rotation matrix
581 +  
582 +  q0Sqr = the_q[0] * the_q[0];
583 +  q1Sqr = the_q[1] * the_q[1];
584 +  q2Sqr = the_q[2] * the_q[2];
585 +  q3Sqr = the_q[3] * the_q[3];
586 +  
587 +  A[0][0] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
588 +  A[0][1] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
589 +  A[0][2] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
590 +  
591 +  A[1][0]  = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
592 +  A[1][1] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
593 +  A[1][2] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
594 +  
595 +  A[2][0] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
596 +  A[2][1] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
597 +  A[2][2] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
598 +
599 +  // perform the rotation
600 +  
601 +  for( i=0; i<3; i++ ){
602 +    u[i] = 0.0;
603 +    
604 +    for( j=0; j<3; j++ ){
605 +      u[i] += A[j][i] * rb[j];
606 +    }
607 +  }
608 + }
609 +
610 + void normalizeQ( double q[4] ){
611 +  
612 +  double qSqr, qL;
613 +  int i;
614 +  
615 +  qSqr = 0.0;
616 +  for(i=0;i<4;i++)
617 +    qSqr += q[i] * q[i];
618 +  
619 +  for(i=0;i<4;i++)
620 +    q[i] /= qSqr;
621 + }
622 +
623 + double scanSmallestZ(int startFrame){
624 +  double smallest;
625 +  int i;
626 +  
627 +  smallest = posArray[startFrame].Hmat[2][2];
628 +  for(i=startFrame;i<nFrames;i++)
629 +    smallest =
630 +      (smallest<posArray[i].Hmat[2][2])? smallest : posArray[i].Hmat[2][2];
631 +
632 +  return smallest;
633 + }
634 +  
635 + double calcAverageXY( double startTime ){
636 +  
637 +  int i, count, startFrame, startFound;
638 +  double avg;
639 +
640 +  startFound = 0;
641 +  startFrame = -1;
642 +  while( !startFound ){
643 +
644 +    startFrame++;
645 +
646 +    if(startFrame >= nFrames){
647 +      
648 +      fprintf( stderr,
649 +               "Start Time, %G, was not found in the dump file.\n",
650 +               startTime );
651 +      exit(0);
652 +    }
653 +    
654 +    if(startTime <= frameTimes[startFrame])
655 +      startFound = 1;
656 +
657 +
658 +  }
659 +
660 +  avg = 0.0;
661 +  count = 0;
662 +  for(i=startFrame;i<nFrames;i++){
663 +    avg += posArray[i].Hmat[0][0] * posArray[i].Hmat[1][1];
664 +    count++;
665 +  }
666 +  avg /= (double)count;
667 +
668 +  return avg;
669 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines