#ifdef _CVI_
# include <ansi_c.h>
# else _CVI_
# include <stdlib.h>
# include <stdio.h>
# include <string.h>
#endif _CVI_
#include "H3D.h"
#define H3DMAX 500
H3D *h3[H3DMAX];
//int Printf(char *format, ...);
int H3DClear(int h3d) {
if (!h3[h3d]) return -1;
memset(h3
[h3d
]->data
, 0,h3
[h3d
]->size
); h3[h3d]->min=0;
h3[h3d]->max=0;
h3[h3d]->nentries=0;
return 0;
}
int H3DPrint(int h3d) {
if (!h3[h3d]) return -1;
//Printf("PrintH3D nx=%d minx=%f stepx=%f ny=%d miny=%f stepy=%f size=%d\n", h3[h3d]->nx, h3[h3d]->minx, h3[h3d]->stepx, h3[h3d]->ny, h3[h3d]->miny, h3[h3d]->stepy, h3[h3d]->size ) ;
return 0;
}
int H3DExist(int h) {
if (h3[h]) return 1;
else return 0;
}
int H3DCalculateBin(int h, int axis, double value) {
int nx=1,xmin=0,dx=0;
int bin;
switch (axis) {
case 0:
nx = H3DGetNbinsX(h);
xmin= H3DGetMinX(h);
dx = H3DGetStepX(h);
break;
case 1:
nx = H3DGetNbinsY(h);
xmin= H3DGetMinY(h);
dx = H3DGetStepY(h);
break;
case 2:
nx = H3DGetNbinsZ(h);
xmin= H3DGetMinZ(h);
dx = H3DGetStepZ(h);
break;
}
if (dx==0) return -1;
if (value<xmin) return -1;
bin = (int)((value-xmin)/dx);
if (bin>=nx) return -1;
else return bin;
}
int H3DFill(int h3d,double x, double y, double z, double val) {
int idx, ix,iy, iz;
if (!h3[h3d]) return -1;
ix = H3DCalculateBin(h3d,0,x);
if (ix<0) return ix;
iy = H3DCalculateBin(h3d,1,y);
if (iy<0) return iy;
iz = H3DCalculateBin(h3d,2,z);
if (iz<0) return iz;
H3DFillBin(h3d, ix, iy, iz, val);
return 0;
}
int H3DGetBin(int h3d,int x, int y, int z) {
return x+h3[h3d]->nx *(y+z*h3[h3d]->ny);
}
int H3DFillBin(int h3d,int x, int y, int z, double val) {
int idx;
if (!h3[h3d]) {
//Printf("FillH3D error h3d is not initialized\n");
return -1;
}
idx = H3DGetBin(h3d,x,y,z);
h3[h3d]->data[idx]+=val;
//Printf("%d %d data %f %f\n",x,y,val, h3[h3d]->data[idx]);
if (h3[h3d]->data[idx]>h3[h3d]->max) h3[h3d]->max= h3[h3d]->data[idx];
if (h3[h3d]->data[idx]<h3[h3d]->min) h3[h3d]->min= h3[h3d]->data[idx];
h3[h3d]->nentries++;
return 0;
}
double H3DGetBinContent(int h3d,int atx, int aty, int atz) {
int idx;
if (!h3[h3d]) return 0;
if (h3[h3d]->nx <= atx) return 0;
if (h3[h3d]->ny <= aty) return 0;
if (h3[h3d]->nz <= atz) return 0;
idx = H3DGetBin(h3d,atx,aty,atz);
if (idx*sizeof(double) < h3[h3d]->size ) return h3[h3d]->data[idx];
return 0;
}
int H3DInit(int h3d,char *name, char *title,int nx, double minx, double stepx, int ny, double miny, double stepy, int nz, double minz, double stepz) {
if (h3[h3d]) {
h3[h3d] = NULL;
}
//printf("InitH3D hID=%d\n",h3d);
h3
[h3d
] = (H3D
*) malloc(sizeof(H3D
));//h2 =h3d;
H3DSetTitle(h3d,title);
H3DSetName(h3d,name);
h3[h3d]->id=H3D_ID;
h3[h3d]->nx = nx;
h3[h3d]->ny = ny;
h3[h3d]->nz = nz;
h3[h3d]->minx = minx;
h3[h3d]->miny = miny;
h3[h3d]->minz = minz;
h3[h3d]->stepx = stepx;
h3[h3d]->stepy = stepy;
h3[h3d]->stepz = stepz;
h3[h3d]->size = h3[h3d]->nx*h3[h3d]->ny*h3[h3d]->nz*sizeof(double);
h3
[h3d
]->data
= (double *) malloc(h3
[h3d
]->size
); h3[h3d]->len=sizeof(H3D)-sizeof(double *)+h3[h3d]->size;
H3DClear(h3d);
H3DPrint(h3d);
//Printf("InitH3D 0x%x\n", h3d );
return 0;
}
double H3DGetXBinCenter(int h3d,int xbin) {
return h3[h3d]->minx+xbin*h3[h3d]->stepx;
}
double H3DGetYBinCenter(int h3d,int ybin) {
return h3[h3d]->miny+ybin*h3[h3d]->stepy;
}
double H3DGetZBinCenter(int h3d,int zbin) {
return h3[h3d]->minz+zbin*h3[h3d]->stepz;
}
int H3DWrite2File(int h3d,FILE *fp) {
if (!fp) return -1;
//printf("H3D sizeof(H3D)=%lu len-datasize=%d len=%lu datasize=%d\t",sizeof(H3D)-sizeof(double *),h3[h3d]->len-h3[h3d]->size,h3[h3d]->len,h3[h3d]->size);
//printf("H3D sz=%d %d\n",sizeof(double),sizeof(double *));
fwrite (h3
[h3d
], 1, sizeof(H3D
)-sizeof(double *), fp
); fwrite (h3
[h3d
]->data
, 1, h3
[h3d
]->size
, fp
); return 0;
}
int H3DWrite(int h3d,const char *fname,const char *opt) {
FILE
*fp
=fopen(fname
,opt
); H3DWrite2File(h3d,fp);
return 0;
}
int H3DSetTitle(int h3d,char *title) {
sprintf(h3
[h3d
]->title
,"%s",title
); return 0;
}
int H3DSetTitleX(int h3d,char *title) {
sprintf(h3
[h3d
]->titlex
,"%s",title
); return 0;
}
int H3DSetTitleY(int h3d,char *title) {
sprintf(h3
[h3d
]->titley
,"%s",title
); return 0;
}
int H3DSetTitleZ(int h3d,char *title) {
sprintf(h3
[h3d
]->titlez
,"%s",title
); return 0;
}
int H3DSetName(int h3d,char *title) {
return 0;
}
int H3DGetNbinsY(int h) {
if (h3[h]) return h3[h]->ny;
else return 0;
}
int H3DGetNbinsX(int h) {
if (h3[h]) return h3[h]->nx;
else return 0;
}
int H3DGetNbinsZ(int h) {
if (h3[h]) return h3[h]->nz;
else return 0;
}
double H3DGetMinX(int h) {
if (h3[h]) return h3[h]->minx;
else return 0;
}
double H3DGetMinY(int h) {
if (h3[h]) return h3[h]->miny;
else return 0;
}
double H3DGetMinZ(int h) {
if (h3[h]) return h3[h]->minz;
else return 0;
}
double H3DGetStepX(int h) {
if (h3[h]) return h3[h]->stepx;
else return 0;
}
double H3DGetStepY(int h) {
if (h3[h]) return h3[h]->stepy;
else return 0;
}
double H3DGetStepZ(int h) {
if (h3[h]) return h3[h]->stepz;
else return 0;
}
double H3DGetMin(int h) {
if (h3[h]) return h3[h]->min;
else return 0;
}
double H3DGetMax(int h) {
if (h3[h]) return h3[h]->max;
else return 0;
}
double * H3DGetSliceXYData(int h, int slice) {
if (h3[h]) return (h3[h]->data+slice*h3[h]->nx*h3[h]->ny);
else return NULL;
}
double * H3DGetData(int h) {
if (h3[h]) return h3[h]->data;
else return NULL;
}