/*
// isocheck.c 
// created by F.Mizutani  960328
// modify 1996.04.23 add program code
// modify 1996.05.01 unix, dos support
*/

/* unix header file */
#include <stdio.h>
#include <string.h>

#define         DOSCOD  0x1111
#define         UNXCOD  0x2222
#define		PCODE	11002
#ifdef  DOS
#define		DCODE	DOSCOD
#else
#define		DCODE   UNXCOD
#endif
#define		OBTMAX	256
#define		PMAX	32
#define		FALSE   0
#define		TRUE    !FALSE

float    xmax, ymax, zmax, datamax;
float    xmin, ymin, zmin, datamin;

typedef struct {
	long    mn;
	float	mw;
	float	x;
	float	y;
	float	z;
} MOL;

long    idsz,   /* Mesh data size. */
        idmx,   /* Mesh size for x direction. */
        idmy,   /* Mesh size for y direction. */
        idmz;   /* Mesh size for z direction. */
float   mbgx, mbgy, mbgz;
float   mstx, msty, mstz;
float   mdvx, mdvy, mdvz;
long    dtype;
long	oct;
long	obt[OBTMAX];
MOL	*mol;
long	molct;
int	sfg;

void Swapword(s)
unsigned long *s;
{
	unsigned long   s1, s2, s3;

	s1 = (*s & 0xff000000) >> 24;
	s2 = (*s & 0x00ff0000) >>  8;
	s3 = (*s & 0x0000ff00) <<  8;
	*s &= 0x000000ff;
	*s = *s << 24;
	*s |= s1;
	*s |= s2;
	*s |= s3;
}

void InputData(fp)
FILE  *fp;
{
        long           i, j, k, n, dcode, pcode, *lpt;
        float          x, y, z, data, *dt;
	int	       st;	
        static float   smax =   (float)(1.0e30);
        static float   smin = - (float)(1.0e30);

	fread(&dcode, sizeof(long), 1, fp);
	if (dcode != DOSCOD && dcode != UNXCOD){
	    fclose(fp);
	    printf("This data is not made by isoconv.\n");
	    exit(1);
	}
	if (dcode == DCODE) sfg = FALSE;
	else                sfg = TRUE;
	fread(&pcode, sizeof(long), 1, fp);
	if (sfg) Swapword((unsigned long *) &pcode);
	if (pcode != PCODE){
	    fclose(fp);
	    printf("This data is not made by isoconv.\n");
	    exit(1);
	}
	if (dcode == DOSCOD){
	    printf("This data is DOS format.\n");
	}
	else {
	    printf("This data is not DOS format.\n");
	}
        fread(&dtype, sizeof(long), 1, fp);
	fread(&datamax, sizeof(float), 1, fp);
	fread(&datamin, sizeof(float), 1, fp);
	if (sfg) Swapword((unsigned long *) &dtype);
	if (sfg) Swapword((unsigned long *) &datamax);
	if (sfg) Swapword((unsigned long *) &datamin);
        fread(&idmx,  sizeof(long), 1, fp);
        fread(&idmy,  sizeof(long), 1, fp);
        fread(&idmz,  sizeof(long), 1, fp);
	if (sfg) Swapword((unsigned long *) &idmx);
	if (sfg) Swapword((unsigned long *) &idmy);
	if (sfg) Swapword((unsigned long *) &idmz);
        idsz = idmx * idmy * idmz;
	printf("Data        Min=%+e  Max=%+e\n", datamin, datamax);
        printf("Resolution  X=%d Y=%d Z=%d TOTAL=%d\n",idmx,idmy,idmz,idsz);
        fread(&mbgx, sizeof(float), 1, fp);
        fread(&mbgy, sizeof(float), 1, fp);
        fread(&mbgz, sizeof(float), 1, fp);
        fread(&mstx, sizeof(float), 1, fp);
        fread(&msty, sizeof(float), 1, fp);
        fread(&mstz, sizeof(float), 1, fp);
	if (sfg) Swapword((unsigned long *) &mbgx);
	if (sfg) Swapword((unsigned long *) &mbgy);
	if (sfg) Swapword((unsigned long *) &mbgz);
	if (sfg) Swapword((unsigned long *) &mstx);
	if (sfg) Swapword((unsigned long *) &msty);
	if (sfg) Swapword((unsigned long *) &mstz);
        printf("Begin AXIS  X=%+f Y=%+f Z=%+f\n",mbgx,mbgy,mbgz);
        printf("Step  AXIS  X=%+f Y=%+f Z=%+f\n",mstx,msty,mstz);

	fread(&molct, sizeof(long), 1, fp);
	if (sfg) Swapword((unsigned long *) &molct);
	printf("Moleculer count is %d\n", molct);
	mol = (MOL *)malloc(molct * sizeof(MOL));
	fread(mol, sizeof(MOL), molct, fp);
	if (sfg){
	    lpt = (unsigned long *)mol;
	    st  = sizeof(MOL) / sizeof(long);
	    for (i = 0; i < molct; i++){
		for (j = 0; j < st; j++){
		    Swapword(lpt + i * st + j);
		}
	    }
	}
	for (i = 0; i < molct; i++){
	    printf("%3d: No.%3d  W=%+f  X=%+f Y=%+f Z=%+f\n", i+1, 
		mol[i].mn, mol[i].mw, mol[i].x, mol[i].y, mol[i].z);
	}

	if (dtype == 1){
	    fread(&oct, sizeof(long), 1,   fp);
	    fread(obt,  sizeof(long), oct, fp);
	    if (sfg){
		Swapword((unsigned long *) &oct);
	        lpt = (unsigned long *)obt;
	        for (i = 0; i < oct; i++){
		    Swapword(lpt + i);
	        }
	    }
	    printf("Orbital data  COUNT=%d\n  Orbital number ", oct);
	    for (i = 0; i < oct; i++){
	        printf("No.%d ", obt[i]);
	    }
	    printf("\n");
	}
	else printf("Density data\n");

        dt = (float *)malloc(idmz * sizeof(float));

	printf("Data sample \n");
        for(i = 0; i < 1; i++){
            for (j = 0; j < 1; j++){
                if (fread(dt, sizeof(float), idmz, fp) != idmz){
	            if (sfg){
	                lpt = (unsigned long *)dt;
	                for (k = 0; k < idmz; k++){
		            Swapword(lpt + k);
	                }
	            }
                    printf("data read error ");
                    exit(1);
		}
                for (k = 0; k < 8; k++){
		    printf("  D%d=%e\n", k + 1, dt[k]);
                }
            }
	}
}

main(argc, argv)
int	argc;
char	**argv;
{
    int     db, ac;
    char    filein[1024], fileout[1024];
    FILE    *fi, *fo;
    float   dt[18];
    long    pl[PMAX];
    float   pf[PMAX];

    printf("\n");
    printf("ISOCHECK  Copyright(C) 1996  Computer Center of IMS\n\n");
    db = 0;
    if (argc != 2) {
        printf("USAGE: ISOCHECK infile\n");
        exit(1);
    }
    strcpy(filein,  argv[1]);
    if((fi = fopen(filein, "rb"))==NULL){
            perror("fopen error ");
            exit(1);
    }
    InputData(fi);
    fclose(fi);
}