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 1053 by mmeineke, Sat Feb 14 02:48:52 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;
67    fpos_t *currPT;
68    char readBuffer[BUFFER_SIZE];
69 +  char trash[BUFFER_SIZE];
70    char* foo;
71    int lineNum = 0;
72    
73    
74 <
54 <  inFile = fopen(inName);
55 <  if(inFile ==NULL){
74 >  if( !fileOpen ){
75      fprintf(stderr,
76 <            "Error opening file \"%s\"\n",
58 <            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  
100 <    currPT = (fpos_t *)malloc(sizeof(fpos_t));
101 <    fgetpos(inFile, currPT);
102 <    
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 ) ){
# Line 76 | Line 114 | int setFrames( char* inFile ){
114        exit(0);
115      }
116      
79    i = atoi(readBuffer);
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;
85  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