#include <sys/types.h>
#include <sys/uio.h>
#include <stdio.h>
#include <fcntl.h>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

// define SWAP_ENDIAN on Mac
// #define SWAP_ENDIAN

const int nx=256;
int np=8;
int nzp=nx/np;

float min, max;
double sum, sq;
long numz;

float v[nx][nx][nx];

typedef unsigned char BYTE;

inline void swap(BYTE *b) //swap endian little <-> big
{
  BYTE tmp;
  tmp=*b;*b=*(b+3);*(b+3)=tmp;
  tmp=*(b+1);*(b+1)=*(b+2);*(b+2)=tmp;
}

void err()
{
  fprintf(stderr,"file not found\n");
  exit(1);
}


int main(int argc, char** argv)
{
  int x,y,z,xc,yc,zc,xc1,yc1,zc1,zn,yn,res;
  float v1,v2;
  FILE *f1,*f1i;
  BYTE *b,*b0,*b1;
  char filename[]="R000498QXP000";
  int file1;
  int i,j;
  long seekp;
  int ii;
  int step;

  switch(argc)
  {
  case 1: step=4799;np=32;nzp=nx/np;break;
  case 2: sscanf(argv[1],"%d",&step);
  if(step<2||step>200000) {fprintf(stderr,"incorrect step\n");exit(1);};
  np=32;nzp=nx/np;break;
  case 3: sscanf(argv[1],"%d",&step);
  if(step<2||step>200000) {fprintf(stderr,"incorrect step\n");exit(1);};
  sscanf(argv[2],"%d",&np);
  if(np==1||np==2||np==4||np==8||np==16||np==32||np==64||np==128){;}
  else {fprintf(stderr,"incorrect np\n");exit(1);}
  nzp=nx/np; break;
  default: fprintf(stderr,"usage: test [nstep [np]]\n"); exit(1);
  }

  sprintf(filename,"R%06dQXP000",step);  

for(ii=0;ii<8;ii++)
  {
  for(i=0;i<np;i++){
   filename[12]=0x30+i%10;filename[11]=0x30+i/10;
   file1=open(filename, O_RDONLY);
   if(file1==-1) err();
   seekp=ii*(long(nx)*nx*nzp*4+8)+4; // 512x512x8 (nx.ny.nzp), 4 bytes before and after (fortran!)
   lseek(file1,seekp,SEEK_SET);
   read(file1, (BYTE*)v+i*(long(nx)*nx*nzp*4), long(nx)*nx*nzp*4);
   close(file1);
  }

#ifdef SWAP_ENDIAN
  b0=(BYTE*)v;b1=(BYTE*)v+long(nx)*nx*nx*4;
  for(b=b0;b<b1;b+=4) swap(b);
#endif  

  min=1000.0;max=-1000.0;
  sum=0.0;sq=0.0;numz=0;
  for(z=0;z<nx;z++)
  for(y=0;y<nx;y++)
  for(x=0;x<nx;x++)
  {
    v1=v[z][y][x];
    if(v1>max) max=v1;
    if(v1<min) min=v1;
    if(v1==0) numz++;
    sum+=v1;
    sq+=v1*v1;
  }
  sum/=float(nx)*nx*nx;
  sq/=float(nx)*nx*nx;
  printf("%11.5f %11.5f %11.5f %11.5f %d\n", min, max, sum, sqrt(sq-sum*sum), numz);

}
  return 0;
}