/*
// d2h.c
// created by F.Mizutani  960411
// modify 1996.05.01 unix, dos support
*/

/* unix header file */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>

#define  PCODE   11002    /* program code */
#define  DOSCOD  0x1111
#define  UNXCOD  0x2222
#ifdef  DOS
#define  DCODE   DOSCOD
#else
#define  DCODE   UNXCOD
#endif
#define  PMAX    32       /* parameters   */
#define  DIV     50       /* find neighbor polygon  < step/DIV */
#define	 OBTMAX  64
#define  HIST	 128
#define  FALSE   0
#define  TRUE    !FALSE

float    datamax;
float    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;

float   hst[HIST];
long    hct[HIST];

int     sfg;

void swapper(p, n)
char    *p;
int     n;
{
    char d[4];
    int  i, j;

    for (j = 0; j < n; j++){
        for (i = 0; i < 4; i++) d[3 - i] = p[i + 4 * j];
        for (i = 0; i < 4; i++) p[i + 4 * j] = d[i];
    }
}

/* 
 * read mesh data from input file
 * idsz :  Mesh data size
 * idmx :  Mesh size for x direction
 * idmy :  Mesh size for y direction
 * idmz :  Mesh size for z direction
 * ids  :  Mesh struct data
 */
void InputData(fp)
FILE  *fp;
{
        long           i, j, k, l, n, dcode, pcode;
        float          x, y, z, data, *dt;
	float          min, max, step;

        printf("ReadFile ... ");
	fread(&dcode, sizeof(long), 1, fp);
	if (dcode != DOSCOD && dcode != UNXCOD){
	    fclose(fp);
	    printf("This data is not made by isoconv.\n");
	    exit(1);
	}
	printf("Done\n");
	if (dcode == DCODE) sfg = FALSE;
	else                sfg = TRUE;
	fread(&pcode, sizeof(long), 1, fp);
	if (sfg) swapper((char *)&pcode, 1);
	if (pcode != PCODE){
	    fclose(fp);
	    printf("This data is not made by isoconv.\n");
	    exit(1);
	}
        fread(&dtype, sizeof(long), 1, fp);
	fread(&datamax, sizeof(float), 1, fp);
	fread(&datamin, sizeof(float), 1, fp);
	if (sfg) swapper((char *)&dtype, 1);
	if (sfg) swapper((char *)&datamax, 1);
	if (sfg) swapper((char *)&datamin, 1);
        fread(&idmx,  sizeof(long), 1, fp);
        fread(&idmy,  sizeof(long), 1, fp);
        fread(&idmz,  sizeof(long), 1, fp);
	if (sfg) swapper((char *)&idmx, 1);
	if (sfg) swapper((char *)&idmy, 1);
	if (sfg) swapper((char *)&idmz, 1);
        idsz = idmx * idmy * idmz;
	printf("Data        Max=%+e  Min=%+e\n", datamax, datamin);
        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) swapper((char *)&mbgx, 1);
	if (sfg) swapper((char *)&mbgy, 1);
	if (sfg) swapper((char *)&mbgz, 1);
	if (sfg) swapper((char *)&mstx, 1);
	if (sfg) swapper((char *)&msty, 1);
	if (sfg) swapper((char *)&mstz, 1);
        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);
        mdvx = mstx / DIV; mdvy = msty / DIV; mdvz = mstz / DIV;

	fread(&molct, sizeof(long), 1, fp);
	if (sfg) swapper((char *)&molct, 1);
	mol = (MOL *)malloc(molct * sizeof(MOL));
	fread(mol, sizeof(MOL), molct, fp);
	if (sfg) swapper((char *)mol, sizeof(MOL) / sizeof(long));
	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){
	    printf("Orbital data\n");
	    fread(&oct, sizeof(long), 1,   fp);
	    if (sfg) swapper((char *)&oct, 1);
	    fread(obt,  sizeof(long), oct, fp);
	    if (sfg) swapper((char *)obt, oct);
	}
	else printf("Density data\n");

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

	data = datamax;
	max = log10(datamax);
	min = log10(datamin);
	step = (max - min) / HIST;
	for (i = 0; i < HIST; i++){
	    hst[i] = (float)pow((double)10.0, (double)(max - i * step));
	}

	printf("X loop ");
        for (i = 0; i < idmx; i++){
            printf("[%d]", i + 1);
            for (j = 0; j < idmy; j++){
                if (fread(dt, sizeof(float), idmz, fp) != idmz){
                    printf("data read error ");
                    exit(1);
		}
	        if (sfg) swapper((char *)dt, idmz);
                for (k = 0; k < idmz; k++){
		    if (dtype == 1) data = (float)fabs((double)dt[k]);
		    else            data = dt[k];
		    for (l = 1; l < HIST; l++){
			if (data > hst[l]){ hct[l - 1]++; break; }
                    }
		    if (l == HIST) hct[HIST - 1]++;
                }
            }
        }
	printf("\n");
}



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

    db = 0;
    printf("\n");
    printf("D2H  Copyright(C) 1996  Computer Center of IMS\n\n");
    if (argc != 2 && argc != 3 && argc != 4) {
        printf("USAGE: D2H infile\n");
        printf("       output is change extention to .his\n");
        exit(1);
    }
    strcpy(filein,  argv[1]);
    strcpy(fileout, filein);
    for (i = strlen(fileout) - 1; i >= 0; i--){
	if (fileout[i] == '.') break;
    }
    if (i > 0){
	fileout[i] = 0x00;
	strcat(fileout, ".his");
    }
    ac = 3;
    while (1){
        if (ac >= argc) break;
        if (argv[ac][0] == '-' || argv[ac][0] == '/'){
            if (toupper(argv[ac][1]) == 'D') db    = 1;
	    if (toupper(argv[ac][1]) == 'O') db    = -1;
        }
        ac++;
    }
    printf("argc= %d infile = %s outfile = %s\n",
            argc, filein, fileout);

    if((fi = fopen(filein, "rb"))==NULL){
            perror("fopen error ");
            exit(1);
    }
    InputData(fi);
    fclose(fi);

    if ((fo = fopen(fileout, "wt")) != NULL){
	for (i = 0; i < HIST; i++) fprintf(fo, "%e %ld\n", hst[i], hct[i]);
        fclose(fo);
    }
    printf("Terminate\n");
}