supickamp.c

su 的源代码库

	dfprint(dt_resamp);
	dfprint(dt);
	dfprint(tmin_in); */


	val_resamp=ealloc1float(nts);
	t_resamp=ealloc1float(nts);
   
	/* initialize val_resamp */
	for(i=0;i<nts;i++)
	   val_resamp[i]=0.0;
	
	/* initialize trace picks */
	if(! initpick(pick1) )
	   err("ERROR initializing local Pick structure");

	/* initialize global picks */
	if(! initpick(gpick) )
	   err("ERROR initializing global Pick structure");
	
	trnum=0;
	
	/* Loop over traces */
	do {
	
		register int itime;

		pick1->iabs=trnum;
		
		trnum++;
		
		/* Get value for second dimension */
		if(d2 == 0)
		{
		   /* get value of key and convert to float */
		   gethval(&tr, index, &val);
		   fval = vtof(type,val);
		} else 
		   fval=x2beg + (trnum-1)*d2;
 
		/* Get t_min in window */
		if(nx==0)
		   tmin=tmin_in;
		else

		  /* linearly interpolate between input values 
			can be taken out of loop for efficiency */
		  intlin(nx,x_above,t_above,t_above[0],t_above[nx-1],
			1,&fval,&tmin);
 
		t1=tmin_in + (nt-1)*dt;

		/* Get t_max in window */	
		if(bfilep!= NULL)
		   intlin(nx,x_below,t_below,t_below[0],t_below[nx-1],
			1,&fval,&tmax);
		else if(wl != 0)
		   tmax=tmin+wl;
		else 
		   tmax=t1;

		/* check for bounds */
		if(tmin <tmin_in)
		   tmin=tmin_in;
		if(tmax >t1)
		   tmax=t1;

		/* debug output */
		if(verbose == 1)
		fprintf(outparfp,"trace %i x=%f t_above=%f t_below=%f \n",
			trnum-1,fval,tmin,tmax);

		/* Compute output times */
		{ 
		   register float tvalue;
		   for (itime=0,tvalue=tmin; tvalue<=tmax; itime++,
			tvalue+=dt_resamp)
			t_resamp[itime] = tvalue;
		}

		/* sinc interpolate new data */
		ints8r(nt, dt, tmin_in, tr.data, 
			0.0, 0.0, itime+1,t_resamp, val_resamp); 
		
	  
		/* search for max,min,abs,energy,tmax,tmin,tabs */
		if(! tracepick(pick1,t_resamp,val_resamp,itime) )
		   err(" ERROR in tracepick ");

		/* normalize energy */
		pick1->energy *=dt_resamp;
		
		if(! globalpick(pick1,gpick) )
		   err(" ERROR in subroutine globalpick");
		
		/* debug output */
		if(verbose == 1)
		{
		fprintf(outparfp,"trace #=%i \t   x2=%e \n",
			pick1->iabs,fval);
		fprintf(outparfp,"max=%e \t tmax=%e \n",
			pick1->max,pick1->tmax);
		fprintf(outparfp,"min=%e \t tmin=%e \n",
			pick1->min,pick1->tmin);		   
		fprintf(outparfp,"abs=%e \t tabs=%e \n",
			pick1->abs,pick1->tabs);		   
		fprintf(outparfp,"energy in pick window =%e \n \n",
			pick1->energy);
		}

		/* column output */
		if(verbose == 2)
		fprintf(outparfp,
			"%i\t %f\t %f\t %f\t"
			"%i\t %e\t "
			"%e\t %e\t "
			"%e\t %e\t "
			"%e\t %e\t "
			"%e\n",
			trnum-1,fval,tmin,tmax,
			pick1->iabs,fval,
			pick1->max,pick1->tmax,
			pick1->min,pick1->tmin,
			pick1->abs,pick1->tabs,
			pick1->energy) ;

		if(! getarg(&t1,mode1,pick1,fval,trnum-1) )
		   err(" ERROR in subroutine <getarg>, argument 1");

		if(! getarg(&x1,mode2,pick1,fval,trnum-1) )
		   err(" ERROR in subroutine <getarg>, argument 2");

		if(modef==1)
		{
		   fprintf(stdout,"%e \t %e\n",t1,x1);
		} else {
		  /* Output raw BINARY pairs to stdout */  
			efwrite(&x1, FSIZE, 1, stdout);
			efwrite(&t1, FSIZE, 1, stdout);
		}
	} while (gettr(&tr));

	/* Output global information */
	fprintf(outparfp,"\nglobal max=%e \t #=%i \t tmax=%e \n",
	   gpick->max,gpick->imax,gpick->tmax);
	fprintf(outparfp,"global min=%e \t #=%i \t tmin=%e \n",
	   gpick->min,gpick->imin,gpick->tmin);
	fprintf(outparfp,"global abs=%e \t #=%i \t tabs=%e \n",
	   gpick->abs,gpick->iabs,gpick->tabs);
	fprintf(outparfp,"max energy in pick window =%e \t #=%i \n \n",
	   gpick->energy,gpick->ienergy);

	if (bfile_set) fclose(bfilep);
	if (afile_set) fclose(afilep);

	/* If output=binary, print number of pairs to outpar */
	if (modef==2)
		fprintf(outparfp, "\n npairs = %d\n", trnum);

	return(CWP_Exit());
}


int tracepick(Pick *pick1, float *t, float *val, int nt)
{
   register  int it;
   float min,max,abs,tmax,tmin,tabs,energy;

   min = FLT_MAX;
   max = -FLT_MAX;
  
   
   tmin = tmax = tabs = t[0];
   abs = 0.0;
   energy=0.0;
   
   for (it=0; it <nt;it++)
   {
	

		if (val[it] > max) { /* Find max */
				max = val[it];
				tmax = t[it];
		}

		if (val[it] < min) { /* Find min */
				min = val[it];
				tmin = t[it];
		}

		if (ABS(val[it]) > abs) {/* Find absmax */
				abs = ABS(val[it]);
				tabs = t[it];
		}

		energy += val[it]*val[it];
   
			
   
		val[it] = 0.0; /* zero data values */

   }
   
   pick1->min=min;
   pick1->max=max;
   pick1->abs=abs;
   pick1->tmin=tmin;
   pick1->tmax=tmax;
   pick1->tabs=tabs;
   pick1->energy=energy;
   
   return 1;
}


int initpick(Pick *pick1)
{
   pick1->min=FLT_MAX;
   pick1->max=FLT_MIN;
   pick1->abs=-.1;
   pick1->tmin=0.;
   pick1->tmax=0.;
   pick1->tabs=0.;
   pick1->energy=-.1;
   pick1->imax=0;
   pick1->imin=0;
   pick1->iabs=0;
   pick1->ienergy=0;  
 
   return 1;
}


int globalpick(Pick *pick1, Pick *gpick)
{
   		if (pick1->max > gpick->max) { /* Find global max */
			gpick->max = pick1->max;
			gpick->tmax = pick1->tmax;
			gpick->imax  = pick1->iabs;
		}

		if (pick1->min < gpick->min) { /* Find global min */
			gpick->min = pick1->min;
			gpick->tmin = pick1->tmin;
			gpick->imin  = pick1->iabs;
		}
		if (pick1->abs > gpick->abs) { /* Find global abs */
			gpick->abs = pick1->abs;
			gpick->tabs = pick1->tabs;
			gpick->iabs  = pick1->iabs;
		}
		if (pick1->energy > gpick->energy) { /* Find global energy */
			gpick->energy = pick1->energy;
			gpick->ienergy  = pick1->iabs;
		}

		return 1;
		
}
int get_mode(char *arg1, int *mode)
{
   if( STREQ(arg1, "max") )
   {
	*mode=1; return 1;
   }
   else if( STREQ(arg1, "tmax") )
   {
	*mode=2; return 1;
   }
   else if( STREQ(arg1, "min") )
   {
	*mode=3; return 1;
   }
   else if( STREQ(arg1, "tmin") )
   {
	*mode=4; return 1;
   }
   else if( STREQ(arg1, "abs") )
   {
	*mode=5; return 1;
   }
   else if( STREQ(arg1, "tabs") )
   {
	*mode=6; return 1;
   }
   else if( STREQ(arg1, "energy") )
   {
	*mode=7; return 1;
   }
   else if( STREQ(arg1, "x2") )
   {
	*mode=8; return 1;
   }	
   else if( STREQ(arg1, "i2") )
   {
	*mode=9; return 1;
   }
   else
	return 0 ;
}

int getarg(float *arg, int mode, Pick *pick, float x2, int trnum)
{
   switch(mode)
   {
	case 1:
		*arg = pick->max; break;
	case 2:
		*arg = pick->tmax; break;
	case 3:
		*arg = pick->min; break;
	case 4:
		*arg = pick->tmin; break;
	case 5:
		*arg = pick->abs; break;
	case 6:
		*arg = pick->tabs; break;
	case 7:
		*arg = pick->energy; break;
	case 8:
		*arg = x2; break;
	case 9:
		*arg = (float) trnum; break;
   }
   return 1;
}
int getformat(char *format, int *modef)	
{
   if (STREQ(format, "ascii") )
   {
	
	*modef=1; return 1;
   }
   else if(STREQ(format, "binary") )
   {
	
	*modef=2; return 1;
   }
   else return 0;
}