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

License is hereby granted for anyone to make any use of this software and
its associated source and compiled binary codes, including free re-distribution,
and incorporation into commercial products, provided this copyright notice
is retained within all files, versions, and distribution media substantially
incorporating John Louie's original work.

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.*;
import java.util.*;

class testsegy
	{
	public static void main(String args[])
		{
		SEGYFileHead sfhd = new SEGYFileHead();
		FltVol fvol;
		String path = "/crunch/s0/louie/d.reno/d.taber/seatoun-refr/refl1001-1029.sgy";
		sfhd.printheader = true;
		fvol = sfhd.importSEGYFile(path, true, false, 1028, 1029, 0);
		fvol.viewinFrame("Seatoun Reflection", true, 0.04F, RgCtab.rwb, false, -1F);
		fvol.f.Errout(sfhd.toString());
		}
	}

public class SEGYFileHead
	{
	public static int SEGY_THEAD_SIZE=3600;	/* segy tape header blocks total in bytes */
	public static int SEGY_CARD_SIZE=3200;	/* segy tape header EBCDIC card block size in bytes */
	public static int SEGY_HEAD_SIZE=240;	/* segy header size in bytes */
	public Vector addresses, headers;
	public RgFileHead fhd;
	public boolean printheader;
	public boolean geogr;

	public byte cards[];	/* SEG-Y Reel Header Part 1: 3200-byte EBCDIC card image block */
	public int jobid;	/* Job identification number */
	public int lineno;	/* * Line number */
	public int reelno;	/* * Tape reel number */
	public short ntrace;	/* * no. of data traces per record */
	public short naux;	/* * no. of auxiliary traces per record */
	public short mudt;	/* * sample interval of these data, microseconds */
	public short omudt;	/* sample interval of original field data, microseconds */
	public short nt;	/* * no. of samples per data trace here; 32767 max - see num_samps below */
	public short ont;	/* no. of samples per original field data trace */
	public short samptype;	/* * data sample format code: 1= 4-byte float; 2= 4-byte fixed-point
				     3= 2-byte fixed-point; 4= 4-byte fixed point with gain code
				     Adding: 0= 4-byte IEEE/Java/Sun/SGI/HP float */
	public short fold;	/* * expected CDP fold */
	public short sort;	/* Trace sorting code: 1= as recorded; 2= CDP ensemble,
				   3= single-fold coninuous profile, 4= horizontally stacked */
	public short sums;	/* Vertical sum code (how many summed) */
	public short sweep0;	/* Sweep frequency start, Hz */
	public short sweepf;	/* Sweep frequency end, Hz */
	public short sweepl;	/* Sweep length, milliseconds */
	public short sweept;	/* Sweep type code: 1=linear; 2=parabolic; 3=exponential; 4=other */
	public short sweeptr;	/* Trace number of sweep channel */
	public short taper0;	/* Sweep trace taper length at start, milliseconds */
	public short taperf;	/* Sweep trace taper length at end, milliseconds */
	public short tapert;	/* Sweep taper type: 1=linear; 2=cos^2; 3=other */
	public short cor;	/* Correlated data traces: 1=no; 2=yes */
	public short gain;	/* Binary gain recovered: 1=yes; 2=no */
	public short recov;	/* Gain recovery method: 1=none; 2= spherical divergence; 3=AGC;
				   4=other */
	public short units;	/* * Measurement system units: 1=meters; 2=feet */
	public short pol;	/* Impulse signal polarity: 1= Increase in pressure or upward
				   geophone case movement gives negative sample value; or
				   2= positive sample value */
	public short vibpol;	/* Vibratory polarity code (seismic signal lags pilot signal by):
				   1=337.5-22.5 deg; 2; 3; 4; 5; 6; 7; 8=292.5-337.5 deg */
	public byte unas[];	/* unassigned, optional; here we use 336 bytes */
	public int num_samps;	/* last officially unassigned 4 bytes go to
						    hold trace length larger than 32767 */

	
	/* Constructor */
	public SEGYFileHead()
		{
		cards = new byte[SEGY_CARD_SIZE];
		jobid = 1;
		lineno = 1;
		reelno = 1;
		samptype = (short)(1);
		fold = (short)(1);
		sort = (short)(1);
		sums = (short)(1);
		cor = (short)(1);
		gain = (short)(1);
		recov = (short)(1);
		units = (short)(1);
		pol = (short)(2);
		unas = new byte[336];
		addresses = null;
		headers = null;
		fhd = null;
		printheader = false;
		geogr = false;
		}

	public int getnt()
		{
		if (nt == (short)(-1))
			return -1;
		else if ((nt < 1 || nt >= (short)(32767))
			&& num_samps > 0 && num_samps < 2000000000)
			return num_samps;
		else
			return nt;
		}
	public void setnt(int n)
		{
		if (n == -1)
			nt = (short)(-1);
		else if (n >= 32767)
			nt = (short)(32767);
		else
			nt = (short)(n);
		num_samps = n;
		}
	public void copyRgFileHead(RgFileHead fhd)
		{
		this.ntrace = (short)(fhd.ntrace);
		/*this.nt = (short)(fhd.nt);*/
		this.setnt(fhd.nt);
		this.mudt = (short)(1e6F*fhd.dt);
		if (this.units != (short)(2))
			this.units = (short)(1);
		}
	public String toString()
		{
		return "ntrace=" + ntrace
			+ " nt=" + nt
			+ " num_samps=" + num_samps
			+ " lineno=" + lineno
			+ " reelno=" + reelno
			+ " naux=" + naux
			+ " mudt(microsec)=" + mudt
			+ " samptype=" + samptype
			+ " fold=" + fold
			+ " units=" + units;
		}
	public boolean isValid()
		{
		if (ntrace < (short)(1) || ntrace > (short)(32000)
		/*|| nt < (short)(1) || nt > (short)(32000)*/
		|| getnt() < 1 || getnt() > 2000000000
		|| mudt < (short)(1) || mudt > (short)(32000) || mudt == (short)(0)
		|| units < (short)(0) || units > (short)(3))
			return false;
		return true;
		}
	public int scanSEGYFile(String path, boolean tapehead, boolean swabhead,
		int firstrec, int lastrec, int samptype)
		{
		return scanSEGYFile(path, tapehead, swabhead, firstrec, lastrec,	
			samptype, (ViewmatFrame)(null));
		}
	public int scanSEGYFile(String path, boolean tapehead, boolean swabhead,
		int firstrec, int lastrec, int samptype, ViewmatFrame win)
		{
		int countrec, i;
		long bytes;
		boolean eof;
		RandomAccessFile raf;
		SEGYTraceHead thd = new SEGYTraceHead();
		fhd = new RgFileHead();
		RgHead rgthd, temphd;
		rgthd = new RgHead();
		if (firstrec > -1 && lastrec > -1)
			System.out.println("Searching SEG-Y file " + path
			+ "\nfor record numbers " + firstrec + " through " + lastrec);
		else if (firstrec < 0 && lastrec > -1)
			System.out.println("Searching SEG-Y file " + path
			+ "\nfor record numbers <= " + lastrec);
		else if (lastrec < 0 && firstrec > -1)
			System.out.println("Searching SEG-Y file " + path
			+ "\nfor record numbers >= " + firstrec);
		else
			System.out.println("Searching SEG-Y file " + path
			+ "\nfor all record numbers.");
		if (tapehead)
			{
			readSEGYFileHead(path, swabhead);
			regularize();
			if (!isValid())
				{
				System.out.println("File " + path
				+ "\ndoes not begin with the 3600-byte SEG-Y reel header,"
				+ "\nor the reel header could not be read.");
				return 0;
				}
			addresses = new Vector(ntrace, ntrace);
			headers = new Vector(ntrace, ntrace);
			fhd.copySEGYFileHead(this);
			}
		else
			{
			fhd.ntrace = -1;
			fhd.nt = -1;
			fhd.dt = -1F;
			addresses = new Vector(48, 48);
			headers = new Vector(48, 48);
			}
		fhd.starttime = -9e25F;
		try
			{
			raf = new RandomAccessFile(path, "r");
			}
		catch(IOException ioe)
			{
			System.out.println("Cannot open file " + path
				+ "\nfor random read access: " + ioe);
			return 0;
			}
		catch(SecurityException se)
			{
			System.out.println("Read access to file " + path
				+ "\nnot permitted: " + se);
			return 0;
			}
		bytes = (long)(0);
		/* Skip over SEG-Y tape header blocks if indicated */
		if (tapehead)
			{
			try
				{
				raf.seek((long)(SEGY_THEAD_SIZE));
				}
			catch(IOException ioe)
				{
				System.out.println("Cannot seek in file " + path
					+ "\npast SEG-Y reel header: " + ioe);
				return 0;
				}
			bytes += SEGY_THEAD_SIZE;
			}
		eof = false;
		/* Read SEG-Y trace headers to EOF */
		for (i=0; !eof && i<1000000; i++)
			{
			try
				{
				thd.read(raf);
				}
			catch(EOFException eofe)
				{
				eof = true;
				break;
				}
			catch(IOException ioe)
				{
				System.out.println("IO Exception reading file " + path
					+ "\nfor SEG-Y trace header, i=" + i + ": " + ioe);
				return 0;
				}
			bytes += SEGY_HEAD_SIZE;
			/* Convert locations, using coordinate scaler
			   from header, and convert to meters if feet */
			rgthd.copySEGYTraceHead(thd, this);
			this.geogr = thd.isGeogr();
			if ((firstrec < 0 && lastrec < 0)
				|| (firstrec < 0 && rgthd.rec <= lastrec)
				|| (lastrec < 0 && rgthd.rec >= firstrec)
				|| (rgthd.rec >= firstrec && rgthd.rec <= lastrec))
				{
				if (fhd.starttime < -1e25F)
					fhd.starttime = thd.delay/1e3F;
				if (fhd.nt < 1)
					{
					fhd.nt = thd.getnt();
					fhd.dt = thd.dt/1e6F;
					this.copyRgFileHead(fhd);
					}
				else if ((short)(fhd.starttime*1e3F) != thd.delay
				|| this.mudt != thd.dt
				|| this.getnt() != thd.getnt())
					continue;
				addresses.addElement(new Long(bytes));
				headers.addElement(new RgHead(rgthd));
				if (printheader)
					System.out.println("header " + i + ": " + thd);
				}
			/* Skip over trace */
			bytes += thd.getnt()*(this.samptype==3 ? 2 : 4);
			try
				{
				raf.seek(bytes);
				}
			catch(IOException ioe)
				{
				System.out.println("Cannot seek in file " + path
					+ "\npast SEG-Y reel header to "
					+ bytes + " bytes: " + ioe);
				return 0;
				}
			if (i==0 || i==1 || i==9 || i==99 || i%1000 == 999)
				System.out.println("Read SEG-Y trace header " + (i+1));
			if (printheader)
				System.out.println("header " + i + ": " + thd);
			}
		System.out.println("" + i + " trace headers read from file " + path);
		if (addresses.size() < 1)
			{
			if (firstrec > -1 && lastrec > -1)
				System.out.println("No traces found for record numbers "
				+ firstrec + " through " + lastrec);
			else if (firstrec < 0 && lastrec > -1)
				System.out.println("No traces found for record numbers <= "
				+ lastrec);
			else if (lastrec < 0 && firstrec > -1)
				System.out.println("No traces found for record numbers >= "
				+ firstrec);
			else
				System.out.println("No traces found for any record number.");
			return 0;
			}
		if (firstrec > -1 && lastrec > -1)
			System.out.println("" + addresses.size()
			+ " traces found for record numbers " + firstrec
			+ " through " + lastrec);
		else if (firstrec < 0 && lastrec > -1)
			System.out.println("" + addresses.size()
			+ " traces found for record numbers <= " + lastrec);
		else if (lastrec < 0 && firstrec > -1)
			System.out.println("" + addresses.size()
			+ " traces found for record numbers >= " + firstrec);
		else
			System.out.println("" + addresses.size()
			+ " traces found (using any record number).");
		boolean getall;
		if (lastrec < 0 && firstrec < 0)
			getall = true;
		else
			getall = false;
		addresses.trimToSize();
		headers.trimToSize();
		try
			{
			raf.close();
			}
		catch(IOException ioe)
			{
			System.out.println("Could not close file " + path
				+ "\nafter random read access.");
			}

		rgthd = (RgHead)(headers.elementAt(0));
		countrec = 0;
		for (i=0; i<headers.size(); i++)
			{
			temphd = (RgHead)(headers.elementAt(i));
			if (temphd.rec == rgthd.rec)
				countrec++;
			}
		if (this.naux < 0)
			this.naux = (short)(0);
		if (countrec != (this.ntrace + this.naux))
			{
			if (getall)
				{
				this.ntrace = (short)(headers.size());
				this.naux = (short)(0);
				}
			else
				{
				if (tapehead)
					System.out.println("SEG-Y reel header lists "
					+ this.ntrace + " data traces + " + this.naux
					+ " auxiliary traces per record; first record in file has "
					+ countrec + " traces.");
				this.ntrace = (short)(countrec);
				this.naux = (short)(0);
				}
			}
		if (!tapehead)
			{
			lineno = 1;
			reelno = 1;
			naux = (short)(0);
			if (samptype < 0 || samptype > 4)
				this.samptype = (short)(1);
			else if (samptype == 2)
				this.samptype = (short)(2);
			else if (samptype == 3)
				this.samptype = (short)(3);
			else
				this.samptype = (short)(0);
			fold = (short)(1);
			}
		else if (samptype == 0)
				this.samptype = (short)(0);
		/* Now count records like the first requested */
		countrec = 1;
		lastrec = rgthd.rec;
		/*temphd = rgthd;*/
		for (i=0; i<headers.size(); i++)
			{
			temphd = (RgHead)(headers.elementAt(i));
			if (temphd.rec != lastrec)
				{
				lastrec = temphd.rec;
				countrec++;
				}
			}
		if (headers.size()/countrec != (this.ntrace+this.naux))
			{
			if (getall)
				{
				countrec = 1;
				fhd.ntrace = headers.size();
				}
			else
				{
				System.out.println("Found " + countrec
				+ " contiguous field records within " + headers.size()
				+ " traces similar to the first requested\n"
				+ " but expected " + (headers.size()/countrec)
				+ " records of " + (this.ntrace+this.naux)
				+ " traces each.");
				countrec = headers.size()/(this.ntrace+this.naux);
				String advice="This SEG-Y file does not appear to be sorted by\n"
				+ "field record number or FFID. Extracting " + countrec
				+ "\nrecords of " + (this.ntrace+this.naux)
				+ " traces each, but traces from different FFIDs will\n"
				+ "be mixed together and may require sorting for coherent\n"
				+ "display, or to apply geometry information. Requesting\n"
				+ "one field record at a time will extract just the traces\n"
				+ "from that record. Picking of records not sorted by FFID\n"
				+ "may be possible if geometries are correct.";
				if (win == null)
					System.out.println(advice);
				else
					win.Errout(advice);
				fhd.ntrace = (this.ntrace + this.naux)*countrec;
				}
			}
		else
			fhd.ntrace = (this.ntrace + this.naux)*countrec;
		return countrec;
		}
	public void regularize()
		{
		if (units <= 0)
			units = 1;
		}
	public void readSEGYFileHead(String name, boolean swabhead)
		{
		FileInputStream fis = null;
		BufferedInputStream bis = null;
		DataInputStream dis = null;
		System.out.println("readSEGYFileHead: reading reel header from file " + name + " ... ");
		try {fis = new FileInputStream(name); }
		catch(FileNotFoundException fnf)
			{
			System.out.println(fnf.toString());
			return;
			}
		bis = new BufferedInputStream(fis);
		dis = new DataInputStream(bis);
		try {read(dis); }
		catch(IOException ioe)
			{
			System.out.println("Could not read SEG-Y Reel header from file " + name);
			return;
			}
		try {fis.close(); }
		catch(IOException ioe)
			{
			System.out.println(ioe.toString());
			return;
			}
		regularize();
		System.out.println("readSEGYFileHead done.");
		}
	public void read(DataInputStream dis) throws IOException
		{
		int i;
		for (i=0; i<SEGY_CARD_SIZE; i++)
			cards[i] = dis.readByte();
		jobid = dis.readInt();
		lineno = dis.readInt();
		reelno = dis.readInt();
		ntrace = dis.readShort();
		naux = dis.readShort();
		mudt = dis.readShort();
		omudt = dis.readShort();
		nt = dis.readShort();
		ont = dis.readShort();
		samptype = dis.readShort();
		fold = dis.readShort();
		sort = dis.readShort();
		sums = dis.readShort();
		sweep0 = dis.readShort();
		sweepf = dis.readShort();
		sweepl = dis.readShort();
		sweept = dis.readShort();
		sweeptr = dis.readShort();
		taper0 = dis.readShort();
		taperf = dis.readShort();
		tapert = dis.readShort();
		cor = dis.readShort();
		gain = dis.readShort();
		recov = dis.readShort();
		units = dis.readShort();
		pol = dis.readShort();
		vibpol = dis.readShort();
		for (i=0; i<336; i++)
			unas[i] = dis.readByte();
		num_samps = dis.readInt();
		}
	public void write(DataOutputStream dos) throws IOException
		{
		int i;
		for (i=0; i<SEGY_CARD_SIZE; i++)
			 dos.writeByte(cards[i]);
		dos.writeInt(jobid);
		dos.writeInt(lineno);
		dos.writeInt(reelno);
		dos.writeShort(ntrace);
		dos.writeShort(naux);
		dos.writeShort(mudt);
		dos.writeShort(omudt);
		dos.writeShort(nt);
		dos.writeShort(ont);
		dos.writeShort(samptype);
		dos.writeShort(fold);
		dos.writeShort(sort);
		dos.writeShort(sums);
		dos.writeShort(sweep0);
		dos.writeShort(sweepf);
		dos.writeShort(sweepl);
		dos.writeShort(sweept);
		dos.writeShort(sweeptr);
		dos.writeShort(taper0);
		dos.writeShort(taperf);
		dos.writeShort(tapert);
		dos.writeShort(cor);
		dos.writeShort(gain);
		dos.writeShort(recov);
		dos.writeShort(units);
		dos.writeShort(pol);
		dos.writeShort(vibpol);
		for (i=0; i<336; i++)
			 dos.writeByte(unas[i]);
		dos.writeInt(num_samps);
		}

	public FltVol importSEGYFile(String path, boolean tapehead, boolean swabhead,
		int firstrec, int lastrec, int samptype)
		{
		int iv0, iv1, iv2, iv3, temp;
		SEGYFileHead sfhd = this;
		if (firstrec > -1 && lastrec > -1 && firstrec > lastrec)
			{
			temp = firstrec;
			firstrec = lastrec;
			lastrec = temp;
			}
		iv3 = this.scanSEGYFile(path, tapehead, swabhead, firstrec, lastrec, samptype);
		if (sfhd.isValid() && iv3 > 0)
			{
			iv1 = sfhd.getnt();
			iv2 = sfhd.ntrace;
			}
		else if (firstrec < 0 && lastrec < 0)
			{
			System.out.println("File " + path
			+ "\ndoes not appear to be a SEG-Y file:\n"
			+ sfhd);
			return null;
			}	
		else if (firstrec < 0 && lastrec > -1)
			{
			System.out.println("File " + path
			+ "\ndoes not appear to be a SEG-Y file,\n"
			+ "or does not include record numbers <= " + lastrec + ":\n"
			+ sfhd);
			return null;
			}	
		else if (firstrec > -1 && lastrec < 0)
			{
			System.out.println("File " + path
			+ "\ndoes not appear to be a SEG-Y file,\n"
			+ "or does not include record numbers >= " + firstrec + ":\n"
			+ sfhd);
			return null;
			}	
		else
			{
			System.out.println("File " + path
			+ "\ndoes not appear to be a SEG-Y file,\n"
			+ "or does not include record numbers " + firstrec
			+ " through " + lastrec + ":\n"
			+ sfhd);
			return null;
			}
		long addrs[];
		addrs = new long[sfhd.addresses.size()];
		Long Ltemp;
		for (int i=0; i<addrs.length; i++)
			{
			Ltemp = (Long)(sfhd.addresses.elementAt(i));
			addrs[i] = Ltemp.longValue();
			}
		RgHead rghd[] = new RgHead[sfhd.headers.size()];
		for (int i=0; i<rghd.length; i++)
			rghd[i] = new RgHead((RgHead)(sfhd.headers.elementAt(i)));
		FltVol fvol = new FltVol(iv3, (int)(sfhd.ntrace), sfhd.getnt());
		fvol.headers = rghd;
		fvol.geogr = sfhd.geogr;
		fvol.setUnit0ofDim(0, sfhd.fhd.starttime);
		fvol.setDUnitofDim(0, sfhd.fhd.dt);
		fvol.setUnitsofDim(0, "Time, sec");
		fvol.setUnit0ofDim(1, 0);
		fvol.setUnitsofDim(1, "Trace Sequence");
		fvol.setDUnitofDim(1, 1);
		fvol.setUnit0ofDim(2, 0);
		fvol.setUnitsofDim(2, "Record Index");
		fvol.setDUnitofDim(2, 1);
		fvol.readSEGYFile(path, addrs, samptype);

		return fvol;
		}
		
	}