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 1055 by mmeineke, Mon Feb 16 21:50:34 2004 UTC vs.
Revision 1080 by mmeineke, Wed Mar 3 15:16:15 2004 UTC

# Line 13 | Line 13 | int isScanned;
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  
18 –
22   struct linkedPos{
23    
24    fpos_t *myPos;
# Line 27 | Line 30 | struct staticPos{
30   struct staticPos{
31      
32    fpos_t *myPos;
30 –  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  
43 –  int nFrames = 0;
64    int i,j,k;
65    struct linkedPos* headPos;
66    struct linkedPos* currPos;
# Line 51 | Line 71 | int setFrames( char* inFile ){
71    int lineNum = 0;
72    
73    
74 <
55 <  inFile = fopen(inName);
56 <  if(inFile ==NULL){
74 >  if( !fileOpen ){
75      fprintf(stderr,
76 <            "Error opening file \"%s\"\n",
59 <            inName);
76 >            "Attempt to scan the file without first opening the file.\n" );
77      exit(0);
78    }
79  
80 <  
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 +  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  
68 –    currPT = (fpos_t *)malloc(sizeof(fpos_t));
69 –    fgetpos(inFile, currPT);
70 –    
71 –    fgets( readBuffer, sizeof( readBuffer ), inFile );
72 –    lineNum++;
73 –    if( feof( inFile ) ){
74 –      fprintf( stderr,
75 –               "File \"%s\" ended unexpectedly at line %d\n",
76 –               inName,
77 –               lineNum );
78 –      exit(0);
79 –    }
80 –    
100      currPos->next = (struct linkedPos*)malloc(sizeof(struct linkedPos));
101      currPos = currPos->next;
102      currPos->myPos = currPT;      
# Line 211 | Line 230 | int setFrames( char* inFile ){
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 +      

Diff Legend

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