/*<
Stack

 DESCRIPTION
	Stack of common-midpoint gathers

 USAGE 
	Stack norm=yes < in.H > out.H

 INPUT PARAMETERS
	norm - char*    yes  normalize stack by number of non-zero samples
	                no   no normalization
	n1	 - integer  trace length
	n2	 - integer  number of offsets
  n3   - integer  number of midpoints
  o3,d3- real     first midpoint and sampling of midpoints  

CATEGORY
Seis:Filter

COMPILE LEVEL
DISTR
>*/
/*
 KEYWORDS    stack common-midpoint data

 SEE ALSO

 WHERE
	./cube/seis/Stack.c

*/
/*
 * Keyword: stack common-midpoint data
 
  Modified:
  Bob -Nov'98  Changed stackf to C
 */

#include <sep.main>

#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int  stackf(float*,float*,float*,int*,int,int,int,int,char*,char*);
_XFUNCPROTOEND
#else
int  stackf();
#endif

char intag[]="in";
char outtag[]="out";

MAIN()
  {
	int one=1,ny,nh,nt,ntbytes,nhbytes,itemp;
	float *p,*q;
	float *hnorm;
	int *d, dlen,ierr;
	char norm[10];
	float d3,o3;
	char label3[256];

	/* obtain parameters and open files */
	if(!fetch("n1","d",&nt))
		if(!fetch("nt","d",&nt)) seperr ("Can't find nt or n1\n");

	if(!fetch("n2","d",&nh))
		if(!fetch("nh","d",&nh)) seperr ("Can't find nh or n2\n");

	if(!fetch("n3","d",&ny))
		if(!fetch("ny","d",&ny)) seperr ("Can't find ny or n3\n");

	if(!fetch("d3","f",&d3))
		if(!fetch("dy","f",&d3)) d3 = 1.0;

	if(!fetch("o3","f",&o3))
		if(!fetch("y0","f",&o3)) o3 = 0.0;

	strcpy(label3," ");			/* default for label3 */
	if(!fetch("label3","s",label3))
		fetch("ylabel","s",label3);

	strcpy(norm,"yes"); 	getch("norm","s",norm);

	/* output parameters  */
	putch("n1","d",&nt);
	putch("n2","d",&ny);   /* ny becomes n2 */
	putch("d2","f",&d3);   /* new d2 is old d3 */
	putch("o2","f",&o3);   /* new o2 is old o3 */
	putch("label2","s",label3); /* ditto */
	putch("n3","d",&one);  /* n3 becomes 1  */
	d3=1.0; o3=0.0; strcpy(label3," ");
	putch("d3","f",&d3);
	putch("o3","f",&o3);
	putch("label3","s",label3);
	putch("#norm","s",norm);
	hclose();

	/* obtain space for arrays and initialize array pointers */
	ntbytes = nt*sizeof(float);
	dlen = nt*sizeof(int);
	nhbytes = (++nh)*sizeof(float);
	hnorm = (float *) alloc( dlen + 2*ntbytes + nhbytes);
	p = hnorm + nh;
	q = p + nt;
	d = (int *) (q + nt);
	itemp = norm[0] == 'y';

        ierr=stackf(p,q,hnorm,d,nh,ny,nt,itemp,intag,outtag);
	return(0);
  }

#if NeedFunctionPrototypes
_XFUNCPROTOBEGIN
int  stackf(float *p, float *q, float *h, int *d, int nhp1, int ny, int nt, int norm, char *intag, char *outtag)
_XFUNCPROTOEND
#else
int  stackf(p,q,h,d,nhp1,ny,nt,norm,intag,outtag)
float *p, *q, *h; /* p-input, q-output, h     */
int *d, nhp1,  ny, nt, norm; /*d- , nhp1-noffsets, ny-nmidpts, nt-nt, norm-normalize*/
char *intag, *outtag; /*input tag,output tag */
#endif

{
int  ih,iy,it,ntbyte,count;

  ntbyte=nt*4;
      if(norm != 0) {
          h[0] = 1.0;
					for(ih=1; ih < nhp1; ih++) h[ih] = 1.0/ih;
					for(iy=0; iy < ny; iy++) {
						for(it=0; it < nt; it++) {
							q[it] = 0.;
							d[it] = 0;
						}
						for(ih=1; ih < nhp1; ih++){
              count=sreed(intag,p,ntbyte);
      				if(count < ntbyte) seperr("Stack: unexpected end of data\\n");
							for(it=0; it < nt; it++) {
              	q[it] = q[it]+p[it];
                if(p[it] != 0.0) d[it] = d[it]+1;
              }
           	}
						for(it=0; it< nt; it++) q[it]=q[it]*h[d[it]];
						if(ntbyte!=srite(outtag,q,ntbyte))
							seperr("trouble writing out \n");
					}
				}
				else{
					for(iy=0; iy < ny; iy++) {
          count = sreed(intag,q,ntbyte);
        	if(count < ntbyte)
        		seperr("Stack: unexpected end of data\n");
							for(ih=2; ih < nhp1; ih++){
                count = sreed(intag,p,ntbyte);
      					if(count < ntbyte) 
									seperr("Stack: unexpected end of data\n");
								for(it=0; it < nt; it++) q[it] = q[it]+p[it];
              }
              if(ntbyte!=srite(outtag,q,ntbyte))
								seperr("trouble writing out data \n");
          }
        }
		}