/* JRG, the Resource Geology Seismic Processing System for Java,
   and Viewmat for Java are Copyright 1998-2000 by John N. Louie

The software and methods here are the subject of academic research,
not commercial products. I would like to know what use you make of
my methods, and have your feedback on their success or failure.
Also, by letting me know who you are, I can inform you if bugs or
errors are discovered.
 
Please send me an email message with: your name and email address;
whether you are a student or faculty member, consultant or employee;
the name of your university or company, and department; and a sentence
or two describing what use you intend to make of these methods.
Send your message to louie@seismo.unr.edu */


import java.lang.Math;
import java.util.*;
import java.io.*;

public class CMPStack
	{
	static float dt;
	boolean geomexists;
	RgHead headers[], midpoints[];
	static float off0, doff, s0, ds;
	static int ntau;
	static float tau0;
	static int VP0;
	static float VP0projdist, mangle, m0;
	static int nm;
	static float dm, arc, v0;
	static int nv;
	static float dv;
	FltVol fvol;
	String picktext;
	VInfo velinfo[];
	boolean dixvel;
	static String pickformat = "Use the format\n"
		+ "below, which is what would result from saving\n"
		+ "picks made from a cvstack display.\n"
		+ "The amplitude and indices columns must exist,\n"
		+ "although their values will be ignored.\n\n"
		+ "Amp\titau\tTau(s)\tim\tMdist(m)\tiv\tVnmo(m/s)\n"
		+ "0\t0\t0.0\t\t0\t0.0\t\t0\t1500.0\n\n"
		+ "This example will apply a constant NMO velocity\n"
		+ "of 1500 m/s.";
	
	/* Constructor - sets class instance data as default parameters for stack */
	CMPStack()
		{
		this.setNull();
		}
	void setNull()
		{
		dt = 1F;
		geomexists = false;
		headers = null;
		midpoints = null;
		off0 = 0F;
		doff= 1F;
		s0 = 0F;
		ds = 1F;
		ntau = 0;
		tau0 = 0F;
		VP0 = 0;
		VP0projdist = 0F;
		mangle = 90F;
		m0 = 0F;
		nm = 0;
		dm = 1F;
		arc = 5F*dm;
		v0 = 1500F;
		nv = 1;
		dv = 200F;
		fvol = null;
		picktext = pickformat;
		velinfo = null;
		dixvel = false;
		}
	/* Construct CMPStack parameters appropriate for a FltVol of records */
	CMPStack(FltVol fvol)
		{
		this.fvol = fvol;
		if (fvol == null)
			this.setNull();
		else
			{
			if (fvol.getDUnitofDim(0) != 1F)
				dt = fvol.getDUnitofDim(0);
			headers = fvol.headers;
			if (headers == null)
				geomexists = false;
			else
				geomexists = true;
			if (geomexists)
				{
				off0 = 0F;
				doff= 1F;
				s0 = 0F;
				ds = 1F;
				midpoints = RgHead.getMidpointArray(headers);
				}
			else
				{
				if (fvol.getUnit0ofDim(1) != 0F)
					off0 = fvol.getUnit0ofDim(1);
				if (fvol.getDUnitofDim(1) != 1F)
					doff= fvol.getDUnitofDim(1);
				if (fvol.getUnit0ofDim(2) != 0F)
					s0 = fvol.getUnit0ofDim(2);
				if (fvol.getDUnitofDim(2) != 1F)
					ds = fvol.getDUnitofDim(2);
				midpoints = null;
				}
			if (ntau > fvol.getNMembofDim(0) || ntau < fvol.getNMembofDim(0)/20)
				ntau = fvol.getNMembofDim(0);
			if (fvol.getUnit0ofDim(0) != 0 || tau0 >= fvol.getNMembofDim(0)*dt)
				tau0 = fvol.getUnit0ofDim(0);
			/*m0 = 0F;*/
			if (fvol.getDUnitofDim(1) != 1F && dm < 1.1F*fvol.getDUnitofDim(1))
				dm = fvol.getDUnitofDim(1);
			if (dm < 0F)
				dm = -dm;
			if (geomexists)
				{
				RgHead projection = new RgHead();
				/* Projection info is packed into the RgHead object */
				projection.getProjection(midpoints);
				VP0 = projection.svp;
				mangle = (float)(projection.gy);
				nm = (int)(projection.off/dm);
				VP0projdist = (float)(projection.gz);
				if (dm > projection.gz)
					arc = 2F*dm;
				else
					arc = 5F*(float)(projection.gz);				
				}
			else
				{
				if (s0 > (s0 + ds*fvol.getNMembofDim(2)))
					VP0 = fvol.getNMembofDim(2) - 1;
				else
					VP0 = 0;
				mangle = 90F;
				nm = (int)(ds*fvol.getNMembofDim(2)/dm);
				VP0projdist = 0F;
				arc = 5F*dm;
				}
			if (nm < 0)
				nm = -nm;
			if (nm == 0)
				nm = 1;
			/*
			v0 = 1500F;
			nv = 1;
			dv = 200F;
			*/
			dixvel = false;
			picktext = "Copy and paste your complete set of NMO\n"
				+ "velocity picks into this window.\n" + pickformat;
			}
		}
	
	/* Produce interpolated velocity picks on stack grid */
	public FltVol makevels()
		{
		int it, im, taddr, saddr, iptdrop, ip, ix, iv;
		int ps, prs, t, pptsp, rs;
		double vel, amp, linedist, arcdist, rat, datz=0;
   		double mindist, tmangle, angledif;
   		int vp0 = VP0;
   		FltVol pvol = new FltVol(nv, nm, ntau);
		float ppts[] = pvol.vec;
		if (fvol == null)
			fvol = pvol;
		/* Parse the pick text into the velinfo array of VInfo objects */
		parseVelinfo();
		FltVol stack = pvol;
		if (dixvel)
			stack.AmpUnits = new String("Dix Interval Velocity");
		else
			stack.AmpUnits = new String("Interpolated Velocity");
		stack.Amp0 = 0F;
		stack.AmpFac = 1F;
		stack.setUnitsofDim(0, "Tau, sec");
		stack.setUnit0ofDim(0, tau0);
		stack.setDUnitofDim(0, dt);
		if (geomexists)
			stack.setUnitsofDim(1, "Dist. N"
			+ Viewmat.getSignif("" + mangle, 4) + "E of VP" + VP0 + ", m");
		else
			stack.setUnitsofDim(1, "Dist. N"
			+ Viewmat.getSignif("" + mangle, 4) + "E of X="
			+ Viewmat.getSignif("" + (s0+ds*VP0), 6) + ", m");
		stack.setUnit0ofDim(1, m0);
		stack.setDUnitofDim(1, dm);
		if (picktext == null)
			{
			stack.setUnitsofDim(2, "m/s NMO Velocity");
			stack.setUnit0ofDim(2, v0);
			stack.setDUnitofDim(2, dv);
			}
		else
			{
			stack.setUnitsofDim(2, "m/s Vnmo Adjustment");
			stack.setUnit0ofDim(2, -((nv-1)/2*dv));
			stack.setDUnitofDim(2, dv);
			}
		/* Set up midpoint line with geometry */
		stack.headers = new RgHead[stack.vecs];
		RgHead vpsp, stktrhdp;
		/* put VP0 coords into vpsp */
		if (geomexists)
			vpsp = getvp(vp0);
		else
			{
			vpsp = new RgHead();
			vpsp.off = 0;
			vpsp.gvp = vpsp.svp = VP0;
			vpsp.gx = vpsp.sx = s0 + ds*VP0;
			vpsp.gy = vpsp.sy = 0;
			vpsp.gz = vpsp.sz = 0;
			vpsp.allset = true;
			}
		if (geomexists)
			{
			for (iv=0; iv<nv; iv++)
				{
				for (ix=0; ix<nm; ix++)
					{
					stack.headers[iv*nm + ix] = new RgHead();
					stktrhdp = stack.headers[iv*nm + ix];
					stktrhdp.off = ix*dm + m0;
					/*vpsp = getvp(vp0);*/
					if (vpsp == null)
						break;
					/* svp and coords in stack are for VP0 */
					stktrhdp.svp = vp0;
					stktrhdp.sx = vpsp.sx;
					stktrhdp.sy = vpsp.sy;
					stktrhdp.sz = vpsp.sz;
					/* receiver coords in stack are on midpoint line */
					stktrhdp.gx = vpsp.sx + stktrhdp.off * Math.sin(mangle);
					stktrhdp.gy = vpsp.sy + stktrhdp.off * Math.cos(mangle);
					/* give each midpoint trace elevation of VP0 datum elevation */
					datz = stktrhdp.gz = vpsp.sz;
					/* find VP closest to midpoint for gvp */
					RgHead vps = fvol.headers[0];
					stktrhdp.gvp = vps.svp;
					mindist = (stktrhdp.gx - vps.sx)*(stktrhdp.gx - vps.sx)
						+ (stktrhdp.gy - vps.sy)*(stktrhdp.gy - vps.sy);
					for (it=0; it<fvol.headers.length; it++)
						{
						vps = fvol.headers[it];
						linedist = (stktrhdp.gx - vpsp.sx)*(stktrhdp.gx - vpsp.sx)
							+ (stktrhdp.gy - vpsp.sy)*(stktrhdp.gy - vpsp.sy);
						if (linedist < mindist)
							{
							mindist = linedist;
							stktrhdp.gvp = vps.svp;
							}
						}
					stktrhdp.allset = true;
					}
				}
			}
		else
			{
			for (iv=0; iv<nv; iv++)
				{
				for (ix=0; ix<nm; ix++)
					{
					/* insert a simple geometry into the stack headers
					where there are no data trace headers */
					stack.headers[iv*nm + ix] = new RgHead();
					stktrhdp = stack.headers[iv*nm + ix];
					stktrhdp.off = ix*dm + m0;
					stktrhdp.svp = VP0;
					stktrhdp.gvp = ix + 1;
					stktrhdp.sx = s0 + ds*VP0;
					stktrhdp.sy = 0;
					stktrhdp.sz = 0;
					stktrhdp.gx = stktrhdp.sx + stktrhdp.off * Math.sin(mangle);
					stktrhdp.gy = stktrhdp.off * Math.cos(mangle);
					stktrhdp.gz = 0;
					stktrhdp.allset = true;
					}
				}
			}
		/* find v increments */
		for (im=0; im<nm; im++)
			{
			/* for CV stacks */
			if (picktext == null)
				{
				for (ip = 0; ip < nv; ip++)
				   {
				   for (it=0; it<ntau; it++)
					   {
					   vel = v0 + ip*dv;
					   ppts[ip*nm*ntau + im*ntau + it] = (float)(vel);
					   }
				   }
				}
			/* for stacks with velocity picks */
			else
				{
				for (ip = 0; ip < nv; ip++)
				   {
				   for (it=0; it<ntau; it++)
					   {
					   /* use velocity picks from nearest distances
					   along midpoint line */
					   vel = ip*dv - (nv-1)/2*dv
						 + velocity(stack.headers[ip*nm + im], tau0+dt*it);
					   ppts[ip*nm*ntau + im*ntau + it] = (float)(vel);
					   }
				   }
				}
			}
   		return pvol;
		}

	/* Constant-Velocity Stack Tests, do not use picks */
	public FltVol cvstack()
		{
		picktext = null;
		return stack();
		}
	/* Stacking from NMO-Velocity Picks */
	public FltVol cmpstack()
		{
		nv = 1;
		return stack();
		}
	/* Main stacking routine, velocity variation tests possible */
	public FltVol stack()
		{
		int nt, noff, plane, nplanes;
		nt = fvol.ne;
		noff = fvol.nv;
		nplanes = fvol.planes;
		if (geomexists)
			{
			fvol.Noticeln("stack: stacking " + nplanes + " labeled records of " 
			+ noff + " traces of " + nt
  			+ " samples each with dt=" + dt);
  			fvol.Noticeln("\t" + fvol.getOffsets()
  			+ " ntau=" + ntau
			+ " tau0=" + tau0);
			}
		else
  			{
  			fvol.Noticeln("stack: stacking " + nplanes + " records of " 
			+ noff + " traces of " + nt
  			+ " samples each with dt=" + dt);
			fvol.Noticeln("\tgeometry from plot parameters:\n\toff0=" + off0
			+ " doff=" + doff
			+ " s0=" + s0
			+ " ds=" + ds
			+ " ntau=" + ntau
			+ " tau0=" + tau0);
			}
		fvol.Noticeln("\tVP0=" + VP0
		+ " mangle=" + mangle
		+ " m0=" + m0
		+ " nm=" + nm
		+ " dm=" + dm
		+ " arc=" + arc);
		if (picktext == null)
			fvol.Noticeln("\tv0=" + v0
			+ " nv=" + nv
			+ " dv=" + dv
			+ "... ");
		else
			fvol.Noticeln("\tUsing " + picktext.length()
			+ " characters of velocity pick text...");
		FltVol stack;
		/*if (nv > 1)*/
			stack = new FltVol(nv, nm, ntau);
		/*else
			stack = new FltPlane(nm, ntau);*/
		stack.AmpUnits = new String("Stacked " + fvol.AmpUnits);
		stack.Amp0 = fvol.Amp0;
		stack.AmpFac = fvol.AmpFac;
		stack.setUnitsofDim(0, "Time, sec");
		stack.setUnit0ofDim(0, tau0);
		stack.setDUnitofDim(0, dt);
		if (geomexists)
			stack.setUnitsofDim(1, "Dist. N"
			+ Viewmat.getSignif("" + mangle, 4) + "E of VP" + VP0 + ", m");
		else
			stack.setUnitsofDim(1, "Dist. N"
			+ Viewmat.getSignif("" + mangle, 4) + "E of X="
			+ Viewmat.getSignif("" + (s0+ds*VP0), 6) + ", m");
		stack.setUnit0ofDim(1, m0);
		stack.setDUnitofDim(1, dm);
		/*if (nv > 1)
			{*/
			if (picktext == null)
				{
				stack.setUnitsofDim(2, "m/s NMO Velocity");
				stack.setUnit0ofDim(2, v0);
				stack.setDUnitofDim(2, dv);
				}
			else
				{
				stack.setUnitsofDim(2, "m/s Vnmo Adjustment");
				stack.setUnit0ofDim(2, -((nv-1)/2*dv));
				stack.setDUnitofDim(2, dv);
				}
			/*}*/
		/* Parse the pick text into the velinfo array of VInfo objects */
		parseVelinfo();
		int it, im, taddr, saddr, iptdrop, ip, ix, iv;
		int ps, prs, t, pptsp, rs;
		double vel, amp, linedist, arcdist, rat;
   		double mindist, tmangle, angledif;
   		int istarttime, halflen, ilen;
   		double drop, aoffset, a, mx, my, mz, sz, gz, datz=0, elcor;
   		double dhld, thld, sqhld, doffset, ddt, dtau0;
   		int vp0 = VP0;
   		mangle *= (float)(Math.PI/180);
   		float ppts[] = new float[nm*nv*ntau];
   		float [] tr = fvol.vec;
		float [] s = stack.vec;
		/* Set up midpoint line with geometry */
		stack.headers = new RgHead[stack.vecs];
		RgHead vpsp, stktrhdp;
		/* put VP0 coords into vpsp */
		if (geomexists)
			vpsp = getvp(vp0);
		else
			{
			vpsp = new RgHead();
			vpsp.off = 0;
			vpsp.gvp = vpsp.svp = VP0;
			vpsp.gx = vpsp.sx = s0 + ds*VP0;
			vpsp.gy = vpsp.sy = 0;
			vpsp.gz = vpsp.sz = 0;
			vpsp.allset = true;
			}
		if (geomexists)
			{
			for (iv=0; iv<nv; iv++)
				{
				for (ix=0; ix<nm; ix++)
					{
					stack.headers[iv*nm + ix] = new RgHead();
					stktrhdp = stack.headers[iv*nm + ix];
					stktrhdp.off = ix*dm + m0;
					/*vpsp = getvp(vp0);*/
					if (vpsp == null)
						break;
					/* svp and coords in stack are for VP0 */
					stktrhdp.svp = vp0;
					stktrhdp.sx = vpsp.sx;
					stktrhdp.sy = vpsp.sy;
					stktrhdp.sz = vpsp.sz;
					/* receiver coords in stack are on midpoint line */
					stktrhdp.gx = vpsp.sx + stktrhdp.off * Math.sin(mangle);
					stktrhdp.gy = vpsp.sy + stktrhdp.off * Math.cos(mangle);
					/* give each midpoint trace elevation of VP0 datum elevation */
					datz = stktrhdp.gz = vpsp.sz;
					/* find VP closest to midpoint for gvp */
					RgHead vps = fvol.headers[0];
					stktrhdp.gvp = vps.svp;
					mindist = (stktrhdp.gx - vps.sx)*(stktrhdp.gx - vps.sx)
						+ (stktrhdp.gy - vps.sy)*(stktrhdp.gy - vps.sy);
					for (it=0; it<fvol.headers.length; it++)
						{
						vps = fvol.headers[it];
						linedist = (stktrhdp.gx - vpsp.sx)*(stktrhdp.gx - vpsp.sx)
							+ (stktrhdp.gy - vpsp.sy)*(stktrhdp.gy - vpsp.sy);
						if (linedist < mindist)
							{
							mindist = linedist;
							stktrhdp.gvp = vps.svp;
							}
						}
					stktrhdp.allset = true;
					}
				}
			}
		else
			{
			for (iv=0; iv<nv; iv++)
				{
				for (ix=0; ix<nm; ix++)
					{
					/* insert a simple geometry into the stack headers
					where there are no data trace headers */
					stack.headers[iv*nm + ix] = new RgHead();
					stktrhdp = stack.headers[iv*nm + ix];
					stktrhdp.off = ix*dm + m0;
					stktrhdp.svp = VP0;
					stktrhdp.gvp = ix + 1;
					stktrhdp.sx = s0 + ds*VP0;
					stktrhdp.sy = 0;
					stktrhdp.sz = 0;
					stktrhdp.gx = stktrhdp.sx + stktrhdp.off * Math.sin(mangle);
					stktrhdp.gy = stktrhdp.off * Math.cos(mangle);
					stktrhdp.gz = 0;
					stktrhdp.allset = true;
					}
				}
			}
		ddt = (double)(dt);
		dtau0 = (double)(tau0);
		float starttime = fvol.getUnit0ofDim(0);
		istarttime = (int)(starttime/dt);
		/* pre-compute velocities */
		pptsp = 0;
		for (im=0; im<nm; im++)
			{
			/* find v increments */
			/* for CV stacks */
			if (picktext == null)
				{
				for (ip = 0; ip < nv; ip++)
				   {
				   for (it=0; it<ntau; it++)
					   {
					   vel = v0 + ip*dv;
					   ppts[pptsp++] = (float)(vel*vel);
					   }
				   }
				}
			/* for stacks with velocity picks */
			else
				{
				for (ip = 0; ip < nv; ip++)
				   {
				   for (it=0; it<ntau; it++)
					   {
					   /* use velocity picks from nearest distances
					   along midpoint line */
					   vel = ip*dv - (nv-1)/2*dv
						 + velocity(stack.headers[ip*nm + im], tau0+dt*it);
					   ppts[pptsp++] = (float)(vel*vel);
					   }
				   }
				}
			}
		/*fvol.Noticeln("velocities assigned");*/

		/* Loop over record planes */
		RgHead tracehd = null;
		int irec, ivol;
		for(irec=0; irec<nplanes; irec++)
		{
		fvol.Noticeln("stack: stacking record " + (irec+1) + " of " + nplanes);
		/* Loop over traces in record */
		for(ix=0; ix<noff; ix++)
			{
			if (ix%20 == 0 || ix == 1)
				fvol.Noticeln("\tstacking trace " + ix);
			/*fvol.Noticeln("getting trace pointers");*/
			ivol = irec*noff + ix;
			t = ivol*nt;
			/*fvol.Noticeln("looking for midpoint");*/
			/* Find midpoint of this trace */
			if (geomexists)
				{
				tracehd = fvol.headers[ivol];
				mx = (tracehd.gx + tracehd.sx)/2;
				my = (tracehd.gy + tracehd.sy)/2;
				mz = (tracehd.gz + tracehd.sz)/2;
				sz = tracehd.sz;
				gz = tracehd.gz;
				}
			else
				{
				mx = s0+ds*irec + (off0+doff*ix)/2;
				my = 0;
				mz = 0;
				sz = 0;
				gz = 0;
				}
			/*fvol.Noticeln("looking for midpoint bin");*/
			/* find bin holding this midpoint */
			linedist = Math.sqrt((mx - vpsp.sx) *(mx - vpsp.sx) +
				 (my - vpsp.sy)*(my - vpsp.sy));
			if (linedist < m0-dm/2 || linedist > m0+(nm-1)*dm+dm/2)
				continue;	/* go to next trace if midpoint is too close
					or too far from VP0 */
			else if (linedist == 0)
				tmangle = -1;
			else if ((mx-vpsp.sx) >= 0 && (my-vpsp.sy) >= 0)
				{
				rat = (mx - vpsp.sx)/linedist;
				if (rat > 1.0)
					rat = 1.0;
				if (rat < -1.0)
					rat = -1.0;
				tmangle = Math.asin(rat);
				}
			else if ((mx-vpsp.sx) >= 0 && (my-vpsp.sy) <= 0)
				{
				rat = (mx - vpsp.sx)/linedist;
				if (rat > 1.0)
					rat = 1.0;
				if (rat < -1.0)
					rat = -1.0;
				tmangle = Math.PI/2 + Math.acos(rat);
				}
			else if ((mx-vpsp.sx) <= 0 && (my-vpsp.sy) <= 0)
				{
				rat = (vpsp.sx-mx)/linedist;
				if (rat > 1.0)
					rat = 1.0;
				if (rat < -1.0)
					rat = -1.0;
				tmangle = Math.PI + Math.asin(rat);
				}
			else if ((mx-vpsp.sx) <= 0 && (my-vpsp.sy) >= 0)
				{
				rat = (my-vpsp.sy)/linedist;
				if (rat > 1.0)
					rat = 1.0;
				if (rat < -1.0)
					rat = -1.0;
				tmangle = 1.5*Math.PI + Math.asin(rat);
				}
			else
				{
				fvol.Noticeln("stack: ignoring trace "
					+ (ix+1) + ": cannot locate midpoint x="
					+ mx + " y=" + my);
				continue;
				}
		
			if (tmangle == -1)
				arcdist = 0;
			else
				{
				angledif = tmangle - mangle;
				angledif = angledif>0 ? angledif : -angledif;
				angledif = angledif>Math.PI ? 2*Math.PI-angledif : angledif;
				arcdist = angledif*linedist;
				}
			if (arcdist > arc/2)
				continue;	/* go to next trace if midpoint arc to defined stacking
					   line is too long (midpoint is outside sectoral bin) */ 
		
			im = (int)((linedist - m0)/dm + 0.5);
			if (im < 0 || im > nm-1)
				{
				fvol.Noticeln("stack: ignoring trace " + (ix+1) + ", get im="
				+ im + " from linedist=" + linedist + ", m0="
				+ m0 + ", dm=" + dm + ", nm=" + nm);
				continue;
				}
			/* hyperbolic stack, use linear interpolation */
			if (geomexists)
				{
				doffset = tracehd.off * tracehd.off;
				aoffset = tracehd.off>0.0 ? tracehd.off : -tracehd.off;
				}
			else
				{
				aoffset = off0 + doff*ix;
				aoffset = aoffset > 0 ? aoffset : -aoffset;
				doffset = aoffset*aoffset;
				}
			/*fvol.Noticeln("midpoint binned");*/
			/*pptsp = 0;*/
			for(ip=0; ip<nv; ip++)
				{
				/* Elevation correction, assuming vertical rays and
				uppermost velocity */
				elcor = (mz - datz)/Math.sqrt(ppts[im*nv*ntau + ip*ntau + 0]);
				ps = ip*nm*ntau + im*ntau;
				/*fvol.Noticeln("starting to stack into trace " + im + " iv " + ip);*/
				for(it=0; it<ntau; it++)
			         	{
					dhld = dtau0 + it*ddt;
					/* get velocity from model and find nmo */
					vel = ppts[im*nv*ntau + ip*ntau + it];
					if (vel == 0)
						{
						ps++;
						continue;
						}
					thld = doffset/vel;
					sqhld = dhld*dhld + thld;
					if (sqhld < 0)
						{
						ps++;
						continue;
						}
					drop = Math.sqrt(sqhld) - dhld + elcor;
				 	iptdrop = (int) (drop/dt);
				 	a = (drop - iptdrop*dt)/dt;
					taddr = (int)(dhld/ddt) + iptdrop - istarttime;
					if (taddr > nt-2 || taddr < 0)
						{
						ps++;
						continue;
						}
					amp = tr[t+taddr] * (1-a) + tr[t+taddr+1] * a;
					s[ps] += amp;
					ps++;
					} /* end of loop over stack trace samples */
				} /* end of loop over velocities */
			} /* end of loop over traces in each record*/
		} /* end of loop over all records */
		fvol.Noticeln("stack: done.");
		return stack;
		}

	/* method to return the coordinates of the given receiver gvp or
	source svp number, copied to all the sx, sy, sz and gx, gy, gz data
	fields of the returned RgHead object,
	from within the RgHead[] fvol.headers list of the data gathers */
	RgHead getvp(int vp)
		{
		RgHead copyloc;
		if (fvol.headers == null)
			return null;
		for (int i=0; i<fvol.headers.length; i++)
			{
			if (fvol.headers[i].svp == vp)
				{
				copyloc = new RgHead(fvol.headers[i]);
				copyloc.gvp = copyloc.svp;
				copyloc.gx = copyloc.sx;
				copyloc.gy = copyloc.sy;
				copyloc.gz = copyloc.sz;
				return copyloc;
				}
			else if (fvol.headers[i].gvp == vp)
				{
				copyloc = new RgHead(fvol.headers[i]);
				copyloc.svp = copyloc.gvp;
				copyloc.sx = copyloc.gx;
				copyloc.sy = copyloc.gy;
				copyloc.sz = copyloc.gz;
				return copyloc;
				}
			}
		return null;
		}

	/*	CMPStack.velocity()
		Routine to return the linearly interpolated velocity between the
		points read into the velinfo array the from picktext. 
		Where the supplied midpoint distance or time is
		outside the limits of those supplied by the velinfo, the velocity
		returned will be that nearest in distance and time. */
	public double velocity(RgHead cmp, double t)
		{
		int m, modelp;
		double vel, thisvel, d1, d2;
		vel = v0;
		if (velinfo == null)
			return vel;
		if (cmp.off <= velinfo[0].dist)
			vel = vpick(t, 0);
		else if (cmp.off >= velinfo[velinfo.length-1].dist)
			vel = vpick(t, velinfo.length-1);
		else
			{
			modelp = 0;
			for (m=0; m < velinfo.length-1; m++)
				{
				if (cmp.off > velinfo[modelp].dist
				&& cmp.off < velinfo[modelp+1].dist)
					{
					thisvel = vpick(t, modelp);
					if (dixvel)
						{
						d1 = cmp.off - velinfo[modelp].dist;
						d1 = d1 < 0 ? -d1 : d1;
						d2 = velinfo[modelp+1].dist - cmp.off;
						d2 = d2 < 0 ? -d2 : d2;
						if (d1 < d2)
							vel = thisvel;
						else
							vel = vpick(t, modelp+1);
						}
					else
						vel = thisvel + (cmp.off - velinfo[modelp].dist)/
							(velinfo[modelp+1].dist - velinfo[modelp].dist)
							*(vpick(t, modelp+1) - thisvel);
					break;
					}
				else if (cmp.off == velinfo[modelp].dist)
					{
					vel = vpick(t, modelp);
					break;
					}
				modelp++;
				}
			}
		return(vel);
		}
	
	/*	Vpick.c
		Routine to return the linearly interpolated velocity pick between the
		points read into a velinfo array, when supplied with the time and
		the index into the velinfo array nearest the midpoint distance. The
		5 nearest picks are averaged according to the inverse of the distance
		and time-distance squared */
	public double vpick(double t, int modp)
		{
		int p, nnear, range0, rangef, pmaxtau, pmintau;
		nnear = 7/*5*/;
		double temp, dist, vel, maxtau, mintau;
		if (t == velinfo[modp].tau && !dixvel)
			return velinfo[modp].vel;
		if (velinfo.length <= nnear)
			{
			nnear = velinfo.length;
			range0 = 0;
			rangef = velinfo.length - 1;
			}
		else if (modp < nnear/2)
			{
			range0 = 0;
			rangef = nnear - 1;
			}
		else if (modp > velinfo.length - nnear/2 - 1)
			{
			range0 = velinfo.length - nnear;
			rangef = velinfo.length - 1;
			}
		else
			{
			range0 = modp - nnear/2;
			rangef = modp + nnear/2;
			}
		pmaxtau = range0;
		pmintau = range0;
		maxtau = velinfo[pmaxtau].tau;
		mintau = velinfo[pmintau].tau;
		for (p=range0; p<=rangef; p++)
			{
			if (velinfo[p].tau > maxtau)
				{
				maxtau = velinfo[p].tau;
				pmaxtau = p;
				}
			if (velinfo[p].tau < mintau)
				{
				mintau = velinfo[p].tau;
				pmintau = p;
				}
			}
		if (t <= mintau && dixvel)
			return velinfo[pmintau].vel;
		if (t >= maxtau)
			return velinfo[pmaxtau].vel;
		if (dixvel)
			{
			maxtau = maxtau - t;
			mintau = t - mintau;
			for (p=range0; p<=rangef; p++)
				{
				if (t > velinfo[p].tau && (t - velinfo[p].tau) < mintau)
					{
					mintau = t - velinfo[p].tau;
					pmintau = p;
					}
				if (t <= velinfo[p].tau && (velinfo[p].tau - t) < maxtau)
					{
					maxtau = velinfo[p].tau - t;
					pmaxtau = p;
					}
				}
			vel = dixeqn(pmintau, pmaxtau);
			}
		else
			{
			vel = dist = 0;
			for (p=range0; p<=rangef; p++)
				{
				/*
				if (t >= maxtau && velinfo[p].vel < velinfo[pmaxtau].vel)
					continue;
				*/
				temp = (velinfo[modp].dist - velinfo[p].dist)
					*(velinfo[modp].dist - velinfo[p].dist)
					+ ((t - velinfo[p].tau)*(t - velinfo[p].tau)
					*(velinfo[p].vel*velinfo[p].vel));
				if (temp == 0)
					return velinfo[p].vel;
				vel += velinfo[p].vel/temp;
				dist += 1/temp;
				}
			vel /= dist;
			}
		return vel;
		}
	
	public double dixeqn(int i1, int i2)
		{
		int tmp;
		double rad;
		if (velinfo[i1].tau > velinfo[i2].tau)
			{
			tmp = i1;
			i1 = i2;
			i2 = tmp;
			}
		else if (velinfo[i1].tau == velinfo[i2].tau)
			{
			return 0;
			}
		rad = (velinfo[i2].vel*velinfo[i2].vel*velinfo[i2].tau
			-  velinfo[i1].vel*velinfo[i1].vel*velinfo[i1].tau)
			/(velinfo[i2].tau - velinfo[i1].tau);
		if (rad <= 0)
			return 0;
		return Math.sqrt(rad);
		}
	
	public void parseVelinfo()
		{
		StringBufferInputStream sbis;
		InputStreamReader isr;
		BufferedReader br;
		StringTokenizer st;
		VInfo temp;
		String line;
		int npicks;
		if (picktext == null)
			return;
		npicks = 0;
		/* Count readable lines in picktext */
		sbis = new StringBufferInputStream(picktext);
		isr = new InputStreamReader(sbis);
		br = new BufferedReader(isr);
		try
			{
			while ((line = br.readLine()) != null)
				{
				if (line.length() < 13)
					continue;
				st = new StringTokenizer(line);
				temp = parseVInfoLine(st);
				if (temp == null)
					continue;
				npicks++;
				}
			}
		catch (IOException ioe)
			{
			fvol.Errout("Internal IOException reading velocity pick text;\n"
				+ "attempting to continue:\n" + ioe);
			}
		if (npicks < 1)
			{
			fvol.Errout("No valid lines of velocity picks could be read,\n"
				+ "using a constant velocity of " + v0 + " m/s.\n\n"
				+ "Copy and paste your complete set of NMO\n"
				+ "velocity picks into the pick text window.\n"
				+ pickformat);
			picktext = null;
			return;
			}
		/* dimension array for VInfo objects */
		velinfo = new VInfo[npicks];
		/* read the picktext and parse into array */
		sbis = new StringBufferInputStream(picktext);
		isr = new InputStreamReader(sbis);
		br = new BufferedReader(isr);
		try
			{
			int i=0;
			while ((line = br.readLine()) != null)
				{
				if (line.length() < 13)
					continue;
				st = new StringTokenizer(line);
				temp = parseVInfoLine(st);
				if (temp == null)
					continue;
				velinfo[i++] = new VInfo(temp);
				}
			}
		catch (IOException ioe)
			{
			fvol.Errout("Internal IOException reading velocity pick text;\n"
				+ "attempting to continue:\n" + ioe);
			}
		/* Now sort the velinfo array, first by distance, second by tau */
		sortVelinfo(0, velinfo.length - 1, true);
		sortVelinfo(0, velinfo.length - 1, false);
		}
	/* standard quicksort for velinfo arrays - if sorttau==true, sorts on tau
	field; if false, sorts on dist field */
	public void sortVelinfo(int left, int right, boolean sorttau)
		{
		int i, last;
		if (left >= right)
			return;
		swap(left, (left+right)/2);
		last = left;
		for (i=left+1; i<=right; i++)
			{
			if (sorttau)
				{
				if (velinfo[i].tau < velinfo[left].tau)
					swap(++last, i);
				}
			else
				{
				if (velinfo[i].dist < velinfo[left].dist)
					swap(++last, i);
				}
			}
		swap(last, left);
		sortVelinfo(left, last-1, sorttau);
		sortVelinfo(last+1, right, sorttau);
		}
	/* swap elements in a velinfo array */
	public void swap(int d1, int d2)
		{
		VInfo tmp = velinfo[d1];
		velinfo[d1] = velinfo[d2];
		velinfo[d2] = tmp;
		}
	/* parse a line of pick text and return VInfo object if valid pick,
	return null if not a valid pick line */
	public VInfo parseVInfoLine(StringTokenizer st)
	throws IOException
		{
		double dist, tau, vel, temp;
		int itemp;
		if (st.hasMoreTokens()) /* Amp column */
			{
			try {temp = Double.valueOf(st.nextToken()).doubleValue(); }
			catch (NumberFormatException nfe) {return null; }
			}
		else
			return null;
		if (st.hasMoreTokens()) /* tau index column */
			{
			try {itemp = Integer.valueOf(st.nextToken()).intValue(); }
			catch (NumberFormatException nfe) {return null; }
			}
		else
			return null;
		if (st.hasMoreTokens()) /* tau column */
			{
			try {tau = Double.valueOf(st.nextToken()).doubleValue(); }
			catch (NumberFormatException nfe) {return null; }
			}
		else
			return null;
		if (st.hasMoreTokens()) /* Mdist index column */
			{
			try {itemp = Integer.valueOf(st.nextToken()).intValue(); }
			catch (NumberFormatException nfe) {return null; }
			}
		else
			return null;
		if (st.hasMoreTokens()) /* Mdist column */
			{
			try {dist = Double.valueOf(st.nextToken()).doubleValue(); }
			catch (NumberFormatException nfe) {return null; }
			}
		else
			return null;
		if (st.hasMoreTokens()) /* Vel index column */
			{
			try {itemp = Integer.valueOf(st.nextToken()).intValue(); }
			catch (NumberFormatException nfe) {return null; }
			}
		else
			return null;
		if (st.hasMoreTokens()) /* Vel column */
			{
			try {vel = Double.valueOf(st.nextToken()).doubleValue(); }
			catch (NumberFormatException nfe) {return null; }
			}
		else
			return null;
		if (vel <= 0.0)
			return null;
		VInfo vi = new VInfo();
		vi.set(dist, tau, vel);
		return vi;
		}

	/* End of class CMPStack */
	}

class VInfo
	{
	double dist;
	double tau;
	double vel;
	
	/* Constructor */
	VInfo()
		{
		this.dist = 0;
		this.tau = 0;
		this.vel = 0;
		}
	/* Construct a copy */
	VInfo(VInfo vi)
		{
		this.dist = vi.dist;
		this.tau = vi.tau;
		this.vel = vi.vel;
		}
	void set(double dist, double tau, double vel)
		{
		this.dist = dist;
		this.tau = tau;
		this.vel = vel;
		}
	}