/*  from Kevin Cowtan   */
/*  Taken from xtalview */
#include <stdio.h>
#include <string.h>
#include <math.h>

#define ROUND(a) ((a)>0 ? (int)(a+0.5) : (int)(a-0.5))

typedef struct {
        int nx,ny,nz;
        float a,b,c,alpha,beta,gamma;
        double scale; /* scale on density */
        int nsym;
        } MAPHEADER;





void main(argc, argv)
     int argc;
     char **argv;
{
  MAPHEADER mh;
  FILE *file_ip,*file_op;
  float *map,f;
  int x,y,z;

  mh.scale=100.0;
  mh.nsym=1;
  printf("enter grid\n");
  scanf("%i %i %i",&(mh.nx),&(mh.ny),&(mh.nz));
  printf("enter cell\n");
  scanf("%f %f %f %f %f %f",
	&(mh.a),&(mh.b),&(mh.c),&(mh.alpha),&(mh.beta),&(mh.gamma));

  printf("%i %i %i\n",mh.nx,mh.ny,mh.nz);
  printf("%f %f %f\n",mh.a,mh.b,mh.c);
  printf("%f %f %f\n",mh.alpha,mh.beta,mh.gamma);

  map=malloc(mh.nx*mh.ny*mh.nz*sizeof(float));

  file_ip=fopen(argv[1],"r");
  for (x=0;x<mh.nx;x++)
    for (y=0;y<mh.ny;y++)
      for (z=0;z<mh.nz;z++) {
	fscanf(file_ip,"%*f %*f %*f %f\n",&f);
	map[x+mh.nx*(y+mh.ny*z)]=100.0*f;
      }
  fclose(file_ip);

  file_op=fopen(argv[2],"w");
  fswrite(map,file_op,&mh,0,0);
  fclose(file_op);
}


int 
fswrite(rho,fp,mh,norn,ncent)
float * rho; /* rho[nx*ny*nz] is sorted x fast, y medium, z slow in C rho[nz][ny][nx] */
FILE *fp; /* the file to output to must be opened in calling routine */
MAPHEADER *mh; /* Kevin */
int norn; /* 0 = y planes; 1 = x planes; 2 = zplanes */
int ncent; /* 0 for acentric; 1 for centric (i.e. Patterson) */
{
	int irow,icol,iplane,ii,nrh,ix,iy,iz,i;
	int j;
	int nrow, ncol, nplane, n;
	int nsym, index;
	int header, trailer;
	int sint=sizeof(int), irho;
	int buf[256];
	int nx,ny,nz;
	float para[6];
	int mnx[6];
	float scale;
	float rspmx = 0.0;
	int nbit = 0;
	int latt = 1; /* fsfour code for p lattice */
	int npic = 10;
	int noset =0;
	float zero = 0.0;
	char itle[80];

	for(i=0;i<80;i++)itle[i]=' ';
	strcpy(itle,"Fsfour format map from map2fs");
	scale = mh->scale;
	nsym = mh->nsym;
	nx = mh->nx;
	ny = mh->ny;
	nz = mh->nz;
	para[0]=mh->a;
	para[1]=mh->b;
	para[2]=mh->c;
	para[3]=mh->alpha;
	para[4]=mh->beta;
	para[5]=mh->gamma;
	for(i=0;i<6;i++)mnx[i]=0;
	if( norn < 0 || norn > 2){
		norn = 0;
	}
	if(ncent <0 || ncent > 1)ncent = 0;
	if ( norn == 0 ){
		nrow = nx;
		ncol = nz;
		nplane = ny;
	}
	if ( norn == 1 ){
		nrow = ny;
		ncol = nz;
		nplane = nx;
	}
	if ( norn == 2 ){
		nrow = nx;
		ncol = ny;
		nplane = nz;
	}
	header = (20 + 1 + 6 + 4*1)*4;
	fwrite(&header,sint,1,fp);
	fwrite(itle, 1, 80,fp);
	fwrite(&noset, sint,1,fp);
	fwrite(para,sint, 6, fp);
	fwrite(&nsym,sint,1,fp);
	fwrite(&ncent,sint,1,fp);
	fwrite(&latt,sint,1,fp);
	fwrite(&npic,sint,1,fp);
	fwrite(&header,sint,1,fp);
	
	/* there are nsym records with 9 numbers next */
	header = 9 * 4;
	for( i = 0; i<nsym; i++){
		fwrite(&header, sint, 1,fp);
		for(j=0;j<9;j++)fwrite(&zero,sint,1,fp);
		fwrite(&header, sint, 1,fp);
	}
	header = 4 * ( 7 + 6 );
	fwrite(&header,sint,1,fp);
	fwrite(&nrow, sizeof(int), 1, fp);
	fwrite(&ncol, sizeof(int), 1, fp);
	fwrite(&nplane, sizeof(int), 1, fp);
	fwrite(&scale, sizeof(int), 1, fp); 
	fwrite(&norn, sizeof(int), 1, fp);
	fwrite(mnx, sint, 6, fp);
	fwrite(&rspmx,sint,1,fp);
	fwrite(&nbit,sint,1,fp);
	fwrite(&header,sint,1,fp);

	ii = 0;
	header = 4 * nrow;
	for(iplane=0;iplane<nplane;iplane++){
		for(icol=0;icol<ncol;icol++){
			fwrite(&header, sint, 1, fp); 
			for(irow=0;irow<nrow;irow++){
				ii++;
				if(norn==0){
					ix= irow;
					iy= iplane;
					iz= icol;}
				if(norn==1){
					ix= iplane;
					iy= irow;
					iz= icol;}
				if(norn==2){
					ix= irow;
					iy= icol;
					iz= iplane;}
				index = nx*(ny*iz+iy) + ix;
				irho = ROUND(rho[index]);
				fwrite(&irho,sint,1,fp);
			}
			fwrite(&header, sint, 1, fp); 
		}
	}
	return(ii);
}