grids.c~

seismic software,very useful

#include "gridhd.h"
#include "ghdr.h"
#include "grid.h"
#include "par.h"

int fgetghdr(FILE  *gfp, ghed *gh) {
/*  read in a grid header from a file pointer
  input:
	gfp 	--	file pointer of disk grid file
  output:
	gh	--	grid header 
  return:
	0	--	grid header found
	1	--	not a standard grid file (no header found)
  author:
	Zhiming Li, 	8/17/92		      

  notes:
	the grid file consists of n1*n2*n3*n4*n5 grid data values followed 
        by NGHDRS (25) floating-point header numbers 
*/
	
	int ierr;
	float tmp; 
	int nread, idtype;
	long long int n1, n2, n3, n4, n5;
	float scale;

	filetype ftype;

	long long isize, jsize, iseek;
	

	ierr = 0;
	/* check to see if gfp is pipe */
	/*
	ftype = filestat(fileno(gfp));
	*/
	
	/*
	if(ftype == DISK) {
	*/
		/* retain the pointer position */
		file2g(gfp);
		iseek = ftell2g(gfp);
		/* check grid file size */ 
		fseek2g(gfp,0L,2);
		jsize = ftell2g(gfp);

		if(jsize <  GHDRBYTES ) {
			ierr = 1;
			warn(" warning: grid file size less than %d bytes",
				GHDRBYTES);
			return ierr;
		}

		fseek2g(gfp,-GHDRBYTES,2);

	/*
	}
	*/

	/* read in header */
	nread = fread( (char *)gh,1,GHDRBYTES,gfp);
	if(nread!=GHDRBYTES) {
		warn(" warning: grid header read error \n");
		ierr = 1;
		return ierr; 
	} 
	/* check parameters */
	scale = gh->scale;
	if(scale==0.) {
		warn(" warning: gh.scale = 0.; 1 used \n");
		scale = 1.;
	}
	tmp  = gh->n1 / scale + 0.5;
	n1 = (int) tmp;
	if(n1==0) n1 = 1; 
	tmp  = gh->n2 / scale + 0.5;
	n2 = (int) tmp;
	if(n2==0) n2 = 1;
	tmp  = gh->n3 / scale + 0.5;
	n3 = (int) tmp;
	if(n3==0) n3 = 1;
	tmp  = gh->n4 / scale + 0.5;
	n4 = (int) tmp;
	if(n4==0) n4=1;
	tmp  = gh->n5 / scale + 0.5;
	n5 = (int) tmp;
	if(n5==0) n5 = 1;

	/* get grid data type */
	tmp = gh->dtype / scale + 0.5;
	idtype = (int) tmp;
	if(idtype!=1 && idtype!=2 && idtype!=4 && idtype!=8) {
		warn(" warning: invalid dtype=%d; dtype=4 used \n",idtype);
		idtype = 4;
	} 

/*
	if(ftype==DISK) {
*/
		isize = n1*n2*n3;
		isize = isize*n4*n5*idtype + GHDRBYTES;
		if(jsize != isize ) {
			ierr = 1;
			warn(
	" warning: grid file size: %g bytes not consistant with grid header\n",
				(float)jsize);
			return ierr;
		}
		/* return to the original pointer position */
		fseek2g(gfp,iseek,0);
/*
	}
*/


	return ierr;
}
	  
	
int fputghdr(FILE  *gfp, ghed *gh) {
/*  write a grid header to a grid file
  input:
        gfp     --      file pointer of disk grid file
        gh      --      grid header
  return:
        0       --      grid header written/updated
        1       --      grid size error
  author:
        Zhiming Li,     8/17/92               

  notes:
        the grid file consists of n1*n2*n3*n4*n5 grid data values followed
        by NGHDRS (25) floating-point header numbers
*/
	int ierr;
        float tmp;
        int nwrite, idtype;
        int n1, n2, n3, n4, n5;
        float scale;

	filetype ftype;
	long long isize, jsize;

	ierr = 0;

	/* check if gfp is pipe */
	/*
	ftype = filestat(fileno(gfp));
	*/


/*
	if(ftype==DISK) {
*/
		/* flush gfp buffer, in case call program forgets to do so */
			file2g(gfp);
        	fflush(gfp);
        	/* check grid file size */
        	fseek2g(gfp,0L,2);
        	jsize = ftell2g(gfp);
			/*
        	if(jsize <  GHDRBYTES ) {
                	ierr = 1;
                	warn(" warning: grid file size less than %d bytes",
				GHDRBYTES);
                	return ierr;
        	}
			*/
/*
	}
*/

	/* check parameters */
	/* make sure n1,n2,n3,n4,n5 are not zero */
        scale = gh->scale;
        if(scale==0.) {
                warn(" warning: gh.scale = 0.; 1 used \n");
                scale = 1.;
        }
        tmp  = gh->n1 / scale + 0.5;
        n1 = (int) tmp;
        if(n1==0) n1 = 1;
        tmp  = gh->n2 / scale + 0.5;
        n2 = (int) tmp;
        if(n2==0) n2 = 1;
        tmp  = gh->n3 / scale + 0.5;
        n3 = (int) tmp;
        if(n3==0) n3 = 1;
        tmp  = gh->n4 / scale + 0.5;
        n4 = (int) tmp;
        if(n4==0) n4=1;
        tmp  = gh->n5 / scale + 0.5;
        n5 = (int) tmp;
        if(n5==0) n5 = 1;

   	/* get grid data type */
        tmp = gh->dtype / scale + 0.5;
        idtype = (int) tmp;
	if(idtype!=1 && idtype!=2 && idtype!=4 && idtype!=8) {
		warn(" warning: invalid dtype=%d; dtype=4 used \n",idtype);
		idtype = 4;
	} 


	gh->scale = scale;
	gh->dtype = (float)idtype * scale;
	gh->n1 = (float)n1 * scale;
	gh->n2 = (float)n2 * scale;
	gh->n3 = (float)n3 * scale;
	gh->n4 = (float)n4 * scale;
	gh->n5 = (float)n5 * scale;

	/* position gfp for grid header */
	if(ftype==DISK) {
		isize = n1*n2*n3;
		isize = isize*n4*n5*idtype + GHDRBYTES;
        	if(jsize==isize || jsize==isize-GHDRBYTES) {
			fseek2g(gfp,isize-GHDRBYTES,0);
		} else {
                	ierr = 1;
        		warn(
	" warning: grid file size: %g bytes not consistant with grid header\n",
                	(float)jsize);
			/* position to the end of the file */
			fseek2g(gfp,0L,2);
        	}
	}
	nwrite = efwrite((char*)gh,1,GHDRBYTES,gfp); 
	if(nwrite!=GHDRBYTES) {
		ierr=1;
		warn(" warning: grid header write error \n");
	}

	return ierr;
}

void fromghdr(ghed *gh, float *scale, int *dtype, int *n1, int *n2, int *n3,
	int *n4, int *n5, float *d1, float *d2, float *d3, float *d4,
	float *d5, float *o1, float *o2, float *o3, float *o4, float *o5,
	float *dcdp2, float *dline3, float *ocdp2, float *oline3,
	float *gmin, float *gmax, int *orient, int *gtype) {
/* from grid header gh, find the values of the header parameters */ 
/* return values with scale factor removed */ 

	float tmp, sc;
	
	sc = gh->scale;
	if(sc==0.) err(" grid file header value: gh.scale = 0. ");
	*scale = sc;
	tmp = gh->dtype / sc + 0.5;
	*dtype = (int) tmp;
	tmp = gh->n1 / sc + 0.5;
	*n1 = (int) tmp;
	tmp = gh->n2 / sc + 0.5;
	*n2 = (int) tmp;
	tmp = gh->n3 / sc + 0.5;
	*n3 = (int) tmp;
	tmp = gh->n4 / sc + 0.5;
	*n4 = (int) tmp;
	tmp = gh->n5 / sc + 0.5;
	*n5 = (int) tmp;
	*d1 = gh->d1 / sc; 	
	*d2 = gh->d2 / sc; 	
	*d3 = gh->d3 / sc; 	
	*d4 = gh->d4 / sc; 	
	*d5 = gh->d5 / sc; 	
	*o1 = gh->o1 / sc; 	
	*o2 = gh->o2 / sc; 	
	*o3 = gh->o3 / sc; 	
	*o4 = gh->o4 / sc; 	
	*o5 = gh->o5 / sc; 	
	*dcdp2 = gh->dcdp2 / sc; 	
	*dline3 = gh->dline3 / sc; 	
	*ocdp2 = gh->ocdp2 / sc; 	
	*oline3 = gh->oline3 / sc; 	
	*gmin = gh->gmin / sc; 	
	*gmax = gh->gmax / sc; 	
	tmp = gh->orient / sc + 0.5; 	
	*orient = tmp; 
	tmp = gh->gtype / sc + 0.5; 	
	*gtype = tmp; 
}
		
void toghdr(ghed *gh, float *scale, int *dtype, int *n1, int *n2, int *n3,
	int *n4, int *n5, float *d1, float *d2, float *d3, float *d4,
	float *d5, float *o1, float *o2, float *o3, float *o4, float *o5,
	float *dcdp2, float *dline3, float *ocdp2, float *oline3,
	float *gmin, float *gmax, int *orient, int *gtype) {
/* update the grid header gh, using the values of the header parameters */ 
/* header values will be updated with scale factor applied */ 

	float sc;
	
	sc = *scale;
	if(sc==0.) err(" input scale factor: scale = 0. ");
	gh->scale = sc;
	gh->dtype = *dtype * sc;
	gh->n1 = *n1 * sc;
	gh->n2 = *n2 * sc;
	gh->n3 = *n3 * sc;
	gh->n4 = *n4 * sc;
	gh->n5 = *n5 * sc;
	gh->d1 = *d1 * sc;
	gh->d2 = *d2 * sc;
	gh->d3 = *d3 * sc;
	gh->d4 = *d4 * sc;
	gh->d5 = *d5 * sc;
	gh->o1 = *o1 * sc;
	gh->o2 = *o2 * sc;
	gh->o3 = *o3 * sc;
	gh->o4 = *o4 * sc;
	gh->o5 = *o5 * sc;
	gh->dcdp2 = *dcdp2 * sc;
	gh->dline3 = *dline3 * sc;
	gh->ocdp2 = *ocdp2 * sc;
	gh->oline3 = *oline3 * sc;
	gh->gmin = *gmin * sc;
	gh->gmax = *gmax * sc;
	gh->orient = *orient * sc;
	gh->gtype = *gtype * sc;
}

/* get header value from key */
void getgval(ghed *gh, char *key, float *val) { 

	int index = getgindex(key);
	char *tp = (char*) gh;

	float scale = *( (float*) (tp + ghdr[0].offs) );
	*val = *( (float*) (tp + ghdr[index].offs) );
	*val = *val / scale;
}

/* put header value from key */
void putgval(ghed *gh, char *key, float val) { 

	int index = getgindex(key);
	char *tp = (char*) gh;
	float scale = *( (float*) (tp + ghdr[0].offs) );
	if ( -1 == (index))
                err("%s: key word not in ghdr.h: '%s'", __FILE__, key);


	if(scale==0.) scale = 1.;
	*( (float*) (tp + ghdr[index].offs) ) = val * scale;
}

/* get grid header index by key */

int getgindex(char *key)        /* get index for this key */
{
        register int i;

        for (i = 0; i < NGHDRS; i++)
                if (STREQ(ghdr[i].key, key))
                        return i;       /* key found */

        /* not found */
        return -1;
}

/* get key from grid header index  */
char *getgkey(int index)
{
        return (index < NGHDRS && index >= 0) ? ghdr[index].key : NULL;
}

/* get grid header value type from key */
char *ghdtype(char *key)
{
        int index = getgindex(key);

        if (-1 == (index))
                err("%s: key word not in ghdr.h: '%s'", __FILE__, key);

        return ghdr[index].type;
}


/* find min and max absolute values of a short int array */
void sminmax(short *sa, int n, float *gmin, float *gmax) {
	int i;
	float norm;
	norm = (float) sa[0];
	*gmin = *gmax = norm;
	for(i=0;i<n;i++) {
		norm  = (float) sa[i];
		if(*gmin>norm) *gmin  = norm;
		if(*gmax<norm) *gmax  = norm;
	}
}

/* find min and max values of a float array */
void fminmax(float *fa, int n, float *gmin, float *gmax) {
	int i;
	*gmin = *gmax = fa[0];
	for(i=0;i<n;i++) {
		if(*gmin>fa[i]) *gmin  = fa[i];
		if(*gmax<fa[i]) *gmax  = fa[i];
	}
}

/* find min and max absolute values of a complex array */
void cminmax(complex *ca, int n, float *gmin, float *gmax) {
	int i;
	float norm;
	*gmin = *gmax = sqrt(ca[0].r*ca[0].r+ca[0].i*ca[0].i);
	for(i=0;i<n;i++) {
		norm  = sqrt(ca[i].r*ca[i].r+ca[i].i*ca[i].i);
		if(*gmin>norm) *gmin  = norm;
		if(*gmax<norm) *gmax  = norm;
	}
}