PhiMap file (.phi) files

  Spock's PhiMap format is based on that of Grasp and DelPhi for maximum compatibility. However, those programs are limited to a 65 cubed grid size, while spock supports grid sizes ranging from 3 cubed to 257 cubed, which requires some trickery in the file format. The file is a fortran unformatted data file, which may be read by the following routine taken from the spock source (a similar routine may be used to create PhiMap files that spock can read):

      subroutine readphimap(phi,N,mid,scale,fname)
c/*-
c     read a phimap from the specified file.
c     return 0 in scale if error
c-*/
      real*4 phi(N,N,N)
      real*4 mid(3)
      real*4 scale
      character*(*) fname

      character*20 uplbl
      character*10 nxtlbl
      character*60 toplbl
      character*16 botlbl

      open(3,file=fname,form='unformatted',err=999)
      read(3,err=999) uplbl
      read(3,err=999) nxtlbl,toplbl

      read(3,err=999) phi

      read(3,err=999) botlbl
      read(3,err=999) scale,mid
      close(3)
      return
 999  call perror("spock")
      write(*,*) 'Error reading phimap'
      scale=0
      close(3)
      return
      end

The trickery is in the PHI variable's dimensions, here specified as N. In DelPhi/Grasp, N is always 65, while in spock N may vary. When writing a PhiMap file N is known in advance, so there is no problem. However, when reading a PhiMap file, there's no way to specify the size of the map and still maintain compatibility. However, the size of the file itself can serve as an indicator of how much data is in it, as the only varying parameter is the grid size. You'll have to use dynamic memory allocation (or reserve a chunk of memory large enough to hold the largest possible grid size -- 257 cubed). Here's a C function to determine the size of the grid in a PhiMap file based on its file size:

static int GetPhiMapDimension(char *fname)
{
  int cubes[8]={27,125,729,4913,35937,274625,2146689,16974593};
  int dims[8]={3,5,9,17,33,65,129,257};
  int dim,i,fsize;
  int extra=106 /* characters */ + 4*sizeof(float) /* floats */;
  FILE *fil;

  /* get dimension of file */
  if (!(fil=fopen(fname,"r"))) {
    printf("Couldn't get file %s\n",fname);
    return 0;
  }
  fseek(fil,0,SEEK_END);
  fsize=ftell(fil)-extra;
  fclose(fil);

  dim=0;
  for(i=0;i<7;i++) {
    if(fsize>=sizeof(float)*cubes[i] && fsize<sizeof(float)*cubes[i+1]) {
      dim=dims[i];
      break;
    }
  }

  if (dim==0)
    printf("Couldn't determine dimension of phimap in file\n");
  else 
    printf("This is a %d cubed PhiMap\n",dim);

  return dim;
}

So, if you want to read in a phimap of unknown dimension, call GetPhiMapDimension() which will return the dimension of the map. Then allocate enough memory to hold the map and call readphimap_().

Note: I use fortran as little as possible in spock, but it's necessary for the unformatted I/O. That's why GetPhiMapDimension is in C and readphimap is in fortran. I won't be offended if you write a fortran version of get GetPhiMapDimension, but it's harder to do dynamic memory allocation portably in fortran. You're welcome to use these routines in your own programs. If you write a program to store some other interesting property in a phimap file, I'd like to know about it.