/* 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.io.*;

public class RgHead
	{
	public int rec;				/* field record number */
	public int svp;				/* source location number */
	public int gvp;				/* receiver location number */
	public double sx;		/* source x-coordinate, m */
	public double sy;		/* source y-coordinate, m */
	public double sz;		/* source z-coordinate, m */
	public double gx;		/* receiver x-coordinate, m */
	public double gy;		/* receiver y-coordinate, m */
	public double gz;		/* receiver z-coordinate, m */
	public double off;			/* offset, m, may be signed */
	public byte unas[];			/* 3 unassigned bytes to keep 4-byte alignment */
	public boolean allset;		/* false if not all values set */

	public static final int RGHEADBYTES = 18*4;
	public static double RAD=0.0174532;
	public static double DAR=57.29578;
	public static double COR=0.993306;	/* COR= (1-e)**2, where e= ecentricity */
						/* COR= 0.993306 for e=1.0/289.78 */
	public static double TOL=0.001;
	public static double RTOL=0.001;
	public static double SQRT2=1.414213562;
	public static double PI=3.1415926;
	public static double FOURPI=12.56637062;
	public static double MPDEG=111190.0;
	
	/* Constructors */
	RgHead()
		{
		unas = new byte[3];
		allset = false;
		}
	/* Construct a copy */
	RgHead(RgHead hd)
		{
		this.rec = hd.rec;
		this.svp = hd.svp;
		this.gvp = hd.gvp;
		this.sx = hd.sx;
		this.sy = hd.sy;
		this.sz = hd.sz;
		this.gx = hd.gx;
		this.gy = hd.gy;
		this.gz = hd.gz;
		this.off = hd.off;
		this.unas = new byte[3];
		this.unas[0] = hd.unas[0];
		this.unas[1] = hd.unas[1];
		this.unas[2] = hd.unas[2];
		this.allset = hd.allset;
		}
	/* Construct a header from FltPlane dimensional attributes */
	RgHead(int iv, FltPlane fp)
		{
		int np=1, nv, ip;
		nv = fp.vecs;
		if (fp.getNMembofDim(2) > 1)
			{
			np = fp.getNMembofDim(2);
			nv = fp.vecs/np;
			if (np*nv != fp.vecs)
				{
				np = 1;
				nv = fp.vecs;
				}
			}
		ip = iv/nv;
		this.rec = ip;
		this.svp = ip;
		this.gvp = this.svp + iv - ip*nv;
		this.sx = fp.getUnit0ofDim(2) + ip*fp.getDUnitofDim(2);
		this.sy = 0.0;
		this.sz = 0.0;
		this.off = fp.getUnit0ofDim(1) + (iv - (ip*nv))*fp.getDUnitofDim(1);
		this.gx = this.sx + this.off;
		this.gy = 0.0;
		this.gz = 0.0;
		this.unas = new byte[3];
		this.allset = false;
		}

	/* Find a projection line through source-point coordinates or midpoints.
	   Returned svp and gvp set to location closest to origin;
	   sx, sy, sz give origin coordinates;
	   gx is Slope on map in meters North/meters East;
	   gy is Direction in degrees East of North;
	   gz is the distance in meters from the origin to the closest station
	   off is the length of the projection */
	public void getProjection(RgHead[] headers)
		{
		/* set defaults for a projection */
		this.rec = 0;
		this.allset = true;
		/* Compute projection line from sx, sy in headers */
		double sx = 0;
		double sx2 = 0;
		double sy = 0;
		double sy2 = 0;
		double sxy = 0;
		double m, b, ang, temp;
		RgHead[]  coord = headers;
		int ncoord = coord.length;
		double westx = coord[0].sx;
		double westy = coord[0].sy;
		for (int icoord=0; icoord<ncoord; icoord++)
			{
			sx += coord[icoord].sx;
			sx2 += coord[icoord].sx * coord[icoord].sx;
			sy += coord[icoord].sy;
			sy2 += coord[icoord].sy * coord[icoord].sy;
			sxy += coord[icoord].sx * coord[icoord].sy;
			if (coord[icoord].sx < westx)
				{
				westx = coord[icoord].sx;
				westy = coord[icoord].sy;
				}
			}
		if (sx2 - sx*sx/ncoord == 0)
			{
			/* S-N Projection Line */
			m = 0;
			b = 0;
			ang = 0;
			this.sx = sx/ncoord;
			this.sy = coord[0].sy;
			this.off = coord[0].sy;
			for (int i=1; i<ncoord; i++)
				{
				if (coord[i].sy < this.sy)
					this.sy = coord[i].sy;
				if (coord[i].sy > this.off)
					this.off = coord[i].sy;
				}
			this.off -= this.sy;
			}
		else
			{
			m = (sxy - sx*sy/ncoord)/(sx2 - sx*sx/ncoord);
			b = sy/ncoord - m*sx/ncoord;
			this.off = coord[0].sx;
			for (int i=1; i<ncoord; i++)
				{
				if (coord[i].sx > this.off)
					this.off = coord[i].sx;
				}
			if (m == 0.0)
				{
				/* W-E Projection Line */
				ang = 90;
				this.sy = b;
				this.sx = westx;
				this.off -= this.sx;
				}
			else
				{
				temp = (b + westx/m -westy)/(1.0/m - m);
				this.sx = temp;
				this.sy = m*temp + b;
				ang = Math.atan2(1.0, m)/RAD;
				temp = m*this.off + b;
				this.off = Math.sqrt((this.off - this.sx)*(this.off - this.sx)
					+ (temp - this.sy)*(temp - this.sy));
				}
			}
		this.gx = m;
		this.gy = ang;
		/* find point in header closest to projection origin */
		this.svp = coord[0].svp;
		this.sz = coord[0].sz;
		this.gz = (this.sx - coord[0].sx)*(this.sx - coord[0].sx)
			+ (this.sy - coord[0].sy)*(this.sy - coord[0].sy);
		for (int i=1; i<ncoord; i++)
			{
			temp = (this.sx - coord[i].sx)*(this.sx - coord[i].sx)
				+ (this.sy - coord[i].sy)*(this.sy - coord[i].sy);
			if (temp < this.gz)
				{
				this.svp = coord[i].svp;
				this.gz = temp;
				this.sz = coord[i].sz;
				}
			}
		this.gz = Math.sqrt(this.gz);
		this.gvp = this.svp;
		}

	/* Fill an RgHead array with midpoint locations for
	   both source and receiver coordinates */
	public static RgHead[] getMidpointArray(RgHead[] headers)
		{
		RgHead midpoints[] = new RgHead[headers.length];
		for (int i=0; i<headers.length; i++)
			midpoints[i] = headers[i].getMidpoint();
		return midpoints;
		}
	public RgHead getMidpoint()
		{
		RgHead midpoint = new RgHead(this);
		midpoint.gx = midpoint.sx = (this.sx + this.gx)/2;
		midpoint.gy = midpoint.sy = (this.sy + this.gy)/2;
		midpoint.gz = midpoint.sz = (this.sz + this.gz)/2;
		return midpoint;
		}
	/* Recompute source-receiver offsets */
	public double calcOffset(boolean geogr)
		{
		double off;
		if (geogr)
			off = MPDEG*getDegreesDist(sy, sx, gy, gx);
		else
			off = Math.sqrt((sx-gx)*(sx-gx)
				+ (sy-gy)*(sy-gy) + (sz-gz)*(sz-gz));
		return off;
		}
	/* distaz originally by Charles Ammon, 1988, Penn State */
	public static double getDegreesDist(double lat1, double lon1,
		double lat2, double lon2)
		{
		double delta, az, test, azp;
		/*
		lat1= atan( tan( RAD*lat1) /0.992277 );
		lat2= atan( tan( RAD*lat2) /0.992277 );
		*/
		lat1 *= RAD;
		lat2 *= RAD;
		lon1 *= RAD;
		lon2 *= RAD;
	
		delta = Math.acos(Math.sin(lat1)*Math.sin(lat2) 
			+ Math.cos(lat1)*Math.cos(lat2)*Math.cos(lon2-lon1));
	
		test= Math.sin(lon2 - lon1);
		azp= DAR*Math.acos( (Math.sin(lat2) -Math.sin(lat1)*Math.cos( delta))
			/( Math.sin( delta)*Math.cos(lat1)));
	
		az= (test < 0.0 ? 360.0- azp : azp);
		delta *= DAR;
		return delta;
		}
	public void set(int rec, int svp, int gvp, double sx, double sy, double sz,
		double gx, double gy, double gz, double off)
		{
		this.rec = rec;
		this.svp = svp;
		this.gvp = gvp;
		this.sx = sx;
		this.sy = sy;
		this.sz = sz;
		this.gx = gx;
		this.gy = gy;
		this.gz = gz;
		this.off = off;
		this.allset = true;
		}
	public void copy(RgHead hd)
		{
		this.rec = hd.rec;
		this.svp = hd.svp;
		this.gvp = hd.gvp;
		this.sx = hd.sx;
		this.sy = hd.sy;
		this.sz = hd.sz;
		this.gx = hd.gx;
		this.gy = hd.gy;
		this.gz = hd.gz;
		this.off = hd.off;
		this.unas[0] = hd.unas[0];
		this.unas[1] = hd.unas[1];
		this.unas[2] = hd.unas[2];
		this.allset = hd.allset;
		}
	public void copySEGYTraceHead(SEGYTraceHead thd, SEGYFileHead fhd)
		{
		float lscl, coscl, metscl;
		thd.regularize();
		fhd.regularize();
		if (thd.lscaler < (short)(0))
			lscl = 1F/(float)(-thd.lscaler);
		else
			lscl = (float)(thd.lscaler);
		if (thd.coscaler < (short)(0))
			coscl = 1F/(float)(-thd.coscaler);
		else
			coscl = (float)(thd.coscaler);
		if (fhd.units == (short)(2) && thd.counit != (short)(2))
			metscl = 0.3048F;
		else
			metscl = 1F;
		this.rec = thd.frecno;
		this.sx = thd.sx*coscl*metscl;
		this.sy = thd.sy*coscl*metscl;
		this.sz = lscl*(thd.szsurf - thd.sdepth)*metscl;
		this.gx = thd.gx*coscl*metscl;
		this.gy = thd.gy*coscl*metscl;
		this.gz = thd.gz*lscl*metscl;
		this.off = thd.offset;
		if (fhd.units == (short)(2))
			this.off *= 0.3048;
		this.svp = thd.sp;
		this.gvp = thd.frtr1gno + thd.trnof - 1;
		this.allset = true;
		}
	public void copyOldRG(OldRGTraceHead rgthd)
		{
		this.off = rgthd.offset;		/* source-receiver offset, m */
		this.gx = rgthd.gx;		/* receiver x-coordinate, m */
		this.gy = rgthd.gy;		/* receiver y-coordinate, m */
		this.gz = rgthd.gz;		/* receiver z-coordinate, m */
		this.svp = rgthd.sp;			/* source number */
		}
	public void read(DataInputStream dis) throws IOException
		{
		this.rec = dis.readInt();
		this.svp = dis.readInt();
		this.gvp = dis.readInt();
		this.sx = dis.readDouble();
		this.sy = dis.readDouble();
		this.sz = dis.readDouble();
		this.gx = dis.readDouble();
		this.gy = dis.readDouble();
		this.gz = dis.readDouble();
		this.off = dis.readDouble();
		this.unas[0] = dis.readByte();
		this.unas[1] = dis.readByte();
		this.unas[2] = dis.readByte();
		this.allset = dis.readBoolean();
		}
	public void write(DataOutputStream dos) throws IOException
		{
		dos.writeInt(rec);
		dos.writeInt(svp);
		dos.writeInt(gvp);
		dos.writeDouble(sx);
		dos.writeDouble(sy);
		dos.writeDouble(sz);
		dos.writeDouble(gx);
		dos.writeDouble(gy);
		dos.writeDouble(gz);
		dos.writeDouble(off);
		dos.writeByte(unas[0]);
		dos.writeByte(unas[1]);
		dos.writeByte(unas[2]);
		dos.writeBoolean(allset);
		}
	public void print(PrintStream p)
		{
		p.println(toString());
		}
	public void printLine(PrintStream p)
		{
		p.println(toLine());
		}
	public String toString()
		{
		return "rec=" + rec + " svp=" + svp + " gvp=" + gvp
			+ " sx=" + Viewmat.getSignif("" + sx)
			+ " sy=" + Viewmat.getSignif("" + sy)
			+ " sz=" + Viewmat.getSignif("" + sz)
			+ " gx=" + Viewmat.getSignif("" + gx)
			+ " gy=" + Viewmat.getSignif("" + gy)
			+ " gz=" + Viewmat.getSignif("" + gz)
			+ " off=" + Viewmat.getSignif("" + off)
			+ " allset=" + allset;
		}
	public String toLine()
		{
		return "" + rec + " " + svp + " " + gvp
			+ " " + Viewmat.getSignif("" + sx, 12)
			+ " " + Viewmat.getSignif("" + sy, 12)
			+ " " + Viewmat.getSignif("" + sz, 12)
			+ " " + Viewmat.getSignif("" + gx, 12)
			+ " " + Viewmat.getSignif("" + gy, 12)
			+ " " + Viewmat.getSignif("" + gz, 12)
			+ " " + Viewmat.getSignif("" + off, 12);
		}
	}