#include <fstream>
#include <stdio.h>

#include <netinet/in.h>
/****************************************************************************/
/*                                                                          */
/* This program is used to read segb data files produced by David Abensours */
/* data capture program.                                                    */
/*                                                                          */
/* Plans for this program.                                                  */
/* Support the following output format:                                     */
/*                                      SEGY disk files                     */
/*                                      SEGY disk image                     */
/*                                      SEGY tape write                     */
/*                                      SU (seismic unix) data files        */
/*                                      ASCII columns for octave, gnuplot,  */
/*                                            matlab, GMT etc.              */
/*                                                                          */
/*    COPYRIGHT 1999 Řistein Aanensen                                       */
/*                                                                          */
/****************************************************************************/

unsigned int power_of (unsigned short gr_tall, short potens);
int int_to_4bcd(int input);


struct segy {	/* segy - trace identification header */

	int tracl;	/* trace sequence number within line */

	int tracr;	/* trace sequence number within reel */

	int fldr;	/* field record number */

	int tracf;	/* trace number within field record */

	int ep;	/* energy source point number */

	int cdp;	/* CDP ensemble number */

	int cdpt;	/* trace number within CDP ensemble */

	short trid;	
	short nvs;	
	short nhs;	
	short duse;	
	int offset;	/* distance from source point to receiver
			   group (negative if opposite to direction
			   in which the line was shot) */

	int gelev;	/* receiver group elevation from sea level
			   (above sea level is positive) */

	int selev;	/* source elevation from sea level
			   (above sea level is positive) */

	int sdepth;	/* source depth (positive) */

	int gdel;	/* datum elevation at receiver group */

	int sdel;	/* datum elevation at source */

	int swdep;	/* water depth at source */

	int gwdep;	/* water depth at receiver group */

	short scalel;	
	short scalco;	
	int  sx;	/* X source coordinate */

	int  sy;	/* Y source coordinate */

	int  gx;	/* X group coordinate */

	int  gy;	/* Y group coordinate */

	short counit;	
	short wevel;	/* weathering velocity */

	short swevel;	/* subweathering velocity */

	short sut;	/* uphole time at source */

	short gut;	/* uphole time at receiver group */

	short sstat;	/* source static correction */

	short gstat;	/* group static correction */

	short tstat;	/* total static applied */

	short laga;	
	short lagb;	
	short delrt;	
	short muts;	/* mute time--start */

	short mute;	/* mute time--end */

	unsigned short ns;	/* number of samples in this trace */

	unsigned short dt;	/* sample interval; in micro-seconds */

	short gain;	
	short igc;	/* instrument gain constant */

	short igi;	/* instrument early or initial gain */

	short corr;	/* correlated:
				1 = no
				2 = yes */

	short sfs;	/* sweep frequency at start */

	short sfe;	/* sweep frequency at end */

	short slen;	/* sweep length in ms */

	short styp;	
	short stas;	
	short stae;	
	short tatyp;	
	short afilf;	
	short afils;	
	short nofilf;	
	short nofils;	
	short lcf;	
	short hcf;	
	short lcs;	
	short hcs;	
	short year;	
	short day;	
	short hour;	
	short minute;	
	short sec;	
	short timbas;	
	short trwf;	
	short grnors;	
	short grnofr;	
	short grnlof;	
	short gaps;	
	short otrav;	
	float d1;	
	float f1;	
	float d2;	
	float f2;	
	float ungpow;	
	float unscale;	
	int ntr; 	
	short mark;	
        short shortpad; 
	short unass[14];	
	float  data[6144];

};

int main(int args, char* argv[]) {
  int input;
  int test;
  FILE *fil;
  float data;
  unsigned short byte_pr_sample, nch, rec_length;
  segy tr;
  unsigned short gain[120], ifp[4];
 
  if (fread(&input,2,1,stdin))
    tr.fldr=int_to_4bcd(input);
  if (fread(&input,8,1,stdin));
  
  if (fread(&input,2,1,stdin));
  unsigned short mask = 0xF0;
  test = (input & mask) >> 4;
  byte_pr_sample = test*100;
  mask = 0xF;
  test = (input & mask);
  byte_pr_sample += test*10;
  mask = 0xF000;
  test = (input & mask) >> 12;
  byte_pr_sample += test;             /* Final byte per sample */
  nch = (((byte_pr_sample/2)-7)*4)/5; /* Number of channels */
  mask = 0xF00;
  tr.dt = (unsigned short) (input & mask) >> 8; /* Sample intervall */
   
   if (fread(&input,4,1,stdin));
   if (fread(&input,1,1,stdin));
   mask = 0xF0;
   test = (input & mask) >> 4;
   rec_length = test*10;
   mask = 0xF;
   test = (input & mask);
   rec_length += test;                /* Record length */
   tr.ns = (unsigned short) (rec_length*1024)/tr.dt;         /* Number of samples */
   tr.dt=tr.dt*1000;
   float outp[nch][tr.ns];
   
   if (fread(&input,7,1,stdin));
   for (int i=0; i < nch; i++) {    /* Loop through n/2 to get constant gain */
     if (fread(&input,2,1,stdin));
     mask = 0xF;
     test = (input & mask) >> 1;
     gain[i]=power_of(2,test);
   }
   
   if (fread(&input,12,1,stdin));
   int iter=1;
   for (int k=0; k < tr.ns; k++) {
     if (fread(&input,14,1,stdin));
     for (int i=0; i < nch/4; i++) {    /* Loop through traces, get gain codes, and trace values */
       if (fread(&input,2,1,stdin));
       mask = 0xF0;
       test = (input & mask) >> 5;
       ifp[0]=power_of(4,test);
       mask = 0xF;
       test = (input & mask) >> 1;
       ifp[1]=power_of(4,test);
       mask = 0xF000;
       test = (input & mask) >> 13;
       ifp[2]=power_of(4,test);
       mask = 0xF00;
       test = (input & mask) >> 9;
       ifp[3]=power_of(4,test);
       
       /* Four data channels */
       for (int j=0; j < 4; j++) {
	 iter ++;
	 if (fread(&input,2,1,stdin));
	 input = htons(input);
	 mask = 0x8000;
	 test = (input & mask) >> 15;
	 mask = 0x6FF7;
	 data = (input & mask) >> 1;
	 if (test == 1) {
	   outp[j+(i*4)][k]=-1*(data/2)/ifp[j]/gain[j*(i+1)];
	 } else {
	   outp[j+(i*4)][k]=(data/2)/ifp[j]/gain[j*(i+1)];
	 }
       }
     }
   }
   
   //fil=fopen("test.su","w");
   
   for (int i=0; i<nch; i++) {
     tr.data=outp[i];
     tr.tracf=i+1;
     fwrite(&tr,240+(4*tr.ns),1,stdout);
   }
   //fclose(fil);
   exit(0);
}

int int_to_4bcd (int input) {
  int out;
  short test;
  unsigned short mask;
  mask = 0xF0;
  test = (input & mask) >> 4;
  out=test*1000;
  mask = 0xF;
  test = (input & mask);
  out += test*100;
  mask = 0xF000;
  test = (input & mask) >> 12;
  out += test*10;
  mask = 0xF00;
  test = (input & mask) >> 8;
  out += test;
  return out;
}

unsigned int power_of (unsigned short gr_tall, short potens) {
  unsigned int out=gr_tall; 
  for (short i=1; i<potens; i++)
    out *= gr_tall;
  return out;
}