--- trunk/tcProps/readWrite.c 2004/02/14 02:48:52 1053 +++ trunk/tcProps/readWrite.c 2004/02/18 21:20:51 1058 @@ -13,6 +13,9 @@ FILE* inFile; #define BUFFER_SIZE 2000 int isScanned; +int fileOpen; +int nFrames; +double *frameTimes; FILE* inFile; @@ -27,40 +30,57 @@ struct staticPos{ struct staticPos{ fpos_t *myPos; - double timeStamp; double Hmat[3][3]; }; struct staticPos* posArray; +void parseDumpLine( char readBuffer[BUFFER_SIZE], int index, + struct atomCoord* atoms ); +void setU( double q[4], double u[3] ); +void normalizeQ( double q[4] ); +void openFile( char* inName ){ + + inFile = fopen(inName); + if(inFile ==NULL){ + fprintf(stderr, + "Error opening file \"%s\"\n", + inName); + exit(0); + } + + fileOpen = 1; +} + void closeFile( void ){ fclose( inFile ); + + fileOpen = 0; } -int setFrames( char* inFile ){ +int setFrames( void ){ - int nFrames = 0; int i,j,k; struct linkedPos* headPos; struct linkedPos* currPos; fpos_t *currPT; char readBuffer[BUFFER_SIZE]; + char trash[BUFFER_SIZE]; char* foo; int lineNum = 0; - - inFile = fopen(inName); - if(inFile ==NULL){ + if( !fileOpen ){ fprintf(stderr, - "Error opening file \"%s\"\n", - inName); + "Attempt to scan the file without first opening the file.\n" ); exit(0); } - + + nFrames = 0; headPos = (struct linkedPos*)malloc(sizeof(struct linkedPos)); + currPos = headPos; while( !feof( inFile ) ){ currPT = (fpos_t *)malloc(sizeof(fpos_t)); @@ -76,11 +96,512 @@ int setFrames( char* inFile ){ exit(0); } + currPos->next = (struct linkedPos*)malloc(sizeof(struct linkedPos)); + currPos = currPos->next; + currPos->myPos = currPT; + nFrames++; + i = atoi(readBuffer); + + fgets( readBuffer, sizeof( readBuffer ), inFile ); + lineNum++; + if( feof( inFile ) ){ + fprintf( stderr, + "File \"%s\" ended unexpectedly at line %d\n", + inName, + lineNum ); + exit(0); + } + // parse the comment line + + foo = strtok(readBuffer, " ,;\t"); + + if(foo == NULL){ + fprintf(stderr, + "error in reading time from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->timeStamp = atof( foo ); + + // get the Hx vector + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hx[0] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[0][0] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hx[1] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[1][0] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hx[2] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[2][0] = atof( foo ); + + // get the Hy vector + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hy[0] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[0][1] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hy[1] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[1][1] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hy[2] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[2][1] = atof( foo ); + + // get the Hz vector + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hz[0] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[0][2] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hz[1] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[1][2] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf(stderr, + "error in reading Hz[2] from %s at line %d\n", + inName, lineNum ); + exit(0); + } + currPos->Hmat[2][2] = atof( foo ); + + // skip the atoms + + for(j=0; jnext; + while( currPos != NULL ){ + + if(i>=nFrames){ + + fprintf( stderr, + "The number of frame pointers exceeded the number of frames.\n" + " OOPS!\n" ); + exit(0); + } + + frameTimes[i] = currPos->timeStamp; + posArray[i].myPos = currPos->myPos; + for(j=0;j<3;j++){ + for(k=0;k<3;k++){ + posArray[i].Hmat[j][k] = currPos->Hmat[j][k]; + } + } + + i++; + currPos = currPos->next; + } + + // clean up the linked list + + currPos = headPos; + while( currPos != NULL ){ + headPos = currPos; + currPos = headPos->next; + free(headPos); + } + + isScanned = 1; - return nFrames; } + + +void readFrame(int theFrame, struct atomCoord* atoms, double Hmat[3][3] ){ + + // list of 'a priori' constants + + const int nLipAtoms = NL_ATOMS; + const int nBonds = NBONDS; + const int nLipids = NLIPIDS; + const int nSSD = NSSD; + const int nAtoms = nLipAtoms * nLipids + nSSD; + + fpos_t *framePos; + char readBuffer[BUFFER_SIZE]; + + int i, j, k; + + + // quick checks + + if (theFrame >= nFrames){ + + fprintf(stderr, + "frame %d, is out of range\n", + theFrame ); + exit(0); + } + + if( !fileOpen ){ + + fprintf( stderr, + "File is closed, cannot read frame %d.\n" + theFrame ); + exit(0); + } + + // access the frame + + framePos = posArray[theFrame].myPos; + fsetPos( inFile, framePos ); + + for(j=0;j<3;j++){ + for(k=0;k<3;k++){ + Hmat[j][k] = posArray[theFrame].Hmat[j][k]; + } + } + + fgets( readBuffer, sizeof( readBuffer ), inFile ); + if( feof( inFile ) ){ + fprintf( stderr, + "File \"%s\" ended unexpectedly in Frame %d\n", + inName, + theFrame ); + exit(0); + } + + i = atoi( readBuffer ); + if( i != nAtoms ){ + fprintf( stderr, + "Error in frame %d: frame atoms, %d, does not params atoms, %d\n", + theFrame, + i, + nAtoms ); + exit(0); + } + + // read and toss the comment line + + fgets( readBuffer, sizeof( readBuffer ), inFile ); + if( feof( inFile ) ){ + fprintf( stderr, + "File \"%s\" ended unexpectedly in Frame %d\n", + inName, + theFrame ); + exit(0); + } + + // read the atoms + + for( i=0; i < nAtoms; i++){ + + fgets( readBuffer, sizeof( readBuffer ), inFile ); + if( feof( inFile ) ){ + fprintf( stderr, + "File \"%s\" ended unexpectedly in Frame %d at atom %d\n", + inName, + theFrame, + i ); + exit(0); + } + + + parseDumpLine( readBuffer, i, atoms ); + } +} + +void parseDumpLine( char readLine[BUFFER_SIZE], int i, + struct atomCoord* atoms ){ + char* foo; + double q[4]; + + // set the string tokenizer + + foo = strtok(readLine, " ,;\t"); + + // check the atom ID + + switch( atoms[i].type ){ + + case HEAD: + if( strcmp(foo, "HEAD") ){ + fprintf( stderr, + "Atom %s does not match master array type \"HEAD\".\n" + foo ); + exit(0); + } + break; + + case CH: + if( strcmp(foo, "CH") ){ + fprintf( stderr, + "Atom %s does not match master array type \"CH\".\n" + foo ); + exit(0); + } + break; + + case CH2: + if( strcmp(foo, "CH2") ){ + fprintf( stderr, + "Atom %s does not match master array type \"CH2\".\n" + foo ); + exit(0); + } + break; + + case CH3: + if( strcmp(foo, "CH3") ){ + fprintf( stderr, + "Atom %s does not match master array type \"CH3\".\n" + foo ); + exit(0); + } + break; + + case SSD: + if( strcmp(foo, "SSD") ){ + fprintf( stderr, + "Atom %s does not match master array type \"SSD\".\n" + foo ); + exit(0); + } + break; + + default: + fprintf(stderr, + "Atom %d is an unrecognized enumeration\n", + i ); + exit(0); + } + + // get the positions + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading postition x from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + atoms[i].pos[0] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading postition y from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + atoms[i].pos[1] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading postition z from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + atoms[i].pos[2] = atof( foo ); + + + // get the velocities + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading velocity x from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + atoms[i].vel[0] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading velocity y from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + atoms[i].vel[1] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading velocity z from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + atoms[i].vel[2] = atof( foo ); + + + // get the quaternions + + if( (atoms[i].type == HEAD) || (atoms[i].type == SSD) ){ + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + sprintf(painCave.errMsg, + "error in reading quaternion 0 from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + q[0] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading quaternion 1 from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + q[1] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading quaternion 2 from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + q[2] = atof( foo ); + + foo = strtok(NULL, " ,;\t"); + if(foo == NULL){ + fprintf( stderr, + "error in reading quaternion 3 from %s\n" + "natoms = %d, i = %d\n", + inName, nAtoms, i ); + exit(0); + } + q[3] = atof( foo ); + + normalizeQ( q ); + + setU( q, atoms[i].u ); + } +} + + +void setU( double the_q[4], double u[3] ){ + + double q0Sqr, q1Sqr, q2Sqr, q3Sqr; + double A[3][3]; + double rb[3]; + + // initialize the axis + + rb[0] = 0.0; + rb[1] = 0.0; + rb[2] = 1.0; + + // set the rotation matrix + + q0Sqr = the_q[0] * the_q[0]; + q1Sqr = the_q[1] * the_q[1]; + q2Sqr = the_q[2] * the_q[2]; + q3Sqr = the_q[3] * the_q[3]; + + A[0][0] = q0Sqr + q1Sqr - q2Sqr - q3Sqr; + A[0][1] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] ); + A[0][2] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] ); + + A[1][0] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] ); + A[1][1] = q0Sqr - q1Sqr + q2Sqr - q3Sqr; + A[1][2] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] ); + + A[2][0] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] ); + A[2][1] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] ); + A[2][2] = q0Sqr - q1Sqr -q2Sqr +q3Sqr; + + // perform the rotation + + for( i=0; i<3; i++ ){ + u[i] = 0.0; + + for( j=0; j<3; j++ ){ + u[i] += A[j][i] * rb[j]; + } + } +} + +void normalizeQ( double q[4] ){ + + double qSqr, qL; + int i; + + qSqr = 0.0; + for(i=0;i<4;i++) + qSqr += q[i] * q[i]; + + for(i=0;i<4;i++) + q[i] /= qSqr; +} +