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 1058 by mmeineke, Wed Feb 18 21:20:51 2004 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines