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 1060 by mmeineke, Thu Feb 19 21:10:06 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 +      

Diff Legend

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