|
|
/*ScianFiles.c
File reading routines for scian
Eric Pepke
August 17, 1990
12/4/91 Fixed bug with colored polygons
12/4/91 Made default file format work on command line
*/
#include "Scian.h"
#include "ScianTypes.h"
#include "ScianArrays.h"
#include "ScianIcons.h"
#include "ScianWindows.h"
#include "ScianObjWindows.h"
#include "ScianVisWindows.h"
#include "ScianVisObjects.h"
#include "ScianControls.h"
#include "ScianColors.h"
#include "ScianDialogs.h"
#include "ScianFiles.h"
#include "ScianFileSystem.h"
#include "ScianLists.h"
#include "ScianPictures.h"
#include "ScianErrors.h"
#include "ScianTimers.h"
#include "ScianDatasets.h"
#include "ScianFilters.h"
#include "ScianTextBoxes.h"
#include "ScianTitleBoxes.h"
#include "ScianButtons.h"
#include "ScianSliders.h"
#include "ScianScripts.h"
#include "ScianIDs.h"
#include "ScianStyle.h"
#include "ScianMethods.h"
#include "ScianObjFunctions.h"
#include "ScianHwuFiles.h"
#include "ScianSciences.h"
#include "ScianVisObjects.h"
#include "ScianVisSticks.h"
#include "ScianWhisselFiles.h"
#include "ScianAxes.h"
#include "ScianTemplates.h"
#include "ScianTemplateHelper.h"
#include "ScianDatabase.h"
#include "ScianSymbols.h"
#include "ScianSnap.h"
#ifdef HDF31
#define DFSDgetrange DFSDgetmaxmin
#endif
char *curFileName = 0; /*Current file name*/
#define SHMDIDDLE /*Diddle using shared memory technique*/
#if 1
#define CARBONALPHA
#endif
ObjPtr fileClass; /*Class of files*/
ObjPtr fileReaderClass;
#ifdef SHMDIDDLE
#include
#include
real *sharedSegment = 0; /*Attached segment*/
#endif
#ifdef PROTO
ObjPtr MakeDatasetName(char *);
#else
ObjPtr MakeDatasetName();
#endif
void SkipBlanks(file)
FILE *file;
/*Skips blanks in the file*/
{
int c;
while ((c = getc(file)) == ' ' || c == '\t');
ungetc(c, file);
}
void SkipBlanksAndCommas(file)
FILE *file;
/*Skips blanks in the file*/
{
int c;
while ((c = getc(file)) == ' ' || c == '\t' || c == ',');
ungetc(c, file);
}
void SkipNonBlanks(file)
FILE *file;
/*Skips non-blanks in the file*/
{
int c;
while ((c = getc(file)) != ' ' && c >= 0 && c != '\n');
ungetc(c, file);
}
void ReadLn(file)
FILE *file;
/*Reads to next line in the file*/
{
int c;
while ((c = getc(file)) >= 0 && c != '\n');
}
#ifdef PROTO
void FileFormatError(char *routine, char *what)
#else
void FileFormatError(routine, what)
char *routine, *what;
#endif
/*Reports a file format error*/
{
ReportError(routine, what);
}
char *ShortNameOf(n)
char *n;
/*Returns the short file name of n*/
{
char *t;
/*Make name last file name in argument*/
for (t = n; *t; ++t)
{
if (*t == '/') n = t + 1;
}
return n;
}
ObjPtr MakeDatasetName(s)
char *s;
/*Makes and returns a dataset name for file name s*/
{
Bool hasExtension = false;
Bool hasNumber = false;
s = ShortNameOf(s);
strcpy(tempStr, s);
s = tempStr;
/*Trim off anything non-numeric after a period*/
while (*s)
{
if (*s == '.')
{
hasNumber = true;
}
else if (hasNumber)
{
if (*s <= '0' || *s >= '9')
{
hasNumber = false;
hasExtension = true;
}
}
++s;
}
if (hasExtension)
{
while (*s != '.') --s;
*s = 0;
}
return NewString(tempStr);
}
#define ADDCOLORS 100
#ifdef PROTO
static VertexPtr ReadNFFVertices(FILE *inFile, int *nVertices)
#else
static VertexPtr ReadNFFVertices(inFile, nVertices)
FILE *inFile;
int *nVertices;
#endif
/*Reads a bunch of NFF vertices from inFile. Puts the number of vertices in
nVertices.*/
{
VertexPtr vertices = 0;
int nAllocVertices = 0;
char lookahead;
float x, y, z, nx, ny, nz;
int nRead;
int k;
real vec1[3], vec2[3];
real normalLength;
Bool anyNonZero = false;
*nVertices = 0;
for (;;)
{
SkipBlanks(inFile);
lookahead = getc(inFile);
if (lookahead == '#')
{
/*It's a comment. Skip over it*/
ReadLn(inFile);
}
else if (lookahead == '.' || lookahead == '+' || lookahead == '-' || isdigit(lookahead))
{
/*It's a possible vertex line*/
ungetc(lookahead, inFile);
fgets(tempStr, TEMPSTRSIZE, inFile);
nRead = sscanf(tempStr, " %g %g %g %g %g %g",
&x, &y, &z, &nx, &ny, &nz);
if (nRead != 3 && nRead != 6)
{
FileFormatError("ReadNFFVertices", "Badly formed vertex");
continue;
}
if (*nVertices >= nAllocVertices)
{
/*Must reallocate*/
if (nAllocVertices)
{
nAllocVertices += 25;
vertices = Realloc(vertices, nAllocVertices * sizeof(Vertex));
}
else
{
nAllocVertices = 25;
vertices = Alloc(nAllocVertices * sizeof(Vertex));
}
}
vertices[*nVertices] . position[0] = x;
vertices[*nVertices] . position[1] = y;
vertices[*nVertices] . position[2] = z;
if (nRead == 6)
{
vertices[*nVertices] . normal[0] = nx;
vertices[*nVertices] . normal[1] = ny;
vertices[*nVertices] . normal[2] = nz;
}
else
{
vertices[*nVertices] . normal[0] = 0.0;
vertices[*nVertices] . normal[1] = 0.0;
vertices[*nVertices] . normal[2] = 0.0;
}
++*nVertices;
}
else
{
/*End of the vertices*/
ungetc(lookahead, inFile);
break;
}
}
/*Go through and remake normals, if it's necessary*/
for (k = 1; k < *nVertices - 1; ++k)
{
if (vertices[k] . normal[0] == 0.0 &&
vertices[k] . normal[1] == 0.0 &&
vertices[k] . normal[2] == 0.0)
{
vec1[0] = vertices[k + 1] . position[0] - vertices[k] . position[0];
vec1[1] = vertices[k + 1] . position[1] - vertices[k] . position[1];
vec1[2] = vertices[k + 1] . position[2] - vertices[k] . position[2];
vec2[0] = vertices[k - 1] . position[0] - vertices[k] . position[0];
vec2[1] = vertices[k - 1] . position[1] - vertices[k] . position[1];
vec2[2] = vertices[k - 1] . position[2] - vertices[k] . position[2];
CROSS(vec1, vec2, vertices[k] . normal);
normalLength = sqrt(SQUARE(vertices[k] . normal[0]) +
SQUARE(vertices[k] . normal[1]) +
SQUARE(vertices[k] . normal[2]));
if (normalLength > 0.0)
{
vertices[k] . normal[0] /= normalLength;
vertices[k] . normal[1] /= normalLength;
vertices[k] . normal[2] /= normalLength;
anyNonZero = true;
}
}
else
{
anyNonZero = true;
}
}
/*If none nonzero, set the first one arbitrarily*/
if (!anyNonZero)
{
vertices[0] . normal[2] = 1.0;
}
/*Copy previous nonzero normals going up*/
for (k = 1; k < *nVertices; ++k)
{
if (vertices[k] . normal[0] == 0.0 &&
vertices[k] . normal[1] == 0.0 &&
vertices[k] . normal[2] == 0.0)
{
vertices[k] . normal[0] = vertices[k - 1] . normal[0];
vertices[k] . normal[1] = vertices[k - 1] . normal[1];
vertices[k] . normal[2] = vertices[k - 1] . normal[2];
}
}
/*Copy next nonzero normals going down*/
for (k = (*nVertices) - 2; k >= 0; --k)
{
if (vertices[k] . normal[0] == 0.0 &&
vertices[k] . normal[1] == 0.0 &&
vertices[k] . normal[2] == 0.0)
{
vertices[k] . normal[0] = vertices[k + 1] . normal[0];
vertices[k] . normal[1] = vertices[k + 1] . normal[1];
vertices[k] . normal[2] = vertices[k + 1] . normal[2];
}
}
return vertices;
}
static ObjPtr ReadNFFXFile(name)
char *name;
/*Reads an extended NFF file from name*/
{
ObjPtr picture, retVal;
int c;
real curTime = 0.0;
ObjPtr timeSteps = NULLOBJ, timeData = NULLOBJ;
Bool timeRead = false;
ObjPtr eternalPicture = NULLOBJ;
ObjPtr timedObj;
char cmdStr[256];
FILE *inFile;
int whichStep = 0;
real ox = 0.0, oy = 0.0, oz = 0.0;
int curColorIndex = -2;
int nextColorIndex = 2;
long nColorsAllocated = 0;
short3 *curColors;
Bool *colorAllocated; /*Flag to see if this color is allocated*/
int k; /*Counter*/
int movingStart = -1; /*Start of moving colors*/
char *s; /*String pointer for testing*/
ObjPtr palette;
inFile = fopen(name, "r");
if (!inFile)
{
Error("ReadNFFXFile", OPENFILEERROR, name);
return NULLOBJ;
}
nColorsAllocated = ADDCOLORS;
curColors = (short3 *) Alloc(nColorsAllocated * sizeof(short3));
curColors[0][0] = 255;
curColors[0][1] = 255;
curColors[0][2] = 255;
curColors[1][0] = 0;
curColors[1][1] = 0;
curColors[1][2] = 0;
colorAllocated = (Bool *) Alloc(nColorsAllocated * sizeof(Bool));
for (k = 0; k < nColorsAllocated; ++k)
{
colorAllocated[k] = false;
}
picture = NewPicture();
/*Read the picture in*/
for (SkipBlanks(inFile); (c = getc(inFile)) >= 0 && c < 255 && !feof(inFile); SkipBlanks(inFile))
{
ungetc(c, inFile);
fgets(tempStr, 255, inFile);
if (1 == sscanf(tempStr, "%s", cmdStr))
{
if (*cmdStr == '#')
{
/*It's a comment*/
}
else if (0 == strcmp(cmdStr, "t"))
{
/*It's a time marker*/
real nextTime;
if (1 != sscanf(tempStr, "t %g", &nextTime))
{
FileFormatError("ReadNFFXFile", "Bad time step");
continue;
}
if (!timeRead)
{
timeRead = true;
timeSteps = NewList();
timeData = NewList();
eternalPicture = picture;
picture = NewPicture();
}
else if (picture)
{
/*Got to spit out the old picture*/
PostfixList(timeSteps, NewReal(curTime));
PostfixList(timeData, picture);
picture = NewPicture();
}
if (movingStart < 0)
{
/*First encounter with a time object*/
movingStart = nextColorIndex;
if (movingStart < 0) movingStart = 0;
}
else
{
/*Not first encounter, erase colorAllocated*/
if (nextColorIndex >= nColorsAllocated - 1)
{
nColorsAllocated += ADDCOLORS;
curColors = (short3 *) Realloc(curColors,
nColorsAllocated * sizeof(short3));
colorAllocated = (Bool *) Realloc(colorAllocated,
nColorsAllocated * sizeof(Bool));
for (k = nextColorIndex; k < nColorsAllocated; ++k)
{
curColors[k][0] = 0;
curColors[k][1] = 0;
curColors[k][2] = 0;
colorAllocated[k] = false;
}
}
for (k = movingStart; k < nextColorIndex; ++k)
{
colorAllocated[k] = false;
}
}
curTime = nextTime;
SkipBlanks(inFile);
continue;
}
else if (0 == strcmp(cmdStr, "f"))
{
float r, g, b, d;
if (8 == sscanf(tempStr, "f %g %g %g %g %g %g %g %g \n",
&r, &g, &b, &d, &d, &d, &d, &d))
{
short rs, gs, bs;
/*It's a valid color.*/
rs = r * 255.0;
gs = g * 255.0;
bs = b * 255.0;
if (movingStart >= 0)
{
/*Already in moving portion. Search for colors*/
for (k = movingStart; k < nextColorIndex; ++k)
{
if (!colorAllocated[k] &&
curColors[k][0] == rs &&
curColors[k][1] == gs &&
curColors[k][2] == bs)
{
curColorIndex = k - 2;
colorAllocated[k] = true;
goto foundColor;
}
}
}
if (nextColorIndex >= nColorsAllocated)
{
nColorsAllocated += ADDCOLORS;
curColors = (short3 *) Realloc(curColors,
nColorsAllocated * sizeof(short3));
colorAllocated = (Bool *) Realloc(colorAllocated,
nColorsAllocated * sizeof(Bool));
for (k = nextColorIndex; k < nColorsAllocated; ++k)
{
colorAllocated[k] = false;
}
}
curColorIndex = nextColorIndex - 2;
curColors[nextColorIndex][0] = rs;
curColors[nextColorIndex][1] = gs;
curColors[nextColorIndex][2] = bs;
++nextColorIndex;
}
else
{
FileFormatError("ReadNFFXFile", "Bad color format statement");
}
foundColor:;
}
else if (0 == strcmp(cmdStr, "s"))
{
/*Sphere*/
float center[3];
float radius;
if (4 == sscanf(tempStr, "s %g %g %g %g",
&(center[0]), &(center[1]), &(center[2]), &radius))
{
center[0] += ox;
center[1] += oy;
center[2] += oz;
AppendSphereToPicture(picture, center, radius, curColorIndex);
}
else
{
FileFormatError("ReadNFFXFile", "Badly formatted sphere");
}
}
else if (0 == strcmp(cmdStr, "o"))
{
/*Origin*/
sscanf(tempStr, "o %g %g %g", &ox, &oy, &oz);
}
else if (0 == strcmp(cmdStr, "c"))
{
float end1[3], end2[3], rad1, rad2;
/*It's a conical frustum*/
if (4 != fscanf(inFile, "%g %g %g %g\n",
&(end1[0]), &(end1[1]), &(end1[2]), &rad1))
{
FileFormatError("ReadNFFXFile", "Badly formed frustum");
continue;
}
end1[0] += ox;
end1[1] += oy;
end1[2] += oz;
if (4 != fscanf(inFile, "%g %g %g %g\n",
&(end2[0]), &(end2[1]), &(end2[2]), &rad2))
{
FileFormatError("ReadNFFXFile", "Badly formed frustum");
continue;
}
end2[0] += ox;
end2[1] += oy;
end2[2] += oz;
AppendFrustumToPicture(picture, end1, rad1, end2, rad2, curColorIndex);
}
else if ((0 == strcmp(cmdStr, "p")) ||
(0 == strcmp(cmdStr, "pp")))
{
int nVertices, nReadVertices;
Vertex *vertices;
int k;
/*Polygon*/
s = tempStr;
SKIPBLANKS(s);
while (!isspace(*s)) ++s;
SKIPBLANKS(s);
if (*s == '*')
{
nVertices = -1;
}
else
{
if (1 != sscanf(s, "%d", &nVertices))
{
FileFormatError("ReadNFFFile", "No vertices specified");
continue;
}
}
vertices = ReadNFFVertices(inFile, &nReadVertices);
if (nVertices > 0 && (nVertices != nReadVertices))
{
FileFormatError("ReadNFFFile", "Vertices number mismatch");
}
if (vertices)
{
for (k = 0; k < nReadVertices; ++k)
{
vertices[k] . colorIndex = curColorIndex;
}
AppendPolyToPicture(picture, nReadVertices, vertices);
Free(vertices);
}
}
else if (0 == strcmp(cmdStr, "pl"))
{
int nVertices, nReadVertices;
Vertex *vertices;
int k;
/*Polyline*/
s = tempStr;
SKIPBLANKS(s);
while (!isspace(*s)) ++s;
SKIPBLANKS(s);
if (*s == '*')
{
nVertices = -1;
}
else
{
if (1 != sscanf(s, "%d", &nVertices))
{
FileFormatError("ReadNFFFile", "No vertices specified");
continue;
}
}
vertices = ReadNFFVertices(inFile, &nReadVertices);
if (nVertices > 0 && (nVertices != nReadVertices))
{
FileFormatError("ReadNFFFile", "Vertices number mismatch");
}
if (vertices)
{
for (k = 0; k < nReadVertices; ++k)
{
vertices[k] . colorIndex = curColorIndex;
}
AppendPolylineToPicture(picture, 1, 0, nReadVertices, vertices);
Free(vertices);
}
}
else
{
char t, err[200];
t = *tempStr;
sprintf(err, "Bad object type: %c (%x)", t, t);
FileFormatError("ReadNFFXFile", err);
}
}
else
{
/*Blank line in file*/
}
}
/*Now make the object*/
retVal = NewObject(geometryClass, 0);
if (picture)
{
if (timeRead)
{
/*Got to spit out the old picture into the time var*/
PostfixList(timeSteps, NewReal(curTime));
PostfixList(timeData, picture);
}
else
{
/*Just make it eternal*/
eternalPicture = picture;
}
}
if (timeRead)
{
timedObj = NewTimedObject(timeSteps, timeData);
SetVar(retVal, DATA, timedObj);
SetVar(retVal, ETERNALPART, eternalPicture);
}
else
{
SetVar(retVal, DATA, eternalPicture);
}
SetVar(retVal, NAME, MakeDatasetName(name));
/*Create a new palette*/
if (nextColorIndex > 2)
{
curColors[nextColorIndex][0] = 255;
curColors[nextColorIndex][1] = 255;
curColors[nextColorIndex][2] = 255;
++nextColorIndex;
palette = NewPalette(nextColorIndex);
FieldPaletteName(palette, retVal);
CopyColorsToPalette(palette, curColors);
SetVar(retVal, CPALETTE, palette);
SetVar(retVal, COLORBYSELF, ObjTrue);
}
Free(curColors);
Free(colorAllocated);
fclose(inFile);
RegisterDataset(retVal);
return NULLOBJ;
}
static ObjPtr ReadJAKFile(name)
char *name;
/*Reads one of jayakumar's files*/
{
ObjPtr timeSteps, timeData;
ObjPtr timedObj;
FILE *inFile;
int whichFrame;
real dims[2]; /*The dimensions of the data, temp.*/
real bounds[6]; /*The bounds of the data form*/
long xSize, ySize;
real dataMin, dataMax;
real tempMin, tempMax;
ObjPtr dataForm;
ObjPtr dimsArray, boundsArray;
ObjPtr retVal;
inFile = fopen(name, "r");
if (!inFile)
{
Error("ReadJAKFile", OPENFILEERROR, name);
return NULLOBJ;
}
timeSteps = NewList();
timeData = NewList();
/*Read in the frames*/
whichFrame = 0;
while (fscanf(inFile, " %ld %ld", &xSize, &ySize) == 2)
{
ObjPtr data;
int i, j;
real *dataPtr;
if (fscanf(inFile, " %g %g", &tempMin, &tempMax) != 2)
{
char err[200];
sprintf(err, "Error reading minmax in frame %d", whichFrame);
FileFormatError("ReadJAKFile", err);
fclose(inFile);
return NULLOBJ;
}
if (whichFrame == 0)
{
/*If it's the first frame, set sizes and minmax*/
dims[0] = xSize;
dims[1] = ySize;
dataMin = tempMin;
dataMax = tempMax;
}
else
{
/*Enlarge minmax*/
if (tempMin < dataMin) dataMin = tempMin;
if (tempMax > dataMax) dataMax = tempMax;
}
data = NewRealArray(2, xSize, ySize);
dataPtr = ArrayMeat(data) + (ySize - 1) * xSize;
for (j = 0; j < ySize; ++j)
{
for (i = 0; i < xSize; ++i)
{
if (1 != fscanf(inFile, " %g", dataPtr))
{
char err[200];
sprintf(err, "Error in frame %d at %d %d", whichFrame, i, j);
FileFormatError("ReadJAKFile", err);
fclose(inFile);
return NULLOBJ;
}
++dataPtr;
}
dataPtr -= 2 * xSize;
}
PostfixList(timeSteps, NewReal((real) whichFrame));
PostfixList(timeData, data);
++whichFrame;
}
fclose(inFile);
/*Create the data form*/
dataForm = NewObject(dataFormClass, 0);
dimsArray = NewRealArray(1, (long) 2);
/*Put in some dimensions*/
CArray2Array(dimsArray, dims);
SetVar(dataForm, DIMENSIONS, dimsArray);
/*Put in the bounds*/
bounds[0] = 0.;
bounds[1] = 1000.0;
bounds[2] = 0.;
bounds[3] = 1000.0;
bounds[4] = dataMin;
bounds[5] = dataMax;
boundsArray = NewRealArray(1, (long) 6);
CArray2Array(boundsArray, bounds);
SetVar(dataForm, BOUNDS, boundsArray);
/*Create the field*/
retVal = NewObject(data2DScalar, 0);
SetVar(retVal, DATA, NewTimedObject(timeSteps, timeData));
SetVar(retVal, DATAFORM, dataForm);
SetVar(retVal, NAME, MakeDatasetName(name));
RegisterDataset(retVal);
return NULLOBJ;
}
#define DS(k) \
if (1 != fread(&dataSize, sizeof(long), 1, inFile) || \
dataSize != k) \
{ \
FileFormatError("ReadJAKFile", "Data length error"); \
fclose(inFile); \
return NULLOBJ; \
}
static ObjPtr ReadJAKBFile(name)
char *name;
/*Reads one of jay a kumar's binary files*/
{
ObjPtr timeSteps, timeData;
ObjPtr timedObj;
FILE *inFile;
long whichFrame;
real dims[2]; /*The dimensions of the data, temp.*/
real bounds[6]; /*The bounds of the data form*/
long xSize, ySize;
long dataSize; /*Dummy to hold size of data*/
real dataMin, dataMax;
ObjPtr dataForm;
ObjPtr dimsArray, boundsArray;
ObjPtr retVal;
long nFrames;
inFile = fopen(name, "r");
if (!inFile)
{
Error("ReadJAKFile", OPENFILEERROR, name);
return NULLOBJ;
}
timeSteps = NewList();
timeData = NewList();
whichFrame = 0;
DS(4);
if (1 != fread(&nFrames, sizeof(int), 1, inFile))
{
FileFormatError("ReadJAKFile", "Error reading nframes");
fclose(inFile);
return NULLOBJ;
}
DS(4);
DS(8);
if (1 != fread(&xSize, sizeof(int), 1, inFile))
{
FileFormatError("ReadJAKFile", "Error reading size");
fclose(inFile);
return NULLOBJ;
}
if (1 != fread(&ySize, sizeof(int), 1, inFile))
{
FileFormatError("ReadJAKFile", "Error reading size");
fclose(inFile);
return NULLOBJ;
}
DS(8);
dims[0] = xSize;
dims[1] = ySize;
DS(8);
if (1 != fread(&dataMin, sizeof(real), 1, inFile))
{
FileFormatError("ReadJAKFile", "Error reading size");
fclose(inFile);
return NULLOBJ;
}
if (1 != fread(&dataMax, sizeof(real), 1, inFile))
{
FileFormatError("ReadJAKFile", "Error reading size");
fclose(inFile);
return NULLOBJ;
}
DS(8);
FileFormatError("ReadJAKFile", "Error reading size");
/*Read in the frames*/
while (1 == fread(&dataSize, sizeof(long), 1, inFile))
{
ObjPtr data;
int i, j;
real *dataPtr;
if (dataSize != 4)
{
FileFormatError("ReadJAKFile", "Data length error");
fclose(inFile);
return NULLOBJ;
}
if (1 != fread(&whichFrame, sizeof(long), 1, inFile))
{
char err[256];
sprintf(err, "Error reading frame #\n");
FileFormatError("ReadJAKFile", err);
fclose(inFile);
return NULLOBJ;
}
DS(4);
data = NewRealArray(2, xSize, ySize);
dataPtr = ArrayMeat(data) + (ySize - 1) * xSize;
for (j = 0; j < ySize; ++j)
{
DS(xSize * sizeof(real));
if (xSize != fread(dataPtr, sizeof(real), xSize, inFile))
{
char err[200];
sprintf(err, "Error reading frame %d", whichFrame);
FileFormatError("ReadJAKBFile", err);
fclose(inFile);
return NULLOBJ;
}
DS(xSize * sizeof(real));
dataPtr -= xSize;
}
PostfixList(timeSteps, NewReal((real) whichFrame));
PostfixList(timeData, data);
}
fclose(inFile);
/*Create the data form*/
dataForm = NewObject(dataFormClass, 0);
dimsArray = NewRealArray(1, (long) 2);
/*Put in some dimensions*/
CArray2Array(dimsArray, dims);
SetVar(dataForm, DIMENSIONS, dimsArray);
/*Put in the bounds*/
bounds[0] = 0.;
bounds[1] = 1000.0;
bounds[2] = 0.;
bounds[3] = 1000.0;
bounds[4] = dataMin;
bounds[5] = dataMax;
boundsArray = NewRealArray(1, (long) 6);
CArray2Array(boundsArray, bounds);
SetVar(dataForm, BOUNDS, boundsArray);
/*Create the field*/
retVal = NewObject(data2DScalar, 0);
SetVar(retVal, DATA, NewTimedObject(timeSteps, timeData));
{
ObjPtr minMaxArray;
real minMax[2];
minMax[0] = dataMin;
minMax[1] = dataMax;
minMaxArray = NewRealArray(1, 2L);
CArray2Array(minMaxArray, minMax);
SetVar(retVal, MINMAX, minMaxArray);
}
SetVar(retVal, DATAFORM, dataForm);
SetVar(retVal, NAME, MakeDatasetName(name));
RegisterDataset(retVal);
return NULLOBJ;
}
static ObjPtr ReadSYFile(name)
char *name;
/*Reads one of Saul Youssef's fields from file name*/
{
float dx, dy, dz;
int lx, ly, lz;
float minX, minY, minZ;
float maxX, maxY, maxZ;
ObjPtr dataForm; /*The grid holding the data*/
ObjPtr boundsArray; /*Array containing bounds of grid*/
ObjPtr dimsArray; /*Array containing dimensions*/
ObjPtr fieldData; /*The array of data in the field*/
real dims[3]; /*The dimensions of the data, temp.*/
real bounds[6]; /*The bounds of the data form*/
real *dataPtr; /*The meat of the array*/
int i, j, k; /*Index into the grid*/
ObjPtr retVal; /*The field to return*/
FILE *inFile;
inFile = fopen(name, "r");
if (!inFile)
{
Error("ReadSYFile", OPENFILEERROR, name);
return NULLOBJ;
}
/*Read header*/
if (fscanf(inFile, " %d %d %d %g %g %g %g %g %g",
&lx, &ly, &lz, &minX, &maxX, &minY, &maxY, &minZ, &maxZ)
!= 9)
{
FileFormatError("ReadSYFile", "Error reading header");
}
/*Create the data form*/
dataForm = NewObject(dataFormClass, 0);
dimsArray = NewRealArray(1, (long) 3);
/*Put in some dimensions*/
dims[0] = (real) lx;
dims[1] = (real) ly;
dims[2] = (real) lz;
CArray2Array(dimsArray, dims);
SetVar(dataForm, DIMENSIONS, dimsArray);
/*Put in the bounds*/
bounds[0] = minX;
bounds[1] = maxX;
bounds[2] = minY;
bounds[3] = maxY;
bounds[4] = minZ;
bounds[5] = maxZ;
boundsArray = NewRealArray(1, (long) 6);
CArray2Array(boundsArray, bounds);
SetVar(dataForm, BOUNDS, boundsArray);
/*Read the data*/
fieldData = NewRealArray(3, (long) lx, (long) ly, (long) lz);
dataPtr = ArrayMeat(fieldData);
for (k = 0; k < lz; ++k)
{
for (j = 0; j < ly; ++j)
{
for (i = 0; i < lx; ++i)
{
real curVal;
if (1 != fscanf(inFile, " %g", &curVal))
{
FileFormatError("ReadSYFile", "Error reading data");
fclose(inFile);
return NULLOBJ;
}
*(dataPtr + i * ly * lz + j * lz + k) = curVal;
}
}
}
/*Create the field*/
retVal = NewObject(data3DScalar, 0);
SetVar(retVal, DATA, fieldData);
SetVar(retVal, DATAFORM, dataForm);
SetVar(retVal, NAME, MakeDatasetName(name));
fclose(inFile);
RegisterDataset(retVal);
return NULLOBJ;
}
#ifdef LORENZ
#define LRZNMAX 10
#ifdef FORTRAN_
extern void odeint_();
#else
extern void odeint();
#endif
static ObjPtr ReadLRZFile(name)
char *name;
/*Reads a Lorenz attractor initial value file*/
{
float y[LRZNMAX]; /*3 initial values*/
float dydx[LRZNMAX];
int nstep, nskip;
int nvar;
float xstart;
float eps, stpmin, stpsize, sig, b, r;
real *dataPtr;
ObjPtr theData, theClass, lrzObj;
FILE *inFile;
inFile = fopen(name, "r");
if (!inFile)
{
Error("ReadLRZFile", OPENFILEERROR, name);
return NULLOBJ;
}
/*Read the file*/
if (11 != fscanf(inFile, " %f %f %f %f %d %d %f %f %f %f %f",
&(y[0]), &(y[1]), &(y[2]), &stpsize, &nstep,
&nskip, &eps, &stpmin, &sig, &b, &r))
{
FileFormatError("ReadLRZFile", "Cannot read parameter file");
fclose(inFile);
return NULLOBJ;
}
/*Now we know how big it is, make some data*/
theData = NewRealArray(2, (long) nstep, (long) LRZNMAX);
if (!theData)
{
fclose(inFile);
return NULLOBJ;
}
dataPtr = ArrayMeat(theData);
nvar = 3;
xstart = 0.0;
#ifdef FORTRAN_
odeint_
#else
odeint
#endif
(y, &nvar, &xstart, &nskip, &nstep, dataPtr,
&eps, &stpsize, &stpmin, &sig, &b, &r);
/*Now create the data object*/
theClass = NewObject(data1DVector, 0);
lrzObj = NewObject(theClass, 0);
SetVar(lrzObj, DATA, theData);
SetVar(lrzObj, NAME, MakeDatasetName(name));
fclose(inFile);
RegisterDataset(lrzObj);
return NULLOBJ;
}
#endif
static ObjPtr ReadDDFile(name)
char *name;
/*Reads a Dennis Duke format file*/
{
long idim, jdim;
int c, k;
ObjPtr theData, theClass, eegData;
FILE *inFile;
long index[1];
inFile = fopen(name, "r");
if (!inFile)
{
Error("ReadDDFile", OPENFILEERROR, name);
return NULLOBJ;
}
/*First determine idim and jdim*/
idim = 0; jdim = 0;
SkipBlanks(inFile);
while ((c = getc(inFile)) >= 0)
{
ungetc(c, inFile);
++jdim;
k = 0;
while ((c = getc(inFile)) >= 0 && c != '\n')
{
ungetc(c, inFile);
SkipNonBlanks(inFile);
SkipBlanks(inFile);
++k;
}
if (idim > 0)
{
if (k != idim)
{
char err[200];
sprintf(err,
"Error in %s line %d: bad number of items\n",
name, jdim);
FileFormatError("ReadDDFile", err);
fclose(inFile);
return NULLOBJ;
}
else
{
idim = k;
}
}
else
{
idim = k;
}
getc(inFile);
SkipBlanks(inFile);
}
/*Now we know how big it is, make some data*/
theData = NewArray(AT_OBJECT, 1, &idim);
for (k = 0; k < idim; ++k)
{
((ObjPtr *) ELEMENTS(theData))[k] = NewRealArray(1, (long) jdim);
}
/*Now create the data object*/
theClass = NewObject(data1DVector, 0);
eegData = NewObject(theClass, 0);
SetVar(eegData, DATA, theData);
SetVar(eegData, NAME, MakeDatasetName(name));
if (!SetCurField(FIELD1, eegData)) return NULLOBJ;
/*Read in the data*/
rewind(inFile);
SkipBlanks(inFile);
index[0] = 0;
while ((c = getc(inFile)) >= 0)
{
ungetc(c, inFile);
k = 0;
while ((c = getc(inFile)) >= 0 && c != '\n')
{
real data;
ungetc(c, inFile);
if (1 != fscanf(inFile, "%g", &data))
{
FileFormatError("ReadDDFile", "Bad number");
fclose(inFile);
return NULLOBJ;
}
PutFieldComponent(FIELD1, k, index, data);
++k;
if (k >= idim)
{
k = 0;
++index[0];
}
SkipBlanksAndCommas(inFile);
}
getc(inFile);
SkipBlanksAndCommas(inFile);
}
fclose(inFile);
SetVar(eegData, NCOMPONENTS, NewInt(idim));
RegisterDataset(eegData);
return NULLOBJ;
}
static ObjPtr brainForm = 0;
long numNodes; /*Number of nodes in brain form*/
/*States for data form reading machine*/
#define S_READNODES 1
#define S_READEDGES 2
#define S_READCELLS 3
#define S_READPOLATIONS 4
#define DD2SECS 8.0
#define DD2SAMPLES 1024
#define DD2POINTS 19
#define NPOLATIONS 100
int pDerived[NPOLATIONS];
int pNear[NPOLATIONS];
int pFar[NPOLATIONS];
real pFac[NPOLATIONS];
int nPolations = 0;
static ObjPtr ReadDD2File(name)
char *name;
/*Reads a brain surface EEG file*/
{
ObjPtr retVal;
ObjPtr timeSteps;
ObjPtr timeData;
FILE *inFile;
int sample;
long k;
real *data, dummy;
real minMax[2];
ObjPtr minMaxArray;
if (!brainForm)
{
/*Must make the brain data form*/
threeArray *nodeBuf;
long nodeBufSize;
TwoReals *edgeBuf;
long edgeBufSize;
long numEdges;
long numCells;
ObjPtr cellList;
ObjPtr nodeData;
ObjPtr edgeArray;
int state;
char line[257];
char token[257];
char *s;
real val;
int lineNum;
char pTypeChar;
long tempDims[3];
real bounds[6];
/*Open the file*/
lineNum = 0;
state = 0;
inFile = fopen("brainform", "r");
if (!inFile)
{
FileFormatError("ReadDD2File", "File brainform was not found");
return NULLOBJ;
}
/*Create the buffers. First the nodes*/
numNodes = 0;
nodeBufSize = 1000;
nodeBuf = (threeArray *) Alloc(nodeBufSize * 3 * sizeof(real));
if (!nodeBuf)
{
OMErr();
return NULLOBJ;
}
/*Then the edges*/
numEdges = 0;
edgeBufSize = 1000;
edgeBuf = (TwoReals *) Alloc(edgeBufSize * 2 * sizeof(real));
if (!edgeBuf)
{
OMErr();
Free(nodeBuf);
return NULLOBJ;
}
/*Now the cells*/
numCells = 0;
cellList = NewList();
if (!cellList)
{
Free(edgeBuf);
Free(nodeBuf);
return NULLOBJ;
}
/*Now read the file*/
while (++lineNum, fgets(line, 256, inFile))
{
s = &(line[0]);
/*Get the first token*/
SHIFTNUM(token, s);
if (0 == strcmp2(token, "NODES"))
{
state = S_READNODES;
}
else if (0 == strcmp2(token, "EDGES"))
{
state = S_READEDGES;
}
else if (0 == strcmp2(token, "CELLS"))
{
state = S_READCELLS;
}
else if (0 == strcmp2(token, "POLATIONS"))
{
state = S_READPOLATIONS;
}
else if (1 == sscanf(token, "%g", &val))
{
/*It's a number. Do something based on the state*/
switch(state)
{
case S_READNODES:
/*It must be the x of a node*/
{
real x, y, z;
if (numNodes > nodeBufSize)
{
/*Must expand the node buffer*/
nodeBufSize += 200;
nodeBuf = (threeArray *) Realloc(nodeBuf, nodeBufSize * sizeof(threeArray));
if (!nodeBuf)
{
OMErr();
Free(edgeBuf);
return NULLOBJ;
}
}
x = val;
SHIFTNUM(token, s);
if (1 != sscanf(token, "%g", &y))
{
char err[200];
sprintf(err, "Error in line %d of brainform: Missing y\n", lineNum);
FileFormatError("ReadDD2File", err);
break;
}
SHIFTNUM(token, s);
if (1 != sscanf(token, "%g", &z))
{
char err[200];
sprintf(err, "Error in line %d of brainform: Missing z\n", lineNum);
FileFormatError("ReadDD2File", err);
break;
}
nodeBuf[numNodes][0] = x;
nodeBuf[numNodes][1] = y;
nodeBuf[numNodes][2] = z;
++numNodes;
}
break;
case S_READEDGES:
/*It must be the first node of an edge*/
{
real n1, n2;
if (numEdges > edgeBufSize)
{
/*Must expand the node buffer*/
edgeBufSize += 200;
edgeBuf = (TwoReals *) Realloc(edgeBuf, edgeBufSize * sizeof(TwoReals));
if (!nodeBuf)
{
OMErr();
Free(nodeBuf);
return NULLOBJ;
}
}
n1 = val;
SHIFTNUM(token, s);
if (1 != sscanf(token, "%g", &n2))
{
char err[200];
sprintf(err, "Error in line %d of brainform: Missing second node\n", lineNum);
FileFormatError("ReadDD2File", err);
break;
}
edgeBuf[numEdges][0] = n1;
edgeBuf[numEdges][1] = n2;
++numEdges;
}
break;
case S_READCELLS:
/*It must be the first node of a cell*/
{
real cellBuf[600];
ObjPtr cellArray;
long c;
c = 0;
do
{
cellBuf[c++] = val;
SHIFTNUM(token, s);
} while (1 == sscanf(token, "%g", &val));
cellArray = NewRealArray(1, c);
CArray2Array(cellArray, cellBuf);
PostfixList(cellList, cellArray);
++numCells;
}
break;
default:
FileFormatError("ReadDD2File", "Error in brainform");
}
}
else if (1 == sscanf(token, "%c", &pTypeChar))
{
switch (state)
{
case S_READPOLATIONS:
SHIFTNUM(token, s);
if (1 != sscanf(token, "%d", &(pDerived[nPolations])))
{
FileFormatError("ReadDD2File", "Error in brainform");
break;
}
SHIFTNUM(token, s);
if (1 != sscanf(token, "%d", &(pNear[nPolations])))
{
FileFormatError("ReadDD2File", "Error in brainform");
break;
}
SHIFTNUM(token, s);
if (1 != sscanf(token, "%d", &(pFar[nPolations])))
{
FileFormatError("ReadDD2File", "Error in brainform");
break;
}
/*Calculate multiplying factor*/
{
int n, f, d;
float d1, d2;
n = pNear[nPolations];
f = pFar[nPolations];
d = pDerived[nPolations];
d1 = sqrt(SQUARE(nodeBuf[f][0] - nodeBuf[n][0]) +
SQUARE(nodeBuf[f][1] - nodeBuf[n][1]) +
SQUARE(nodeBuf[f][2] - nodeBuf[n][2]));
d2 = sqrt(SQUARE(nodeBuf[d][0] - nodeBuf[n][0]) +
SQUARE(nodeBuf[d][1] - nodeBuf[n][1]) +
SQUARE(nodeBuf[d][2] - nodeBuf[n][2]));
if (pTypeChar == 'e' || pTypeChar == 'E')
{
pFac[nPolations] = - d2 / d1;
}
else
{
pFac[nPolations] = d2 / d1;
}
}
++nPolations;
break;
default:
FileFormatError("ReadDD2File", "Error in brainform");
}
}
else
{
FileFormatError("ReadDD2File", "Error in brainform");
}
}
fclose(inFile);
/*Everything's been read.*/
/*Make the vector dataset containing the nodes*/
tempDims[0] = numNodes;
tempDims[1] = numEdges;
tempDims[2] = numCells;
nodeData = NewDataset("brainform", 1, tempDims, 3);
SetCurField(FIELD1, nodeData);
bounds[0] = bounds[2] = bounds[4] = PLUSINF;
bounds[1] = bounds[3] = bounds[5] = MINUSINF;
for (k = 0; k < numNodes; ++k)
{
if (nodeBuf[k][0] < bounds[0]) bounds[0] = nodeBuf[k][0];
if (nodeBuf[k][0] > bounds[1]) bounds[1] = nodeBuf[k][0];
PutFieldComponent(FIELD1, 0, &k, nodeBuf[k][0]);
if (nodeBuf[k][0] < bounds[2]) bounds[2] = nodeBuf[k][0];
if (nodeBuf[k][0] > bounds[3]) bounds[3] = nodeBuf[k][0];
PutFieldComponent(FIELD1, 1, &k, nodeBuf[k][1]);
if (nodeBuf[k][0] < bounds[4]) bounds[4] = nodeBuf[k][0];
if (nodeBuf[k][0] > bounds[5]) bounds[5] = nodeBuf[k][0];
PutFieldComponent(FIELD1, 2, &k, nodeBuf[k][2]);
}
brainForm = NewUnstructuredDataForm("brainform", 2, tempDims, bounds, nodeData);
edgeArray = NewRealArray(2, numEdges, 2L);
if (!edgeArray)
{
Free(nodeBuf);
Free(edgeBuf);
return NULLOBJ;
}
CArray2Array(edgeArray, edgeBuf);
SetVar(brainForm, EDGES, edgeArray);
SetVar(brainForm, CELLS, cellList);
Free(nodeBuf);
Free(edgeBuf);
}
/*Get ready to read*/
inFile = fopen(name, "r");
if (!inFile)
{
Error("ReadDD2File", OPENFILEERROR, name);
return NULLOBJ;
}
timeSteps = NewList();
timeData = NewList();
data = (real *) Alloc(sizeof(real) * numNodes);
for (sample = 0; sample < DD2SAMPLES; ++sample)
{
ObjPtr sampleArray;
if (21 != fscanf(inFile, " %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g",
&data[0],
&dummy,
&data[1],
&data[2],
&data[3],
&data[4],
&data[5],
&data[6],
&data[7],
&data[8],
&data[9],
&data[10],
&data[11],
&data[12],
&data[13],
&data[14],
&data[15],
&data[16],
&data[17],
&dummy,
&data[18]
))
{
FileFormatError("ReadDD2File", "Not enough data");
fclose(inFile);
return NULLOBJ;
}
/*Fill in the polated data*/
for (k = 0; k < nPolations; ++k)
{
data[pDerived[k]] = data[pNear[k]] +
(data[pFar[k]] - data[pNear[k]]) * pFac[k];
}
sampleArray = NewRealArray(1, (long) numNodes);
CArray2Array(sampleArray, data);
PostfixList(timeData, sampleArray);
PostfixList(timeSteps, NewReal(DD2SECS * ((real) sample) / ((real) DD2SAMPLES)));
}
Free(data);
fclose(inFile);
retVal = NewObject(data3DUnstructSurface, 0);
SetVar(retVal, DATA, NewTimedObject(timeSteps, timeData));
SetVar(retVal, DATAFORM, brainForm);
SetVar(retVal, NAME, MakeDatasetName(name));
minMax[0] = -60.0;
minMax[1] = 60.0;
minMaxArray = NewRealArray(1, 2L);
CArray2Array(minMaxArray, minMax);
SetVar(retVal, MINMAX, minMaxArray);
RegisterDataset(retVal);
return NULLOBJ;
}
#define MAX_HDF_RANK 100 /*Maximum rank of HDF array*/
/*Values for out of bounds data*/
#define OB_MISSING 0 /*Treat as missing*/
#define OB_CLIP 1 /*Clip to bounds*/
#define OB_USE 2 /*Just use as is*/
#ifdef HDFDEF
ObjPtr ReadHDFFile(reader, fileName)
ObjPtr reader;
char *fileName;
/*Reads an HDF scientific data set (SDS) file.*/
{
int nDataSets; /*Number of data sets*/
int whichSet; /*Index into current data set*/
ObjPtr var; /*Random variable*/
int outOfBoundsAction; /*Action on out of bounds, based on
out of bounds radio button*/
/*Step 1: Look for scientific data sets*/
/*Find out the number of scientific datasets in the file*/
nDataSets = DFSDnumber(fileName);
if (nDataSets == -1)
{
char errMes[255];
sprintf(errMes, "Cannot read HDF file %s. The file is either missing or empty.\n", fileName);
FileFormatError("ReadHDFFile", errMes);
return NULLOBJ;
}
/*Determine action for out-of-bounds*/
if ((var = GetIntVar("ReadHDFFile", reader, OUTOFBOUNDS)) == NULLOBJ)
{
outOfBoundsAction = OB_MISSING;
}
else
{
outOfBoundsAction = GetInt(var);
}
/*Go through the datasets, reading them in*/
for (whichSet = 0; whichSet < nDataSets; ++whichSet)
{
/*Dataset characteristics from file*/
int initRank; /*Rank of the dataset*/
long initDimSizes[MAX_HDF_RANK];/*Size of each dimension*/
int adjustedRank; /*Adjusted rank*/
long *adjustedDimSizes; /*Adjusted dim sizes*/
int initDimPermutations[MAX_HDF_RANK];
/*Permutations of dimensions, i.e.,
for each dimension, which spatial
dim it corresponds to*/
int *adjustedDimPermutations; /*Adjusted for vector*/
Bool dimUsed[MAX_HDF_RANK]; /*Flag to see if dim used*/
real **initScales = 0; /*Array for initial scales*/
real **adjustedScales; /*Scales adjusted for vector*/
real **permutedScales = 0; /*Permuted scales*/
int nSpatialDimensions; /*# of spatial dimensions used*/
int *hdfToTopological = 0; /*HDF dimension to topological dimension*/
long *permutedDimSizes = 0;
float min = PLUSINF; /*Current minimum and maximum*/
float max = MINUSINF;
char dimsLabel[256]; /*Label of a dimension, used to determine which one*/
char dimsUnit[256]; /*Units of a dimension, not used*/
char dimsFormat[256]; /*Format of a dimension, not used*/
char datasetName[256]; /*Name of a dataset*/
char dataUnit[256];
char dataFormat[256];
char dataCoord[256];
real cTable[256]; /*Compression table*/
char *blanklessName; /*Dataset name w/o leading blanks*/
int k; /*Random counter*/
int vectorDim; /*Dimension of vector*/
int nComponents; /*# of components*/
ObjPtr xName = NULLOBJ;
ObjPtr yName = NULLOBJ;
ObjPtr zName = NULLOBJ;
Bool isLeftHanded = false; /*True iff dataset is likely to be left handed*/
/*Variables for the dataset*/
ObjPtr curField; /*Field being assembled*/
ObjPtr dataForm; /*Data form*/
int whichDim; /*Which dimension working on*/
int i; /*Random counter*/
long arraySize; /*Total size of the array*/
real intMissing = missingData; /*Internal missing data*/
Bool minMaxSet = false; /*True iff min and max have been set*/
Bool readVector = false; /*True iff reading vector*/
Bool wrapPolar = false; /*True iff wrapping polar*/
Bool compressData = true; /*True iff compress data*/
Bool botherWithAxes = true; /*True iff compress data*/
char dummy[256]; /*Dummy value for ignored strings*/
float *tempData;
long dataIndex;
unsigned char cd;
long *index, *permutedIndex;
long componentSize; /*Size of each component*/
/*Start off with null field*/
curField = 0;
/*Assume it will be x, y, z*/
xName = NewString("X");
yName = NewString("Y");
zName = NewString("Z");
/*Get dimensions of dataset*/
if (-1 == DFSDgetdims(fileName, &initRank, initDimSizes, (int) MAX_HDF_RANK))
{
FileFormatError("ReadHDFFile", "Unable to get dimensions.");
return NULLOBJ;
}
/*See if there's a min and max defined in the dataset*/
if (0 == DFSDgetrange(&max, &min))
{
if (min > max)
{
real temp;
temp = max;
max = min;
min = temp;
}
minMaxSet = true;
}
else
{
min = MINUSINF;
max = PLUSINF;
minMaxSet = false;
}
/*Temporary holding place for scales*/
if (!(initScales = Alloc(sizeof(real *) * initRank)))
{
OMErr();
return NULLOBJ;
}
for (whichDim = 0; whichDim < initRank; ++whichDim)
{
initScales[whichDim] = 0;
}
/*See if there's any possibility of reading vector data*/
readVector = GetPredicate(reader, READVECTOR);
/*See if there's any possibility of wrapping polar*/
wrapPolar = GetPredicate(reader, WRAPPOLAR);
/*See if we should bother with axes*/
botherWithAxes = GetPredicate(reader, BOTHERWITHAXES);
/*See if we want to compress data*/
compressData = GetPredicate(reader, COMPRESSDATA);
/*See if it's a scalar or vector dataset*/
if (readVector && (initDimSizes[0] == 2 || initDimSizes[0] == 3))
{
/*It's vector, with the vector dimension at 0*/
vectorDim = 0;
nComponents = initDimSizes[0];
adjustedRank = initRank - 1;
adjustedDimSizes = &(initDimSizes[1]);
adjustedScales = &(initScales[1]);
adjustedDimPermutations = &(initDimPermutations[1]);
}
else if (readVector && (initDimSizes[initRank - 1] == 2 || initDimSizes[initRank - 1] == 3))
{
/*It's vector, with the vector dimension at initRank - 1*/
vectorDim = initRank - 1;
nComponents = initDimSizes[initRank - 1];
adjustedRank = initRank - 1;
adjustedDimSizes = initDimSizes;
adjustedScales = initScales;
adjustedDimPermutations = initDimPermutations;
}
else
{
/*It's scalar*/
vectorDim = -1;
adjustedRank = initRank;
adjustedDimSizes = initDimSizes;
adjustedScales = initScales;
adjustedDimPermutations = initDimPermutations;
}
/*Get the name of the dataset*/
datasetName[0] = '\0';
if ((-1 == DFSDgetdatastrs(datasetName, dataUnit, dataFormat, dataCoord))
|| (datasetName[0] == '\0'))
{
strcpy(datasetName, fileName);
}
else
{
/*See about the coordinate system*/
if ((0 == strcmp2(dataCoord, "polar")) ||
(0 == strcmp2(dataCoord, "spherical")))
{
/*It's polar coordinates*/
}
else
{
/*It's not polar coordinates*/
wrapPolar = false;
}
}
/*Get the dimension information*/
for (whichDim = 0; whichDim < initRank; ++whichDim)
{
float *tempScale;
/*Get the strings for this dimension.*/
initDimPermutations[whichDim] = -1;
switch (whichDim)
{
case 0:
strcpy(dimsLabel, "X");
break;
case 1:
strcpy(dimsLabel, "Y");
break;
case 2:
strcpy(dimsLabel, "Z");
break;
default:
strcpy(dimsLabel, "");
break;
}
if (botherWithAxes)
{
if (DFSDgetdimstrs(whichDim + 1, dimsLabel, dimsUnit, dimsFormat) != -1)
{
/*See if dimsLabel corresponds to a dimension*/
initDimPermutations[whichDim] = FindAxisDimension(dimsLabel);
if (initDimPermutations[whichDim] < 0)
{
/*Not found*/
sprintf(tempStr, "Warning: Invalid axis", dimsLabel);
FileFormatError("ReadHDFFile", tempStr);
}
else
{
switch(initDimPermutations[whichDim])
{
case 0:
xName = NewString(dimsLabel);
break;
case 1:
yName = NewString(dimsLabel);
break;
case 2:
zName = NewString(dimsLabel);
break;
default:
break;
}
}
}
}
/*Set up a place for the scales*/
initScales[whichDim] = (real *) Alloc(initDimSizes[whichDim] * sizeof(real));
if (!initScales[whichDim]) {OMErr(); return NULLOBJ;}
/*And a temporary scale*/
tempScale = (float *) Alloc(initDimSizes[whichDim] * sizeof(float));
if (!tempScale) {OMErr(); return NULLOBJ;}
/*Get the scale from the file*/
if (DFSDgetdimscale(whichDim + 1, initDimSizes[whichDim], tempScale) != -1)
{
double circ = 0.0;
/*Just copy the read scale*/
for (i = 0; i < initDimSizes[whichDim]; ++i)
{
initScales[whichDim][i] = tempScale[i];
}
/*See if we need to wrap*/
if (wrapPolar)
{
if ((strcmp(dimsUnit, "degrees") == 0)) circ = 360.0;
else if ((strcmp(dimsUnit, "radians") == 0)) circ = 2.0 * M_PI;
else if ((strcmp(dimsUnit, "gradians") == 0)) circ = 400.0;
if (circ > 0.0)
{
/*Let's wrap*/
if (initScales[whichDim][1] < initScales[whichDim][0])
{
/*Down*/
for (i = 2; i < initDimSizes[whichDim]; ++i)
{
while (initScales[whichDim][i] > initScales[whichDim][i - 1])
{
initScales[whichDim][i] -= circ;
}
}
}
else
{
/*Up*/
for (i = 2; i < initDimSizes[whichDim]; ++i)
{
while (initScales[whichDim][i] < initScales[whichDim][i - 1])
{
initScales[whichDim][i] += circ;
}
}
}
}
}
}
else
{
/*Guess the scale*/
if (whichDim != vectorDim)
{
sprintf(tempStr, "Guessing scale %d from %d to %d", whichDim, 0, initDimSizes[whichDim] - 1);
FileFormatError("ReadHDFFile", tempStr);
}
/*Fill in the scale based on the guess*/
for (i = 0; i < initDimSizes[whichDim]; ++i)
{
initScales[whichDim][i] = (real) i;
}
}
Free(tempScale);
}
/*Fill in unpermuted axes*/
for (k = 0; k < MAX_HDF_RANK; ++k)
{
dimUsed[k] = false;
}
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
if (adjustedDimPermutations[whichDim] >= 0)
{
if (dimUsed[adjustedDimPermutations[whichDim]])
{
sprintf(tempStr, "Dimension %d is used more than once", whichDim);
FileFormatError("ReadHDFFile", tempStr);
adjustedDimPermutations[whichDim] = -1;
}
else
{
dimUsed[adjustedDimPermutations[whichDim]] = true;
}
}
}
k = 0;
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
if (adjustedDimPermutations[whichDim] < 0)
{
while (dimUsed[k]) ++k;
adjustedDimPermutations[whichDim] = k;
dimUsed[k] = true;
++k;
}
}
/*Find # of spatial dimensions used*/
nSpatialDimensions = 0;
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
nSpatialDimensions = MAX(nSpatialDimensions, adjustedDimPermutations[whichDim] + 1);
}
/*Make HDF to topological dimension mapping*/
hdfToTopological = (int *) Alloc(sizeof(int) * adjustedRank);
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
hdfToTopological[whichDim] = 0;
for (k = 0; k < adjustedDimPermutations[whichDim]; ++k)
{
if (dimUsed[k])
{
++hdfToTopological[whichDim];
}
}
}
/*Make permuted dim sizes*/
permutedDimSizes = (long *) Alloc(sizeof(long) * adjustedRank);
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
permutedDimSizes[hdfToTopological[whichDim]] = adjustedDimSizes[whichDim];
}
/*Make permuted scales*/
permutedScales = (real **) Alloc(sizeof(real *) * adjustedRank);
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
permutedScales[hdfToTopological[whichDim]] = adjustedScales[whichDim];
}
/*See if it's left-handed or not*/
isLeftHanded = false;
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
if (adjustedScales[whichDim][1] < adjustedScales[whichDim][0])
{
isLeftHanded = isLeftHanded ? false : true;
}
}
/*Remove leading blanks on label*/
blanklessName = datasetName;
while (*blanklessName && *blanklessName == ' ') ++blanklessName;
/*Make the new data form*/
dataForm = NewSeparableDataForm(blanklessName, adjustedRank, permutedDimSizes, permutedScales);
/*Calculate the size of the array we'll need for the dataset*/
arraySize = 1;
for (i = 0; i < initRank; ++i)
{
arraySize *= initDimSizes[i];
}
/*Make the temp data*/
tempData = (float *) Alloc(arraySize * sizeof(float));
/*Calculate the size of one component*/
if (vectorDim < 0)
{
componentSize = arraySize;
}
else
{
componentSize = arraySize / nComponents;
}
/*Read in the data*/
if (-1 == DFSDgetdata(fileName, initRank, initDimSizes, tempData))
{
FileFormatError("ReadHDFFile", "Can't read field data.");
return NULLOBJ;
}
if (minMaxSet && (outOfBoundsAction != OB_USE))
{
/*Search for out of bounds data*/
for (dataIndex = 0; dataIndex < arraySize; ++dataIndex)
{
if (tempData[dataIndex] < min)
{
tempData[dataIndex] = (outOfBoundsAction == OB_CLIP) ? min : missingData;
}
else
{
if (tempData[dataIndex] > max)
{
tempData[dataIndex] = (outOfBoundsAction == OB_CLIP) ? max : missingData;
}
}
}
}
if (!minMaxSet)
{
/*Go through and look for min and max*/
for (dataIndex = 0; dataIndex < arraySize; ++dataIndex)
{
if (tempData[dataIndex] != missingData)
{
if (minMaxSet)
{
min = MIN(min, tempData[dataIndex]);
max = MAX(max, tempData[dataIndex]);
}
else
{
min = max = tempData[dataIndex];
minMaxSet = true;
}
}
}
}
/*Make compression table*/
if (compressData)
{
cTable[0] = missingData;
for (k = 1; k < 256; ++k)
{
cTable[k] = min + (max - min) * (k - 1) / 254.0;
}
}
/*Make the new dataset*/
if (compressData)
{
curField = NewCompressedDataset(blanklessName, adjustedRank, permutedDimSizes, vectorDim > 0 ? nComponents : 0, cTable);
}
else
{
curField = NewDataset(blanklessName, adjustedRank, permutedDimSizes, vectorDim > 0 ? nComponents : 0);
}
if (isLeftHanded)
{
SetVar(curField, ISLEFTHANDED, ObjTrue);
}
/*Put in axis names*/
if (xName)
{
SetVar(curField, XNAME, xName);
}
if (yName)
{
SetVar(curField, YNAME, yName);
}
if (zName)
{
SetVar(curField, ZNAME, zName);
}
/*Link the data form with the field*/
SetVar(curField, DATAFORM, dataForm);
/*Set the current field for reading data in*/
SetCurField(FIELD1, curField);
/*Set up the indices*/
index = (long *) Alloc(sizeof(long) * adjustedRank);
permutedIndex = (long *) Alloc(sizeof(long) * adjustedRank);
if (vectorDim < 0)
{
/*Scalar dataset*/
/*Zero the indices*/
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
index[whichDim] = 0;
}
for (dataIndex = 0; dataIndex < componentSize; ++dataIndex)
{
/*Permute the index*/
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
permutedIndex[hdfToTopological[whichDim]] = index[whichDim];
}
if (compressData)
{
if (tempData[dataIndex] == missingData)
{
cd = 0;
}
else
{
cd = 1.5 + (tempData[dataIndex] - min) / (max - min) * 254.0;
if (cd < 1) cd = 1;
if (cd > 255) cd = 255;
}
PutCompressedFieldComponent(FIELD1, 0, permutedIndex, cd);
}
else
{
PutFieldComponent(FIELD1, 0, permutedIndex, tempData[dataIndex]);
}
/*Advance to next index*/
for (whichDim = adjustedRank - 1; whichDim >= 0; --whichDim)
{
if ((++index[whichDim]) >= adjustedDimSizes[whichDim])
{
index[whichDim] = 0;
}
else
{
break;
}
}
}
}
else if (vectorDim == 0)
{
float *dataChunk;
/*Vector index on 0*/
for (k = 0; k < nComponents; ++k)
{
/*Zero the indices*/
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
index[whichDim] = 0;
}
dataChunk = tempData + k * componentSize;
for (dataIndex = 0; dataIndex < componentSize; ++dataIndex)
{
/*Permute the index*/
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
permutedIndex[hdfToTopological[whichDim]] = index[whichDim];
}
if (compressData)
{
if (dataChunk[dataIndex] == missingData)
{
cd = 0;
}
else
{
cd = 1.5 + (dataChunk[dataIndex] - min) / (max - min) * 254.0;
if (cd < 1) cd = 1;
if (cd > 255) cd = 255;
}
PutCompressedFieldComponent(FIELD1, k, permutedIndex, cd);
}
else
{
PutFieldComponent(FIELD1, k, permutedIndex, dataChunk[dataIndex]);
}
/*Advance to next index*/
for (whichDim = adjustedRank - 1; whichDim >= 0; --whichDim)
{
if ((++index[whichDim]) >= adjustedDimSizes[whichDim])
{
index[whichDim] = 0;
}
else
{
break;
}
}
}
}
}
else
{
/*Vector index on last dim*/
/*Zero the indices*/
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
index[whichDim] = 0;
}
for (dataIndex = 0; dataIndex < componentSize * nComponents; dataIndex += nComponents)
{
/*Permute the index*/
for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
{
permutedIndex[hdfToTopological[whichDim]] = index[whichDim];
}
for (k = 0; k < nComponents; ++k)
{
if (compressData)
{
if (tempData[dataIndex + k] == missingData)
{
cd = 0;
}
else
{
cd = 1.5 + (tempData[dataIndex + k] - min) / (max - min) * 254.0;
if (cd < 1) cd = 1;
if (cd > 255) cd = 255;
}
PutCompressedFieldComponent(FIELD1, k, permutedIndex, cd);
}
else
{
PutFieldComponent(FIELD1, k, permutedIndex, tempData[dataIndex + k]);
}
}
/*Advance to next index*/
for (whichDim = adjustedRank - 1; whichDim >= 0; --whichDim)
{
if ((++index[whichDim]) >= adjustedDimSizes[whichDim])
{
index[whichDim] = 0;
}
else
{
break;
}
}
}
}
SAFEFREE(index);
SAFEFREE(permutedIndex);
/*Set the min and max*/
if (minMaxSet)
{
if (vectorDim < 0)
{
/*Scalar field, just do as normal*/
real minMax[2];
ObjPtr minMaxArray;
minMax[0] = min;
minMax[1] = max;
minMaxArray = NewRealArray(1, 2L);
CArray2Array(minMaxArray, minMax);
SetVar(curField, MINMAX, minMaxArray);
}
else
{
/*Vector field, adjust it a bit*/
real minMax[2];
ObjPtr minMaxArray;
minMax[0] = 0.0;
minMax[1] = MAX(max, -min);
minMaxArray = NewRealArray(1, 2L);
CArray2Array(minMaxArray, minMax);
SetVar(curField, MINMAX, minMaxArray);
}
}
/*Free the temp data*/
SAFEFREE(tempData);
SAFEFREE(permutedDimSizes);
SAFEFREE(hdfToTopological);
/*Free scales*/
for (i = 0; i < initRank; ++i)
{
if (initScales[i]) Free(initScales[i]);
}
SAFEFREE(initScales);
SAFEFREE(permutedScales);
/*Register the dataset*/
if (curField)
{
RegisterDataset(curField);
}
}
/*Make the file go back to the beginning*/
DFSDrestart();
return NULLOBJ;
} /* ReadHDFFile */
#endif
#ifdef PROTO
ObjPtr NewFileReader(char *name)
#else
ObjPtr NewFileReader(name)
char *name;
#endif
{
ObjPtr fileReader;
ThingListPtr runner;
ObjPtr nameVar;
fileReader = NewObject(fileReaderClass, 0);
SetVar(fileReader, NAME, NewString(name));
AddObjToDatabase(fileReader);
return fileReader;
}
void DefineFormat(name, extension, reader, writer)
char *name;
char *extension;
ObjPtr (*reader)();
Bool (*writer)();
/*Defines a format with name, reader, and writer*/
{
ObjPtr fileReader;
fileReader = NewFileReader(name);
SetVar(fileReader, EXTENSION, NewString(extension));
SetMethod(fileReader, OLDREAD, reader);
ApplySavedSettings(fileReader);
}
ObjPtr FindFormatReader(format)
char *format;
/*Finds the file reader named format*/
{
ObjPtr keyList;
ThingListPtr runner;
ObjPtr fileReaders;
/*Search database for datasets*/
keyList = NewList();
PostfixList(keyList, NewSymbol(CLASSID));
PostfixList(keyList, NewInt(CLASS_FILEREADER));
PostfixList(keyList, NewSymbol(NAME));
PostfixList(keyList, NewString(format));
fileReaders = SearchDatabase(keyList);
if (fileReaders)
{
runner = LISTOF(fileReaders);
if (runner)
{
if (runner -> next)
{
ReportError("FindFormatReader", "Duplicate file reader");
}
return runner -> thing;
}
}
return NULLOBJ;
}
ObjPtr FindExtensionReader(ext)
char *ext;
/*Finds the file reader for extension ext*/
{
ObjPtr keyList;
ThingListPtr runner;
ObjPtr fileReaders;
/*Search database for datasets*/
keyList = NewList();
PostfixList(keyList, NewSymbol(CLASSID));
PostfixList(keyList, NewInt(CLASS_FILEREADER));
PostfixList(keyList, NewSymbol(EXTENSION));
PostfixList(keyList, NewString(ext));
fileReaders = SearchDatabase(keyList);
if (fileReaders)
{
runner = LISTOF(fileReaders);
if (runner)
{
if (runner -> next)
{
ReportError("FindExtensionReader", "Duplicate file reader");
}
return runner -> thing;
}
}
return NULLOBJ;
}
ObjPtr ReadFormattedFile(reader, name)
char *name;
ObjPtr reader;
/*Reads file name with reader.*/
{
FuncTyp method;
ObjPtr retVal;
if (!reader) return ObjFalse;
timedDatasets = GetPredicate(reader, TIMEDDATASETS) ? true : false;
method = GetMethodSurely("ReadFormattedFile", reader, READALL);
if (!method)
{
return ObjFalse;
}
/*Read the file*/
curFileName = name;
retVal = (*method)(reader, name);
curFileName = 0;
return retVal;
}
void ReadFile(name, format)
char *name;
char *format;
/*Reads a file from the current directory with name and format format.
If format points to a null string, looks for a format*/
{
ObjPtr reader;
ObjPtr allFileReaders;
if (format[0])
{
/*There's a format specified*/
reader = FindFormatReader(format);
if (!reader)
{
WinInfoPtr errWindow;
char err[200];
sprintf(err, "There is no format named %s\n", format);
errWindow = AlertUser(UIERRORALERT, (WinInfoPtr) 0, err, 0, 0, "OK");
SetVar((ObjPtr) errWindow, HELPSTRING,
NewString("You have specified on the SciAn command line a format \
that does not exist, or you left out the format string entirely. This has \
prevented SciAn from reading the file. You can read it now using the \
file browser, which you can show by choosing New File Browser from the File menu."));
SetVar((ObjPtr) errWindow, INHIBITLOGGING, ObjTrue);
return;
}
}
else
{
char *s;
/*Must find a format*/
s = name;
while (*s) ++s;
while (s >= name && *s != '.') --s;
if (*s == '.')
{
reader = FindExtensionReader(s + 1);
}
else
{
reader = NULLOBJ;
}
if (!reader)
{
WinInfoPtr errWindow;
char err[200];
sprintf(err, "There is no default format for file %s\n", name);
errWindow = AlertUser(UIERRORALERT, (WinInfoPtr) 0, err, 0, 0, "yeah?");
SetVar((ObjPtr) errWindow, HELPSTRING,
NewString("You have specified on the SciAn command line a file \
whose format cannot be deduced from the name. This has \
prevented SciAn from reading the file. You can read it now using the \
file browser, which you can show by choosing New File Browser from the File menu."));
SetVar((ObjPtr) errWindow, INHIBITLOGGING, ObjTrue);
return;
}
}
ReadFormattedFile(reader, name);
}
static ObjPtr HideFileReadersWindow(window)
ObjPtr window;
/*Hide a file readers window, just return OK*/
{
return ObjTrue;
}
WinInfoPtr NewFileReadersWindow(void)
/*Create a new file readers window*/
{
WinInfoPtr objWin;
ThingListPtr runner;
ObjPtr panel, contents, corral, button;
int bw;
int l, r, b, t;
ObjPtr allFileReaders;
ObjPtr keyList;
/*Create the window*/
objWin = NewObjWindow(NULLOBJ, "File Readers", WINDBUF, FRWINWIDTH, FRWINHEIGHT, SCRWIDTH, SCRHEIGHT);
/*Set a null but successful HIDE routine*/
SetMethod((ObjPtr) objWin, HIDE, HideFileReadersWindow);
/*Add a help string*/
SetVar((ObjPtr) objWin, HELPSTRING,
NewString("This window shows all the file readers in Scian."));
/*Put in a panel*/
panel = NewPanel(greyPanelClass, 0, FRWINWIDTH, 0, FRWINHEIGHT);
SetVar(panel, STICKINESS, NewInt(STICKYLEFT + STICKYRIGHT +
STICKYBOTTOM + STICKYTOP));
contents = GetVar((ObjPtr) objWin, CONTENTS);
PrefixList(contents, panel);
SetVar(panel, PARENT, (ObjPtr) objWin);
l = 0; r = FRWINWIDTH; b = 0; t = FRWINHEIGHT;
/*Put in buttons and an icon corral*/
contents = GetListVar("NewDatasetsWindow", panel, CONTENTS);
if (!contents)
{
return 0;
}
/*Make an icon corral*/
corral = NewIconCorral(NULLOBJ, l + MINORBORDER, r - MINORBORDER, b + 2 * MINORBORDER + BUTTONHEIGHT, t - MINORBORDER, BARRIGHT + BARBOTTOM);
SetVar(corral, STICKINESS, NewInt(STICKYLEFT + STICKYRIGHT +
STICKYBOTTOM + STICKYTOP));
SetVar(corral, TOPDOWN, ObjTrue);
SetVar(corral, NAME, NewString("File Readers Corral"));
SetVar(corral, HELPSTRING,
NewString("This corral contains icons for all the file readers in \
SciAn. You can show the controls of any file readers by selecting \
the icons and pressing the Show Controls button at the bottom of the window."));
PrefixList(contents, corral);
SetVar(corral, PARENT, panel);
SetVar((ObjPtr) objWin, CORRAL, corral);
l += MINORBORDER;
r -= MINORBORDER;
b += MINORBORDER;
t = b + BUTTONHEIGHT;
bw = (r - l - MINORBORDER) / 2;
/*Make a show info button*/
button = NewFunctionButton(objWin,
l, l + bw,
b, b + BUTTONHEIGHT, OF_SHOW_CONTROLS);
if (button)
{
SetVar(button, PARENT, panel);
SetVar(button, STICKINESS, NewInt(STICKYBOTTOM + STICKYLEFT + FLOATINGRIGHT));
PrefixList(contents, button);
}
/*Drop all the file readers into the window*/
keyList = NewList();
PostfixList(keyList, NewSymbol(CLASSID));
PostfixList(keyList, NewInt(CLASS_FILEREADER));
allFileReaders = SearchDatabase(keyList);
runner = LISTOF(SortListByStringVar(allFileReaders, NAME, true));
while(runner)
{
ObjPtr icon, name;
FuncTyp method;
icon = GetVar(runner -> thing, DEFAULTICON);
if (!icon)
{
icon = NewIcon(0, 0, ICONFILEREADER, "?");
}
else
{
icon = NewObject(icon, 0L);
}
name = GetVar(runner -> thing, NAME);
SetVar(icon, NAME, name);
SetVar(icon, ICONLOC, NULLOBJ);
SetVar(icon, REPOBJ, runner -> thing);
SetVar(icon, CORRAL, corral);
DropIconSeriesInCorral(corral, icon);
runner = runner -> next;
}
return objWin;
}
WinInfoPtr FileReadersWindow()
/*Returns or creates a file readers window*/
{
WinInfoPtr retVal;
retVal = GetWinFromTitle("File Readers");
if (!retVal)
{
retVal = NewFileReadersWindow();
}
return retVal;
}
static ObjPtr ReadOldFile(fileReader, name)
ObjPtr fileReader;
char *name;
/*Method to read a file using an old style file reader*/
{
FuncTyp method;
method = GetMethodSurely("ReadOldFile", fileReader, OLDREAD);
if (!method)
{
return ObjFalse;
}
(*method)(name);
return ObjTrue;
}
static ObjPtr ShowFileReaderControls(fileReader, windowName)
ObjPtr fileReader;
char *windowName;
/*Makes a new control window to control a file reader*/
{
WinInfoPtr controlWindow;
ObjPtr var;
ObjPtr panel;
ObjPtr contents;
WinInfoPtr dialogExists;
dialogExists = DialogExists((WinInfoPtr) fileReader, NewString("Controls"));
if (!dialogExists)
{
ObjPtr textBox, checkBox;
int contentsRight, contentsTop;
int left, right, bottom, top, mid;
int l, r, b, t;
ObjPtr list;
ThingListPtr runner;
FuncTyp method;
ObjPtr button;
/*Make the panel*/
panel = NewPanel(greyPanelClass, 0, 10, 0, 10);
if (!panel)
{
return ObjFalse;
}
contents = GetVar(panel, CONTENTS);
SetVar(contents, PARENT, panel);
/*Give file reader chance to add controls*/
method = GetMethod(fileReader, ADDCONTROLS);
if (method)
{
(*method)(fileReader, contents);
}
/*Calculate the bounds*/
contentsRight = contentsTop = 0;
runner = LISTOF(contents);
while (runner)
{
Get2DIntBounds(runner -> thing, &left, &right, &bottom, &top);
contentsRight = MAX(contentsRight, right);
contentsTop = MAX(contentsTop, top);
runner = runner -> next;
}
bottom = contentsTop + MAJORBORDER;
left = MAJORBORDER;
/*Create the check box for time dependency*/
top = bottom + CHECKBOXHEIGHT;
right = left + FRTIMEWIDTH;
checkBox = NewCheckBox(left, right, bottom, top,
"Use \"field@time\" format for time samples",
GetPredicate(fileReader, TIMEDDATASETS));
SetVar(checkBox, PARENT, panel);
PrefixList(contents, checkBox);
bottom = top + MINORBORDER;
if (!GetVar(fileReader, TIMEDDATASETS))
{
SetVar(fileReader, TIMEDDATASETS, ObjTrue);
}
AssocDirectControlWithVar(checkBox, fileReader, TIMEDDATASETS);
SetVar(checkBox, HELPSTRING,
NewString("If this check box is checked, any dataset name of the \
form \"name@time\" will be interpreted as a sample of field name at time time. \
The time string can be a single number, in which case it will be interpreted \
as seconds or timesteps, or it can be a string in the form hh:mm:ss, or clock \
format. For example, a dataset named \"Z@12:04:00\" is a sample for the file \
Z taken at 12:04:00. This feature provides an easy way to get time-dependent \
data from file formats that are only suited to static data.\n\
\n\
If this check box is not checked, dataset names will be used as they appear, \
even if they contain the \"@\" character. Any inherent capability of the file \
reader to read time-dependent data will continue to work, of course."));
/*Create the default file extension box*/
top = bottom + EDITBOXHEIGHT;
right = left + FREXTTITLEWIDTH;
mid = (bottom + top) / 2;
textBox = NewTextBox(left, right,
mid - TEXTBOXHEIGHT / 2 + EDITBOXDOWN,
mid + TEXTBOXHEIGHT / 2 + EDITBOXDOWN,
0, "Default Extension Text", "Default Extension:");
SetVar(textBox, PARENT, panel);
PrefixList(contents, textBox);
/*Create the editable box*/
left = right + MAJORBORDER;
right = left + FREXTTEXTWIDTH;
var = GetVar(fileReader, EXTENSION);
if (!var)
{
var = NewString("");
SetVar(fileReader, EXTENSION, var);
}
textBox = NewTextBox(left, right,
mid - EDITBOXHEIGHT / 2,
mid + EDITBOXHEIGHT / 2,
EDITABLE + WITH_PIT + ONE_LINE, "Default Extension", var ? GetString(var) : "");
SetVar(textBox, PARENT, panel);
PrefixList(contents, textBox);
AssocDirectControlWithVar(textBox, fileReader, EXTENSION);
SetVar(textBox, HELPSTRING,
NewString("This text box shows the default file extension, or the portion \
of the file name after the last period, for files \
which use this format. If you change it to another extension, the new extension \
will be used when you open new file windows. The period is implied; do not include \
it in the new file extension."));
top += MINORBORDER;
/*Now put in a Save All Settings button*/
GetTemplateBounds(SettingsTemplate, "Save All Settings", &l, &r, &b, &t);
button = NewButton(l, r, top, top + (t - b), "Save All Settings");
top = top + (t - b) + MINORBORDER;
PrefixList(contents, button);
SetVar(button, PARENT, panel);
SetVar(button, REPOBJ, fileReader);
SetVar(button, HELPSTRING, NewString("Press this button to save all \
the settings for this file reader to a disk file.\n\n\
Settings are saved in a subdirectory named .scianSettings in your home directory. \
If that directory does not exist, it will automatically be created. Settings are \
saved in a file with a name of the form .filrdr where is the name of \
the object (with spaces converted to underscores)."));
ActivateButton(button, false);
list = AssembleSnapVars(fileReader);
if (list)
{
ThingListPtr runner;
runner = LISTOF(list);
while (runner)
{
NameTyp symbolID;
symbolID = GetSymbolID(runner -> thing);
DeclareIndirectDependency(button, APPEARANCE, REPOBJ, symbolID);
runner = runner -> next;
}
}
SetVar(button, APPEARANCE, ObjTrue);
SetMethod(button, APPEARANCE, MakeSaveButtonAppearance);
SetMethod(button, CHANGEDVALUE, SaveButtonChanged);
left = bottom = 0;
right = MAX(MAJORBORDER + FROVERALLWIDTH, contentsRight) + MAJORBORDER;
/*Finally, make the window around the panel*/
Set2DIntBounds(panel, left, right, bottom, top);
controlWindow = GetDialog((WinInfoPtr) fileReader, NewString("Controls"), windowName,
right, top, right, top, WINDBUF + WINFIXEDSIZE);
contents = GetVar((ObjPtr) controlWindow, CONTENTS);
PrefixList(contents, panel);
SetVar(panel, PARENT, (ObjPtr) controlWindow);
SetVar((ObjPtr) controlWindow, HELPSTRING,
NewString("This window shows controls for an individual file reader. \
Use Help in Context to get help on each control."));
}
else
{
PopWindow(dialogExists);
}
}
ObjPtr MakeReaderIconHelp(readerIcon, class)
ObjPtr readerIcon, class;
/*Makes help for a file reader icon*/
{
ObjPtr retVal, additionalHelp;
ObjPtr repObj;
ObjPtr var;
repObj = GetVar(readerIcon, REPOBJ);
if (!repObj)
{
return NULLOBJ;
}
var = GetVar(repObj, NAME);
sprintf(tempStr, "This icon represents the %s file reader. ",
var ? GetString(var) : "");
retVal = NewString(tempStr);
additionalHelp = GetVar(repObj, HELPSTRING);
if (!additionalHelp)
{
additionalHelp = NewString("For more information on this file reader, see the User Manual.");
}
retVal = ConcatStrings(retVal, additionalHelp);
SetVar(class, HELPSTRING, retVal);
return retVal;
}
static ObjPtr AddHDFControls(fileReader, panelContents)
ObjPtr fileReader, panelContents;
/*Adds controls appropriate to an HDF fileReader*/
{
ObjPtr titleBox, radio, checkBox, button;
int left, right, bottom, top;
ObjPtr var;
left = MAJORBORDER;
bottom = MAJORBORDER;
/*Do the read vectors check box*/
right = left + HDFVECTORWIDTH;
top = bottom + CHECKBOXHEIGHT;
checkBox = NewCheckBox(left, right, bottom, top, "Read 2- and 3-Vector Data",
GetPredicate(fileReader, READVECTOR));
PrefixList(panelContents, checkBox);
SetVar(checkBox, PARENT, panelContents);
SetVar(checkBox, HELPSTRING, NewString("This controls whether the HDF reader \
will automatically read vector data with 2 and 3 dimensions. If it is checked, \
whenever the first dimension of the dataset is either 2 or 3, it will \
be interpreted as containing 2- or 3-vector data. The first component is in the X \
direction, the second in the Y direction, and the third in the Z direction. \
If the box is not checked, the data will be read in as scalar data.\n\
\n\
For example, say you have a 3-D velocity field (40 by 50 by 60), where the velocity \
has three cartesian components. Store it as an HDF dataset 40 by 50 by 60 by 3. \
To store the velocity at (23, 37, 42), for example, put the X component at (1, 23, \
37, 42), the Y component at (2, 23, 37, 42), and the Z component at (3, 23, 37, 42). \
(This example uses FORTRAN conventions; a C example would go from 0 to 2 instead.)\n\
\n\
You may, if you like, use the last dimension instead of the first to specify the components \
of the vector."));
AssocDirectControlWithVar(checkBox, fileReader, READVECTOR);
bottom = top + MINORBORDER;
top = bottom + CHECKBOXHEIGHT;
checkBox = NewCheckBox(left, right, bottom, top, "Wrap Polar Coordinates",
GetPredicate(fileReader, WRAPPOLAR));
if (!GetVar(fileReader, WRAPPOLAR)) SetVar(fileReader, WRAPPOLAR, ObjTrue);
PrefixList(panelContents, checkBox);
SetVar(checkBox, PARENT, panelContents);
SetVar(checkBox, HELPSTRING, NewString("This controls whether the HDF reader \
will automatically wrap scales in polar or spherical coordinates when reading the dataset. When it \
is off, the scales along the axes will be used as is. We recommend that you always \
leave this check box on. \n\n\
This option is needed because most SciAn visualization objects expect the data to \
be defined over a topologically continuous grid and convert spherical coordinates \
to Mercator projection by default. If SciAn sees a sequence along a scale \
such as 178, 179, -180, -179, there is a discontinuity between 179 and -180. \
Even though this sequence is continuous in spherical coordinates, it is discontinuous \
in the Mercator projection. When the check box is on, a sequence like this will \
automatically be converted to 178, 179, 180, 181. When SciAn is enhanced to do spherical \
coordinates natively, both ways will work.\n\n\
In order for this check box to have any effect, the coordinate system in the DFSDsetdatastrs \
call must be spherical or polar and the units along each axis as specified in the \
DFSDsetdimstrs call must be degrees or units. Upper or lower case doesn't matter."));
AssocDirectControlWithVar(checkBox, fileReader, WRAPPOLAR);
bottom = top + MINORBORDER;
top = bottom + CHECKBOXHEIGHT;
checkBox = NewCheckBox(left, right, bottom, top, "Compress Data",
GetPredicate(fileReader, COMPRESSDATA));
if (!GetVar(fileReader, COMPRESSDATA)) SetVar(fileReader, COMPRESSDATA, ObjFalse);
PrefixList(panelContents, checkBox);
SetVar(checkBox, PARENT, panelContents);
SetVar(checkBox, HELPSTRING, NewString("This controls whether the HDF reader \
will store the data in compressed form. In compressed form, only one byte is used \
to store each numerical value. Compressed form takes only 1/4 as much space \
as normal, with an associated loss in precision. If the minimum and maximum \
are set within the file, they are used to determine the comrpession mapping. \
Otherwise, the minimum and maximum are determined from looking at the data."));
AssocDirectControlWithVar(checkBox, fileReader, COMPRESSDATA);
bottom = top + MINORBORDER;
top = bottom + CHECKBOXHEIGHT;
checkBox = NewCheckBox(left, right, bottom, top, "Read Axis Names",
GetPredicate(fileReader, BOTHERWITHAXES));
if (!GetVar(fileReader, BOTHERWITHAXES)) SetVar(fileReader, BOTHERWITHAXES, ObjFalse);
PrefixList(panelContents, checkBox);
SetVar(checkBox, PARENT, panelContents);
SetVar(checkBox, HELPSTRING, NewString("This controls whether the HDF reader \
will read the names of the axes from the HDF file. If it is checked, the names \
of the axes as set using the DFSDsetdimstrs will be read and interpreted as axis names. \
If it is not checked, the first axis will be assumed to be X, the second Y, and the third Z. \
If you don't use the DFSDsetdimstrs routine, you should keep this check box off, or else \
you will see warning messages and the occasional crash of SciAn due to a bug \
in NCSA's HDF routines.\n\n\
The axis names X, Y, and Z are recognized automatically. You can set up additional axes by putting a \
file named .scianAxes in your home directory. Upper or lower case doesn't matter. Here is an example .scianAxes file \
containing definitions for Longitude and Latitude:\n\
\n#Axes file containing supplemental axes to SciAn\n\
Longitude = X\n\
Latitude = Y"));
AssocDirectControlWithVar(checkBox, fileReader, BOTHERWITHAXES);
bottom = top + MINORBORDER;
/*Do the title box around the radio group*/
right = left + 2 * MINORBORDER + HDFRADIOWIDTH;
top = bottom + 2 * MINORBORDER + 2 * CHECKBOXSPACING + 3 * CHECKBOXHEIGHT + TITLEBOXTOP;
titleBox = NewTitleBox(left, right, bottom, top, "Handle Out-of-Bounds Data");
PrefixList(panelContents, titleBox);
SetVar(titleBox, PARENT, panelContents);
/*Make the radio buttons*/
radio = NewRadioButtonGroup("Out of Bounds Radio");
SetVar(radio, PARENT, panelContents);
SetVar(radio, REPOBJ, fileReader);
PrefixList(panelContents, radio);
left += MINORBORDER;
right -= MINORBORDER;
top -= TITLEBOXTOP + MINORBORDER;
bottom = top - CHECKBOXHEIGHT;
button = NewRadioButton(left, right, bottom, top, "Treat as missing");
AddRadioButton(radio, button);
SetVar(button, HELPSTRING, NewString("When this button is down, values in \
the dataset which are out of bounds are treated as missing data."));
top = bottom - CHECKBOXSPACING;
bottom = top - CHECKBOXHEIGHT;
button = NewRadioButton(left, right, bottom, top, "Clip to bounds");
AddRadioButton(radio, button);
SetVar(button, HELPSTRING, NewString("When this button is down, values in \
the dataset which are out of bounds are clipped to the bounds."));
top = bottom - CHECKBOXSPACING;
bottom = top - CHECKBOXHEIGHT;
button = NewRadioButton(left, right, bottom, top, "Use as is");
AddRadioButton(radio, button);
SetVar(button, HELPSTRING, NewString("When this button is down, values in \
the dataset which are out of bounds are used as ordinary data values."));
var = GetIntVar("AddHDFControls", fileReader, OUTOFBOUNDS);
if (!var) SetVar(fileReader, OUTOFBOUNDS, NewInt(0));
AssocDirectControlWithVar(radio, fileReader, OUTOFBOUNDS);
SetVar(radio, HELPSTRING, NewString("This radio button group controls how values \
in the dataset which are out of bounds are treated. A data value is considered \
out of bounds if it is outside the range given by the minimum and maximum, \
as specified by the DFSDSetMinMax HDF call. If no minimum and maximum have \
been set, this radio button group has no effect."));
return ObjTrue;
}
#ifdef PROTO
void ReadObjectControls(ObjPtr object, ObjPtr directory, Bool cpNeeded)
#else
void ReadObjectControls(object, directory, cpNeeded)
ObjPtr object, directory;
Bool cpNeeded;
#endif
/*Sets an object's controls from a possibly existing file within directory
cpNeeded true iff a control panel is needed*/
{
char oldDir[401];
ObjPtr name;
char *runner;
ObjPtr extension;
FILE *temp;
getcwd(oldDir, 400);
if (directory)
{
chdir(GetString(directory));
}
name = GetVar(object, NAME);
if (!name)
{
return;
}
strcpy(tempStr, GetString(name));
/*Modify file name*/
runner = tempStr;
while (*runner)
{
if (*runner == ' ') *runner = '_';
++runner;
}
extension = GetVar(object, SAVEEXTENSION);
if (extension)
{
strcat(runner, ".");
strcat(runner, GetString(extension));
}
SetVar(object, DIRECTORY, directory);
if (temp = fopen(tempStr, "r"))
{
/*Read the file into the object's control panel*/
Bool oldShowControlPanels;
Bool oldSettingUp;
ObjPtr controlWindow;
WinInfoPtr oldScriptWindow;
WinInfoPtr oldSelWinInfo;
fclose(temp);
oldSelWinInfo = selWinInfo;
InhibitLogging(true);
if (cpNeeded)
{
oldShowControlPanels = showControlPanels;
oldScriptWindow = scriptWindow;
showControlPanels = false;
/*Show and create the control panel*/
scriptWindow = (WinInfoPtr) NewControlWindow(object);
}
oldSettingUp = settingUp;
settingUp = true;
BeginScript(tempStr);
while(ReadScriptLine());
settingUp = oldSettingUp;
if (cpNeeded)
{
showControlPanels = oldShowControlPanels;
scriptWindow = oldScriptWindow;
selWinInfo = oldSelWinInfo;
}
InhibitLogging(false);
}
chdir(oldDir);
}
void SaveObjectControls(object, directory)
ObjPtr object, directory;
/*Sets an object's controls from a possibly existing file within directory,
or in the home directory if directory is null*/
{
char oldDir[401], fileName[401];
ObjPtr name;
char *runner;
ObjPtr extension;
getcwd(oldDir, 400);
if (directory)
{
chdir(GetString(directory));
}
name = GetVar(object, NAME);
if (!name)
{
return;
}
strcpy(fileName, GetString(name));
/*Modify file name*/
runner = fileName;
while (*runner)
{
if (*runner == ' ') *runner = '_';
++runner;
}
extension = GetVar(object, SAVEEXTENSION);
if (extension)
{
strcat(runner, ".");
strcat(runner, GetString(extension));
}
InhibitLogging(true);
if (OpenLogFile(fileName))
{
/*Read the file into the object's control panel*/
FuncTyp method;
Bool oldShowControlPanels, oldSettingUp;
LogNoWindow("# Saved settings for ");
LogNoWindow(GetString(name));
LogNoWindow("\n");
oldSettingUp = settingUp;
oldShowControlPanels = showControlPanels;
settingUp = true;
showControlPanels = false;
method = GetMethod(object, SAVEALLCONTROLS);
if (method)
{
(*method)(object);
}
showControlPanels = oldShowControlPanels;
settingUp = oldSettingUp;
CloseLogFile();
}
InhibitLogging(false);
chdir(oldDir);
}
void DoSaveObject()
/*Saves the object whose controls are on the top window*/
{
if (selWinInfo)
{
ObjPtr repObj;
repObj = GetVar((ObjPtr) selWinInfo, REPOBJ);
if (repObj)
{
Log("save controls");
DeferMessage(repObj, SAVECPANEL);
}
}
}
/*Fixed format data reading helpers*/
#ifdef PROTO
Bool GetFString(char *d, char *s, int beg, int end)
#else
Bool GetFString(d, s, beg, end)
char *d;
char *s;
int beg;
int end;
#endif
/*Puts into d a string from beg to end in s. Beg and end start at 1, not 0*/
{
int k;
for (k = beg - 1; k <= end - 1; ++k)
{
*d++ = s[k];
}
*d = 0;
return true;
}
#ifdef PROTO
Bool GetFLong(long *d, char *s, int beg, int end)
#else
Bool GetFLong(d, s, beg, end)
long *d;
char *s;
int beg;
int end;
#endif
{
int k;
char buffer[256];
for (k = 0; k <= end - beg; ++k)
{
buffer[k] = s[beg + k - 1];
}
buffer[k] = '\0';
return sscanf(buffer, "%ld", d) ? true : false;
}
#ifdef PROTO
Bool GetFReal(real *d, char *s, int beg, int end)
#else
Bool GetFReal(d, s, beg, end)
real *d;
char *s;
int beg;
int end;
#endif
{
int k;
char buffer[256];
for (k = 0; k <= end - beg; ++k)
{
buffer[k] = s[beg + k - 1];
}
buffer[k] = '\0';
return sscanf(buffer, "%g", d) ? true : false;
}
TextBuf *NewTextBuf()
/*Returns a new text buffer*/
{
TextBuf *retVal;
retVal = (TextBuf *) Alloc(sizeof(TextBuf));
retVal -> builtTextLength = ALLOCTEXTBUF;
retVal -> builtText = Alloc(ALLOCTEXTBUF);
retVal -> builtText[0] = 0;
retVal -> lastChar = 0;
return retVal;
}
#ifdef PROTO
void AppendTextChar(TextBuf *buf, char c)
#else
void AppendTextChar(TextBuf *buf, char c)
TextBuf *buf;
char c;
#endif
/*Appends a single char to buf*/
{
buf -> builtText[buf -> lastChar] = c;
++(buf -> lastChar);
if (buf -> lastChar >= buf -> builtTextLength)
{
buf -> builtTextLength += ALLOCTEXTBUF;
buf -> builtText = Realloc(buf -> builtText, buf -> builtTextLength);
}
buf -> builtText[buf -> lastChar] = 0;
}
#ifdef PROTO
void SmartAppendTextLine(TextBuf *buf, char *line)
#else
void SmartAppendTextLine(buf, line)
TextBuf *buf;
char *line;
#endif
/*Appends line to the text in buf. Maybe should do some thinking to justify
the word "smart." Basically, line is supposed to be a line formatted for
a teletype, and this is supposed to figure out some way to deal with the
problem of formatting it into a text box.*/
{
int k;
/*Find the end of the string and trim it down to the last non-space char*/
k = strlen(line);
if (k) --k;
while (k && isspace(line[k])) --k;
if (!k)
{
/*Just an end of line*/
}
else
{
int i;
/*There is some text.*/
for (i = 0; i <= k; ++i)
{
AppendTextChar(buf, line[i]);
}
}
AppendTextChar(buf, '\n');
}
#ifdef PROTO
char *GetBuiltText(TextBuf *buf)
#else
char *GetBuiltText(buf)
TextBuf *buf;
#endif
/*Gets the built text from the buffer*/
{
return buf -> builtText;
}
#ifdef PROTO
void DisposeTextBuf(TextBuf *buf)
#else
void DisposeTextBuf(buf)
TextBuf *buf;
#endif
{
if (!buf) return;
Free(buf -> builtText);
Free(buf);
}
#define ALLOCATOMS 100 /*Number of atoms to allocate at one time*/
#ifdef PROTO
void Titleize(char *s)
#else
void Titleize(s)
char *s;
#endif
/*Capitalizes a character string, sort of like a title*/
{
Bool readyForCapital = true;
while (*s)
{
if (isalpha(*s))
{
if (readyForCapital)
{
*s = toupper(*s);
readyForCapital = false;
}
else
{
*s = tolower(*s);
}
}
else if (!isdigit(*s))
{
readyForCapital = true;
}
++s;
}
}
#ifndef RELEASE
static ObjPtr FiddleData(dataset)
ObjPtr dataset;
/*Fiddles the data in the dataset*/
{
SetVar(dataset, DATA, GetVar(dataset, DATA));
DeferMessage(dataset, MARKTIME);
DatasetChanged(dataset);
return ObjTrue;
}
static ObjPtr DiddleData(dataset)
ObjPtr dataset;
/*Diddles the data in the dataset based on shared memory*/
{
real *address;
ObjPtr x, y, z;
real *xel, *yel, *zel;
int k;
ObjPtr var;
long dim;
var = GetIntVar("DiddleData", dataset, SHMADDRESS);
if (!var) return ObjFalse;
address = (real *) GetInt(var);
var = GetRealVar("DiddleData", dataset, SHMSEMAPHORE);
if (!var) return ObjFalse;
DeferMessage(dataset, MARKTIME);
if (GetReal(var) == *address)
{
/*Hasn't changed*/
return ObjTrue;
}
SetVar(dataset, SHMSEMAPHORE, NewReal(*address));
var = GetVar(dataset, DATA);
if (!var)
{
return ObjFalse;
}
x = ((ObjPtr *) ELEMENTS(var))[0];
y = ((ObjPtr *) ELEMENTS(var))[1];
z = ((ObjPtr *) ELEMENTS(var))[2];
xel = ELEMENTS(x);
yel = ELEMENTS(y);
zel = ELEMENTS(z);
dim = DIMS(x)[0];
for (k = 0; k < dim; ++k)
{
xel[k] = address[1 + k * 3];
yel[k] = address[2 + k * 3];
zel[k] = address[3 + k * 3];
}
SetVar(dataset, DATA, var);
DatasetChanged(dataset);
return ObjTrue;
}
static ObjPtr DiddleTimestep(dataset)
ObjPtr dataset;
/*Diddles the timestep in the dataset based on shared memory*/
{
real *address;
ObjPtr x, y, z;
real *xel, *yel, *zel;
int k;
ObjPtr var;
long dim;
var = GetIntVar("DiddleData", dataset, SHMADDRESS);
if (!var) return ObjFalse;
address = (real *) GetInt(var);
var = GetRealVar("DiddleData", dataset, SHMSEMAPHORE);
if (!var) return ObjFalse;
DeferMessage(dataset, MARKTIME);
if (GetReal(var) == *address)
{
/*Hasn't changed*/
return ObjTrue;
}
SetVar(dataset, SHMSEMAPHORE, NewReal(*address));
var = GetVar(dataset, DATA);
if (!var)
{
return ObjFalse;
}
((real *) ELEMENTS(var))[0] = *address;
SetVar(dataset, DATA, var);
DatasetChanged(dataset);
return ObjTrue;
}
static ObjPtr MakeWigglyMoleculeCurData(dataset)
ObjPtr dataset;
/*Makes the curdata of a wiggly molecule*/
{
register ObjPtr oldData, newData;
register ObjPtr *components;
register real *ox, *oy, *oz, *nx, *ny, *nz;
register long k;
register long dim;
register real r;
long temp;
/*Get the old data*/
oldData = GetVar(dataset, DATA);
components = ELEMENTS(oldData);
dim = DIMS(components[0])[0];
ox = ELEMENTS(components[0]);
oy = ELEMENTS(components[1]);
oz = ELEMENTS(components[2]);
temp = 3;
newData = NewArray(AT_OBJECT, 1, &temp);
components = ELEMENTS(newData);
components[0] = NewRealArray(1, dim);
components[1] = NewRealArray(1, dim);
components[2] = NewRealArray(1, dim);
nx = ELEMENTS(components[0]);
ny = ELEMENTS(components[1]);
nz = ELEMENTS(components[2]);
for (k = 0; k < dim; ++k)
{
r = ((real) (random() % 100 - 50)) / 150.0;
nx[k] = ox[k] + r;
r = ((real) (random() % 100 - 50)) / 150.0;
ny[k] = oy[k] + r;
r = ((real) (random() % 100 - 50)) / 150.0;
nz[k] = oz[k] + r;
}
SetVar(dataset, DATA, newData);
return ObjTrue;
}
#endif
/***************************begin getAtomNum***********************************/
/*
Tzong-Yow Hwu, 6/10/1993
7/09/1993, revised by distinguishing among
Regular Atom: atoms belong to ATOM record
this is further distinguished between
1. amino acid residue atom,
2. nucleic acid residue atom, and
3. others.
Nonstandard Atom: atoms of HETATM record
parseAtomName returns the atomic number of atomName and store the
remoteness indicator in remoteness. It returns 0 if unable to
recognize the atomName ,or for atoms of ambiguity (whose symbol is A).
It stores '\0' in remoteness if unable to recognize the
remoteness indicator or if no remoteness indicator found in atomName.
Legitimate remoteness indicators are A, B, G, D, E, Z, H, T, U, F, X, Q,
and W.
atomName string format:
ATOM record
1. amino acid residue atom:
standard: the atomName string is four character long.
Characters 1-2 are for chemical symbols(right justified),
character 3 is a remoteness indicator(alphabetic), and
character 4 is a branch designator(numeric.)
exceptions: 1. Atoms for which some ambiguity exists in the
crystallographic results are designated A.
2. The extra oxygen atom of the carboxy terminal amino acid is
designated OXT.
3. When more than one hydrogen atoms are bounded to a single
non-hydrogen atom, branch designator is required to be the
first char. Not necessarily followed, some heurestic rules
used.
2. nucleic acid residue atom:
1. generally no remoteness indicator. (?)
2. char 1-2 are for atom symbols.
char 3-4 for atom position ,and '*' or '''
3. OXT for the oxygen atom at free 5' phosphate terminus.
3. other residue atom:
parse it according to nucleic acid residue atom.
HETATM record
1. Generally, the 2nd char is used for atom symbol.
The 1st char is sometimes used for other purpose such as AO9, PC5,...
etc. in COENZYME A, and AP, AN9, AC2*, ...etc. in NAD group. But
watch out for atom symbols of two char long, such as FE, CA, ZN,
PB, MG, CL, BR, CU, NA...etc.
2. If there is a remoteness indicator, it would appear as the 3rd char.
But the 3rd char is not always a remoteness indicator. Take the 3rd
char as a remoteness indicator if it matches one the remoteness symbols.
The " NA4" and " NA2" of the METHOTREXATE group may be treated wrong
under the assumption. (?)
#include "ScianPeriodicTable.h"
*/
#ifdef PROTO
int getAtomNum(const char *record, char *remoteness);
#else
int getAtomNum();
#endif
#ifdef PROTO
int getAtomNum(const char *record, char *remoteness)
#else
int getAtomNum(record, remoteness)
const char *record;
char *remoteness;
#endif
{
char *aminoAcidResidueAbbre[] =
{
"ABU", "ACD", "ALA", "ALB", "ALI", "ARG", "ARO", "ASN", "ASP", "ASX",
"BAS",
"CYS",
"GLN", "GLU", "GLX", "GLY",
"HIS", "HYP",
"ILE",
"LEU", "LYS",
"MET",
"PCA", "PHE", "PRO",
"SAR", "SER",
"THR", "TRP", "TYR",
"VAL"
};
const int numAminoAbbre = 31;
char atomType[7], atomName[5], residueName[4], atomSym[3];
char s[132];
int atomicNumber;
int i;
char remotenessSyms[] =
{'A', 'B', 'G', 'D', 'E', 'Z', 'H', 'T', 'U', 'F', 'X', 'Q', 'W'};
const int numRemoteSyms = 13;
#undef VERIFYREMOTE
#define VERIFYREMOTE(a)\
{\
int j;\
for (j = 0;\
j < numRemoteSyms && a != remotenessSyms[j];\
j++) ;\
if (j == numRemoteSyms)\
{\
a = '\0';\
}\
}
/* get atomType, atomName, and residueName from record */
for (i = 0; i < 6; i++)
{
atomType[i] = toupper(record[i]);
}
atomType[i] = '\0';
for (i = 0; i < 4; i++)
{
atomName[i] = toupper(record[i + 12]);
}
atomName[i] = '\0';
for (i = 0; i < 3; i++)
{
residueName[i] = toupper(record[i + 17]);
}
residueName[i] = '\0';
if (0 == strcmp(atomType, "ATOM "))
{
/* ATOM record */
int notAminoAcidResidue = 1;
/* see if it's amino acid residue */
for (i = 0; notAminoAcidResidue && i < numAminoAbbre; i++)
{
notAminoAcidResidue = strcmp(residueName, aminoAcidResidueAbbre[i]);
}
if (!notAminoAcidResidue)
{
/* parse as amino acid residue atoms */
if (atomName[0] == ' ')
{
/* then the 2nd char must be an atomic symbol */
if (atomName[1] == 'A')
{
/* it's an atom of ambiguity, assign atomic number 0 */
atomicNumber = 0;
}
else
{
atomSym[0] = atomName[1];
atomSym[1] = '\0';
atomicNumber = AtomNameToNumber(atomSym);
}
if (atomicNumber)
{
/* see if the 3rd char is a legitimate remoteness indicator */
*remoteness = atomName[2];
VERIFYREMOTE(*remoteness);
}
else
{
*remoteness = '\0';
}
}
else if (isdigit(atomName[0]))
{
/* maybe the case of "1HB ", ...etc. */
if (atomName[1] == 'H')
{
atomicNumber = AtomNameToNumber("H");
/* see if the 3rd char is a legitimate remoteness indicator */
*remoteness = atomName[2];
VERIFYREMOTE(*remoteness);
}
else
{
/* illegal atom name of ATOM amino acid group */
atomicNumber = 0;
*remoteness = '\0';
}
}
else if ( atomName[0] == 'H' && isdigit(atomName[2]) &&
isdigit(atomName[3]))
{
/* maybe the case such as "HD22", ...etc. if the 2nd
char is a legitimate remoteness indicator */
*remoteness = atomName[1];
VERIFYREMOTE(*remoteness);
if (*remoteness == '\0')
{
/* not legal */
atomicNumber = 0;
}
else
{
atomicNumber = AtomNameToNumber("H");
}
}
else
{
/* treat the first 2 chars as atomic symbol */
atomSym[0] = atomName[0];
atomSym[1] = atomName[1];
atomSym[2] = '\0';
if (!(atomicNumber = AtomNameToNumber(atomSym)))
{
*remoteness = '\0';
}
else
{
/* see if the 3rd char is a legitimate remoteness indicator */
*remoteness = atomName[2];
VERIFYREMOTE(*remoteness);
}
}
}
else
{
/* assume it's an atom of nucleic acid residue */
/* no remoteness indicator */
*remoteness = '\0';
if (atomName[0] = ' ')
{
atomSym[0] = atomName[1];
atomSym[1] = '\0';
}
else
{
atomSym[0] = atomName[0];
atomSym[1] = atomName[1];
atomSym[2] = '\0';
}
atomicNumber = AtomNameToNumber(atomSym);
}
}
else if (0 == strcmp(atomType, "HETATM"))
{
/* HETATM record */
if (atomName[0] == ' ' || isdigit(atomName[0]))
{
/* then the 2nd char must be an atomic symbol */
if (atomName[1] == 'A')
{
/* atom of ambiquity */
atomicNumber = 0;
}
else
{
atomSym[0] = atomName[1];
atomSym[1] = '\0';
atomicNumber = AtomNameToNumber(atomSym);
}
}
else if (atomName[2] == '*' || atomName[3] == '*')
{
/* atoms of ribose, the 2nd char must be an atomic symbol */
atomSym[0] = atomName[1];
atomSym[1] = '\0';
atomicNumber = AtomNameToNumber(atomSym);
}
else if (atomName[0] == 'A')
{
/* maybe 1. atoms of COENZYME A or NAD, or
2. atoms of two chars starting with A */
/* When the 2nd char is C, N, O, or P, it's most likely case 1 */
if (atomName[1] == 'C' || atomName[1] == 'N' ||
atomName[1] == 'O' || atomName[1] == 'P' )
{
atomSym[0] = atomName[1];
atomSym[1] = '\0';
}
else
{
atomSym[0] = atomName[0];
atomSym[1] = atomName[1];
atomSym[2] = '\0';
}
atomicNumber = AtomNameToNumber(atomSym);
}
else if (atomName[0] == 'P')
{
/* maybe 1. atoms of COENZYME A, or
2. atoms of two chars starting with P */
/* When the 2nd char is C, N, O, or S, it's most likely case 1 */
if (atomName[1] == 'C' || atomName[1] == 'N' ||
atomName[1] == 'O' || atomName[1] == 'S' )
{
atomSym[0] = atomName[1];
atomSym[1] = '\0';
}
else
{
atomSym[0] = atomName[0];
atomSym[1] = atomName[1];
atomSym[2] = '\0';
}
atomicNumber = AtomNameToNumber(atomSym);
}
else if (atomName[0] == 'N')
{
/* maybe 1. atoms of NAD, or
2. atoms of two chars starting with N */
/* When the 2nd char is C, N, O, or P, it's most likely case 1 */
if (atomName[1] == 'C' || atomName[1] == 'N' ||
atomName[1] == 'O' || atomName[1] == 'P' )
{
atomSym[0] = atomName[1];
atomSym[1] = '\0';
}
else
{
atomSym[0] = atomName[0];
atomSym[1] = atomName[1];
atomSym[2] = '\0';
}
atomicNumber = AtomNameToNumber(atomSym);
}
else
{
/* Otherwise, do the best it can to recognize the atom symbol */
atomSym[0] = atomName[0];
atomSym[1] = atomName[1];
atomSym[2] = '\0';
atomicNumber = AtomNameToNumber(atomSym);
if (!atomicNumber)
{
atomSym[0] = atomName[1];
atomSym[1] = '\0';
atomicNumber = AtomNameToNumber(atomSym);
}
}
if (atomicNumber)
{
/* see if the 3rd char is a legitimate remoteness indicator */
*remoteness = atomName[2];
VERIFYREMOTE(*remoteness);
}
else
{
*remoteness = '\0';
}
}
else
{
sprintf(s, "Funny atom type %c%c%c%c%c%c", record[0], record[1],
record[2], record[3], record[4], record[5]);
ReportError("getAtomNum", s);
*remoteness = '\0';
atomicNumber = 0;
return atomicNumber;
}
if (!atomicNumber)
{
sprintf(s, "Funny atom name %c%c%c%c", record[12], record[13],
record[14], record[15]);
ReportError("getAtomNum", s);
}
return atomicNumber;
}
/************* end getAtomNum******************************************************/
/*Atom structure for protein data bank file*/
typedef struct
{
real x;
real y;
real z;
real occupancy;
real tempFactor;
int atomicNumber;
long residue; /*Residue number*/
Bool isCarbonAlpha;
} Atom;
#define MAKEBOND(a, b) if (nBonds >= nBondsAllocated) \
{if (nBondsAllocated) {nBondsAllocated += 200; \
bondsBuffer = (TwoReals *) Realloc(bondsBuffer, \
sizeof(TwoReals) * nBondsAllocated);} \
else {nBondsAllocated = 200; \
bondsBuffer = (TwoReals *) Alloc(sizeof(TwoReals) * nBondsAllocated); }} \
bondsBuffer[nBonds][0] = (real) (a); \
bondsBuffer[nBonds][1] = (real) (b); \
++nBonds;
static ObjPtr ReadPDBFile(reader, fileName)
ObjPtr reader;
char *fileName;
/*Reads a Protein Data Bank file.*/
{
char typeBuf[10];
char noteBuf[100];
char compoundName[100];
Bool nameGiven = false;
FILE *inFile;
/*Go through the file*/
inFile = fopen(fileName, "r");
if (inFile)
{
/*The file exists. Read it*/
Atom *atoms;
long nAtoms;
long nAtomsAllocated;
long k;
real bounds[6];
ObjPtr notes;
TextBuf *remark;
TextBuf *source;
TextBuf *author;
TextBuf *journal;
TextBuf *expData;
int topResidue;
TwoReals *bondsBuffer = 0;
long nBondsAllocated = 0;
long nBonds = 0;
/*Make the text buffers*/
remark = 0;
source = 0;
author = 0;
journal = 0;
expData = 0;
/*Make bounds impossible*/
bounds[0] = bounds[2] = bounds[4] = PLUSINF;
bounds[1] = bounds[3] = bounds[5] = MINUSINF;
/*No atoms to start with*/
nAtomsAllocated = 0;
nAtoms = 0;
topResidue = 0;
while (fgets(tempStr, 100, inFile))
{
/*Determine the type*/
if (GetFString(typeBuf, tempStr, 1, 6))
{
/*First check for sanity*/
for (k = 0; k < 6; ++k)
{
if (!isalpha(typeBuf[k]) && !isspace(typeBuf[k]) && !isdigit(typeBuf[k]))
{
FileFormatError("ReadPDBFile", "This does not appear to be a valid protein data bank file");
break;
}
}
/*Now see what it is*/
if ((0 == strcmp2("ATOM ", typeBuf)) || (0 == strcmp2("HETATM", typeBuf)))
{
/*It's an atom*/
real x, y, z;
/*
char *t;
*/
char atomNameBuf[4];
char atomRemoteness;
int atomicNumber;
Bool isCarbonAlpha;
long atomSerialNumber;
/*Get the atom serial number*/
if (!GetFLong(&atomSerialNumber, tempStr, 7, 11))
{
FileFormatError("ReadPDBFile", "Missing atom serial number");
}
if (atomSerialNumber > nAtoms)
{
nAtoms = atomSerialNumber;
}
/*Make more atoms if need be*/
if (nAtoms >= nAtomsAllocated)
{
if (nAtomsAllocated)
{
nAtomsAllocated += ALLOCATOMS;
atoms = (Atom *) Realloc(atoms, sizeof(Atom) * nAtomsAllocated);
}
else
{
nAtomsAllocated = ALLOCATOMS;
atoms = (Atom *) Alloc(sizeof(Atom) * nAtomsAllocated);
}
for (k = nAtomsAllocated - ALLOCATOMS;
k < nAtomsAllocated; ++k)
{
atoms[k] . x = missingData;
atoms[k] . y = missingData;
atoms[k] . z = missingData;
atoms[k] . occupancy = missingData;
atoms[k] . tempFactor = missingData;
atoms[k] . atomicNumber = 0;
}
}
/*Get the coordinates*/
GetFReal(&(atoms[atomSerialNumber - 1] . x), tempStr, 31, 38);
GetFReal(&(atoms[atomSerialNumber - 1] . y), tempStr, 39, 46);
GetFReal(&(atoms[atomSerialNumber - 1] . z), tempStr, 47, 54);
if (atoms[atomSerialNumber - 1] . x < bounds[0]) bounds[0] = atoms[atomSerialNumber - 1] . x;
if (atoms[atomSerialNumber - 1] . x > bounds[1]) bounds[1] = atoms[atomSerialNumber - 1] . x;
if (atoms[atomSerialNumber - 1] . y < bounds[2]) bounds[2] = atoms[atomSerialNumber - 1] . y;
if (atoms[atomSerialNumber - 1] . y > bounds[3]) bounds[3] = atoms[atomSerialNumber - 1] . y;
if (atoms[atomSerialNumber - 1] . z < bounds[4]) bounds[4] = atoms[atomSerialNumber - 1] . z;
if (atoms[atomSerialNumber - 1] . z > bounds[5]) bounds[5] = atoms[atomSerialNumber - 1] . z;
/*Get the residue number*/
if (GetFLong(&(atoms[atomSerialNumber - 1] . residue), tempStr, 23,26))
{
if (atoms[atomSerialNumber - 1] . residue > topResidue)
{
topResidue = atoms[atomSerialNumber - 1] . residue;
}
if (atoms[atomSerialNumber - 1] . residue < topResidue)
{
/*DIKEO residue out of order, do something?*/
}
}
/*Try to figure out the name of the atom*/
atoms[atomSerialNumber - 1] . atomicNumber =
atomicNumber = getAtomNum(tempStr, &atomRemoteness);
isCarbonAlpha = false;
if (atomicNumber == 6)
{
/*It's carbon, may be an alpha carbon*/
if (atomRemoteness == 'A')
{
isCarbonAlpha = true;
}
}
atoms[atomSerialNumber - 1] . isCarbonAlpha = isCarbonAlpha;
/*Get the occupancy and temperature factor*/
if (!GetFReal(&(atoms[atomSerialNumber - 1] . occupancy), tempStr, 55, 60))
{
atoms[atomSerialNumber - 1] . occupancy = missingData;
}
if (!GetFReal(&(atoms[atomSerialNumber - 1] . tempFactor), tempStr, 61, 66))
{
atoms[atomSerialNumber - 1] . tempFactor = missingData;
}
}
else if (0 == strcmp2("REMARK", typeBuf))
{
if (!remark)
{
remark = NewTextBuf();
}
GetFString(noteBuf, tempStr, 12, 70);
SmartAppendTextLine(remark, noteBuf);
}
else if (0 == strcmp2("SOURCE", typeBuf))
{
if (!source)
{
source = NewTextBuf();
}
GetFString(noteBuf, tempStr, 11, 70);
SmartAppendTextLine(source, noteBuf);
}
else if (0 == strcmp2("AUTHOR", typeBuf))
{
if (!author)
{
author = NewTextBuf();
}
GetFString(noteBuf, tempStr, 11, 70);
SmartAppendTextLine(author, noteBuf);
}
else if (0 == strcmp2("EXPDTA", typeBuf))
{
if (!expData)
{
expData = NewTextBuf();
}
GetFString(noteBuf, tempStr, 11, 70);
SmartAppendTextLine(expData, noteBuf);
}
else if (0 == strcmp2("JRNL ", typeBuf))
{
if (!journal)
{
journal = NewTextBuf();
}
GetFString(noteBuf, tempStr, 11, 70);
SmartAppendTextLine(journal, noteBuf);
}
else if (0 == strcmp2("COMPND", typeBuf) && !nameGiven)
{
/*Extract a reasonable compound name from the line*/
GetFString(compoundName, tempStr, 11, 70);
nameGiven = true;
}
else if (0 == strcmp2("CONECT", typeBuf))
{
/*Connectivity*/
long thisAtom, otherAtom;
/*See what this atom is*/
if (!GetFLong(&thisAtom, tempStr, 7, 11))
{
continue;
}
/*Check all the bonds*/
if (GetFLong(&otherAtom, tempStr, 12, 16))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 17, 21))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 22, 26))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 27, 31))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 32, 36))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 37, 41))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 42, 46))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 47, 51))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 52, 56))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
if (GetFLong(&otherAtom, tempStr, 57, 61))
{
MAKEBOND(thisAtom - 1, otherAtom - 1);
}
}
}
}
fclose(inFile);
/*Make the notes*/
notes = NewList();
if (source)
{
PostfixList(notes, NewString("Source"));
PostfixList(notes, NewString(GetBuiltText(source)));
DisposeTextBuf(source);
}
if (expData)
{
PostfixList(notes, NewString("Exp. Data"));
PostfixList(notes, NewString(GetBuiltText(expData)));
DisposeTextBuf(expData);
}
if (author)
{
PostfixList(notes, NewString("Author"));
PostfixList(notes, NewString(GetBuiltText(author)));
DisposeTextBuf(author);
}
if (journal)
{
PostfixList(notes, NewString("Journal"));
PostfixList(notes, NewString(GetBuiltText(journal)));
DisposeTextBuf(journal);
}
if (remark)
{
PostfixList(notes, NewString("Remarks"));
PostfixList(notes, NewString(GetBuiltText(remark)));
DisposeTextBuf(remark);
}
/*Now maybe there are some atoms*/
if (nAtoms)
{
ObjPtr formVectors, dataForm, occupancyDataset, tempFactorDataset,
atomicNumberDataset, radiusDataset;
char *s, *t;
ObjPtr timeStepDataset;
ObjPtr residues;
long dims[2];
/*Make a reasonable-sounding compound name*/
if (!nameGiven)
{
strcpy(compoundName, GetString(MakeDatasetName(fileName)));
}
s = compoundName;
SKIPBLANKS(s);
if (!s)
{
strcpy(compoundName, GetString(MakeDatasetName(fileName)));
s = compoundName;
}
t = s;
/*Go until end of string or until first (*/
while (*t && *t != '(') ++t;
--t;
/*Back up until last non-blank char*/
while (isspace(*t)) --t;
/*Advance one and plop a 0*/
++t;
*t = 0;
/*Capitalize the compound name like a title*/
Titleize(s);
/*Make a vector dataset for the form*/
formVectors = NewDataset("form vectors", 1, &nAtoms, 3);
/*Make datasets*/
strcpy(t, " Occupancy");
occupancyDataset = NewDataset(s, 1, &nAtoms, 0);
SetVar(occupancyDataset, NOTES, notes);
strcpy(t, " Temp Factor");
tempFactorDataset = NewDataset(s, 1, &nAtoms, 0);
SetVar(tempFactorDataset, NOTES, notes);
strcpy(t, " Atoms");
atomicNumberDataset = NewDataset(s, 1, &nAtoms, 0);
SetVar(atomicNumberDataset, NOTES, notes);
SetVar(atomicNumberDataset, UNITSNAME, NewString("Atomic number"));
strcpy(t, " Radii");
radiusDataset = NewDataset(s, 1, &nAtoms, 0);
SetVar(radiusDataset, NOTES, notes);
/*Fill it all up*/
SetCurField(FORMFIELD, formVectors);
SetCurField(FIELD1, occupancyDataset);
SetCurField(FIELD2, tempFactorDataset);
SetCurField(FIELD3, atomicNumberDataset);
SetCurField(FIELD4, radiusDataset);
/*Fill fields with atoms.*/
for (k = 0; k < nAtoms; ++k)
{
PutFieldComponent(FORMFIELD, 0, &k, atoms[k] . x);
PutFieldComponent(FORMFIELD, 1, &k, atoms[k] . y);
PutFieldComponent(FORMFIELD, 2, &k, atoms[k] . z);
PutFieldComponent(FIELD1, 0, &k, atoms[k] . occupancy);
PutFieldComponent(FIELD2, 0, &k, atoms[k] . tempFactor);
PutFieldComponent(FIELD3, 0, &k,
atoms[k] . atomicNumber ? (real) atoms[k] . atomicNumber : missingData);
PutFieldComponent(FIELD4, 0, &k,
atoms[k] . atomicNumber ? atomInfo[atoms[k] . atomicNumber - 1] . radius : missingData);
}
if (!nBonds)
{
#ifndef RELEASE
/*Get bonds*/
strcpy(tempStr, fileName);
strcat(tempStr, ".con");
if (inFile = fopen(tempStr, "r"))
{
/*Read in bonds*/
int k, l;
nBonds = 0;
nBondsAllocated = 200;
bondsBuffer = (TwoReals *) Alloc(sizeof(TwoReals) * nBondsAllocated);
while (2 == fscanf(inFile, " %d %d", &k, &l))
{
if (nBonds >= nBondsAllocated)
{
nBondsAllocated += 200;
bondsBuffer = (TwoReals *) Realloc(bondsBuffer,
sizeof(TwoReals) * nBondsAllocated);
}
bondsBuffer[nBonds][0] = k - 1;
bondsBuffer[nBonds][1] = l - 1;
++nBonds;
}
}
else
#endif
{
/*Try to guess some bonds*/
printf("Guessing bonds\n");
nBonds = 0;
nBondsAllocated = 200;
bondsBuffer = (TwoReals *) Alloc(sizeof(TwoReals) * nBondsAllocated);
for (k = 0; k < nAtoms; ++k)
{
long l;
real d2, maxd;
maxd = SQUARE(2.0);
/*Only search atoms within 10*/
for (l = k + 1; l < MAX(nAtoms, k + 20); ++l)
{
d2 = SQUARE(atoms[k] . x - atoms[l] . x) +
SQUARE(atoms[k] . y - atoms[l] . y) +
SQUARE(atoms[k] . z - atoms[l] . z);
if (d2 < maxd)
{
if (nBonds >= nBondsAllocated)
{
nBondsAllocated += 200;
bondsBuffer = (TwoReals *) Realloc(bondsBuffer,
sizeof(TwoReals) * nBondsAllocated);
}
bondsBuffer[nBonds][0] = k;
bondsBuffer[nBonds][1] = l;
++nBonds;
}
}
}
}
}
dims[0] = nAtoms;
dims[1] = nBonds;
dataForm = NewUnstructuredDataForm("form", 1, dims, bounds, formVectors);
if (nBonds)
{
/*Set the bonds*/
for (k = 0; k < nBonds; ++k)
{
Put1Edge(dataForm, k,
(long) bondsBuffer[k][0], (long) bondsBuffer[k][1]);
}
}
SetDatasetForm(occupancyDataset, dataForm);
SetDatasetForm(tempFactorDataset, dataForm);
SetDatasetForm(atomicNumberDataset, dataForm);
SetDatasetForm(radiusDataset, dataForm);
/*Make carbon alpha edges*/
residues = NewRealArray(2, topResidue, 2L);
if (residues)
{
TwoReals *residueElements;
int curResidue;
int residueBeginning;
int i;
int alphaCarbon;
residueElements = (TwoReals *) ELEMENTS(residues);
curResidue = 1; /*In file, residues begin at 1*/
alphaCarbon = -1; /*No alpha carbon yet*/
residueBeginning = 0;
for (k = 0; k <= nAtoms; ++k)
{
if (atoms[k] . atomicNumber > 0)
{
if (k == nAtoms || atoms[k] . residue > curResidue)
{
/*Have to finish up this residue*/
residueElements[curResidue - 1][0] = (real) k - 1;
if (alphaCarbon >= 0)
{
/*Set this to the alpha carbon*/
}
else
{
/*Guess an alpha carbon, pick the first carbon*/
for (i = residueBeginning; i < k; ++i)
{
if (atoms[i] . atomicNumber == 6)
{
alphaCarbon = i;
break;
}
}
/*If there Still isn't an alpha carbon, pick middle*/
if (alphaCarbon < 0)
{
alphaCarbon = (i + k - 1) / 2;
}
}
residueElements[curResidue - 1][1] = (real) alphaCarbon;
alphaCarbon = -1;
residueBeginning = k;
if (k < nAtoms)
{
curResidue = atoms[k] . residue;
}
}
if (k < nAtoms)
{
/*Still more to go*/
if (atoms[k] . isCarbonAlpha)
{
alphaCarbon = k;
}
}
}
}
SetVar(dataForm, RESIDUES, residues);
}
SetVar(atomicNumberDataset, SIZEOBJ, radiusDataset);
SetVar(occupancyDataset, SIZEOBJ, radiusDataset);
SetVar(tempFactorDataset, SIZEOBJ, radiusDataset);
SetVar(radiusDataset, SIZEOBJ, radiusDataset);
RegisterNewDataset(atomicNumberDataset);
RegisterNewDataset(occupancyDataset);
RegisterNewDataset(tempFactorDataset);
RegisterNewDataset(radiusDataset);
#ifndef RELEASE
#ifdef WIGGLY
/*Make wiggly bits*/
SetMethod(formVectors, MARKTIME, FiddleData);
DeferMessage(formVectors, MARKTIME);
SetMethod(formVectors, DATA, MakeWigglyMoleculeCurData);
#else
#ifdef SHMDIDDLE
timeStepDataset = NewDataset("Duration (ps)", 0, &nAtoms, 0);
RegisterDataset(timeStepDataset);
/*Set up to diddle with shared memory*/
if (0 == sharedSegment)
{
/*Make a new piece of shared memory*/
key_t key;
long size;
size = 1 + nAtoms * 3;
key = ftok(fileName, 37);
if (key > 0)
{
/*Key OK, make an ID*/
int shmid;
printf("fileName = %s, key = %d, size = %d\n", fileName, key,
sizeof(real) * (1 + nAtoms * 3));
shmid = shmget(key, sizeof(real) * (1 + nAtoms * 3), IPC_CREAT);
if (shmid > 0)
{
/*Key ok, make shared memory segment*/
sharedSegment = (real *) shmat(shmid, (void *) 0, 0);
if (((long) sharedSegment) > 0)
{
/*It's OK, let's do it.*/
*sharedSegment = 0.0;
SetVar(formVectors, SHMSEMAPHORE, NewReal(0.0));
SetVar(formVectors, SHMADDRESS, NewInt((int) sharedSegment));
SetVar(timeStepDataset, SHMSEMAPHORE, NewReal(0.0));
SetVar(timeStepDataset, SHMADDRESS, NewInt((int) sharedSegment));
/*Make initial setup*/
sharedSegment[1] = (real) nAtoms;
for (k = 0; k < nAtoms; ++k)
{
sharedSegment[1 + k * 3] = atoms[k] . x;
sharedSegment[2 + k * 3] = atoms[k] . y;
sharedSegment[3 + k * 3] = atoms[k] . z;
}
SetMethod(formVectors, MARKTIME, DiddleData);
DeferMessage(formVectors, MARKTIME);
SetMethod(timeStepDataset, MARKTIME, DiddleTimestep);
DeferMessage(timeStepDataset, MARKTIME);
}
else
{
perror("can't attach segment");
}
}
else
{
printf("%d is a ", shmid);
perror("bad shared memory id");
}
}
else
{
perror("bad key");
}
}
#endif
#endif
#endif
}
SAFEFREE(bondsBuffer);
SAFEFREE(atoms);
return ObjTrue;
}
else
{
sprintf(tempStr, "File %s cannot be opened.", fileName);
FileFormatError("ReadPDBFile", tempStr);
return ObjFalse;
}
}
static ObjPtr ReadGTFile(fileName)
char *fileName;
/*Reads a geometry test file.*/
{
ObjPtr dataForm = NULLOBJ; /*Data form*/
ObjPtr dataFormData = NULLOBJ; /*Data form*/
ObjPtr normalsData = NULLOBJ; /*Normals*/
ObjPtr geometryDataset = NULLOBJ; /*Geometry dataset*/
ObjPtr scalarDatasets = NULLOBJ;
ThingListPtr runner;
real bounds[6];
long polygonVertices[100];
char keyword[100];
long size = 1;
long nVertices = 0;
long vertex;
long n;
long k;
float x, y, z;
FILE *inFile;
inFile = fopen(fileName, "r");
if (inFile)
{
scalarDatasets = NewList();
while (1 == fscanf(inFile, " %s", keyword))
{
if (0 == strcmp(keyword, "VERTICES"))
{
if (1 != fscanf(inFile, " %ld \n", &size))
{
FileFormatError("ReadGTFile", "# vertices expected");
return NULLOBJ;
}
/*Create a new dataform dataset*/
dataFormData = NewDataset("geometry vertices data", 1, &size, 3);
SetCurField(FIELD1, dataFormData);
bounds[0] = bounds[2] = bounds[4] = plusInf;
bounds[1] = bounds[3] = bounds[5] = minusInf;
for (k = 0; k < size; ++k)
{
if (3 != fscanf(inFile, "%g %g %g \n", &x, &y, &z))
{
FileFormatError("ReadGTFile", "Vertices expected");
return NULLOBJ;
}
PutFieldComponent(FIELD1, 0, &k, x);
PutFieldComponent(FIELD1, 1, &k, y);
PutFieldComponent(FIELD1, 2, &k, z);
bounds[0] = MIN(bounds[0], x);
bounds[1] = MAX(bounds[1], x);
bounds[2] = MIN(bounds[2], y);
bounds[3] = MAX(bounds[3], y);
bounds[4] = MIN(bounds[4], z);
bounds[5] = MAX(bounds[5], z);
}
/*Create the dataform*/
dataForm = NewUnstructuredDataForm("geometry vertices", 1, &size, bounds, dataFormData);
}
else if (0 == strcmp(keyword, "NORMALS"))
{
if (1 != fscanf(inFile, " %ld \n", &size))
{
FileFormatError("ReadGTFile", "# normals expected");
return NULLOBJ;
}
/*Create a new dataform dataset*/
normalsData = NewDataset("geometry vertices data", 1, &size, 3);
SetCurField(FIELD1, normalsData);
bounds[0] = bounds[2] = bounds[4] = plusInf;
bounds[1] = bounds[3] = bounds[5] = minusInf;
for (k = 0; k < size; ++k)
{
if (3 != fscanf(inFile, "%g %g %g \n", &x, &y, &z))
{
FileFormatError("ReadGTFile", "Vertices expected");
return NULLOBJ;
}
PutFieldComponent(FIELD1, 0, &k, x);
PutFieldComponent(FIELD1, 1, &k, y);
PutFieldComponent(FIELD1, 2, &k, z);
bounds[0] = MIN(bounds[0], x);
bounds[1] = MAX(bounds[1], x);
bounds[2] = MIN(bounds[2], y);
bounds[3] = MAX(bounds[3], y);
bounds[4] = MIN(bounds[4], z);
bounds[5] = MAX(bounds[5], z);
}
}
else if (0 == strcmp(keyword, "SCALAR"))
{
ObjPtr tempData;
if (1 != fscanf(inFile, " %s \n", &keyword))
{
FileFormatError("ReadGTFile", "dataset name expected");
return NULLOBJ;
}
if (size < 2)
{
FileFormatError("ReadGTFile", "scalar values only after vertices");
return NULLOBJ;
}
/*Create a new temporary dataset*/
tempData = NewDataset(keyword, 1, &size, 0);
SetCurField(FIELD1, tempData);
for (k = 0; k < size; ++k)
{
if (1 != fscanf(inFile, " %g \n", &x))
{
FileFormatError("ReadGTFile", "Vertices expected");
return NULLOBJ;
}
PutFieldComponent(FIELD1, 0, &k, x);
}
PrefixList(scalarDatasets, tempData);
}
else if (0 == strcmp(keyword, "POLYGON"))
{
if (!geometryDataset)
{
geometryDataset = NewGeometryDataset("geometry");
}
if (1 != fscanf(inFile, " %ld \n", &nVertices))
{
FileFormatError("ReadGTFile", "# poly vertices expected");
return NULLOBJ;
}
for (k = 0; k < nVertices; ++k)
{
if (1 != fscanf(inFile, " %ld \n", &vertex))
{
FileFormatError("ReadGTFile", "Poly vertex expected");
return NULLOBJ;
}
if (k < 100) polygonVertices[k] = vertex;
}
AppendPolygonToDataset(geometryDataset, MIN(nVertices, 100), polygonVertices);
}
else
{
sprintf(tempStr, "Bad keyword: '%s'", keyword);
FileFormatError("ReadGTFile", tempStr);
}
}
if (!feof(inFile))
{
FileFormatError("ReadGTFile", "Not all of the file was read");
}
fclose(inFile);
if (dataForm)
{
if (normalsData)
{
SetDatasetForm(normalsData, dataForm);
if (geometryDataset)
{
SetVar(geometryDataset, NORMALSOBJ, normalsData);
}
}
if (geometryDataset && dataForm)
{
SetDatasetForm(geometryDataset, dataForm);
RegisterNewDataset(geometryDataset);
}
if (scalarDatasets)
{
runner = LISTOF(scalarDatasets);
while (runner)
{
SetDatasetForm(runner -> thing, dataForm);
RegisterNewDataset(runner -> thing);
runner = runner -> next;
}
}
}
}
else
{
sprintf(tempStr, "File %s cannot be opened.", fileName);
FileFormatError("ReadGTFile", tempStr);
return ObjFalse;
}
}
void InitFiles(void)
/*Initializes the file handling routines*/
{
ObjPtr icon, fileReader;
fileReaderClass = NewObject(NULLOBJ, 0);
AddToReferenceList(fileReaderClass);
SetMethod(fileReaderClass, READALL, ReadOldFile);
SetMethod(fileReaderClass, NEWCTLWINDOW, ShowFileReaderControls);
icon = NewIcon(0, 0, ICONFILEREADER, "?");
SetVar(icon, ICONLOC, NULLOBJ);
SetVar(icon, FORMAT, NewString("File Reader"));
SetVar(fileReaderClass, DEFAULTICON, icon);
SetVar(fileReaderClass, SAVEEXTENSION, NewString("filrdr"));
SetVar(fileReaderClass, TIMEDDATASETS, ObjTrue);
SetMethod(fileReaderClass, SHOWCONTROLS, NewControlWindow);
SetVar(fileReaderClass, DOUBLECLICK, NewString(OF_SHOW_CONTROLS));
SetVar(fileReaderClass, CLASSID, NewInt(CLASS_FILEREADER));
SetMethod(icon, MAKE1HELPSTRING, MakeReaderIconHelp);
AddSnapVar(fileReaderClass, EXTENSION);
AddSnapVar(fileReaderClass, TIMEDDATASETS);
SetMethod(fileReaderClass, SAVECPANEL, SaveSnapshotControls);
SetMethod(fileReaderClass, SAVEALLCONTROLS, LogSnapshotControls);
#ifndef RELEASE
DefineFormat("DD", "dd", ReadDDFile, 0);
DefineFormat("DD2", "dd2", ReadDD2File, 0);
#endif
DefineFormat("NFF", "nff", ReadNFFXFile, 0);
DefineFormat("SY", "sy", ReadSYFile, 0);
#ifndef RELEASE
DefineFormat("GT", "gt", ReadGTFile, 0);
#endif
#ifdef LORENZ
DefineFormat("LRZ", "lrz", ReadLRZFile, 0);
#endif
#ifdef HDFDEF
fileReader = NewFileReader("HDF");
SetVar(fileReader, EXTENSION, NewString("hdf"));
SetVar(fileReader, OUTOFBOUNDS, NewInt(0));
AddSnapVar(fileReader, OUTOFBOUNDS);
SetVar(fileReader, HELPSTRING,
NewString("This file reader reads scientific datasets using the HDF data \
format, developed at the National Center for Supercomputing Applications. Libraries \
for reading HDF files and documentation are available via anonymous ftp from \
ftp.ncsa.uiuc.edu."));
SetVar(fileReader, READVECTOR, ObjTrue);
AddSnapVar(fileReader, READVECTOR);
SetVar(fileReader, WRAPPOLAR, ObjTrue);
AddSnapVar(fileReader, WRAPPOLAR);
SetVar(fileReader, COMPRESSDATA, ObjFalse);
AddSnapVar(fileReader, COMPRESSDATA);
SetVar(fileReader, BOTHERWITHAXES, ObjFalse);
AddSnapVar(fileReader, BOTHERWITHAXES);
SetMethod(fileReader, READALL, ReadHDFFile);
SetMethod(fileReader, ADDCONTROLS, AddHDFControls);
ApplySavedSettings(fileReader);
#endif
fileReader = NewFileReader("PDB");
SetVar(fileReader, EXTENSION, NewString("ent"));
SetVar(fileReader, HELPSTRING,
NewString("This file reader reads scientific datasets using the Protein \
Data Bank file format. Files are read as datasets over unstructured grids."));
SetMethod(fileReader, READALL, ReadPDBFile);
ApplySavedSettings(fileReader);
#ifndef RELEASE
DefineFormat("JAK", "", ReadJAKFile, 0);
DefineFormat("JAKB", "", ReadJAKBFile, 0);
#endif
fileClass = NewObject(NULLOBJ, 0);
AddToReferenceList(fileClass);
SetVar(fileClass, CLASSID, NewInt(CLASS_FILE));
InitHwuFiles();
InitWhisselFiles();
#ifdef JERRYFILES
InitJerryFiles();
#endif
InitSIMSFiles();
InitFileSystem();
}
void KillFiles(void)
/*Kills the file handling routines*/
{
int k;
#ifdef SHMDIDDLE
if (sharedSegment)
{
shmdt(sharedSegment);
}
#endif
#ifdef JERRYFILES
KillJerryFiles();
#endif
KillSIMSFiles();
KillWhisselFiles();
KillHwuFiles();
DeleteThing(fileClass);
DeleteThing(fileReaderClass);
}
|
