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 1058 by mmeineke, Wed Feb 18 21:20:51 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  
21  
# 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( 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  
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 <  headPos = (struct linkedPos*)malloc(
80 >
81 >  nFrames = 0;
82 >  headPos = (struct linkedPos*)malloc(sizeof(struct linkedPos));
83 >  currPos = headPos;
84    while( !feof( inFile ) ){
85  
86 <    currPT
86 >    currPT = (fpos_t *)malloc(sizeof(fpos_t));
87 >    fgetpos(inFile, currPT);
88 >    
89 >    fgets( readBuffer, sizeof( readBuffer ), inFile );
90 >    lineNum++;
91 >    if( feof( inFile ) ){
92 >      fprintf( stderr,
93 >               "File \"%s\" ended unexpectedly at line %d\n",
94 >               inName,
95 >               lineNum );
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;
274 <  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