ObjectMap.c

/* 
A* -------------------------------------------------------------------
B* This file contains source code for the PyMOL computer program
C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific. 
D* -------------------------------------------------------------------
E* It is unlawful to modify or remove this copyright notice.
F* -------------------------------------------------------------------
G* Please see the accompanying LICENSE file for further information. 
H* -------------------------------------------------------------------
I* Additional authors of this source file include:
-* 
-* 
-*
Z* -------------------------------------------------------------------
*/

#include"os_predef.h"
#include"os_std.h"
#include"os_gl.h"

#include"OOMac.h"
#include"ObjectMap.h"
#include"Base.h"
#include"MemoryDebug.h"
#include"Map.h"
#include"Parse.h"
#include"Isosurf.h"
#include"Vector.h"
#include"Color.h"
#include"main.h"
#include"Scene.h"
#include"PConv.h"
#include"Word.h"
#include"Vector.h"
#include"PyMOLGlobals.h"
#include"Matrix.h"



#ifndef _PYMOL_NOPY
#ifdef _PYMOL_NUMPY
typedef struct {
  PyObject_HEAD
  char *data;
  int nd;
  int *dimensions, *strides;
  PyObject *base;
  void *descr;
  int flags;
} MyArrayObject;
#endif
#endif

static int ObjectMapIsStateValidActive(ObjectMap *I, int state)
{
  if((state>=0)&&(state<I->NState))
    if(I->State[state].Active)
      return 1;
  return 0;
}

void ObjectMapTransformMatrix(ObjectMap *I, int state, double *matrix)
{
  if(ObjectMapIsStateValidActive(I,state))
    ObjectStateTransformMatrix(&I->State[state].State,matrix);
  ObjectMapUpdateExtents(I);
}

void ObjectMapResetMatrix(ObjectMap *I, int state)
{
  if(ObjectMapIsStateValidActive(I,state))
    ObjectStateResetMatrix(&I->State[state].State);
  ObjectMapUpdateExtents(I);
}
int ObjectMapGetMatrix(ObjectMap *I,int state,double **matrix)
{
  if(ObjectMapIsStateValidActive(I,state)) {
    *matrix = ObjectStateGetMatrix(&I->State[state].State);
    return true;
  }
  return false;
}

int ObjectMapSetMatrix(ObjectMap *I,int state,double *matrix)
{
  if(ObjectMapIsStateValidActive(I,state)) {
    ObjectStateSetMatrix(&I->State[state].State,matrix);
    return true;
  }
  return false;
}

int ObjectMapStateGetExcludedStats(PyMOLGlobals *G,ObjectMapState *ms,float *vert_vla, float beyond,float within, float *level)
{
  double sum=0.0,sumsq=0.0;
  float mean,stdev;  
  int cnt = 0;
  int list_size;
  float cutoff = beyond;
  MapType *voxelmap = NULL;

  if(vert_vla) {
    list_size = VLAGetSize(vert_vla)/3;
  } else {
    list_size = 0;
  }
  if(cutoff<within)
    cutoff = within;

  if(list_size) 
    voxelmap=MapNew(G,-cutoff,vert_vla,list_size,NULL);

  if(voxelmap||(!list_size)) {
    int a,b,c;
    int h,k,l,i,j;
    int *fdim = ms->FDim;
    float *v,f_val;
    int within_flag, within_default=false;
    int beyond_flag;
    
    Isofield *field = ms->Field;
    if(list_size)
      MapSetupExpress(voxelmap);  

    within_flag=true;
    beyond_flag=true;

    if(within<R_SMALL4)
      within_default = true;
    for(c=0;c<fdim[2];c++) {
      for(b=0;b<fdim[1];b++) {
        for(a=0;a<fdim[0];a++) {
          if(list_size) {
            within_flag = within_default;
            beyond_flag = true;

            v = F4Ptr(field->points,a,b,c,0);
            
            MapLocus(voxelmap,v,&h,&k,&l);
            i=*(MapEStart(voxelmap,h,k,l));
            if(i) {
              j=voxelmap->EList[i++];
              while(j>=0) {
                if(!within_flag) {
                  if(within3f(vert_vla+3*j,v,within)) {                  
                    within_flag=true;
                  }
                }
                if(within3f(vert_vla+3*j,v,beyond)) {
                  beyond_flag=false;
                  break;
                }
                j=voxelmap->EList[i++];
              }
            }
          }

          if(within_flag&&beyond_flag) { /* point isn't too close to any vertex */
            f_val = F3(field->data,a,b,c);
            sum+=f_val;
            sumsq+=(f_val*f_val);
            cnt++;
          }
        }
      }
    }
    if(voxelmap) 
      MapFree(voxelmap);
  }
  if(cnt) {
    mean = (float)(sum/cnt);
    stdev = (float)sqrt1d((sumsq - (sum*sum/cnt))/(cnt));
    level[1] = mean;
    level[0] = mean-stdev;
    level[2] = mean+stdev;
  }
  return cnt;
}

int ObjectMapStateGetDataRange(PyMOLGlobals *G,ObjectMapState *ms, float *min, float *max)
{
  float max_val=0.0F, min_val=0.0F;
  CField *data = ms->Field->data;
  int cnt = data->dim[0] * data->dim[1] * data->dim[2];
  float *raw_data = (float*)data->data;
  if(cnt) {
    int a;
    min_val = (max_val = *(raw_data++));
    for(a=1;a<cnt;a++) {
      double f_val = *(raw_data++);
      if(min_val > f_val) min_val = f_val;
      if(max_val < f_val) max_val = f_val;
    }
  }
  *min = min_val;
  *max = max_val;
  return cnt;
}

int ObjectMapInterpolate(ObjectMap *I,int state,float *array,float *result,int *flag,int n)
{
  int ok=false;
  if(state<0) state=0;
  if(state<I->NState)
    if(I->State[state].Active)
      ok = ObjectMapStateInterpolate(&I->State[state],array,result,flag,n);
  return(ok);
}



static int ObjectMapStateTrim(PyMOLGlobals *G, ObjectMapState *ms, 
                                  float *mn, float *mx,int quiet)
{
  int div[3];
  int min[3];
  int max[3];
  int fdim[4];
  int new_min[3], new_max[3], new_fdim[3];
  int a,b,c,d,e,f;
  float *vt,*vr;
  float v[3];
  float grid[3];
  Isofield *field;
  int result = true;
  float orig_size = 1.0F;
  float new_size = 1.0F;
  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:

    {
      float tst[3],frac_tst[3];
      float frac_mn[3];
      float frac_mx[3];
      int hit_flag = false;

      /* compute the limiting box extents in fractional space */

      for(a=0;a<8;a++) {
        tst[0] = (a&0x1) ? mn[0] : mx[0];
        tst[1] = (a&0x2) ? mn[1] : mx[1];
        tst[2] = (a&0x4) ? mn[2] : mx[2];
        transform33f3f(ms->Crystal->RealToFrac, tst,frac_tst);
        if(!a) {
          copy3f(frac_tst,frac_mn);
          copy3f(frac_tst,frac_mx);
        } else {
          for(b=0;b<3;b++) {
            frac_mn[b] = (frac_mn[b]>frac_tst[b]) ? frac_tst[b] : frac_mn[b];
            frac_mx[b] = (frac_mx[b]<frac_tst[b]) ? frac_tst[b] : frac_mx[b];
          }
        }
      }
      for(a=0;a<3;a++) {
        div[a]=ms->Div[a];
        min[a]=ms->Min[a];
        max[a]=ms->Max[a];
        fdim[a]=ms->FDim[a];
      }
      fdim[3]=3;

      {
        int first_flag[3] = {false,false,false};
        for(d=0;d<3;d++) {
          int tst_min,tst_max;
          float v_min, v_max;
          for(c=0;c<(fdim[d]-1);c++) {
            tst_min = c+min[d];
            tst_max = tst_min+1;
            v_min=tst_min/((float)div[d]);
            v_max=tst_max/((float)div[d]);
            if(((v_min>=frac_mn[d]) && (v_min<=frac_mx[d]))||
               ((v_max>=frac_mn[d]) && (v_max<=frac_mx[d]))) {
              if(!first_flag[d]) {
                first_flag[d]=true;
                new_min[d] = tst_min;
                new_max[d] = tst_max;
              } else {
                new_min[d] = (new_min[d] > tst_min) ? tst_min : new_min[d];
                new_max[d] = (new_max[d] < tst_max) ? tst_max : new_max[d];
              }
            }
          }
          new_fdim[d] = (new_max[d] - new_min[d]) + 1;
        }
        hit_flag = first_flag[0] && first_flag[1] & first_flag[2];
      }

      if(hit_flag) 
        hit_flag = (new_fdim[0]!=fdim[0])||(new_fdim[1]!=fdim[1])||
          (new_fdim[2]!=fdim[2]);

      if(hit_flag) {
        orig_size = fdim[0]*fdim[1]*fdim[2];
        new_size = new_fdim[0]*new_fdim[1]*new_fdim[2];
        
        field=IsosurfFieldAlloc(G,new_fdim);
        field->save_points = ms->Field->save_points;


        for(c=0;c<new_fdim[2];c++) {
          f = c+(new_min[2] - min[2]);
          for(b=0;b<new_fdim[1];b++) {
            e = b+(new_min[1] - min[1]);
            for(a=0;a<new_fdim[0];a++) {
              d = a+(new_min[0] - min[0]);
              vt = F4Ptr(field->points,a,b,c,0);
              vr = F4Ptr(ms->Field->points,d,e,f,0);
              copy3f(vr,vt);
              F3(field->data,a,b,c) = F3(ms->Field->data,d,e,f);
            }
          }
        }
        IsosurfFieldFree(G,ms->Field);
        for(a=0;a<3;a++) {
          ms->Min[a]=new_min[a];
          ms->Max[a]=new_max[a];
          ms->FDim[a]=new_fdim[a];
        }
        ms->Field = field;
        
        /* compute new extents */
        v[2]=(ms->Min[2])/((float)ms->Div[2]);
        v[1]=(ms->Min[1])/((float)ms->Div[1]);
        v[0]=(ms->Min[0])/((float)ms->Div[0]);
        
        transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
        
        v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
        v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
        v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
        
        transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
        
        /* new corner */
        {
          float vv[3];
          d = 0;
          for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1)) {
            v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
            for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
              v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
              for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
                v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                transform33f3f(ms->Crystal->FracToReal,v,vv);
                copy3f(vv,ms->Corner+3*d);
                d++;
              }
            }
          }
        }
        result=true;
      }
    }
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:
  case cMapSourceChempyBrick:
    {
      int hit_flag = false;
      float *origin = ms->Origin;

      for(a=0;a<3;a++) {
        min[a]=ms->Min[a];
        grid[a]=ms->Grid[a];
        max[a]=ms->Max[a];
        fdim[a]=ms->FDim[a];
      }
      fdim[3]=3;

      {
        int first_flag[3] = {false,false,false};
        for(d=0;d<3;d++) {
          int tst_min,tst_max;
          float v_min, v_max;
          for(c=0;c<(fdim[d]-1);c++) {
            tst_min = c+min[d];
            tst_max = tst_min+1;
            v_min=origin[d]+grid[d]*(tst_min);
            v_max=origin[d]+grid[d]*(tst_max);
            if(((v_min>=mn[d]) && (v_min<=mx[d]))||
               ((v_max>=mn[d]) && (v_max<=mx[d]))) {
              if(!first_flag[d]) {
                first_flag[d]=true;
                hit_flag=true;
                new_min[d] = tst_min;
                new_max[d] = tst_max;
              } else {
                new_min[d] = (new_min[d] > tst_min) ? tst_min : new_min[d];
                new_max[d] = (new_max[d] < tst_max) ? tst_max : new_max[d];
              }
            }
          }
          new_fdim[d] = (new_max[d] - new_min[d]) + 1;
        }
        hit_flag = first_flag[0] && first_flag[1] & first_flag[2];
      }

      if(hit_flag) 
        hit_flag = (new_fdim[0]!=fdim[0])||(new_fdim[1]!=fdim[1])||
          (new_fdim[2]!=fdim[2]);

      if(hit_flag) {

        orig_size = fdim[0]*fdim[1]*fdim[2];
        new_size = new_fdim[0]*new_fdim[1]*new_fdim[2];

        field=IsosurfFieldAlloc(G,new_fdim);
        field->save_points = ms->Field->save_points;

        for(c=0;c<new_fdim[2];c++) {
          f = c+(new_min[2] - min[2]);
          for(b=0;b<new_fdim[1];b++) {
            e = b+(new_min[1] - min[1]);
            for(a=0;a<new_fdim[0];a++) {
              d = a+(new_min[0] - min[0]);
              vt = F4Ptr(field->points,a,b,c,0);
              vr = F4Ptr(ms->Field->points,d,e,f,0);
              copy3f(vr,vt);
              F3(field->data,a,b,c) = F3(ms->Field->data,d,e,f);
            }
          }
        }
        IsosurfFieldFree(G,ms->Field);
        for(a=0;a<3;a++) {
          ms->Min[a]=new_min[a];
          ms->Max[a]=new_max[a];
          ms->FDim[a]=new_fdim[a];
          if(ms->Dim) 
            ms->Dim[a]=new_fdim[a];
        }

        ms->Field = field;

        for(e=0;e<3;e++) {
          ms->ExtentMin[e] = ms->Origin[e]+ms->Grid[e]*ms->Min[e];
          ms->ExtentMax[e] = ms->Origin[e]+ms->Grid[e]*ms->Max[e];
        }
        
        d=0;
        for(c=0;c<ms->FDim[2];c+=ms->FDim[2]-1) {
          v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
          
          for(b=0;b<ms->FDim[1];b+=ms->FDim[1]-1) {
            v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
            
            for(a=0;a<ms->FDim[0];a+=ms->FDim[0]-1) {
              v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
              copy3f(v,ms->Corner+3*d);
              d++;
            }
          }
        }
        result=true;
      }
    }
    break;
  }
  if(result&&(!quiet)) {
    PRINTFB(G,FB_ObjectMap,FB_Actions)
      " ObjectMap: Map volume reduced by %2.0f%%.\n",
       (100*(orig_size-new_size))/orig_size
      ENDFB(G);
  }
  return result;
}

static int ObjectMapStateDouble(PyMOLGlobals *G,ObjectMapState *ms)
{
  int div[3];
  int min[3];
  int max[3];
  int fdim[4];
  int a,b,c;
  float v[3],vr[3];
  float *vt;
  float x,y,z;
  float grid[3];

  Isofield *field;

  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:
    for(a=0;a<3;a++) {
      div[a]=ms->Div[a]*2;
      min[a]=ms->Min[a]*2;
      max[a]=ms->Max[a]*2;
      fdim[a]=ms->FDim[a]*2-1;
    }
    fdim[3]=3;
    field=IsosurfFieldAlloc(G,fdim);
    field->save_points = ms->Field->save_points;
    for(c=0;c<fdim[2];c++) {
      v[2]=(c+min[2])/((float)div[2]);
      z = (c&0x1) ? 0.5F : 0.0F;
      for(b=0;b<fdim[1];b++) {
        v[1]=(b+min[1])/((float)div[1]);        
        y = (b&0x1) ? 0.5F : 0.0F;
        for(a=0;a<fdim[0];a++) {
          v[0]=(a+min[0])/((float)div[0]);
          transform33f3f(ms->Crystal->FracToReal,v,vr);
          x = (a&0x1) ? 0.5F : 0.0F;
          vt = F4Ptr(field->points,a,b,c,0);
          copy3f(vr,vt);
          if((a&0x1)||(b&0x1)||(c&0x1)) {
            F3(field->data,a,b,c) = FieldInterpolatef(ms->Field->data,
                                    a/2,
                                    b/2,
                                    c/2,x,y,z);
          } else {
            F3(field->data,a,b,c) = F3(ms->Field->data,a/2,b/2,c/2);
          }
        }
      }
    }
    IsosurfFieldFree(G,ms->Field);
    for(a=0;a<3;a++) {
      ms->Min[a]=min[a];
      ms->Max[a]=max[a];
      ms->FDim[a]=fdim[a];
      ms->Div[a]=div[a];
    }

    ms->Field = field;
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:
  case cMapSourceChempyBrick:
  case cMapSourceVMDPlugin:
    for(a=0;a<3;a++) {
      grid[a]=ms->Grid[a]/2.0F;
      min[a]=ms->Min[a]*2;
      max[a]=ms->Max[a]*2;
      fdim[a]=ms->FDim[a]*2-1;
    }
    fdim[3]=3;

    field=IsosurfFieldAlloc(G,fdim);
    field->save_points = ms->Field->save_points;

    for(c=0;c<fdim[2];c++) {
      v[2]=ms->Origin[2]+grid[2]*(c+min[2]);
      z = (c&0x1) ? 0.5F : 0.0F;
      for(b=0;b<fdim[1];b++) {
        v[1]=ms->Origin[1]+grid[1]*(b+min[1]);
        y = (b&0x1) ? 0.5F : 0.0F;
        for(a=0;a<fdim[0];a++) {
          v[0]=ms->Origin[0]+grid[0]*(a+min[0]);
          x = (a&0x1) ? 0.5F : 0.0F;
          vt = F4Ptr(field->points,a,b,c,0);
          copy3f(v,vt);
          if((a&0x1)||(b&0x1)||(c&0x1)) {
            F3(field->data,a,b,c) = FieldInterpolatef(ms->Field->data,
                                    a/2,
                                    b/2,
                                    c/2,x,y,z);
          } else {
            F3(field->data,a,b,c) = F3(ms->Field->data,a/2,b/2,c/2);
          }
        }
      }
    }
    IsosurfFieldFree(G,ms->Field);
    for(a=0;a<3;a++) {
      ms->Min[a]=min[a];
      ms->Max[a]=max[a];
      ms->FDim[a]=fdim[a];
      if(ms->Dim) 
        ms->Dim[a]=fdim[a];
      if(ms->Grid)
        ms->Grid[a]=grid[a];
    }
    ms->Field = field;

    break;
  }
  return 1;
}

static int ObjectMapStateHalve(PyMOLGlobals *G,ObjectMapState *ms,int smooth)
{
  int div[3];
  int min[3];
  int max[3];
  int fdim[4];
  int a,b,c;
  float v[3],vr[3];
  float *vt;
  float x,y,z;
  float grid[3];

  Isofield *field;

  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:
    {
      int *old_div,*old_min,*old_max;
      int a_2, b_2, c_2;

      for(a=0;a<3;a++) {
        div[a]=ms->Div[a]/2; 
        /* note that when Div is odd, a one-cell extrapolation will
           occur in order to preserve map size and get us onto a even
           number of sampling intervals for the cell */

        min[a]=ms->Min[a]/2;
        max[a]=ms->Max[a]/2;
        while((min[a]*2)<ms->Min[a])
          min[a]++;
        while((max[a]*2)>ms->Max[a]) 
          max[a]--;
 
        fdim[a]=(max[a]-min[a])+1;
      }
      fdim[3]=3;
      old_div = ms->Div;
      old_min = ms->Min;
      old_max = ms->Max;
      
      if(smooth) 
        FieldSmooth3f(ms->Field->data);

      field=IsosurfFieldAlloc(G,fdim);
      field->save_points = ms->Field->save_points;
      
      /*
        printf("old_div %d %d %d\n",old_div[0],old_div[1],old_div[2]);
        printf("old_min %d %d %d\n",old_min[0],old_min[1],old_min[2]);
        printf("min %d %d %d\n",min[0],min[1],min[2]);
        printf("%d %d %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]);
      */

      for(c=0;c<fdim[2];c++) {
        v[2]=(c+min[2])/((float)div[2]);
        c_2 = (c+min[2])*2 - old_min[2];
        z = (v[2] - ((c_2 + old_min[2])/(float)old_div[2]))*old_div[2];
        if(c_2>=old_max[2]) {
          c_2 = old_max[2]-1;
          z = (v[2] - ((c_2 + old_min[2])/(float)old_div[2]))*old_div[2];
        }
        for(b=0;b<fdim[1];b++) {
          v[1]=(b+min[1])/((float)div[1]);        
          b_2 = (b+min[1])*2 - old_min[1];
          y = (v[1] - ((b_2 + old_min[1])/(float)old_div[1]))*old_div[1];
          if(b_2>=old_max[1]) {
            b_2 = old_max[1]-1;
            y = (v[1] - ((b_2 + old_min[1])/(float)old_div[1]))*old_div[1];
          }
          for(a=0;a<fdim[0];a++) {
            v[0]=(a+min[0])/((float)div[0]);
            a_2 = (a+min[0])*2 - old_min[0];
            x = (v[0] - ((a_2 + old_min[0])/(float)old_div[0]))*old_div[0];
            if(a_2>=old_max[0]) {
              a_2 = old_max[0]-1;
              x = (v[0] - ((a_2 + old_min[0])/(float)old_div[0]))*old_div[0];
            }
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            vt = F4Ptr(field->points,a,b,c,0);
            copy3f(vr,vt);
            F3(field->data,a,b,c) = FieldInterpolatef(ms->Field->data,
                                                      a_2,
                                                      b_2,
                                                      c_2,x,y,z);
          }
        }
      }
      IsosurfFieldFree(G,ms->Field);
      for(a=0;a<3;a++) {
        ms->Min[a]=min[a];
        ms->Max[a]=max[a];
        ms->FDim[a]=fdim[a];
        ms->Div[a]=div[a];
      }
      
      ms->Field = field;
      
      /* compute new extents */
      v[2]=(ms->Min[2])/((float)ms->Div[2]);
      v[1]=(ms->Min[1])/((float)ms->Div[1]);
      v[0]=(ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
      
      v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
      v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
      v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
      
      /* new corner */
      {
        float vv[3];
        int d=0;
                  
        for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1)) {
          v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
          for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
            v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
            for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
              v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
              transform33f3f(ms->Crystal->FracToReal,v,vv);
              copy3f(vv,ms->Corner+3*d);
              d++;
            }
          }
        }
      }
    }
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:
  case cMapSourceChempyBrick:
  case cMapSourceVMDPlugin:
    for(a=0;a<3;a++) {
      grid[a]=ms->Grid[a]*2.0F;
      min[a]=ms->Min[a]/2;
      max[a]=ms->Max[a]/2;
      fdim[a]=(ms->FDim[a]+1)/2;
    }
    fdim[3]=3;

    field=IsosurfFieldAlloc(G,fdim);
    field->save_points = ms->Field->save_points;

    for(c=0;c<fdim[2];c++) {
      v[2]=ms->Origin[2]+grid[2]*(c+min[2]);
      for(b=0;b<fdim[1];b++) {
        v[1]=ms->Origin[1]+grid[1]*(b+min[1]);
        for(a=0;a<fdim[0];a++) {
          v[0]=ms->Origin[0]+grid[0]*(a+min[0]);
          vt = F4Ptr(field->points,a,b,c,0);
          copy3f(v,vt);
          F3(field->data,a,b,c) = F3(ms->Field->data,a*2,b*2,c*2);
        }
      }
    }
    IsosurfFieldFree(G,ms->Field);
    for(a=0;a<3;a++) {
      ms->Min[a]=min[a];
      ms->Max[a]=max[a];
      ms->FDim[a]=fdim[a];
      if(ms->Dim) 
        ms->Dim[a]=fdim[a];
      if(ms->Grid)
        ms->Grid[a]=grid[a];
    }
    ms->Field = field;

    break;
  }
  return 1;
}

int ObjectMapTrim(ObjectMap *I,int state,float *mn, float *mx, int quiet)
{
  int a;
  int result=true;
  int update=false;

  /* TO DO: convert mn and mx into map local coordinates if map itself is transformed...  */

  if(state<0) {
    for(a=0;a<I->NState;a++) {
      if(I->State[a].Active) {
        if(ObjectMapStateTrim(I->Obj.G,&I->State[a],mn,mx,quiet))
          update=true;
        else
          result=false;
      }
    }
  } else if((state>=0)&&(state<I->NState)&&(I->State[state].Active)) {
    update = result = ObjectMapStateTrim(I->Obj.G,&I->State[state],mn,mx,quiet);
  } else {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors)
      " ObjectMap-Error: invalidate state.\n"
      ENDFB(I->Obj.G);
    result=false;
  }
  if(update)
    ObjectMapUpdateExtents(I);
  return(result);
}

int ObjectMapDouble(ObjectMap *I,int state)
{
  int a;
  int result=true;
  if(state<0) {
    for(a=0;a<I->NState;a++) {
      if(I->State[a].Active)
        result = result && ObjectMapStateDouble(I->Obj.G,&I->State[a]);
    }
  } else if((state>=0)&&(state<I->NState)&&(I->State[state].Active)) {
    ObjectMapStateDouble(I->Obj.G,&I->State[state]);
  } else {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors)
      " ObjectMap-Error: invalidate state.\n"
      ENDFB(I->Obj.G);
    result=false;
  }
  return(result);
}

int ObjectMapHalve(ObjectMap *I,int state,int smooth)
{
  int a;
  int result=true;
  if(state<0) {
    for(a=0;a<I->NState;a++) {
      if(I->State[a].Active)
        result = result && ObjectMapStateHalve(I->Obj.G,&I->State[a],smooth);
    }
    
  } else if((state>=0)&&(state<I->NState)&&(I->State[state].Active)) {
    ObjectMapStateHalve(I->Obj.G,&I->State[state],smooth);
  } else {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors)
      " ObjectMap-Error: invalidate state.\n"
      ENDFB(I->Obj.G);
    result=false;
  }
  ObjectMapUpdateExtents(I);

  return(result);
}

int ObjectMapStateContainsPoint(ObjectMapState *ms,float *point)
{
  register int result=false;
  register float x,y,z;
  register int x_floor, y_floor, z_floor;
  register int x_ceil, y_ceil, z_ceil;

  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:
    {
      float frac[3];

      transform33f3f(ms->Crystal->RealToFrac,point,frac); 
      
      x = (ms->Div[0] * frac[0]);
      y = (ms->Div[1] * frac[1]);
      z = (ms->Div[2] * frac[2]);
      x_floor= floor(x);
      x_ceil = ceil(x);
      y_floor = floor(y);
      y_ceil = ceil(y);
      z_floor = floor(z);
      z_ceil = ceil(z);

      if((x_floor>=ms->Min[0])&&(x_ceil<=ms->Max[0])&&
         (y_floor>=ms->Min[1])&&(y_ceil<=ms->Max[1])&&
         (z_floor>=ms->Min[2])&&(z_ceil<=ms->Max[2]))
        result = true;
    }
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:
  case cMapSourceChempyBrick:
  case cMapSourceVMDPlugin:
    x = (point[0] - ms->Origin[0])/ms->Grid[0];
    y = (point[1] - ms->Origin[1])/ms->Grid[1];
    z = (point[2] - ms->Origin[2])/ms->Grid[2];

    x_floor= floor(x);
    x_ceil = ceil(x);
    y_floor = floor(y);
    y_ceil = ceil(y);
    z_floor = floor(z);
    z_ceil = ceil(z);
    
    if((x_floor>=ms->Min[0])&&(x_ceil<=ms->Max[0])&&
       (y_floor>=ms->Min[1])&&(y_ceil<=ms->Max[1])&&
       (z_floor>=ms->Min[2])&&(z_ceil<=ms->Max[2]))
      result = true;

    if((x>=ms->Min[0])&&(x<=ms->Max[0])&&
       (y>=ms->Min[1])&&(y<=ms->Max[1])&&
       (z>=ms->Min[2])&&(z<=ms->Max[2]))
      result = true;
    break;
    }
  return(result);
}

int ObjectMapStateInterpolate(ObjectMapState *ms,float *array,float *result,int *flag, int n)
{
  int ok=true;
  float *inp;
  float *out;
  int a,b,c;
  float x,y,z;
  inp = array;
  out = result;
  
  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:
    {
      float frac[3];

      while(n--) {
        
        /* get the fractional coordinate */
        transform33f3f(ms->Crystal->RealToFrac,inp,frac); 
        
        inp+=3;

        /* compute the effective lattice offset as a function of cell spacing */

        x = (ms->Div[0] * frac[0]);
        y = (ms->Div[1] * frac[1]);
        z = (ms->Div[2] * frac[2]);
        
        /* now separate the integration and fractional parts for interpolation */
           
        a=(int)floor(x);
        b=(int)floor(y);
        c=(int)floor(z);
        x-=a;
        y-=b;
        z-=c;
        
        if(flag) *flag = 1;
        /* wrap into the map */

        if(a<ms->Min[0]) {
          if(x<0.99F) {
            ok=false;
            if(flag) *flag = 0;
          }
          x=0.0F;
          a=ms->Min[0];
        } else if(a>=ms->FDim[0]+ms->Min[0]-1) {
          if(x>0.01F) {
            ok=false;
            if(flag) *flag = 0;
          }
          x=0.0F;
          a=ms->FDim[0]+ms->Min[0]-1;
        }
        
        if(b<ms->Min[1]) {
          if(y<0.99F) {
            ok=false;
            if(flag) *flag = 0;
          }
          y=0.0F;
          b=ms->Min[1];
        } else if(b>=ms->FDim[1]+ms->Min[1]-1) {
          if(y>0.01F) {
            ok=false;
            if(flag) *flag = 0;
          }
          y=0.0F;
          b=ms->FDim[1]+ms->Min[1]-1;
        }
        
        if(c<ms->Min[2]) {
          if(z<0.99F) {
            ok=false;
            if(flag) *flag = 0;
          }
          z=0.0F;
          c=ms->Min[2];
        } else if(c>=ms->FDim[2]+ms->Min[2]-1) {
          if(z>0.01) {
            ok=false;
            if(flag) *flag = 0;
          }
          z=0.0F;
          c=ms->FDim[2]+ms->Min[2]-1;
        }
        /*      printf("%d %d %d %8.3f %8.3f %8.3f\n",a,b,c,x,y,z);*/
        *(result++)=FieldInterpolatef(ms->Field->data,
                                      a-ms->Min[0],
                                      b-ms->Min[1],
                                      c-ms->Min[2],x,y,z);
        if(flag) flag++;

      }
    }
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:
  case cMapSourceChempyBrick:
  case cMapSourceVMDPlugin:
    while(n--) {

      x = (inp[0] - ms->Origin[0])/ms->Grid[0];
      y = (inp[1] - ms->Origin[1])/ms->Grid[1];
      z = (inp[2] - ms->Origin[2])/ms->Grid[2];
      inp+=3;
      
      a=(int)floor(x);
      b=(int)floor(y);
      c=(int)floor(z);
      x-=a;
      y-=b;
      z-=c;

      if(flag) *flag = 1;
      if(a<ms->Min[0]) {
        x=0.0F;
        a=ms->Min[0];
        ok=false;
        if(flag) *flag = 0;
      } else if(a>=ms->Max[0]) {
        x=1.0F;
        a=ms->Max[0]-1;
        ok=false;
        if(flag) *flag = 0;
      }

      if(b<ms->Min[1]) {
        y=0.0F;
        b=ms->Min[1];
        ok=false;
        if(flag) *flag = 0;
      } else if(b>=ms->Max[1]) {
        y=1.0F;
        b=ms->Max[1]-1;
        ok=false;
        if(flag) *flag = 0;
      }

      if(c<ms->Min[2]) {
        z=0.0F;
        c=ms->Min[2];
        ok=false;
        if(flag) *flag = 0;
      } else if(c>=ms->Max[2]) {
        z=1.0F;
        c=ms->Max[2]-1;
        ok=false;
        if(flag) *flag = 0;
      }
      /*      printf("%d %d %d %8.3f %8.3f %8.3f\n",a,b,c,x,y,z);*/
      *(result++)=FieldInterpolatef(ms->Field->data,
                                    a-ms->Min[0],
                                    b-ms->Min[1],
                                    c-ms->Min[2],x,y,z);
      if(flag) flag++;
    }
    break;
  default:
    ok=false;
    break;
  }
    return(ok);
}

int ObjectMapNumPyArrayToMapState(PyMOLGlobals *G,ObjectMapState *I,PyObject *ary);

void ObjectMapStateRegeneratePoints(ObjectMapState *ms)
{
  int a,b,c,e;
  float v[3],vr[3];
  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:
    for(c=0;c<ms->FDim[2];c++)
      {
        v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
        for(b=0;b<ms->FDim[1];b++) {
          v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
          for(a=0;a<ms->FDim[0];a++) {
            v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            for(e=0;e<3;e++) 
              F4(ms->Field->points,a,b,c,e) = vr[e];
          }
        }
      }
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:
  case cMapSourceChempyBrick:
  case cMapSourceVMDPlugin:
    for(c=0;c<ms->FDim[2];c++) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      for(b=0;b<ms->FDim[1];b++) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        for(a=0;a<ms->FDim[0];a++) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          for(e=0;e<3;e++) {
            F4(ms->Field->points,a,b,c,e) = v[e];
          }
        }
      }
  }
  default: 
    break;
  }
}

#ifndef _PYMOL_NOPY
static PyObject *ObjectMapStateAsPyList(ObjectMapState *I)
{
  PyObject *result = NULL;

  result = PyList_New(16);
  PyList_SetItem(result,0,PyInt_FromLong(I->Active));
  if(I->Crystal) {
    PyList_SetItem(result,1,CrystalAsPyList(I->Crystal));
  } else {
    PyList_SetItem(result,1,PConvAutoNone(Py_None));
  }
  if(I->Origin) {
    PyList_SetItem(result,2,PConvFloatArrayToPyList(I->Origin,3));
  } else {
    PyList_SetItem(result,2,PConvAutoNone(Py_None));
  }
  if(I->Range) {
    PyList_SetItem(result,3,PConvFloatArrayToPyList(I->Range,3));
  } else {
    PyList_SetItem(result,3,PConvAutoNone(Py_None));
  }
  if(I->Dim) {
    PyList_SetItem(result,4,PConvIntArrayToPyList(I->Dim,3));
  } else {
    PyList_SetItem(result,4,PConvAutoNone(Py_None));
  }
  if(I->Grid) {
    PyList_SetItem(result,5,PConvFloatArrayToPyList(I->Grid,3));
  } else {
    PyList_SetItem(result,5,PConvAutoNone(Py_None));
  }
  PyList_SetItem(result,6,PConvFloatArrayToPyList(I->Corner,24));
  PyList_SetItem(result,7,PConvFloatArrayToPyList(I->ExtentMin,3));
  PyList_SetItem(result,8,PConvFloatArrayToPyList(I->ExtentMax,3));
  PyList_SetItem(result,9,PyInt_FromLong(I->MapSource));

  PyList_SetItem(result,10,PConvIntArrayToPyList(I->Div,3));
  PyList_SetItem(result,11,PConvIntArrayToPyList(I->Min,3));
  PyList_SetItem(result,12,PConvIntArrayToPyList(I->Max,3));
  PyList_SetItem(result,13,PConvIntArrayToPyList(I->FDim,4));
  
  PyList_SetItem(result,14,IsosurfAsPyList(I->Field));
  PyList_SetItem(result,15,ObjectStateAsPyList(&I->State));
#if 0
  int Active;
  CCrystal *Crystal;
  int Div[3],Min[3],Max[3],FDim[4];
  Isofield *Field;
  float Corner[24];
  int *Dim;
  float *Origin;
  float *Range;
  float *Grid;
  float ExtentMin[3],ExtentMax[3];
#endif

  return(PConvAutoNone(result));  
}
#endif
#ifndef _PYMOL_NOPY
static PyObject *ObjectMapAllStatesAsPyList(ObjectMap *I)
{
  PyObject *result=NULL;
  int a;
  result = PyList_New(I->NState);
  for(a=0;a<I->NState;a++) {
    if(I->State[a].Active) {
      PyList_SetItem(result,a,ObjectMapStateAsPyList(I->State+a));
    } else {
      PyList_SetItem(result,a,PConvAutoNone(NULL));
    }
  }
  return(PConvAutoNone(result));  

}
#endif
static int ObjectMapStateCopy(PyMOLGlobals *G,ObjectMapState *src,ObjectMapState *I)
{
  int ok=true;
  if(ok) {
    I->Active = src->Active;
    if(I->Active) {

      if(src->Crystal)
        I->Crystal = CrystalCopy(src->Crystal);
      else
        I->Crystal = NULL;
      
      if(src->Origin) {
        I->Origin = Alloc(float,3);
        if(I->Origin) {
          copy3f(src->Origin,I->Origin);
        }
      } else {
        I->Origin = NULL;
      }

      if(src->Range) {
        I->Range = Alloc(float,3);
        if(I->Range) {
          copy3f(src->Range,I->Range);
        }
      } else {
        I->Origin = NULL;
      }

      if(src->Grid) {
        I->Grid = Alloc(float,3);
        if(I->Grid) {
          copy3f(src->Grid,I->Grid);
        }
      } else {
        I->Origin = NULL;
      }

      if(src->Dim) {
        I->Dim = Alloc(int,4);
        if(I->Dim) {
          copy3f(src->Dim,I->Dim);
        }
      } else {
        I->Origin = NULL;
      }
      
      {
        int a;
        for(a=0;a<24;a++)
          I->Corner[a] = src->Corner[a];
      }

      copy3f(src->ExtentMin, I->ExtentMin);
      copy3f(src->ExtentMax, I->ExtentMax);
      
      I->MapSource = src->MapSource;

      copy3f(src->Div, I->Div);
      copy3f(src->Min, I->Min);
      copy3f(src->Max, I->Max);
      copy3f(src->FDim, I->FDim);
      
      I->Field = IsosurfNewCopy(G,src->Field);
      ObjectStateCopy(&I->State, &src->State);
      if(ok) ObjectMapStateRegeneratePoints(I);
    }
  }
  return(ok);
}

#ifndef _PYMOL_NOPY
static int ObjectMapStateFromPyList(PyMOLGlobals *G,ObjectMapState *I,PyObject *list)
{
  int ok=true;
  int ll=0;
  PyObject *tmp;
  if(ok) ok=(list!=NULL);
  if(ok) {
    if(!PyList_Check(list))
      I->Active=false;
    else {
      if(ok) ok=PyList_Check(list);
      if(ok) ll = PyList_Size(list);
      /* TO SUPPORT BACKWARDS COMPATIBILITY...
         Always check ll when adding new PyList_GetItem's */

      if(ok) ok = PConvPyIntToInt(PyList_GetItem(list,0),&I->Active);
      if(ok) {
        tmp = PyList_GetItem(list,1);
        if(tmp == Py_None)
          I->Crystal = NULL;
        else 
          ok = ((I->Crystal=CrystalNewFromPyList(G,tmp))!=NULL);
      }
      if(ok) {
        tmp = PyList_GetItem(list,2);
        if(tmp == Py_None)
          I->Origin = NULL;
        else 
          ok = PConvPyListToFloatArray(tmp,&I->Origin);
      }
      if(ok) {
        tmp = PyList_GetItem(list,3);
        if(tmp == Py_None)
          I->Range = NULL;
        else 
          ok = PConvPyListToFloatArray(tmp,&I->Range);
      }
      if(ok) {
        tmp = PyList_GetItem(list,4);
        if(tmp == Py_None)
          I->Dim = NULL;
        else 
          ok = PConvPyListToIntArray(tmp,&I->Dim);
      }
      if(ok) {
        tmp = PyList_GetItem(list,5);
        if(tmp == Py_None)
          I->Grid = NULL;
        else 
          ok = PConvPyListToFloatArray(tmp,&I->Grid);
      }
      if(ok) ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list,6),I->Corner,24);
      if(ok) ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list,7),I->ExtentMin,3);
      if(ok) ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list,8),I->ExtentMax,3);
      if(ok) ok = PConvPyIntToInt(PyList_GetItem(list,9),&I->MapSource);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,10),I->Div,3);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,11),I->Min,3);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,12),I->Max,3);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,13),I->FDim,4);
      if(ok) ok = ((I->Field=IsosurfNewFromPyList(G,PyList_GetItem(list,14)))!=NULL);
      if(ok&&(ll>15)) ok = ObjectStateFromPyList(G,PyList_GetItem(list,15),&I->State);
      if(ok) ObjectMapStateRegeneratePoints(I);
    }
  }
  return(ok);
}
#endif
#ifndef _PYMOL_NOPY
static int ObjectMapAllStatesFromPyList(ObjectMap *I,PyObject *list)
{
  int ok=true;
  int a;
  VLACheck(I->State,ObjectMapState,I->NState);
  if(ok) ok=PyList_Check(list);
  if(ok) {
    for(a=0;a<I->NState;a++) {
      ok = ObjectMapStateFromPyList(I->Obj.G,I->State+a,PyList_GetItem(list,a));
      if(!ok) break;
    }
  }
  return(ok);
}
#endif

PyObject *ObjectMapAsPyList(ObjectMap *I)
{
#ifdef _PYMOL_NOPY
  return NULL;
#else

  PyObject *result = NULL;

  result = PyList_New(3);
  PyList_SetItem(result,0,ObjectAsPyList(&I->Obj));
  PyList_SetItem(result,1,PyInt_FromLong(I->NState));
  PyList_SetItem(result,2,ObjectMapAllStatesAsPyList(I));

  return(PConvAutoNone(result));  
#endif
}

int ObjectMapNewFromPyList(PyMOLGlobals *G,PyObject *list,ObjectMap **result)
{
#ifdef _PYMOL_NOPY
  return 0;
#else

  int ok = true;
  int ll;
  ObjectMap *I=NULL;
  (*result) = NULL;
  
  if(ok) ok=(list!=NULL);
  if(ok) ok=PyList_Check(list);
  if(ok) ll = PyList_Size(list);
  /* TO SUPPORT BACKWARDS COMPATIBILITY...
   Always check ll when adding new PyList_GetItem's */
  I=ObjectMapNew(G);
  if(ok) ok = (I!=NULL);

  if(ok) ok = ObjectFromPyList(G,PyList_GetItem(list,0),&I->Obj);
  if(ok) ok = PConvPyIntToInt(PyList_GetItem(list,1),&I->NState);
  if(ok) ok = ObjectMapAllStatesFromPyList(I,PyList_GetItem(list,2));
  if(ok) {
    (*result) = I;
    ObjectMapUpdateExtents(I);
  } else {
    /* cleanup? */
  }

  return(ok);
#endif
}

int ObjectMapNewCopy(PyMOLGlobals *G,ObjectMap *src,ObjectMap **result,int source_state, int target_state)
{
  int ok = true;
  ObjectMap *I=NULL;
  I=ObjectMapNew(G);
  if(ok) ok = (I!=NULL);
  if(ok) ok = ObjectCopyHeader(&I->Obj,&src->Obj);
  if(ok) {
    if(source_state==-1) { /* all states */
      int state;
      I->NState = src->NState;
      VLACheck(I->State,ObjectMapState,I->NState);
      for(state=0;state<src->NState;state++) {
        ok = ObjectMapStateCopy(G, src->State + state, I->State + state);
      }
    } else {
      if(target_state<0)
        target_state = 0;
      if(source_state<0)
        source_state = 0;
      VLACheck(I->State,ObjectMapState,target_state);    
      if(source_state<src->NState) {
        ok = ObjectMapStateCopy(G,src->State + source_state, I->State + target_state);    
        if(I->NState<target_state)
          I->NState = target_state;
      } else {
        ok=false;
        /* to do */
      }
    }
  }
  if(ok) 
    *result = I;
  return ok;
}

ObjectMapState *ObjectMapGetState(ObjectMap *I,int state)
{
  ObjectMapState *result = NULL;
  if(state<0)
    state = 0;
  if(state<I->NState)
    result = &I->State[state];
  return(result);
}

ObjectMapState *ObjectMapStatePrime(ObjectMap *I,int state)
{
  ObjectMapState *ms = NULL;
  if(state<0)
    state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);
  return(ms);
}

ObjectMapState *ObjectMapStateGetActive(ObjectMap *I,int state)
{
  ObjectMapState *ms = NULL;
  if(state>=0) {
    if(state<I->NState) {
      ms=&I->State[state];
      if(!ms->Active)
        ms = NULL;
    }
  }
  return(ms);
}

void ObjectMapUpdateExtents(ObjectMap *I)
{
  int a;
  float *min_ext,*max_ext;
  float tr_min[3],tr_max[3];
  I->Obj.ExtentFlag=false;
  
  
  for(a=0;a<I->NState;a++) {
    ObjectMapState *ms = I->State+a;
    if(ms->Active) {
      if(I->State[a].State.Matrix) { 
        transform44d3f(ms->State.Matrix,ms->ExtentMin,tr_min);
        transform44d3f(ms->State.Matrix,ms->ExtentMax,tr_max);
        {
          float tmp;
          int a;
          for(a=0;a<3;a++) 
            if(tr_min[a]>tr_max[a]) { tmp=tr_min[a]; tr_min[a]=tr_max[a]; tr_max[a]=tmp;}
        }
        min_ext = tr_min;
        max_ext = tr_max;
      } else {
        min_ext = ms->ExtentMin;
        max_ext = ms->ExtentMax;
      }
      
      if(!I->Obj.ExtentFlag) {
        copy3f(min_ext,I->Obj.ExtentMin);
        copy3f(max_ext,I->Obj.ExtentMax);
        I->Obj.ExtentFlag=true;
      } else {
        min3f(min_ext,I->Obj.ExtentMin,I->Obj.ExtentMin);
        max3f(max_ext,I->Obj.ExtentMax,I->Obj.ExtentMax);
      }
    }
  }

  if(I->Obj.TTTFlag && I->Obj.ExtentFlag) {
    float *ttt;
    double tttd[16];
    if(ObjectGetTTT(&I->Obj,&ttt,-1)) {
      convertTTTfR44d(ttt,tttd);
      MatrixTransformExtentsR44d3f(tttd,
                                   I->Obj.ExtentMin,I->Obj.ExtentMax,
                                   I->Obj.ExtentMin,I->Obj.ExtentMax);
    }
  }

  PRINTFD(I->Obj.G,FB_ObjectMap)
    " ObjectMapUpdateExtents-DEBUG: ExtentFlag %d\n",I->Obj.ExtentFlag
    ENDFD;
}

void ObjectMapStateClamp(ObjectMapState *I,float clamp_floor, float clamp_ceiling)
{
  int a,b,c;
  float *fp;

  for(a=0;a<I->FDim[0];a++) 
    for(b=0;b<I->FDim[1];b++)
      for(c=0;c<I->FDim[2];c++) {
        fp = F3Ptr(I->Field->data,a,b,c);
        if(*fp<clamp_floor)
          *fp = clamp_floor;
        else if(*fp>clamp_ceiling)
          *fp = clamp_ceiling;
      }
}

int ObjectMapStateSetBorder(ObjectMapState *I,float level)
{
  int result = true;
  int a,b,c;

  c=I->FDim[2]-1;
  for(a=0;a<I->FDim[0];a++) 
    for(b=0;b<I->FDim[1];b++)
      {
        F3(I->Field->data,a,b,0) = level;
        F3(I->Field->data,a,b,c) = level;
      }

  a=I->FDim[0]-1;
  for(b=0;b<I->FDim[1];b++) 
    for(c=0;c<I->FDim[2];c++)
      {
        F3(I->Field->data,0,b,c) = level;
        F3(I->Field->data,a,b,c) = level;
      }

  b=I->FDim[1]-1;
  for(a=0;a<I->FDim[0];a++) 
    for(c=0;c<I->FDim[2];c++)
      {
        F3(I->Field->data,a,0,c) = level;
        F3(I->Field->data,a,b,c) = level;
      }
  return(result);
}

void ObjectMapStatePurge(PyMOLGlobals *G,ObjectMapState *I)
{
  ObjectStatePurge(&I->State);
  if(I->Field) {
    IsosurfFieldFree(G,I->Field);
    I->Field=NULL;
  }
  FreeP(I->Origin);
  FreeP(I->Dim);
  FreeP(I->Range);
  FreeP(I->Grid);
  OOFreeP(I->Crystal);
  I->Active=false;
}

static void ObjectMapFree(ObjectMap *I) {

  int a;
  for(a=0;a<I->NState;a++) {
    if(I->State[a].Active)
      ObjectMapStatePurge(I->Obj.G,I->State+a);
  }
  VLAFreeP(I->State);
  ObjectPurge(&I->Obj);
  OOFreeP(I);
}

static void ObjectMapUpdate(ObjectMap *I) {
  if(!I->Obj.ExtentFlag) {
    ObjectMapUpdateExtents(I);
    if(I->Obj.ExtentFlag)
      SceneInvalidate(I->Obj.G);
  }
}

static void ObjectMapInvalidate(ObjectMap *I,int rep,int level,int state)
{
  if(level>=cRepInvExtents) {
    I->Obj.ExtentFlag=false;
  }
  if((rep<0)||(rep==cRepDot)) {
    int a;
    for(a=0;a<I->NState;a++) {
      if(I->State[a].Active)
        I->State[a].have_range = false;
    }
  }
  SceneInvalidate(I->Obj.G);
}

static void ObjectMapRender(ObjectMap *I,RenderInfo *info)
{
  PyMOLGlobals *G = I->Obj.G;
  int state = info->state;
  CRay *ray = info->ray;
  Picking **pick = info->pick;
  int pass = info->pass;
  ObjectMapState *ms = NULL;
  if(!pass) {
    if(state<I->NState)
      if(I->State[state].Active)
        ms=&I->State[state];
    
    if(ms) {
      float *corner = ms->Corner;
      float tr_corner[24];
      ObjectPrepareContext(&I->Obj,ray);
      
      if(ms->State.Matrix) { /* transform the corners before drawing */
        int a;
        for(a=0;a<8;a++) {
          transform44d3f(ms->State.Matrix,corner+3*a,tr_corner+3*a);
        }
        corner = tr_corner;
      }

      
      if(I->Obj.RepVis[cRepExtent]) {
        if(ray) {
          float *vc;
          float radius = ray->PixelRadius/1.4142F;
          vc = ColorGet(G,I->Obj.Color);
          ray->fColor3fv(ray,vc);
          ray->fSausage3fv(ray,corner+3*0,corner+3*1,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*0,corner+3*2,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*2,corner+3*3,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*1,corner+3*3,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*0,corner+3*4,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*1,corner+3*5,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*2,corner+3*6,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*3,corner+3*7,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*4,corner+3*5,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*4,corner+3*6,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*6,corner+3*7,radius,vc,vc);
          ray->fSausage3fv(ray,corner+3*5,corner+3*7,radius,vc,vc);
        } else if(G->HaveGUI && G->ValidContext) {
          if(pick) {
          } else {
            ObjectUseColor(&I->Obj);
            glDisable(GL_LIGHTING); 
            glBegin(GL_LINES);
            glVertex3fv(corner+3*0);
            glVertex3fv(corner+3*1);
            
            glVertex3fv(corner+3*0);
            glVertex3fv(corner+3*2);
          
            glVertex3fv(corner+3*2);
            glVertex3fv(corner+3*3);
          
            glVertex3fv(corner+3*1);
            glVertex3fv(corner+3*3);
          
            glVertex3fv(corner+3*0);
            glVertex3fv(corner+3*4);
          
            glVertex3fv(corner+3*1);
            glVertex3fv(corner+3*5);
          
            glVertex3fv(corner+3*2);
            glVertex3fv(corner+3*6);
          
            glVertex3fv(corner+3*3);
            glVertex3fv(corner+3*7);
          
            glVertex3fv(corner+3*4);
            glVertex3fv(corner+3*5);
          
            glVertex3fv(corner+3*4);
            glVertex3fv(corner+3*6);
          
            glVertex3fv(corner+3*6);
            glVertex3fv(corner+3*7);
          
            glVertex3fv(corner+3*5);
            glVertex3fv(corner+3*7);
          
            glEnd();
            glEnable(GL_LIGHTING);
          }
        }
      }
      
      if(I->Obj.RepVis[cRepDot]) {
        /* note, the following rep doesn't work with state matrices yet */

        if(!ms->have_range) {
          double sum=0.0,sumsq=0.0;
          CField *data = ms->Field->data;
          int cnt = data->dim[0] * data->dim[1] * data->dim[2];
          float *raw_data = (float*)data->data;
          int a;
          for(a=0;a<cnt;a++) {
            double f_val = *(raw_data++);
            sum+=f_val;
            sumsq+=(f_val*f_val);
          }
          if(cnt) {
            float mean,stdev;
            mean = (float)(sum/cnt);
            stdev = (float)sqrt1d((sumsq - (sum*sum/cnt))/(cnt));
            ms->high_cutoff = mean + stdev;
            ms->low_cutoff = mean - stdev;
            ms->have_range = true;
          }
        }
        if(ms->have_range && SettingGet_b(G,NULL,I->Obj.Setting,cSetting_dot_normals)) {
          IsofieldComputeGradients(G,ms->Field);
        }
        if(ms->have_range) {
          register int a;
          CField *data = ms->Field->data;
          register int cnt = data->dim[0] * data->dim[1] * data->dim[2];
          CField *points = ms->Field->points;
          CField *gradients = NULL;

          if(SettingGet_b(G,NULL,I->Obj.Setting,cSetting_dot_normals)) {
            gradients = ms->Field->gradients;
          }
          if(data && points) {
            register float *raw_data = (float*)data->data;
            register float *raw_point = (float*)points->data;
            register float *raw_gradient = NULL;
            register float high_cut = ms->high_cutoff, low_cut = ms->low_cutoff;
            register float width = SettingGet_f(G,NULL,I->Obj.Setting,cSetting_dot_width);
            
            if(ray) {
              float radius =  ray->PixelRadius * width/1.4142F;
              float vc[3];
              int color = I->Obj.Color;
              int ramped = ColorCheckRamped(G,I->Obj.Color);
              
              {
                float *tmp = ColorGet(G,I->Obj.Color);
                copy3f(tmp,vc);
              }

              for(a=0;a<cnt;a++) {
                register float f_val = *(raw_data++);
                if((f_val >= high_cut) || (f_val <= low_cut)) {
                  if(ramped) {
                    ColorGetRamped(G,color,raw_point,vc,state);                    
                    ray->fColor3fv(ray,vc);
                  }
                  ray->fSphere3fv(ray,raw_point,radius);}
                raw_point+=3;
              }
            } else if(G->HaveGUI && G->ValidContext) {
              if(pick) {
              } else {
                if(gradients) {
                  raw_gradient = (float*)gradients->data;                  
                } else {
                  glDisable(GL_LIGHTING);
                }
                {
                  int ramped = ColorCheckRamped(G,I->Obj.Color);
                  float vc[3];
                  int color = I->Obj.Color;
                  float gt[3];

                  glPointSize(width);
                  glDisable(GL_POINT_SMOOTH);
                  glBegin(GL_POINTS);
                  ObjectUseColor(&I->Obj);
                  for(a=0;a<cnt;a++) {
                    register float f_val = *(raw_data++);
                    if(f_val >= high_cut) {
                      if(raw_gradient) {
                        normalize23f(raw_gradient,gt);
                        invert3f(gt);
                        glNormal3fv(gt);
                      }
                      if(ramped) {
                        ColorGetRamped(G,color,raw_point,vc,state);                    
                        glColor3fv(vc);
                      }
                      glVertex3fv(raw_point);
                    } else if(f_val <= low_cut) {
                      if(raw_gradient) {
                        normalize23f(raw_gradient,gt);
                        glNormal3fv(gt);
                      }
                      if(ramped) {
                        ColorGetRamped(G,color,raw_point,vc,state);                    
                        glColor3fv(vc);
                      }
                      glVertex3fv(raw_point);
                    }
                    if(raw_gradient) 
                      raw_gradient+=3;
                    raw_point+=3;
                  }
                  glEnd();
                }
                glEnable(GL_POINT_SMOOTH);
              }
            }
          }
        }
      }
    }
  }
}

void ObjectMapStateInit(PyMOLGlobals *G,ObjectMapState *I) 
{
  ObjectMapStatePurge(G,I);
  ObjectStateInit(G,&I->State);
  I->Crystal = CrystalNew(G);
  I->Field = NULL;
  I->Origin = NULL;
  I->Dim = NULL;
  I->Range = NULL;
  I->Grid = NULL;
  I->have_range = false;
}
int ObjectMapGetNStates(ObjectMap *I)     
{
  return(I->NState);
}
/*========================================================================*/
ObjectMap *ObjectMapNew(PyMOLGlobals *G)
{
  OOAlloc(G,ObjectMap);

  ObjectInit(G,(CObject*)I);
  I->Obj.type = cObjectMap;

  
  I->NState = 0;
  I->State=VLAMalloc(1,sizeof(ObjectMapState),5,true); /* autozero important */

  {
    int a;
    for(a=0;a<cRepCnt;a++)
      I->Obj.RepVis[a] = false;
    I->Obj.RepVis[cRepExtent]=true; 
  }
  I->Obj.fFree = (void (*)(struct CObject *))ObjectMapFree;
  I->Obj.fUpdate =  (void (*)(struct CObject *)) ObjectMapUpdate;
  I->Obj.fRender =(void (*)(struct CObject *, RenderInfo *))ObjectMapRender;
  I->Obj.fInvalidate =(void (*)(struct CObject *,int,int,int))ObjectMapInvalidate;  
  I->Obj.fGetNFrame = (int (*)(struct CObject *)) ObjectMapGetNStates;

  return(I);
}
/*========================================================================*/
ObjectMapState *ObjectMapNewStateFromDesc(PyMOLGlobals *G,ObjectMap *I,
                                          ObjectMapDesc *inp_md,int state,int quiet)
{
  int ok=true;
  float v[3];
  int a,b,c,d;
  float *fp;
  ObjectMapState *ms = NULL;
  ObjectMapDesc _md,*md;
  ms=ObjectMapStatePrime(I,state);
  
  md = &_md;
  *(md) = *(inp_md);

  if(I) {
    ms->Origin=Alloc(float,3);
    ms->Range=Alloc(float,3);
    ms->Dim=Alloc(int,3);
    ms->Grid=Alloc(float,3);
    ms->MapSource=cMapSourceDesc;
  }
  switch(md->mode) {
  case cObjectMap_OrthoMinMaxGrid: /* Orthorhombic: min, max, spacing, centered over range  */

    subtract3f(md->MaxCorner,md->MinCorner,v);
    for(a=0;a<3;a++) { if(v[a]<0.0) swap1f(md->MaxCorner+a,md->MinCorner+a); };
    subtract3f(md->MaxCorner,md->MinCorner,v);
    for(a=0;a<3;a++) {
      md->Dim[a] = (int)(v[a]/md->Grid[a]);
      if(md->Dim[a]<1) md->Dim[a]=1;
      if((md->Dim[a]*md->Grid[a])<v[a]) md->Dim[a]++;
    }

    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMap: Dim %d %d %d\n",md->Dim[0],md->Dim[1],md->Dim[2]
      ENDFB(I->Obj.G);

    average3f(md->MaxCorner,md->MinCorner,v);
    for(a=0;a<3;a++) { md->MinCorner[a] = v[a]-0.5F*md->Dim[a]*md->Grid[a]; }

    if(Feedback(I->Obj.G,FB_ObjectMap,FB_Blather)) {
      dump3f(md->MinCorner," ObjectMap: MinCorner:");
      dump3f(md->MaxCorner," ObjectMap: MaxCorner:");
      dump3f(md->Grid," ObjectMap: Grid:");
    }
    
    /* now populate the map data structure */

    copy3f(md->MinCorner,ms->Origin);
    copy3f(md->Grid,ms->Grid);
    for(a=0;a<3;a++) ms->Range[a] = md->Grid[a] * (md->Dim[a]-1);

    /* these maps start at zero */
    for(a=0;a<3;a++) {
      ms->Min[a]=0; 
      ms->Max[a]=md->Dim[a]-1;
      ms->Div[a]=ms->Dim[a]-1;
    }
    
    /* define corners */

    for(a=0;a<8;a++) copy3f(ms->Origin,ms->Corner+3*a);

    d = 0;
    for(c=0;c<2;c++) {
      {
        v[2] = (c ? ms->Range[2] : 0.0F);
        for(b=0;b<2;b++) {
          v[1]= (b ? ms->Range[1] : 0.0F);
          for(a=0;a<2;a++) {
            v[0]= (a ? ms->Range[0] : 0.0F);
            add3f(v,ms->Corner+3*d,ms->Corner+3*d);
            d++;
          }
        }
      }
    }
    for(a=0;a<3;a++) ms->FDim[a] = ms->Max[a]-ms->Min[a]+1;
    ms->FDim[3] = 3; 

    ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
    if(!ms->Field) 
      ok=false;
    else {
      for(a=0;a<md->Dim[0];a++) {
        v[0] = md->MinCorner[0] + a * md->Grid[0];
        for(b=0;b<md->Dim[1];b++) {
          v[1] = md->MinCorner[1] + b * md->Grid[1];
          for(c=0;c<md->Dim[2];c++) {
            v[2] = md->MinCorner[2] + c * md->Grid[2];
            fp = F4Ptr(ms->Field->points,a,b,c,0);
            copy3f(v,fp);
          }
        }
      }
    }
    break;
  default:
    ok = false;
  }
  if(ok) {
    switch(md->init_mode) {
    case 0:
      for(a=0;a<md->Dim[0];a++) {      
        for(b=0;b<md->Dim[1];b++) {
          for(c=0;c<md->Dim[2];c++) {
            F3(ms->Field->data,a,b,c)=0.0F;
          }
        }
      }
      break;
    case 1:
      for(a=0;a<md->Dim[0];a++) {      
        for(b=0;b<md->Dim[1];b++) {
          for(c=0;c<md->Dim[2];c++) {
            F3(ms->Field->data,a,b,c)=1.0F;
          }
        }
      }
      break;
    case -2: /* for testing... */
      for(a=0;a<md->Dim[0];a++) {      
        for(b=0;b<md->Dim[1];b++) {
          for(c=0;c<md->Dim[2];c++) {
            F3(ms->Field->data,a,b,c)=(float)sqrt1d(a*a+b*b+c*c);
          }
        }
      }
      break;
    }
  }
  
  if(ok) {
    copy3f(ms->Origin,ms->ExtentMin);
    copy3f(ms->Origin,ms->ExtentMax);
    add3f(ms->Range,ms->ExtentMax,ms->ExtentMax);
    ObjectMapUpdateExtents(I);
  }
  if(!ok) {
    ErrMessage(I->Obj.G,"ObjectMap","Unable to create map");
    ObjectMapFree(I);
    I=NULL;
  } else {
    if(!quiet) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Actions) 
        " ObjectMap: Map created.\n"
        ENDFB(I->Obj.G);
    }
  }
  
  return(ms);
}
/*========================================================================*/
static int ObjectMapCCP4StrToMap(ObjectMap *I,char *CCP4Str,int bytes,int state,int quiet)
{
  char *p;
  int *i;
  unsigned int *u;
  unsigned char *uc,c0,c1,c2,c3;
  float *f;
  float dens;
  int a,b,c,d,e;
  float v[3],vr[3],maxd,mind;
  int ok = true;
  int little_endian = 1,map_endian;
  /* CCP4 named from their docs */
  int nc,nr,ns;
  int map_mode;
  int ncstart,nrstart,nsstart;
  int nx,ny,nz;
  float xlen,ylen,zlen,alpha,beta,gamma;
  int sym_skip;
  int mapc,mapr,maps;
  int cc[3],xref[3];
  int n_pts;
  double sum,sumsq;
  float mean,stdev;
  int normalize;
  ObjectMapState *ms;
  int expectation;

  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);

  normalize=(int)SettingGet(I->Obj.G,cSetting_normalize_ccp4_maps);
  maxd = -FLT_MAX;
  mind = FLT_MAX;
  p=CCP4Str;
  little_endian = *((char*)&little_endian);
  map_endian = (*p||*(p+1));

  if(bytes<256*sizeof(int)) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors)
      " ObjectMapCCP4: Map appears to be truncated -- aborting."
      ENDFB(I->Obj.G);
    return(0);
  }
  if(little_endian!=map_endian) {
    if(!quiet) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
        " ObjectMapCCP4: Map appears to be reverse endian, swapping...\n"
        ENDFB(I->Obj.G);
    }
    c = bytes;
    u = (unsigned int*)p;
    uc = (unsigned char *)u;
    while(c>3) {
      c0 = *(uc++);
      c1 = *(uc++);
      c2 = *(uc++);
      c3 = *(uc++);
      uc = (unsigned char *)u;
      *(uc++) = c3;
      *(uc++) = c2;
      *(uc++) = c1;
      *(uc++) = c0;
      u++;
      c-=4;
    }
  }
  i = (int*)p;
  nc=*(i++); /* columns */
  nr=*(i++); /* rows */
  ns=*(i++); /* sections */
  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: NC %d   NR %d   NS %d\n",
      nc,nr,ns
      ENDFB(I->Obj.G);
  }
  map_mode = *(i++); /* mode */

  
  if(map_mode!=2) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors)
      "ObjectMapCCP4-ERR: Only map mode 2 currently supported (this map is mode %d)",map_mode
      ENDFB(I->Obj.G);
    return(0);
  }

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: Map is mode %d.\n",map_mode
      ENDFB(I->Obj.G);
  }
  ncstart = *(i++);
  nrstart = *(i++);
  nsstart = *(i++);

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: NCSTART %d   NRSTART %d   NSSTART  %d\n",
      ncstart,nrstart,nsstart
      ENDFB(I->Obj.G);
  }

  nx = *(i++);
  ny = *(i++);
  nz = *(i++);

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: NX %d   NY %d   NZ  %d \n",
      nx,ny,nz
      ENDFB(I->Obj.G);
  }

  xlen = *(float*)(i++);
  ylen = *(float*)(i++);
  zlen = *(float*)(i++);

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: X %8.3f   Y %8.3f  Z  %8.3f \n",
      xlen,ylen,zlen
      ENDFB(I->Obj.G);
  }

  alpha = *(float*)(i++);
  beta = *(float*)(i++);
  gamma = *(float*)(i++);

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: alpha %8.3f   beta %8.3f  gamma %8.3f \n",
      alpha,beta,gamma
      ENDFB(I->Obj.G);
  }

  mapc = *(i++);
  mapr = *(i++);
  maps = *(i++);

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: MAPC %d   MAPR %d  MAPS  %d \n",
      mapc,mapr,maps
      ENDFB(I->Obj.G);
  }

  i+=4;
  sym_skip = *(i++);

  {
    int skew = *(i++);

    if(skew) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors)
        "ObjectMapCCP4-ERR: PyMOL doesn't know how to handle skewed maps. Sorry!\n"
        ENDFB(I->Obj.G);
      return(0);
    }
  }

  n_pts = nc*ns*nr;

  /* at least one EM map encountered lacks NZ, so we'll try to guess it */
  
  if((nx==ny)&&(!nz)) {
    nz = nx;

    if(!quiet) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Warnings)
        " ObjectMapCCP4: NZ value is zero, but NX = NY, so we'll guess NZ = NX = NY.\n"
        ENDFB(I->Obj.G);
    }
  }

  expectation = sym_skip + sizeof(int)*(256+n_pts);

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: sym_skip %d bytes %d expectation %d\n",
      sym_skip, bytes,  expectation
      ENDFB(I->Obj.G);
  }

  if(bytes<expectation) {
    if(bytes==(expectation-sym_skip)) { 
      /* accomodate bogus CCP4 map files with bad symmetry length information */
      if(!quiet) {
        PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather)
          " ObjectMapCCP4: Map has invalid symmetry length -- working around.\n"
          ENDFB(I->Obj.G);
      }

      expectation -= sym_skip;
      sym_skip = 0;

    } else {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors)
        " ObjectMapCCP4: Map appears to be truncated -- aborting.\n"
        ENDFB(I->Obj.G);
      return(0);
    }
  }


  if(n_pts>1) {
    f = (float*)(p+(sizeof(int)*256)+sym_skip);
    c = n_pts;
    sum = 0.0;
    sumsq = 0.0;
    while(c--) {
      sumsq+=(*f)*(*f);
      sum+=*f++;
    }
    mean = (float)(sum/n_pts);
    stdev = (float)sqrt1d((sumsq - (sum*sum/n_pts))/(n_pts-1));
    if(stdev<0.000001)
      stdev = 1.0;
    
  } else {
    mean = 1.0;
    stdev = 1.0;
  }
  
  f = (float*)(p+(sizeof(int)*256)+sym_skip);
  mapc--; /* convert to C indexing... */
  mapr--;
  maps--;

  xref[maps]=0;
  xref[mapr]=1;
  xref[maps]=2;

  xref[maps]=0;
  xref[mapr]=1;
  xref[mapc]=2;

  ms->Div[0] = nx;
  ms->Div[1] = ny;
  ms->Div[2] = nz;

  ms->FDim[mapc] = nc;
  ms->Min[mapc] = ncstart;
  ms->Max[mapc] = nc+ncstart-1;

  ms->FDim[mapr] = nr;
  ms->Min[mapr] = nrstart;
  ms->Max[mapr] = nr+nrstart-1;

  ms->FDim[maps] = ns;
  ms->Min[maps] = nsstart;
  ms->Max[maps] = ns+nsstart-1;

  if(!quiet) {
    if(Feedback(I->Obj.G,FB_ObjectMap,FB_Blather)) {
      dump3i(ms->Div," ObjectMapCCP4: ms->Div");
      dump3i(ms->Min," ObjectMapCCP4: ms->Min");
      dump3i(ms->Max," ObjectMapCCP4: ms->Max");
      dump3i(ms->FDim," ObjectMapCCP4: ms->FDim");
    }
  }

  ms->Crystal->Dim[0] = xlen;
  ms->Crystal->Dim[1] = ylen;
  ms->Crystal->Dim[2] = zlen;

  ms->Crystal->Angle[0] = alpha;
  ms->Crystal->Angle[1] = beta;
  ms->Crystal->Angle[2] = gamma;


  ms->FDim[3]=3;
  if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
    ok=false;
  else {
    CrystalUpdate(ms->Crystal);
    /*    CrystalDump(ms->Crystal);*/
    ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
    ms->MapSource = cMapSourceCCP4;
    ms->Field->save_points=false;

    for(cc[maps]=0;cc[maps]<ms->FDim[maps];cc[maps]++)
      {
        v[maps]=(cc[maps]+ms->Min[maps])/((float)ms->Div[maps]);

        for(cc[mapr]=0;cc[mapr]<ms->FDim[mapr];cc[mapr]++) {
          v[mapr]=(cc[mapr]+ms->Min[mapr])/((float)ms->Div[mapr]);

          for(cc[mapc]=0;cc[mapc]<ms->FDim[mapc];cc[mapc]++) {
            v[mapc]=(cc[mapc]+ms->Min[mapc])/((float)ms->Div[mapc]);

            if(normalize) 
              dens = (*f-mean)/stdev;
            else 
              dens = *f;
            F3(ms->Field->data,cc[0],cc[1],cc[2]) = dens;
            if(maxd<*f) maxd = dens;
            if(mind>*f) mind = dens;
            f++;
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            for(e=0;e<3;e++) 
              F4(ms->Field->points,cc[0],cc[1],cc[2],e) = vr[e];
          }
        }
      }
  }
  if(ok) {
    d = 0;
    for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
      {
        v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
        for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
          v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
          for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
            v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            copy3f(vr,ms->Corner+3*d);
            d++;
          }
        }
      }
  }

  if(!quiet) {
    PRINTFB(I->Obj.G,FB_ObjectMap, FB_Details) 
      " ObjectMapCCP4: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
      ENDFB(I->Obj.G);
  }
  
  if(!quiet) {
    if(n_pts>1) {
      if(normalize) {
        PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
          " ObjectMapCCP4: Normalizing with mean = %8.6f and stdev = %8.6f.\n",
          mean,stdev
          ENDFB(I->Obj.G);
      } else {
        
        PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
          " ObjectMapCCP4: Map will not be normalized.\n ObjectMapCCP4: Current mean = %8.6f and stdev = %8.6f.\n",
          mean,stdev
          ENDFB(I->Obj.G);
      }
    }
  }

  if(ok) {
    v[2]=(ms->Min[2])/((float)ms->Div[2]);
    v[1]=(ms->Min[1])/((float)ms->Div[1]);
    v[0]=(ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
    
    v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
    v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
    v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
  }
#ifdef _UNDEFINED
  printf("%d %d %d %d %d %d %d %d %d\n",
         ms->Div[0],
         ms->Min[0],
         ms->Max[0],
         ms->Div[1],
         ms->Min[1],
         ms->Max[1],
         ms->Div[2],
         ms->Min[2],
         ms->Max[2]);
  printf("Okay? %d\n",ok);
  fflush(stdout);
#endif
  if(!ok) {
    ErrMessage(I->Obj.G,"ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    if(!quiet) {
      PRINTFB(I->Obj.G,FB_ObjectMap, FB_Results) 
        " ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd
        ENDFB(I->Obj.G);
    }
  }

  return(ok);
}
/*========================================================================*/
static int ObjectMapPHIStrToMap(ObjectMap *I,char *PHIStr,int bytes,int state) {
  
  char *p;
  float dens,dens_rev;
  int a,b,c,d,e;
  float v[3],maxd,mind;
  int ok = true;
  int little_endian = 1;
  /* PHI named from their docs */
  int map_endian = 0;
  int map_dim;
  int map_bytes;

  ObjectMapState *ms;

  char cc[MAXLINELEN];
  char *rev;

  little_endian = *((char*)&little_endian);


  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);

  maxd = -FLT_MAX;
  mind = FLT_MAX;
  p=PHIStr;

  if(*p)  /* use FORMATTED IO record to determine map endiness */
    map_endian = 1;
  else
    map_endian = 0;

  p+=4;

  ParseNCopy(cc,p,20);
  PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB(I->Obj.G);
  p+=20;
  p+=4;

  p+=4;
  ParseNCopy(cc,p,10);
  PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB(I->Obj.G);
  p+=10;
  ParseNCopy(cc,p,60);
  PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB(I->Obj.G);
  p+=60;
  p+=4;

  rev = (char*)&dens_rev;

  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];
  }

  map_bytes = *((int*)rev);

  map_dim = (int)(pow((map_bytes/4.0),1/3.0)+0.5);

  if((4*map_dim*map_dim*map_dim)!=map_bytes) /* consistency check */
    map_dim = 65;

  PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
      " PHIStrToMap: Map Size %d x %d x %d\n",map_dim,map_dim,map_dim
      ENDFB(I->Obj.G);
  p+=4;

  ms->FDim[0] = map_dim;
  ms->FDim[1] = map_dim;
  ms->FDim[2] = map_dim;
  ms->FDim[3] = 3;

  ms->Div[0] = (map_dim-1)/2;
  ms->Div[1] = (map_dim-1)/2;
  ms->Div[2] = (map_dim-1)/2;
  ms->Min[0] = -ms->Div[0];
  ms->Min[1] = -ms->Div[1];
  ms->Min[2] = -ms->Div[2];
  ms->Max[0] = ms->Div[0];
  ms->Max[1] = ms->Div[1];
  ms->Max[2] = ms->Div[2];

  ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
  ms->MapSource = cMapSourcePHI;
  ms->Field->save_points=false;


  for(c=0;c<ms->FDim[2];c++) { /* z y x ordering into c b a  so that x = a, etc. */
    for(b=0;b<ms->FDim[1];b++) {
      for(a=0;a<ms->FDim[0];a++) {

        if(little_endian!=map_endian) {
          rev[0]=p[3];
          rev[1]=p[2];
          rev[2]=p[1];
          rev[3]=p[0];
        } else {
          rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
          rev[1]=p[1];
          rev[2]=p[2];
          rev[3]=p[3];
        }
        dens = *((float*)rev);
        F3(ms->Field->data,a,b,c) = dens;
        if(maxd<dens) maxd = dens;
        if(mind>dens) mind = dens;
        p+=4;
      }
    }
  }
  p+=4;


  p+=4;
  ParseNCopy(cc,p,16);
  PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB(I->Obj.G);
  p+=16;
  p+=4;

  ms->Grid = Alloc(float,3);
  p+=4;
  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];
  }
  ms->Grid[0] = 1.0F/(*((float*)rev));
  ms->Grid[1] = ms->Grid[0];
  ms->Grid[2] = ms->Grid[0];
  p+=4;

  ms->Origin = Alloc(float,3);
  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];;
  }
  ms->Origin[0] = *((float*)rev);
  p+=4;

  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];;
  }
  ms->Origin[1] = *((float*)rev);
  p+=4;
  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];;
  }
  ms->Origin[2] = *((float*)rev);
  p+=4;

  p+=4;

  if(ok) {
    for(e=0;e<3;e++) {
      ms->ExtentMin[e] = ms->Origin[e]+ms->Grid[e]*ms->Min[e];
      ms->ExtentMax[e] = ms->Origin[e]+ms->Grid[e]*ms->Max[e];
    }
  }

  for(c=0;c<ms->FDim[2];c++) {
    v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
    for(b=0;b<ms->FDim[1];b++) {
      v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
      for(a=0;a<ms->FDim[0];a++) {
        v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
        for(e=0;e<3;e++) {
          F4(ms->Field->points,a,b,c,e) = v[e];
        }
      }
    }
  }

  d=0;
  for(c=0;c<ms->FDim[2];c+=ms->FDim[2]-1) {
    v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);

    for(b=0;b<ms->FDim[1];b+=ms->FDim[1]-1) {
      v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);

      for(a=0;a<ms->FDim[0];a+=ms->FDim[0]-1) {
        v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
        copy3f(v,ms->Corner+3*d);
        d++;
      }
    }
  }

  /* interpolation test code 
  { 
    float test[3];
    float result;
    float cmp;

    for(c=0;c<ms->FDim[0]-1;c++) {
      for(b=0;b<ms->FDim[1]-1;b++) {
        for(a=0;a<ms->FDim[2]-1;a++) {
          for(e=0;e<3;e++) {
            test[e] = (F4(ms->Field->points,a,b,c,e)+
                       F4(ms->Field->points,a,b,c,e))/2.0;
          }
          ObjectMapStateInterpolate(ms,test,&result,1);
          cmp = (F3(ms->Field->data,a,b,c)+
                 F3(ms->Field->data,a,b,c))/2.0;
          if(fabs(cmp-result)>0.001) {
            printf("%d %d %d\n",a,b,c);
            printf("%8.5f %8.5f\n",
                   cmp,
                   result);
          }
        }
      }
    }
  }
  */

  if(!ok) {
    ErrMessage(I->Obj.G,"ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);
  }
  return(ok);
}
/*========================================================================*/
static int ObjectMapXPLORStrToMap(ObjectMap *I,char *XPLORStr,int state,int quiet) {
  
  char *p;
  int a,b,c,d,e;
  float v[3],vr[3],dens,maxd,mind;
  char cc[MAXLINELEN];
  int n;
  int ok = true;
  ObjectMapState *ms;

  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }

  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);

  maxd = -FLT_MAX;
  mind = FLT_MAX;
  p=XPLORStr;

  while(*p) {
    p = ParseNCopy(cc,p,8);
    if(!*cc) 
      p = ParseNextLine(p);
    else if(sscanf(cc,"%i",&n)==1) {
      p=ParseWordCopy(cc,p,MAXLINELEN);
      if(strstr(cc,"!NTITLE")||(!*cc)) {
        p=ParseNextLine(p);
        while(n--) {
          p=ParseNextLine(p);          
        } 
      } else if(strstr(cc,"REMARKS")) {
        p=ParseNextLine(p);          
      } else {
        break;
      }
    }
  }
  if(*p) { /* n contains first dimension */
    ms->Div[0]=n;
    if(sscanf(cc,"%i",&ms->Min[0])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Max[0])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Div[1])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Min[1])!=1) ok=false;    
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Max[1])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Div[2])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Min[2])!=1) ok=false;    
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Max[2])!=1) ok=false;
    p=ParseNextLine(p);
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Dim[0])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Dim[1])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Dim[2])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Angle[0])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Angle[1])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Angle[2])!=1) ok=false;
    p=ParseNextLine(p);
    p = ParseNCopy(cc,p,3);
    if(strcmp(cc,"ZYX")) ok=false;
    p=ParseNextLine(p);
    
  } else {
    ok=false;
  }
  if(ok) {
    
    ms->FDim[0]=ms->Max[0]-ms->Min[0]+1;
    ms->FDim[1]=ms->Max[1]-ms->Min[1]+1;
    ms->FDim[2]=ms->Max[2]-ms->Min[2]+1;
    ms->FDim[3]=3;
    if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
      ok=false;
    else {
      CrystalUpdate(ms->Crystal);
      ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
      ms->MapSource = cMapSourceXPLOR;
      ms->Field->save_points=false;
      for(c=0;c<ms->FDim[2];c++)
        {
          v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
          p=ParseNextLine(p);
          for(b=0;b<ms->FDim[1];b++) {
            v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
            for(a=0;a<ms->FDim[0];a++) {
              v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
              p=ParseNCopy(cc,p,12);
              if(!cc[0]) {
                p=ParseNextLine(p);
                p=ParseNCopy(cc,p,12);                
              }
              if(sscanf(cc,"%f",&dens)!=1) {
                ok=false;
              } else {
                F3(ms->Field->data,a,b,c) = dens;
                if(maxd<dens) maxd = dens;
                if(mind>dens) mind = dens;
              }
              transform33f3f(ms->Crystal->FracToReal,v,vr);
              for(e=0;e<3;e++) {
                F4(ms->Field->points,a,b,c,e) = vr[e];
              }
            }
          }
          p=ParseNextLine(p);
        }
      if(ok) {
        d = 0;
        for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
          {
            v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
            for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
              v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
              for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
                v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                transform33f3f(ms->Crystal->FracToReal,v,vr);
                copy3f(vr,ms->Corner+3*d);
                d++;
              }
            }
          }
      }
    }
  }

  if(ok) {
    v[2]=(ms->Min[2])/((float)ms->Div[2]);
    v[1]=(ms->Min[1])/((float)ms->Div[1]);
    v[0]=(ms->Min[0])/((float)ms->Div[0]);

    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
    
    v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
    v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
    v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);

    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);

  }
#ifdef _UNDEFINED
    printf("%d %d %d %d %d %d %d %d %d\n",
           ms->Div[0],
           ms->Min[0],
           ms->Max[0],
           ms->Div[1],
           ms->Min[1],
           ms->Max[1],
           ms->Div[2],
           ms->Min[2],
           ms->Max[2]);
    printf("Okay? %d\n",ok);
    fflush(stdout);
#endif
  if(!ok) {
    ErrMessage(I->Obj.G,"ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd);
  }
    
  return(ok);
}
/*========================================================================*/
static int ObjectMapFLDStrToMap(ObjectMap *I,char *PHIStr,int bytes,int state) 
{
  char *p;
  float dens,dens_rev;
  int a,b,c,d,e;
  float v[3],maxd,mind;
  int ok = true;
  int little_endian = 1;
  /* PHI named from their docs */
  int map_endian = 1;
  char cc[MAXLINELEN];
  char *rev;
  int ndim=0;
  int veclen=0;
  int nspace=0;
  ObjectMapState *ms;
    int got_ndim=false;
    int got_dim1=false;
    int got_dim2=false;
    int got_dim3=false;
    int got_data=false;
    int got_veclen=false;
    int got_min_ext=false;
    int got_max_ext=false;
    int got_field=false;
    int got_nspace=false;


  little_endian = *((char*)&little_endian);
  map_endian = little_endian;

  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);

  maxd = -FLT_MAX;
  mind = FLT_MAX;

  p = PHIStr;

  while((*p)&&(!(
                 got_ndim&&
                 got_dim1&&
                 got_dim2&&
                 got_dim3&&
                 got_veclen&&
                 got_min_ext&&
                 got_max_ext&&
                 got_field&&
                 got_nspace&&
                 got_data))) {

    if(!got_ndim) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"ndim")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ndim)==1)
          got_ndim=true;
      }
    }
      
    if(!got_dim1) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"dim1")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[0])==1)
          got_dim1=true;
      }
    }
      
    if(!got_dim2) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"dim2")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[1])==1)
          got_dim2=true;
      }
    }
      
    if(!got_dim3) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"dim3")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[2])==1)
          got_dim3=true;
      }
    }
       
    if(!got_nspace) {
      ParseWordCopy(cc,p,6);
      if(!strcmp(cc,"nspace")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&nspace)==1)
          got_nspace=true;
      }
    }

    if(!got_veclen) {
      ParseWordCopy(cc,p,6);
      if(!strcmp(cc,"veclen")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&veclen)==1)
          got_veclen=true;
      }
    }

    if(!got_data) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"data")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(!strcmp(cc,"float"))
          got_data=true;
      }
    }

    if(!got_min_ext) {
      ParseWordCopy(cc,p,7);
      if(!strcmp(cc,"min_ext")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%f",&ms->ExtentMin[0])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%f",&ms->ExtentMin[1])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->ExtentMin[2])==1) {
              got_min_ext=true;
            }
          }
        }
      }
    }

    if(!got_max_ext) {
      ParseWordCopy(cc,p,7);
      if(!strcmp(cc,"max_ext")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%f",&ms->ExtentMax[0])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%f",&ms->ExtentMax[1])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->ExtentMax[2])==1) {
              got_max_ext=true;
            }
          }
        }
      }
    }

    if(!got_field) {
      ParseWordCopy(cc,p,5);
      if(!strcmp(cc,"field")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(!strcmp(cc,"uniform"))
          got_field=true;
      }
    }
    p=ParseNextLine(p);
  }
    
  if(got_ndim&&
     got_dim1&&
     got_dim2&&
     got_dim3&&
     got_veclen&&
     got_min_ext&&
     got_max_ext&&
     got_field&&
     got_nspace&&
     got_data) {

    int pass = 0;

    ms->Origin = Alloc(float,3);
    ms->Range = Alloc(float,3);
    ms->Grid = Alloc(float,3);

    copy3f(ms->ExtentMin,ms->Origin);
    subtract3f(ms->ExtentMax,ms->ExtentMin,ms->Range);
    ms->FDim[3] = 3;

    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
      " FLDStrToMap: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
      ENDFB(I->Obj.G);
      
    for(a=0;a<3;a++) {
      ms->Min[a] = 0;
      ms->Max[a] = ms->FDim[a]-1;
      ms->Div[a] = ms->FDim[a]-1;
        
      if(ms->FDim[a]) 
        ms->Grid[a]=ms->Range[a]/(ms->Max[a]);
      else
        ms->Grid[a]=0.0F;
    }

    ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
    ms->MapSource = cMapSourceFLD;
    ms->Field->save_points=false;

    
    while(1) {
      maxd = -FLT_MAX;
      mind = FLT_MAX;
      p=PHIStr;
      
      while(*p) { /* ^L^L sentinel */
        if((p[0]==12)&&(p[1]==12)) {
          p+=2;
          break;
        }
        p++;
      }
      
      rev = (char*)&dens_rev;
      for(c=0;c<ms->FDim[2];c++) { /* z y x ordering into c b a  so that x = a, etc. */
        for(b=0;b<ms->FDim[1];b++) {
          for(a=0;a<ms->FDim[0];a++) {
            
            if(little_endian!=map_endian) {
              rev[0]=p[3];
              rev[1]=p[2];
              rev[2]=p[1];
              rev[3]=p[0];
            } else {
              rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
              rev[1]=p[1];
              rev[2]=p[2];
              rev[3]=p[3];
            }
            dens = *((float*)rev);
            F3(ms->Field->data,a,b,c) = dens;
            if(maxd<dens) maxd = dens;
            if(mind>dens) mind = dens;
            p+=4;
          }
        }
      }

      /* There's no way to determine the original handedness of input
          field files.  So instead, we simplymake an educated guess about
       whether we're byte-swapped based on the range of the density
       values obtained. */
      
      if(((maxd/FLT_MAX)>0.1F)&&((mind/(-FLT_MAX))>0.1F)) {
        if(pass==0) {
          map_endian = (!map_endian); /* okay, try again swapped */
        } else if(pass==1) {
          /* didn't help, so resort to original order */
          map_endian = (!map_endian); 
        } else {
          break;
        }
      } else {
        break;
      }
      pass++;
    }
    
    for(c=0;c<ms->FDim[2];c++) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      for(b=0;b<ms->FDim[1];b++) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        for(a=0;a<ms->FDim[0];a++) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          for(e=0;e<3;e++) {
            F4(ms->Field->points,a,b,c,e) = v[e];
          }
        }
      }
    }
    
    d=0;
    for(c=0;c<ms->FDim[2];c+=ms->FDim[2]-1) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      if(!c) 
        v[2]=ms->ExtentMin[2];
      else 
        v[2]=ms->ExtentMax[2];
      
      for(b=0;b<ms->FDim[1];b+=ms->FDim[1]-1) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);

      if(!b) 
        v[1]=ms->ExtentMin[1];
      else 
        v[1]=ms->ExtentMax[1];
        
        for(a=0;a<ms->FDim[0];a+=ms->FDim[0]-1) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);

      if(!a) 
        v[0]=ms->ExtentMin[0];
      else 
        v[0]=ms->ExtentMax[0];

          copy3f(v,ms->Corner+3*d);
          d++;
        }
      }
    }
    
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);
      
  } else {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors) 
      " Error: unable to read FLD file.\n"
      ENDFB(I->Obj.G);
            /*  printf(" got_ndim %d\n",got_ndim);
    printf(" got_dim1 %d\n",got_dim1);
    printf(" got_dim2 %d\n",got_dim2);
    printf(" got_dim3 %d\n",got_dim3);
    printf(" got_veclen %d\n",got_veclen);
    printf(" got_min_ext %d\n",got_min_ext);
    printf(" got_max_ext %d\n",got_max_ext);
    printf(" got_field %d\n",got_field);
    printf(" got_nspace %d\n",got_nspace);
    printf(" got_data %d\n",got_data);
            */
  }
    

  return(ok);

}
/*========================================================================*/
static int ObjectMapBRIXStrToMap(ObjectMap *I,char *BRIXStr,int bytes,int state) 
{
  char *p,*pp;
  float dens;
  int a,b,c,d,e;
  float v[3],vr[3],maxd,mind;
  int ok = true;
  char cc[MAXLINELEN];
  ObjectMapState *ms;
  int got_origin=false;
  int got_extent=false;
  int got_grid=false;
  int got_cell=false;
  int got_prod=false;
  int got_plus=false;
  int got_sigma=false;
  int got_smiley=false;
  float prod=1.0F;
  float plus=0.0F;
  float sigma=0.0F;
  float calc_sigma=0.0F;
  float calc_mean=0.0F;
  int normalize;
  int swap_bytes;

  normalize=(int)SettingGet(I->Obj.G,cSetting_normalize_o_maps);
  swap_bytes=(int)SettingGet(I->Obj.G,cSetting_swap_dsn6_bytes);
  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);

  maxd = -FLT_MAX;
  mind = FLT_MAX;
  
  p = BRIXStr;

  ParseNCopy(cc,p,3);
  got_smiley = (strcmp(cc,":-)")==0);

  if(got_smiley) {

    /* ASCII BRIX format header */

    while((*p)&&(!(
                   got_origin&&
                   got_extent&&
                   got_grid&&
                   got_cell&&
                   got_prod&&
                   got_plus&&
                   got_sigma))) {
      
      if(!got_origin) {
        pp=ParseWordCopy(cc,p,6);
        if(WordMatch(I->Obj.G,"origin",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%d",&ms->Min[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Min[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Min[2])==1) {
                got_origin=true;
              }
            }
          }
        }
      }
      
      if(!got_extent) {
        pp=ParseWordCopy(cc,p,6);
        if(WordMatch(I->Obj.G,"extent",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%d",&ms->Max[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Max[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Max[2])==1) {
                got_extent=true;
                ms->Max[0]+=ms->Min[0]-1;
                ms->Max[1]+=ms->Min[1]-1;
                ms->Max[2]+=ms->Min[2]-1;
              }
            }
          }
        }
      }
      
      
      if(!got_grid) {
        pp=ParseWordCopy(cc,p,4);
        if(WordMatch(I->Obj.G,"grid",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%d",&ms->Div[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Div[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Div[2])==1) {
                got_grid=true;
              }
            }
          }
        }
      }
      
      if(!got_cell) {
        pp = ParseWordCopy(cc,p,4);
        if(WordMatch(I->Obj.G,"cell",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          
          if(sscanf(cc,"%f",&ms->Crystal->Dim[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->Crystal->Dim[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%f",&ms->Crystal->Dim[2])==1) {
                p = ParseWordCopy(cc,p,50);
                if(sscanf(cc,"%f",&ms->Crystal->Angle[0])==1) {
                  p = ParseWordCopy(cc,p,50);
                  if(sscanf(cc,"%f",&ms->Crystal->Angle[1])==1) {
                    p = ParseWordCopy(cc,p,50);
                    if(sscanf(cc,"%f",&ms->Crystal->Angle[2])==1) {
                      got_cell=true;
                    }
                  }
                }
              }
            }
          }
        }
      }
      
      if(!got_plus) {
        pp=ParseWordCopy(cc,p,4);
        if(WordMatch(I->Obj.G,"plus",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%f",&plus)==1) {
            got_plus=true;
          }
        }
      }
      
      if(!got_prod) {
        pp=ParseWordCopy(cc,p,4);
        if(WordMatch(I->Obj.G,"prod",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%f",&prod)==1) {
            got_prod=true;
          }
        }
      }
      
      if(!got_sigma) {
        pp = ParseWordCopy(cc,p,5);
        if(WordMatch(I->Obj.G,"sigma",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%f",&sigma)==1) {
            got_sigma=true;
          }
        }
      }
      
      p++;
    }
  } else {
    /* Binary DSN6 format */

    /* first, determine whether or not we need to byte-swap the header... */

    int passed_endian_check = false;
    short int *shint_ptr;
    float scale_factor;
    char tmp_char;
    int a;
    int start_swap_at = 0;
    int stop_swap_at = bytes;

    shint_ptr = (short int*)(p+18*2);
    if(*shint_ptr==100) {
      passed_endian_check = true;
      start_swap_at = 512; /* don't byte-swap header */
    }
    
    if(!swap_bytes) {
      stop_swap_at = 512;
    }
    /* for DSN6, always swap the bricks */

    p = BRIXStr + start_swap_at;
    for( a=start_swap_at; a<stop_swap_at; a+=2) {
      tmp_char = *p;
      *p = *(p+1);
      *(p+1) = tmp_char;
      p+=2;
    }

    p = BRIXStr;

    if(*shint_ptr==100) {
      passed_endian_check = true;
    }
    
    if(!passed_endian_check) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors) 
        " Error: This looks like a DSN6 map file, but I can't match endianness.\n"
        ENDFB(I->Obj.G);
    } else {
      shint_ptr = (short int*)p;
      
      ms->Min[0] = *(shint_ptr++);
      ms->Min[1] = *(shint_ptr++);
      ms->Min[2] = *(shint_ptr++);
      got_origin = true;
      
      ms->Max[0] = *(shint_ptr++);
      ms->Max[1] = *(shint_ptr++);
      ms->Max[2] = *(shint_ptr++);
      got_extent=true;
      ms->Max[0]+=ms->Min[0]-1;
      ms->Max[1]+=ms->Min[1]-1;
      ms->Max[2]+=ms->Min[2]-1;
      
      ms->Div[0] = *(shint_ptr++);
      ms->Div[1] = *(shint_ptr++);
      ms->Div[2] = *(shint_ptr++);
      got_grid = true;
      
      ms->Crystal->Dim[0] = (float)(*(shint_ptr++));
      ms->Crystal->Dim[1] = (float)(*(shint_ptr++));
      ms->Crystal->Dim[2] = (float)(*(shint_ptr++));
      ms->Crystal->Angle[0] = (float)(*(shint_ptr++));
      ms->Crystal->Angle[1] = (float)(*(shint_ptr++));
      ms->Crystal->Angle[2] = (float)(*(shint_ptr++));
      got_cell = true;
      
      prod = (float)(*(shint_ptr++)) / 100.0F;
      got_prod = true;
      
      plus = (float)(*(shint_ptr++));
      got_plus = true;
      
      scale_factor = (float)(*(shint_ptr++));
      
      for(a=0;a<3;a++) {
        ms->Crystal->Dim[a]/=scale_factor;
        ms->Crystal->Angle[a]/= scale_factor;
      }
    }
  }
  
  if(got_origin&&
     got_extent&&
     got_grid&&
     got_cell&&
     got_plus&&
     got_prod) {

    
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Blather) 
      " BRIXStrToMap: Prod = %8.3f, Plus = %8.3f\n",prod,plus
      ENDFB(I->Obj.G);

    ms->FDim[0]=ms->Max[0]-ms->Min[0]+1;
    ms->FDim[1]=ms->Max[1]-ms->Min[1]+1;
    ms->FDim[2]=ms->Max[2]-ms->Min[2]+1;
    ms->FDim[3]=3;
    if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
      ok=false;
    else {
      CrystalUpdate(ms->Crystal);
      ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
      ms->MapSource = cMapSourceBRIX;
      ms->Field->save_points=false;
      
      {
        int block_size = 8;
        int xblocks = ((ms->FDim[0]-1)/block_size)+1;
        int yblocks = ((ms->FDim[1]-1)/block_size)+1;
        int zblocks = ((ms->FDim[2]-1)/block_size)+1;
        int cc,bb,aa,ia,ib,ic,xa,xb,xc;
        double sum = 0.0;
        double sumsq = 0.0;
        int n_pts = 0;

        p=BRIXStr + 512;
        for(cc=0;cc<zblocks;cc++) {
          for(bb=0;bb<yblocks;bb++) {
            for(aa=0;aa<xblocks;aa++) {
              
              for(c=0;c<block_size;c++) {
                xc = c + cc*block_size;
                ic = xc + ms->Min[2];
                v[2]=ic/((float)ms->Div[2]);
                
                for(b=0;b<block_size;b++) {
                  xb = b + bb*block_size;
                  ib = xb + ms->Min[1];
                  v[1]=ib/((float)ms->Div[1]);
                  
                  for(a=0;a<block_size;a++) {
                    xa = a + aa*block_size;
                    ia = xa + ms->Min[0];
                    v[0]=ia/((float)ms->Div[0]);
                    
                    dens = (((float)(*((unsigned char*)(p++)))) - plus) / prod;
                    if((ia<=ms->Max[0])&&(ib<=ms->Max[1])&&(ic<=ms->Max[2])) {
                      F3(ms->Field->data,xa,xb,xc) = dens;
                      sumsq+=dens*dens;
                      sum+=dens;
                      n_pts++;
                      if(maxd<dens) maxd = dens;
                      if(mind>dens) mind = dens;
                      transform33f3f(ms->Crystal->FracToReal,v,vr);
                      for(e=0;e<3;e++) {
                        F4(ms->Field->points,xa,xb,xc,e) = vr[e];
                      }
                      
                    }
                  }
                }
              }
            }
          }
        }

        if(n_pts>1) {
          calc_mean = (float)(sum/n_pts);
          calc_sigma = (float)sqrt1d((sumsq - (sum*sum/n_pts))/(n_pts-1));
        }
        
      }
      
      if(ok) {
        d = 0;
        for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
          {
            v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
            for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
              v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
              for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
                v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                transform33f3f(ms->Crystal->FracToReal,v,vr);
                copy3f(vr,ms->Corner+3*d);
                d++;
              }
            }
          }
      }

      if(ok&&normalize) {

        if(calc_sigma>R_SMALL8) {
          for(c=0;c<ms->FDim[2];c++) {
            for(b=0;b<ms->FDim[1];b++) {
              for(a=0;a<ms->FDim[0];a++) {
                F3(ms->Field->data,a,b,c) =  (F3(ms->Field->data,a,b,c)-calc_mean)/calc_sigma;
              }
            }
          }
        }
      }

      v[2]=(ms->Min[2])/((float)ms->Div[2]);
      v[1]=(ms->Min[1])/((float)ms->Div[1]);
      v[0]=(ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
      
      v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
      v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
      v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
      
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
        " BRIXStrToMap: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
        ENDFB(I->Obj.G);

      if(got_sigma) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
          " BRIXStrToMap: Reported Sigma = %8.3f\n",sigma
          ENDFB(I->Obj.G);
      }

      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
        " BRIXStrToMap: Range = %5.6f to %5.6f\n",mind,maxd
        ENDFB(I->Obj.G);

      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
        " BRIXStrToMap: Calculated Mean = %8.3f, Sigma = %8.3f\n",calc_mean,calc_sigma
        ENDFB(I->Obj.G);

      if(normalize) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
          " BRIXStrToMap: Normalizing...\n"
          ENDFB(I->Obj.G);
      }

      ms->Active=true;
      ObjectMapUpdateExtents(I);

    }
  } else {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors) 
      " Error: unable to read BRIX/DSN6 file.\n"
      ENDFB(I->Obj.G);
  }
  
  
  return(ok);
  
}
/*========================================================================*/
static int ObjectMapGRDStrToMap(ObjectMap *I,char *GRDStr,int bytes,int state) 
{
  /* NOTE: binary GRD reader not yet validated */

  char *p;
  float dens;
  float *f = NULL;
  int a,b,c,d,e;
  float v[3],vr[3],maxd,mind;
  int ok = true;
  char cc[MAXLINELEN];
  ObjectMapState *ms;
  int got_cell=false;
  int got_brick=false;
  int fast_axis=1;
  int got_ranges=false;
  int normalize = false;
  float mean,stdev;
  double sum = 0.0;
  double sumsq = 0.0;
  int n_pts = 0;
  int ascii = true;
  int little_endian = 1,map_endian = 0;

  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);
  normalize=(int)SettingGet(I->Obj.G,cSetting_normalize_grd_maps);
  maxd = -FLT_MAX;
  mind = FLT_MAX;
  
  p = GRDStr;

  if(bytes>4) {
    if((!p[0])||(!p[1])||(!p[2])||(!p[3])) {
      ascii=false;
      
      little_endian = *((char*)&little_endian);
      map_endian = (*p||*(p+1));
           
    }
  }
  /* print map title */

  if(ascii) {
    p = ParseNCopy(cc,p,100);
  } else {
    int header_len;
    char rev[4];
    if(little_endian!=map_endian) {
      rev[0]=p[3];
      rev[1]=p[2];
      rev[2]=p[1];
      rev[3]=p[0];
    } else {
      rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
      rev[1]=p[1];
      rev[2]=p[2];
      rev[3]=p[3];
    }
    header_len = *((int*)rev);
    ParseNCopy(cc,p+4,header_len);
    /* now flip file to correct endianess */
    if(little_endian!=map_endian) {
      char *c = p;
      int cnt = bytes>>2;
      unsigned char c0,c1,c2,c3;
      while(cnt) {
        c0 = c[0];
        c1 = c[1];
        c2 = c[2];
        c3 = c[3];
        c[0]=c3;
        c[1]=c2;
        c[2]=c1;
        c[3]=c0;
        c+=4;
        cnt--;
      }
    }
    p += 4 + header_len + 4;
    f = (float*)p;
  }

  PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
    " ObjectMap: %s\n",cc
    ENDFD;
  
  if(ascii) 
    p = ParseNextLine(p);

  /* skip format -- we're reading float regardless... */
  if(ascii) 
    p = ParseNextLine(p);
  else
    f += 4; /* what are these numbers? */
    
  /* read unit cell */

  if(ok) {
    if(ascii) {
      p = ParseWordCopy(cc,p,50);
      if(sscanf(cc,"%f",&ms->Crystal->Dim[0])==1) {
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%f",&ms->Crystal->Dim[1])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%f",&ms->Crystal->Dim[2])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->Crystal->Angle[0])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%f",&ms->Crystal->Angle[1])==1) {
                p = ParseWordCopy(cc,p,50);
                if(sscanf(cc,"%f",&ms->Crystal->Angle[2])==1) {
                  got_cell=true;
                }
              }
            }
          }
        }
      }
      ok = got_cell;
    } else {
      ms->Crystal->Dim[0] = *(f++);
      ms->Crystal->Dim[1] = *(f++);
      ms->Crystal->Dim[2] = *(f++);
      ms->Crystal->Angle[0] = (*f++);
      ms->Crystal->Angle[1] = (*f++);
      ms->Crystal->Angle[2] = (*f++);
      got_cell = 1;
    }
  }

  if(ascii) 
    p=ParseNextLine(p);

  if(ascii) {
    /* read brick dimensions */
    
    if(ok) {
      p = ParseWordCopy(cc,p,50);
      if(sscanf(cc,"%d",&ms->FDim[0])==1) {
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[1])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%d",&ms->FDim[2])==1) {
            got_brick=true;
            ms->FDim[0]++;
            ms->FDim[1]++;
            ms->FDim[2]++; /* stupid 0-based counters */
          }
        }
      }
      ok = got_brick;
    }
    
    p=ParseNextLine(p);
    
    /* read ranges */
    
    if(ok) {
      p = ParseWordCopy(cc,p,50);
      if(sscanf(cc,"%d",&fast_axis)==1) {
        p = ParseWordCopy(cc,p,50);      
        if(sscanf(cc,"%d",&ms->Min[0])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%d",&ms->Div[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Min[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Div[1])==1) {
                p = ParseWordCopy(cc,p,50);
                if(sscanf(cc,"%d",&ms->Min[2])==1) {
                  p = ParseWordCopy(cc,p,50);
                  if(sscanf(cc,"%d",&ms->Div[2])==1) {
                    got_ranges=true;
                  }
                }
              }
            }
          }
        }
      }
      ok = got_ranges;
    }
  } else {
    float fmin[3];
    float fmax[3];

    fmin[0] = *(f++);
    fmax[0] = *(f++);
    fmin[1] = *(f++);
    fmax[1] = *(f++);
    fmin[2] = *(f++);
    fmax[2] = *(f++);
    ms->FDim[0] = *((int*)f++) + 1;
    ms->FDim[1] = *((int*)f++) + 1;
    ms->FDim[2] = *((int*)f++) + 1;

    ms->Div[0] = pymol_roundf((ms->FDim[0]-1)/(fmax[0]-fmin[0]));
    ms->Div[1] = pymol_roundf((ms->FDim[1]-1)/(fmax[0]-fmin[1]));
    ms->Div[2] = pymol_roundf((ms->FDim[2]-1)/(fmax[0]-fmin[2]));
    
    ms->Min[0] = pymol_roundf(ms->Div[0]*fmin[0]);
    ms->Min[1] = pymol_roundf(ms->Div[1]*fmin[1]);
    ms->Min[2] = pymol_roundf(ms->Div[2]*fmin[2]);

    ms->Max[0] = ms->Min[0]+ms->FDim[0]-1;
    ms->Max[1] = ms->Min[1]+ms->FDim[1]-1;
    ms->Max[2] = ms->Min[2]+ms->FDim[2]-1;
    got_ranges = true;

    /* assumes fast X */

    f+=1; /* advance to data */

  }
  
  if(ok) {

    if(ascii) {
      ms->Div[0] = ms->Div[0] - ms->Min[0];
      ms->Div[1] = ms->Div[1] - ms->Min[1];
      ms->Div[2] = ms->Div[2] - ms->Min[2];
      
      ms->Max[0] = ms->Min[0]+ms->FDim[0]-1;
      ms->Max[1] = ms->Min[1]+ms->FDim[1]-1;
      ms->Max[2] = ms->Min[2]+ms->FDim[2]-1;
    }

    ms->FDim[3]=3;

    if(Feedback(I->Obj.G,FB_ObjectMap,FB_Blather)) {
      dump3i(ms->Div,"ms->Div");
      dump3i(ms->Min,"ms->Min");
      dump3i(ms->Max,"ms->Max");
      dump3i(ms->FDim,"ms->FDim");
    }

    CrystalUpdate(ms->Crystal);
    ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
    ms->MapSource = cMapSourceGRD;
    ms->Field->save_points=false;

    switch(fast_axis) {
    case 3: /* Fast Y - BROKEN! */
    default:
    case 1: /* Fast X */
      for(c=0;c<ms->FDim[2];c++) {
        v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
        for(b=0;b<ms->FDim[1];b++) {
          if(!ascii)
            f++; /* skip block delimiter */
          v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
          for(a=0;a<ms->FDim[0];a++) {
            v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
            if(ascii) {
              p=ParseNextLine(p);
              p=ParseNCopy(cc,p,24);                
              if(sscanf(cc,"%f",&dens)!=1) {
                ok=false;
              }
            } else {
              dens = *(f++);
            }
            if(ok) {
              sumsq+=dens*dens;
              sum+=dens;
              n_pts++;
              F3(ms->Field->data,a,b,c) = dens;
              if(maxd<dens) maxd = dens;
              if(mind>dens) mind = dens;
            }
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            for(e=0;e<3;e++) {
              F4(ms->Field->points,a,b,c,e) = vr[e];
            }
          }
          if(!ascii) f++; /* skip fortran block delimiter */
        }
      }
      break;
    }
  }   

  if(n_pts>1) {
    mean = (float)(sum/n_pts);
    stdev = (float)sqrt1d((sumsq - (sum*sum/n_pts))/(n_pts-1));

    if(normalize) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
        " ObjectMapGRDStrToMap: Normalizing: mean = %8.6f and stdev = %8.6f.\n",
        mean,stdev
        ENDFB(I->Obj.G);
      if(stdev<R_SMALL8)
        stdev = 1.0F;
      for(c=0;c<ms->FDim[2];c++)
        for(b=0;b<ms->FDim[1];b++) {
          for(a=0;a<ms->FDim[0];a++) {
            dens = F3(ms->Field->data,a,b,c);
            F3(ms->Field->data,a,b,c) = (dens-mean)/stdev;
          }
        }
    } else {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
        " ObjectMapGRDStrToMap: Mean = %8.6f and stdev = %8.6f.\n",
        mean,stdev
        ENDFB(I->Obj.G);
    }
  }

  if(ok) {
    d = 0;
    for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
      {
        v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
        for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
          v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
          for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
            v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            copy3f(vr,ms->Corner+3*d);
            d++;
          }
        }
      }
    
    v[2]=(ms->Min[2])/((float)ms->Div[2]);
    v[1]=(ms->Min[1])/((float)ms->Div[1]);
    v[0]=(ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
    
    v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
    v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
    v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
    
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details) 
      " GRDXStrToMap: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
      ENDFB(I->Obj.G);
    
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);

  } else {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Errors) 
      " Error: unable to read GRD file.\n"
      ENDFB(I->Obj.G);
  }
  
  return(ok);
  
}
/*========================================================================*/
static ObjectMap *ObjectMapReadXPLORStr(PyMOLGlobals *G,ObjectMap *I,char *XPLORStr,int state,int quiet)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapXPLORStrToMap(I,XPLORStr,state,quiet);
    
    
    SceneChanged(I->Obj.G);
    SceneCountFrames(I->Obj.G);
  }
  return(I);
}
/*========================================================================*/
static ObjectMap *ObjectMapReadCCP4Str(PyMOLGlobals *G,ObjectMap *I,char *XPLORStr,int bytes,int state,int quiet)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapCCP4StrToMap(I,XPLORStr,bytes,state,quiet);
    SceneChanged(G);
    SceneCountFrames(G);
  }
  return(I);
}
/*========================================================================*/
ObjectMap *ObjectMapLoadCCP4(PyMOLGlobals *G,ObjectMap *obj,char *fname,int state,
                             int is_string,int bytes,int quiet)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f = NULL;
  char *buffer,*p;
  long size;

  if(!is_string) {
    
    f=fopen(fname,"rb");
    if(!f)
      ok=ErrMessage(G,"ObjectMapLoadCCP4File","Unable to open file!");
  } 
  
  if(f || is_string) {
    
    if(!quiet) {
      if((!is_string) && Feedback(G,FB_ObjectMap,FB_Actions)) {
        printf(" ObjectMapLoadCCP4File: Loading from '%s'.\n",fname);
      }
    }
    
    if(!is_string) {
      fseek(f,0,SEEK_END);
      size=ftell(f);
      fseek(f,0,SEEK_SET);
      
            buffer=(char*)mmalloc(size);
            ErrChkPtr(G,buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);
    } else {
      buffer = fname;
      size = (long)bytes;
    }

    I=ObjectMapReadCCP4Str(G,obj,buffer,size,state,quiet);
    
    if(!is_string) 
      mfree(buffer);

    if(!quiet) {
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          CrystalDump(ms->Crystal);
        }
      }
    }
  }
  return(I);
}
/*========================================================================*/
static ObjectMap *ObjectMapReadFLDStr(PyMOLGlobals *G,ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapFLDStrToMap(I,MapStr,bytes,state);
    SceneChanged(G);
    SceneCountFrames(G);
  }
  return(I);
}

/*========================================================================*/
static ObjectMap *ObjectMapReadBRIXStr(PyMOLGlobals *G,ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapBRIXStrToMap(I,MapStr,bytes,state);
    SceneChanged(G);
    SceneCountFrames(G);
  }
  return(I);
}
/*========================================================================*/
static ObjectMap *ObjectMapReadGRDStr(PyMOLGlobals *G,ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapGRDStrToMap(I,MapStr,bytes,state);
    SceneChanged(G);
    SceneCountFrames(G);
  }
  return(I);
}

/*========================================================================*/
static ObjectMap *ObjectMapReadPHIStr(PyMOLGlobals *G,ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapPHIStrToMap(I,MapStr,bytes,state);
    SceneChanged(G);
    SceneCountFrames(G);
  }
  return(I);
}
/*========================================================================*/
ObjectMap *ObjectMapLoadPHIFile(PyMOLGlobals *G,ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
       ok=ErrMessage(G,"ObjectMapLoadPHIFile","Unable to open file!");
  else
       {
            if(Feedback(G,FB_ObjectMap,FB_Actions))
              {
                  printf(" ObjectMapLoadPHIFile: Loading from '%s'.\n",fname);
              }
            
            fseek(f,0,SEEK_END);
      size=ftell(f);
            fseek(f,0,SEEK_SET);

            buffer=(char*)mmalloc(size);
            ErrChkPtr(G,buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);

            I=ObjectMapReadPHIStr(G,obj,buffer,size,state);

            mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}
/*========================================================================*/

static int is_number(char *p)
{
  int result = (*p != 0);
  register char c;
  if(result) 
    while( (c = *(p++)) ) {
      if(! ((c == '.') ||
            (c == '-') ||
            (c == '+') ||
            (c == 'e') ||
            (c == 'E') ||
            ((c>='0')&&(c<='9'))))
        result = false;
      break;
    }
  return result;
}

static int ObjectMapDXStrToMap(ObjectMap *I,char *DXStr,int bytes,int state) {

  int n_items = 0;

  char *p,*pp;
  float dens;
  int a,b,c,d,e;
  float v[3],maxd,mind;
  int ok = true;
  /* DX named from their docs */

  int stage = 0;

  ObjectMapState *ms;

  char cc[MAXLINELEN];


  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(I->Obj.G,ms);

  ms->Origin=Alloc(float,3);
  ms->Grid = Alloc(float,3);

  maxd = -FLT_MAX;
  mind = FLT_MAX;
  p=DXStr;

  /* get the dimensions */

  ms->FDim[3] = 3;

  while(ok&&(*p)&&(stage==0)) {
    pp = p;
    p = ParseNCopy(cc,p,35);
    if((strcmp(cc,"object 1 class gridpositions counts")==0) || is_number(cc)) {
      if(is_number(cc)) 
        p = pp;
      p = ParseWordCopy(cc,p,10);      
      if(sscanf(cc,"%d",&ms->FDim[0])==1) {
        p = ParseWordCopy(cc,p,10);      
        if(sscanf(cc,"%d",&ms->FDim[1])==1) {
          p = ParseWordCopy(cc,p,10);      
          if(sscanf(cc,"%d",&ms->FDim[2])==1) {
            stage = 1;
          }
        }
      }
    }
    p = ParseNextLine(p);
  }

  if(ok&&(stage==1)) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
      " DXStrToMap: Dimensions: %d %d %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
      ENDFB(I->Obj.G);
  }

  /* get the origin */

  while(ok&&(*p)&&(stage==1)) {  
    pp = p;
    p = ParseNCopy(cc,p,6);
    if((strcmp(cc,"origin")==0) || is_number(cc)) {
      if(is_number(cc)) 
        p = pp;
      p = ParseWordCopy(cc,p,20);      
      if(sscanf(cc,"%f",&ms->Origin[0])==1) {
        p = ParseWordCopy(cc,p,20);      
        if(sscanf(cc,"%f",&ms->Origin[1])==1) {
          p = ParseWordCopy(cc,p,20);      
          if(sscanf(cc,"%f",&ms->Origin[2])==1) {
            stage = 2;
          }
        }
      }
    }
    p = ParseNextLine(p);
  }

  if(ok&&(stage==2)) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
      " DXStrToMap: Origin %8.3f %8.3f %8.3f\n",ms->Origin[0],ms->Origin[1],ms->Origin[2]
      ENDFB(I->Obj.G);
  }

  while(ok&&(*p)&&(stage==2)) {  
    pp = p;
    p = ParseNCopy(cc,p,5);
    if(strcmp(cc,"delta")==0) {
      p = ParseWordCopy(cc,p,20);      
      if(sscanf(cc,"%f",&ms->Grid[0])==1) {
        p = ParseNextLine(p);
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        if(sscanf(cc,"%f",&ms->Grid[1])==1) {
        p = ParseNextLine(p);
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          if(sscanf(cc,"%f",&ms->Grid[2])==1) {
            stage = 3;
          }
        }
      }
    } else if(is_number(cc)) {
      p = pp;
      p = ParseWordCopy(cc,p,20);      
      if(sscanf(cc,"%f",&ms->Grid[0])==1) {
        p = ParseNextLine(p);
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        if(sscanf(cc,"%f",&ms->Grid[1])==1) {
          p = ParseNextLine(p);
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          if(sscanf(cc,"%f",&ms->Grid[2])==1) {
            stage = 3;
          }
        }
      }
    }
  }

  if(ok&&(stage==3)) {
    PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
      " DXStrToMap: Grid %8.3f %8.3f %8.3f\n",ms->Grid[0],ms->Grid[1],ms->Grid[2]
      ENDFB(I->Obj.G);
  }

  while(ok&&(*p)&&(stage==3)) {  
    p = ParseNCopy(cc,p,6);
    if(strcmp(cc,"object")==0) {
      p = ParseWordCopy(cc,p,20);      
      p = ParseNTrim(cc,p,29);      
      if(strcmp(cc,"class array type double rank")==0) {
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        if(sscanf(cc,"%d",&n_items)==1) {
          if(n_items == ms->FDim[0]*ms->FDim[1]*ms->FDim[2])
            stage = 4;
        }
      }
    } else if(is_number(cc)) {
      n_items = ms->FDim[0]*ms->FDim[1]*ms->FDim[2];
      stage = 4;
      break;
    }
    p = ParseNextLine(p);
  }

  if(stage == 4) {

    if(ok&&(stage==4)) {
      PRINTFB(I->Obj.G,FB_ObjectMap,FB_Details)
        " DXStrToMap: %d data points.\n",n_items
        ENDFB(I->Obj.G);
    }
    
    ms->Field=IsosurfFieldAlloc(I->Obj.G,ms->FDim);
    ms->MapSource = cMapSourcePHI;
    ms->Field->save_points=false;
    
    for(a=0;a<3;a++) {
      ms->Div[a] = ms->FDim[a] - 1;
      ms->Min[a] = 0;
      ms->Max[a] = ms->FDim[a] - 1;
    }
    
    for(a=0;a<ms->FDim[0];a++) {
      for(b=0;b<ms->FDim[1];b++) {
        for(c=0;c<ms->FDim[2];c++) { 
          
          p = ParseWordCopy(cc,p,20);      
          if(!cc[0]) {
            p = ParseNextLine(p);           
            p = ParseWordCopy(cc,p,20);      
          }
          if(sscanf(cc,"%f",&dens)==1) {
            if(maxd<dens) maxd = dens;
            if(mind>dens) mind = dens;
            F3(ms->Field->data,a,b,c) = dens;
          } else {
            ok=false;
          }
        }
      }

    }
    
    for(e=0;e<3;e++) {
      ms->ExtentMin[e] = ms->Origin[e]+ms->Grid[e]*ms->Min[e];
      ms->ExtentMax[e] = ms->Origin[e]+ms->Grid[e]*ms->Max[e];
    }
    
    for(c=0;c<ms->FDim[2];c++) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      for(b=0;b<ms->FDim[1];b++) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        for(a=0;a<ms->FDim[0];a++) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          for(e=0;e<3;e++) {
            F4(ms->Field->points,a,b,c,e) = v[e];
          }
        }
      }
    }
    
    d=0;
    for(c=0;c<ms->FDim[2];c+=ms->FDim[2]-1) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      
      for(b=0;b<ms->FDim[1];b+=ms->FDim[1]-1) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        
        for(a=0;a<ms->FDim[0];a+=ms->FDim[0]-1) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          copy3f(v,ms->Corner+3*d);
          d++;
        }
      }
    }
    if(ok)
      stage = 5;
  }

  if(stage!=5)
    ok=false;

  if(!ok) {
    ErrMessage(I->Obj.G,"ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);
  }
  return(ok);
}

/*========================================================================*/
static ObjectMap *ObjectMapReadDXStr(PyMOLGlobals *G,ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapDXStrToMap(I,MapStr,bytes,state);
    SceneChanged(G);
    SceneCountFrames(G);
  }
  return(I);
}

/*========================================================================*/
ObjectMap *ObjectMapLoadDXFile(PyMOLGlobals *G,ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
       ok=ErrMessage(G,"ObjectMapLoadDXFile","Unable to open file!");
  else
       {
            if(Feedback(G,FB_ObjectMap,FB_Actions))
              {
                  printf(" ObjectMapLoadDXFile: Loading from '%s'.\n",fname);
              }
            
            fseek(f,0,SEEK_END);
      size=ftell(f);
            fseek(f,0,SEEK_SET);

            buffer=(char*)mmalloc(size);
            ErrChkPtr(G,buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);

            I=ObjectMapReadDXStr(G,obj,buffer,size,state);

            mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}

/*========================================================================*/
ObjectMap *ObjectMapLoadFLDFile(PyMOLGlobals *G,ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
       ok=ErrMessage(G,"ObjectMapLoadFLDFile","Unable to open file!");
  else
       {
            if(Feedback(G,FB_ObjectMap,FB_Actions))
              {
                  printf(" ObjectMapLoadFLDFile: Loading from '%s'.\n",fname);
              }
            
            fseek(f,0,SEEK_END);
      size=ftell(f);
            fseek(f,0,SEEK_SET);

            buffer=(char*)mmalloc(size);
            ErrChkPtr(G,buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);

            I=ObjectMapReadFLDStr(G,obj,buffer,size,state);

            mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}

/*========================================================================*/
ObjectMap *ObjectMapLoadBRIXFile(PyMOLGlobals *G,ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f = NULL;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
    ok=ErrMessage(G,"ObjectMapLoadBRIXFile","Unable to open file!");
  if(ok)
       {
      if(Feedback(G,FB_ObjectMap,FB_Actions))
        {
          printf(" ObjectMapLoadBRIXFile: Loading from '%s'.\n",fname);
        }
      
      fseek(f,0,SEEK_END);
      size=ftell(f);
      fseek(f,0,SEEK_SET);
      
      buffer=(char*)mmalloc(size+255);
      ErrChkPtr(G,buffer);
      p=buffer;
      fseek(f,0,SEEK_SET);
      fread(p,size,1,f);
      p[size]=0;
      fclose(f);
      
            I=ObjectMapReadBRIXStr(G,obj,buffer,size,state);

      mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          CrystalDump(ms->Crystal);
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}
/*========================================================================*/
ObjectMap *ObjectMapLoadGRDFile(PyMOLGlobals *G,ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f = NULL;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
    ok=ErrMessage(G,"ObjectMapLoadGRDFile","Unable to open file!");
  if(ok)
       {
      if(Feedback(G,FB_ObjectMap,FB_Actions))
        {
          printf(" ObjectMapLoadGRDFile: Loading from '%s'.\n",fname);
        }
      
      fseek(f,0,SEEK_END);
      size=ftell(f);
      fseek(f,0,SEEK_SET);
      
      buffer=(char*)mmalloc(size+255);
      ErrChkPtr(G,buffer);
      p=buffer;
      fseek(f,0,SEEK_SET);
      fread(p,size,1,f);
      p[size]=0;
      fclose(f);
      
            I=ObjectMapReadGRDStr(G,obj,buffer,size,state);

      mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          CrystalDump(ms->Crystal);
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}
/*========================================================================*/
ObjectMap *ObjectMapLoadXPLOR(PyMOLGlobals *G,ObjectMap *obj,char *fname,
                                  int state,int is_file,int quiet)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f = NULL;
  long size;
  char *buffer,*p;

  if(is_file) {
    f=fopen(fname,"rb");
    if(!f)
      ok=ErrMessage(G,"ObjectMapLoadXPLOR","Unable to open file!");
  }
  if(ok) {
    if((!quiet) && (Feedback(G,FB_ObjectMap,FB_Actions))) {
      if(is_file) {
        printf(" ObjectMapLoadXPLOR: Loading from '%s'.\n",fname);
      } else {
        printf(" ObjectMapLoadXPLOR: Loading...\n");
      }
    }
       
    if(is_file) {
      fseek(f,0,SEEK_END);
      size=ftell(f);
      fseek(f,0,SEEK_SET);
        
      buffer=(char*)mmalloc(size+255);
      ErrChkPtr(G,buffer);
      p=buffer;
      fseek(f,0,SEEK_SET);
      fread(p,size,1,f);
      p[size]=0;
      fclose(f);
    } else {
      buffer = fname;
    }
      
    I=ObjectMapReadXPLORStr(G,obj,buffer,state,quiet);

    if(is_file) 
      mfree(buffer);
    if(!quiet) {
      if(Feedback(G,FB_ObjectMap, FB_Details)) {
        if(state<0) 
          state=I->NState-1;
          
        if(state<I->NState) {
          ObjectMapState *ms;
          ms = &I->State[state];
          if(ms->Active) {
            CrystalDump(ms->Crystal);
          }
        }
      }
    }
  }
  return(I);

}
/*========================================================================*/
int ObjectMapSetBorder(ObjectMap *I,float level,int state)
{
  int a;
  int result=true;
  if(state==-2)
    state = ObjectGetCurrentState(&I->Obj,false);
  for(a=0;a<I->NState;a++) {
    if((state<0) || (state==a)) {
      if(I->State[a].Active)
        result = result && ObjectMapStateSetBorder(&I->State[a],level);
    }
  }
  return(result);
}
/*========================================================================*/
int ObjectMapNumPyArrayToMapState(PyMOLGlobals *G,ObjectMapState *ms,PyObject *ary) {

#ifdef _PYMOL_NOPY
  return 0;
#else
  
  int a,b,c,d,e;
  float v[3],dens,maxd,mind;
  int ok = true;
#ifdef _PYMOL_NUMPY
  MyArrayObject *pao;
#endif
  
  
#ifdef _PYMOL_NUMPY
  pao = (MyArrayObject*)ary;
#endif
  maxd = -FLT_MAX;
  mind = FLT_MAX;
  if(ok) {
    ms->FDim[0]=ms->Dim[0];
    ms->FDim[1]=ms->Dim[1];
    ms->FDim[2]=ms->Dim[2];
    ms->FDim[3]=3;


    if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
      ok=false;
    else {
      ms->Field=IsosurfFieldAlloc(G,ms->FDim);
      for(c=0;c<ms->FDim[2];c++)
        {
          v[2]=ms->Origin[2]+ms->Grid[2]*c;
          for(b=0;b<ms->FDim[1];b++) {
            v[1]=ms->Origin[1]+ms->Grid[1]*b;
            for(a=0;a<ms->FDim[0];a++) {
              v[0]=ms->Origin[0]+ms->Grid[0]*a;
#ifdef _PYMOL_NUMPY
              dens = (float)(*((double*)
                (pao->data+
                 (pao->strides[0]*a)+
                 (pao->strides[1]*b)+
                 (pao->strides[2]*c))));
#else
              dens = 0.0;
#endif
              F3(ms->Field->data,a,b,c) = dens;
              if(maxd<dens) maxd = dens;
              if(mind>dens) mind = dens;
              for(e=0;e<3;e++) 
                F4(ms->Field->points,a,b,c,e) = v[e];
            }
          }
        }
      d = 0;
      for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
        {
          v[2]=ms->Origin[2]+ms->Grid[2]*c;
          for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
            v[1]=ms->Origin[1]+ms->Grid[1]*b;
            for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
              v[0]=ms->Origin[0]+ms->Grid[0]*a;
              copy3f(v,ms->Corner+3*d);
              d++;
            }
          }
        }
    }
  }
  if(ok) { 
    copy3f(ms->Origin,ms->ExtentMin);
    copy3f(ms->Origin,ms->ExtentMax);
    add3f(ms->Range,ms->ExtentMax,ms->ExtentMax);
  }
  if(!ok) {
    ErrMessage(G,"ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    if(Feedback(G,FB_ObjectMap,FB_Actions)) {
      printf(" ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd);
    }
  }
  return(ok);
#endif
}
/*========================================================================*/
ObjectMap *ObjectMapLoadChemPyBrick(PyMOLGlobals *G,ObjectMap *I,PyObject *Map,
                                           int state,int discrete)
{
#ifdef _PYMOL_NOPY
  return NULL;
#else

  int ok=true;
  int isNew = true;
  PyObject *tmp;
  ObjectMapState *ms;


  if(!I) 
       isNew=true;
  else 
       isNew=false;

  if(ok) {

       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }

    if(state<0) state=I->NState;
    if(I->NState<=state) {
      VLACheck(I->State,ObjectMapState,state);
      I->NState=state+1;
    }
    ms=&I->State[state];
    ObjectMapStateInit(G,ms);

    if(PyObject_HasAttrString(Map,"origin")&&
       PyObject_HasAttrString(Map,"dim")&&
       PyObject_HasAttrString(Map,"range")&&
       PyObject_HasAttrString(Map,"grid")&&
       PyObject_HasAttrString(Map,"lvl"))
      {
        tmp = PyObject_GetAttrString(Map,"origin");
        if(tmp) {
          PConvPyListToFloatArray(tmp,&ms->Origin);
          Py_DECREF(tmp);
        } else 
          ok=ErrMessage(G,"ObjectMap","missing brick origin.");
        tmp = PyObject_GetAttrString(Map,"dim");
        if(tmp) {
          PConvPyListToIntArray(tmp,&ms->Dim);
          Py_DECREF(tmp);
        } else 
          ok=ErrMessage(G,"ObjectMap","missing brick dimension.");
        tmp = PyObject_GetAttrString(Map,"range");
        if(tmp) {
          PConvPyListToFloatArray(tmp,&ms->Range);
          Py_DECREF(tmp);
        } else 
          ok=ErrMessage(G,"ObjectMap","missing brick range.");
        tmp = PyObject_GetAttrString(Map,"grid");
        if(tmp) {
          PConvPyListToFloatArray(tmp,&ms->Grid);
          Py_DECREF(tmp);
        } else
          ok=ErrMessage(G,"ObjectMap","missing brick grid.");
        tmp = PyObject_GetAttrString(Map,"lvl");
        if(tmp) {


          ObjectMapNumPyArrayToMapState(G,ms,tmp);     

          Py_DECREF(tmp);
        } else
          ok=ErrMessage(G,"ObjectMap","missing brick density.");

      }
    SceneChanged(G);
    SceneCountFrames(G);
    if(ok) {
      int a;
      for(a=0;a<3;a++) {
        ms->Min[a]=0; 
        ms->Max[a]=ms->Dim[a]-1;
      }
      
      ms->Active=true;
      ms->MapSource = cMapSourceChempyBrick;
      ObjectMapUpdateExtents(I);
      
    }
  }
  return(I);
#endif
}

/*========================================================================*/
ObjectMap *ObjectMapLoadChemPyMap(PyMOLGlobals *G,ObjectMap *I,PyObject *Map,
                                  int state,int discrete)
{
#ifdef _PYMOL_NOPY
  return NULL;
#else
  int ok=true;
  int isNew = true;
  float *cobj;
  WordType format;
  float v[3],vr[3],dens,maxd,mind;
  int a,b,c,d,e;
  ObjectMapState *ms;

  maxd = -FLT_MAX;
  mind = FLT_MAX;

  if(!I) 
       isNew=true;
  else 
       isNew=false;

  if(ok) {

       if(isNew) {
            I=(ObjectMap*)ObjectMapNew(G);
            isNew = true;
       } else {
            isNew = false;
       }
    
    if(state<0) state=I->NState;
    if(I->NState<=state) {
      VLACheck(I->State,ObjectMapState,state);
      I->NState=state+1;
    }
    ms=&I->State[state];
    ObjectMapStateInit(G,ms);

    if(!PConvAttrToStrMaxLen(Map,"format",format,sizeof(WordType)-1))
      ok=ErrMessage(G,"LoadChemPyMap","bad 'format' parameter.");
    else if(!PConvAttrToFloatArrayInPlace(Map,"cell_dim",ms->Crystal->Dim,3))
      ok=ErrMessage(G,"LoadChemPyMap","bad 'cell_dim' parameter.");
    else if(!PConvAttrToFloatArrayInPlace(Map,"cell_ang",ms->Crystal->Angle,3))
      ok=ErrMessage(G,"LoadChemPyMap","bad 'cell_ang' parameter.");
    else if(!PConvAttrToIntArrayInPlace(Map,"cell_div",ms->Div,3))
      ok=ErrMessage(G,"LoadChemPyMap","bad 'cell_div' parameter.");
    else if(!PConvAttrToIntArrayInPlace(Map,"first",ms->Min,3))
      ok=ErrMessage(G,"LoadChemPyMap","bad 'first' parameter.");
    else if(!PConvAttrToIntArrayInPlace(Map,"last",ms->Max,3))
      ok=ErrMessage(G,"LoadChemPyMap","bad 'last' parameter.");

    if(ok) {
      if (strcmp(format,"CObjectZYXfloat")==0) {
        ok = PConvAttrToPtr(Map,"c_object",(void**)&cobj);
        if(!ok)
          ErrMessage(G,"LoadChemPyMap","CObject unreadable.");        
      } else {
        ok=ErrMessage(G,"LoadChemPyMap","unsupported format.");        
      }
    }
    /* good to go */

    if(ok) {
      if (strcmp(format,"CObjectZYXfloat")==0) {

        ms->FDim[0]=ms->Max[0]-ms->Min[0]+1;
        ms->FDim[1]=ms->Max[1]-ms->Min[1]+1;
        ms->FDim[2]=ms->Max[2]-ms->Min[2]+1;
        if(Feedback(G,FB_ObjectMap,FB_Actions)) {
          printf(" LoadChemPyMap: CObjectZYXdouble %dx%dx%d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]);        
        }
        ms->FDim[3]=3;
        if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
          ok=false;
        else {
          CrystalUpdate(ms->Crystal);
          ms->Field=IsosurfFieldAlloc(G,ms->FDim);
          for(c=0;c<ms->FDim[2];c++)
            {
              v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
              for(b=0;b<ms->FDim[1];b++) {
                v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
                for(a=0;a<ms->FDim[0];a++) {
                  v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                  
                  dens = *(cobj++);

                  F3(ms->Field->data,a,b,c) = dens;
                  if(maxd<dens) maxd = dens;
                  if(mind>dens) mind = dens;
                  transform33f3f(ms->Crystal->FracToReal,v,vr);
                  for(e=0;e<3;e++) 
                    F4(ms->Field->points,a,b,c,e) = vr[e];
                }
              }
            }

          if(ok) {
            d = 0;
            for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
              {
                v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
                for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
                  v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
                  for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
                    v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                    transform33f3f(ms->Crystal->FracToReal,v,vr);
                    copy3f(vr,ms->Corner+3*d);
                    d++;
                  }
                }
              }
          }
        }
      }
    }
    
    if(ok) {
      CrystalDump(ms->Crystal);
      
      v[2]=(ms->Min[2])/((float)ms->Div[2]);
      v[1]=(ms->Min[1])/((float)ms->Div[1]);
      v[0]=(ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
      
      v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
      v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
      v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);

    }

    if(!ok) {
      ErrMessage(G,"ObjectMap","Error reading map");
    } else {
      ms->Active=true;
      ObjectMapUpdateExtents(I);
            if(Feedback(G,FB_ObjectMap,FB_Actions)) {
        printf(" ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd);
      }
    }

    if(ok) {
      SceneChanged(G);
      SceneCountFrames(G);
    }
  }
  return(I);
#endif
}