#include <fstream>
#include <stdio.h>
#include <netinet/in.h>
#include <math.h>

//***************************************************************************
//
//                      s e g c 2 s u . C
//
//
// This program is used to demultiplex SEG-C data files from 9-track tape
// made by DFS-V seismic recording equipment.
//
// Compilation :   g++ -lm -o segc2su segc2su.C
//
// The program reads from stdin and outputs to stdout, and can thus be invoked
// like this:      segc2su < file.segc > file.su
// Input file header parameters are written to STDERR (terminal window)
//
// Output is in Seismic Unix format (modified SEG-Y, omitting reel headers).
// Floats are in native format on the machine where the program is executed.
// Refer to this web page for format description:
// http://www2.geo.uib.no/seismic-unix/
//
// Based on Oistein Aanensen's segb2su program (1999).
//
// University of Bergen, Dept. of Earth Science, Norway
//
// ------------------------------------------------------------------------
// Rev.   Date          By                 Description
// ------------------------------------------------------------------------
// 0.1    Dec 18, 2003  O.Meyer & K.Hiim   First version
//
//***************************************************************************

static int ibm_to_float(int *from, int endian);
unsigned int power_of (unsigned short gr_tall, short potens);
unsigned short bcd_to_ushort(unsigned short *input);
unsigned short bcd_to_char(char input);
int get_nibble(char input, int high);
int is_leapyear(int year);
int day_of_year(int date, int month, int year);
float get_float(char *i);


//   C o n s t a n t   d e c l a r a t i o n

const int FORCE_REC_TIME			= 10;
const short YEAR				= 1982;		// FRAM-IV data
const int HIGH_NIBBLE				= 1;
const int LOW_NIBBLE				= 0;
const int SEG_Y_TRACE_HEADER_SIZE		= 240;		// SEG-Y trace header always this size
const int SEG_Y_DATA_SIZE			= 1024*1023;	// SEG-Y trace data, must be larger or equal to actual size
const int SEG_C_HEADER_SIZE			= 272;		// Size of SEG-C header varies. The first 24 bytes are
												// mandatory, then follows an optional part where the
												// size depends on the number of channels per scan;
												// each channels need four bytes - including any AUX
												// channels in the DFS-V. E.g. 24 seis channels + 6 AUX
												// channels yields a header size of 24 + [(24+6)*4] =
												// 144 bytes

								
//	T y p e   d e c l a r a t i o n 

// SEG-Y trace identification header + data type declaration

struct segy
{
	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];				// unassigned
	int data[SEG_Y_DATA_SIZE];	// data portion
};


//====================================================
//		M A I N   P R O G R A M
//====================================================


int main(int args, char* argv[])
{

	char input[SEG_C_HEADER_SIZE];
	int test;
	FILE *fil;
	float data;
	unsigned short file_nr, byte_pr_scan, sample_interval, record_length, texas, no_of_channels, rec_length, no_of_samples;
	short month, day_of_month, Julian_day;
	segy tr;
 

	// Read complete SEG-C header into memory
	
	test = fread(input, sizeof(input), 1, stdin);
	if (test<1)
	{
		fprintf(stderr, "Error reading SEG-C header.\n");
		exit(1);
	}

	
	// Get INSTRUMENT MANUFACTURER ID, TEXAS INSTRUMENTS = 15 (BCD)
		
	texas = bcd_to_char(input[12]);
	//fprintf(stderr, "Texas Instrument ID: %d\n", texas);
	
	
	// ERROR CORRECTION
	//
	// Some systematic errors in the SEG-C header were detected in FRAM-IV data from 1982 (reel 8).
	// The instrument manufacturer constant parameter in header byte 13 should be 15 (BCD);
	// in some records this parameter was 37, i.e. 22 (BCD) was added. Some neighbouring header
	// bytes were also affected by this. Don't know the reason. Code lines below will clean up the mess.
	// Don't forget, we count header bytes from zero.
	
	if (texas != 15)
	{
		for (int i=0x4; i<0x11; i++)
		{
			if ((i!=10) && (i!=11))			// Byte 10 & 11 were OK, so don't touch them ...
				input[i] -= 0x22;
		}
	}
	
	
	// Get FILE NUMBER
	
	file_nr = bcd_to_ushort((unsigned short *)&input[0]);	// Remember - index from zero!
	fprintf(stderr, "File no.: %d\n", file_nr);
	tr.fldr = file_nr;										// TRACE HEADER: fldr = field record number
	
	
	// GENERAL CONSTANTS : DFS-V FRONT PANEL SWITCHES
	//
	// In the FRAM-IV project the General Constants (header bytes 5 to 10) were used
	// for date and time information. The encoding was (all BCD):
	//
	//		Byte #
	//		  5  = Month
	//		  6	 = Day of month
	//		  7	 = Hour   (GMT)
	//		  8  = Minute (GMT)
	//
	// These parameters are transfererred to Seismic Unix trace header fields:
	//
	//		year (constant declared at beginning of program)
	//		day  (day of year, "Julian day")
	//		hour
	//		minute
	
	tr.year			= YEAR;
	month			= get_nibble(input[4], HIGH_NIBBLE) * 10 + get_nibble(input[4], LOW_NIBBLE);
	day_of_month	= get_nibble(input[5], HIGH_NIBBLE) * 10 + get_nibble(input[5], LOW_NIBBLE);
	tr.day			= day_of_year(day_of_month, month, YEAR);	
	tr.hour			= get_nibble(input[6], HIGH_NIBBLE) * 10 + get_nibble(input[6], LOW_NIBBLE);
	tr.minute		= get_nibble(input[7], HIGH_NIBBLE) * 10 + get_nibble(input[7], LOW_NIBBLE);

	
	// Get BYTES/SCAN, calculate number of channels. A scan includes SYNC GROUP & TIME COUNTER (8 bytes)
	
	byte_pr_scan	= get_nibble(input[10], HIGH_NIBBLE) * 100 + get_nibble(input[10], LOW_NIBBLE) * 10 + get_nibble(input[11], HIGH_NIBBLE);
	fprintf(stderr, "     Bytes per scan: %d\n", byte_pr_scan);
	no_of_channels	= (byte_pr_scan-8)/4;	// 8=size of SYNC GROUP / TIME COUNTER, 4=bytes pr. sample
											// NOTE: Will also include AUX channels
	fprintf(stderr, "     Number of channels: %d\n", no_of_channels);
	
	// Get SAMPLE INTERVAL
	
	sample_interval = get_nibble(input[11], LOW_NIBBLE);
	switch (sample_interval)				// TRACE HEADER: dt = sampling interval in microseconds
	{										// See DFS-V SEG-C format description
		case  1:	tr.dt = 1000; break;
		case  2:	tr.dt = 2000; break;
		case  4:	tr.dt = 4000; break;
		case  8:	tr.dt = 8000; break;
		case  9:	tr.dt =  500; break;
		case 10:	tr.dt = 1000; break;
		case 12:	tr.dt = 2000; break;
		default:	break;
	}
	fprintf(stderr, "     Sample interval: %d microseconds\n", tr.dt);

	
	// Get RECORD LENGTH, calculate number of samples per channel
	
	record_length = bcd_to_char(input[16]);
	//record_length = FORCE_REC_TIME;							// NOTE   FORCE RECORD LENGTH FOR THIS PROJECT
	fprintf(stderr, "     Record length: %d (multiples of 1.024 s)\n", record_length);	// Record length in multiples of 1.024 seconds
	no_of_samples = (record_length*1024)/sample_interval;							// Record length in multiples of 1024 ms
	fprintf(stderr, "     No. of samples: %d\n", no_of_samples);
	tr.ns = no_of_samples;															// TRACE HEADER: ns = number of samples
	
	
	// Get FIXED and EARLY GAIN values from last part of SEG-C header. Put values into two arrays.
	// Also determine number of SEIS channels (future option - in order to get rid of AUX channels)
	
	short fixed_gain[no_of_channels], early_gain[no_of_channels];
	int no_of_seis_ch = 0;
		
	for (int i=0; i < no_of_channels; i++)  
	{
		fixed_gain[i] = get_nibble(input[24+(i*4)], LOW_NIBBLE);
		if (get_nibble(input[24+(i*4)], HIGH_NIBBLE) == 2) 		// Channel ID code, inverse nibble reading
			no_of_seis_ch++;
		early_gain[i] = get_nibble(input[25+(i*4)], LOW_NIBBLE);
		//fprintf(stderr, "FG,EG ch %2d: %d, %d\n", i+1, fixed_gain[i], early_gain[i]);
	}
	fprintf(stderr, "     No. of seis channels: %d\n", no_of_seis_ch);	

	
	// Read data portion of file into array of chars
	
	char raw_data[byte_pr_scan*no_of_samples];					// Declare char array of proper size
	test = fread(raw_data, sizeof(raw_data), 1, stdin);
	if (test<1)
	{
	fprintf(stderr, "     Error reading SEG-C raw data.\n");
	  exit(1);
	}

	// D E M U X I N G
	//
	// First declare 2-dimensional array used in demux process:
	// View no_of_channels as lines, and no_of_samples as columns
	// Then fill in by columns and read by lines later on, when making Seismic Unix file
			
	int demux[no_of_channels][no_of_samples];		// 2-dim array used in demultiplexing
	int k, s, u, v;
	
	for (k=0; k < no_of_samples; k++)				// Loop through all scans
	{
		for (s=0; s < no_of_channels; s++)			// Loop through all channels in a scan set
		{
			//demux[s][k] = get_float(&raw_data[(k*byte_pr_scan)+8+(s*4)]);	// Now points to SEIS CH 1 data in scan k+1
			demux[s][k] = ibm_to_float((int*)&raw_data[(k*byte_pr_scan)+8+(s*4)],0);  // Last parameter = endian 0=little
		}
	}
	fprintf(stderr, "     Scan cnt, ch cnt: %d,%d\n",k, s);
	
	
	
	//ibm_to_float(int from, int endian)
	
	
	
	// W r i t e   S e i s m i c   U n i x   f i l e
	
	// Demux matrix is now filled with values. Extract data for each channel in sequence.
	// Add trace header for each channel (Seismic Unix format, modified SEG-Y)
	
	for (u=0; u < no_of_channels; u++)
	{
		tr.tracf=u+1;								// tracf = Trace number within file
		for (v=0; v < no_of_samples; v++)
		{
			tr.data[v] = demux[u][v];
		}
		tr.igc = fixed_gain[u];			// Trace header igc = INSTRUMENT GAIN CONSTANT, get value for this trace (=channel)
		tr.igi = early_gain[u];			// Trace geader igi = INSTRUMENT INITIAL GAIN, get value for this trace (=channel)
		fwrite(&tr,SEG_Y_TRACE_HEADER_SIZE+(4*tr.ns),1,stdout);
	}
	fprintf(stderr, "     SU write ch cnt: %d\n", u);
	
	//printf("Float: %f\n", get_float(&raw_data[8]));
	

	exit(0);
}

//====================================================
//		E N D   O F   M A I N   P R O G R A M
//====================================================


//*******************************************************
// SEG-C header parameter extraction routines.
//
// Parameters often occupy high or low nibbles of bytes,
// here is code to handle this ...
//*******************************************************

int get_nibble(char input, int high)
{
	unsigned short mask = (high ? 0xF0 : 0xF);
	unsigned short out = input & mask;
	return (high ? out >> 4 : out);
}
 
unsigned short bcd_to_ushort (unsigned short *input)
{
	int out;
	unsigned short test;
	unsigned short i = *input;
	unsigned short mask;
	
	mask = 0xF0;
	test = (i & mask) >> 4;
	out=test*1000;
	mask = 0xF;
	test = (i & mask);
	out += test*100;
	mask = 0xF000;
	test = (i & mask) >> 12;
	out += test*10;
	mask = 0xF00;
	test = (i & mask) >> 8;
	out += test;
	return out;
}

unsigned short bcd_to_char(char input) {
  return (get_nibble(input, 0) + get_nibble(input, 1) * 10);
}

/*
int bcd_to_int (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;
}*/


//****************************************************************
//   I B M   F L O A T I N G   P O I N T   C O N V E R S I O N
//****************************************************************

static int ibm_to_float(int * from, int endian)
/***********************************************************************
ibm_to_float - convert between 32 bit IBM and IEEE floating numbers
************************************************************************
Input::
from            input vector
to              output vector, can be same as input vector
endian          byte order =0 little endian (DEC, PC's)
                            =1 other systems
*************************************************************************
Notes:
Up to 3 bits lost on IEEE -> IBM
                                                                                
Assumes sizeof(int) == 4
                                                                                
IBM -> IEEE may overflow or underflow, taken care of by
substituting large number or zero
                                                                                
Only integer shifting and masking are used.
*************************************************************************
Credits: CWP: Brian Sumner,  c.1985
Modified by O.Meyer, Dec. 18, 2003
*************************************************************************/
{
register int fconv, fmant, i, t;
                                                                             
	fconv = *from;
                                                                                
	/* if little endian, i.e. endian=0 do this */
	if (endian==0) 
		fconv = (fconv<<24) | ((fconv>>24)&0xff) | ((fconv&0xff00)<<8) | ((fconv&0xff0000)>>8);
                                                                                
	if (fconv)
	{
		fmant = 0x00ffffff & fconv;
		/* The next two lines were added by Toralf Foerster */
		/* to trap non-IBM format data i.e. conv=0 data  */
		if (fmant == 0)
		fprintf(stderr, "Data are not in IBM FLOAT format!\n");
		t = (int) ((0x7f000000 & fconv) >> 22) - 130;
		while (!(fmant & 0x00800000)) { --t; fmant <<= 1; }
		if (t > 254) fconv = (0x80000000 & fconv) | 0x7f7fffff;
			else if (t <= 0) fconv = 0;
			else fconv = (0x80000000 & fconv) |(t << 23)|(0x007fffff & fmant);
	}            
    return (fconv);
}



/****************************************************************************
        I B M   F L O A T I N G   P O I N T   C O N V E R S I O N
		
 SEG-C data are stored in IBM floating point format. This routine will convert
 an IBM floating point number to a float with format dictated by the computer
 that the program is running on ...
 
 By K. Hiim.

 IBM floating point format description, quoting from DFS-V SEG-C format
 description:
 
    --------------
       Byte no.
    1   2   3   4
    --------------
    Qs  Qf  Qf  Qf
    Qc  Qf  Qf  0
    Qc  Qf  Qf  0
    Qc  Qf  Qf  0
    Qc  Qf  Qf  0
    Qc  Qf  Qf  0
    Qc  Qf  Qf  0
    Qc  Qf  Qf  0
 
 Qs : Signifies signal polarity
 Qc : Represents a power of 16 in binary excess 64 notation:
        1000000 =  0
	    1000001 = +1
	    0111111 = -1
 Qf : Represents a fractional number with the binary point to the
      left of the most significant bit Qf is always true regardless of
	  polarity. Qf is normalized so that at least one of the four most
	  significant bits of Qf is 1 unless Qf is zero. Signal amplitude =
	  Qf x 16exp(Qc-64) millivolts. Only significant bit positions of
	  Qf for direct (uncompensated) recording are indicated. The least
	  significant bit (bit 7 of the 4th byte) is always a zero. Fixed
	  gain (gain constant) is added to the variable gain when determining
	  Qc and Qf.
****************************************************************************/

float get_float(char *i) {
  unsigned int value = i[0];
  unsigned int tmp;
  int v=10;
//  printf("%d\n", value);
    
  int Qs = ((value & 0x80) ? -1 : 1);
  int Qc = ((value & 0x40) ? -1 : 1) * (value & 0x3f);
  
  tmp = i[1] & 0xff;
//  printf("tmp: %u - value: ", tmp);
  tmp <<= 9;
  value = tmp;
//  printf("%d\n", value);
  tmp = i[2] & 0xff;
//  printf("tmp: %u - value: ", tmp);
  tmp <<= 1;
  value += tmp;
//  printf("%d\n", value);
  tmp = i[3] & 0xff;
//  printf("tmp: %u - value: ", tmp);
  tmp >>= 7;
  value += tmp;
//  printf("%u\n", value);
  
  int Qf = value;
//  printf("Qs: %u\nQc: %u\nQf: %u\n", Qs, Qc, Qf);
  v = 10;
//  printf("%d ", v);
  while ((v-Qf) <= 0) {
//    printf("%d", v);
	v = (v * 10);
  }
//  printf("\n%d\n", v);
  return (Qs * ((float)Qf / (float)v) * pow(16.0, Qc - 64));
}


/******************************************************************************
			J U L I A N   D A Y   C A L C U L A T I O N

 SEG-Y trace headers has parameter for Julian day. The Julian day is really
 day of the year (1 .. 365), not the original 'Julian date' that counts from
 January 1, 4713 BC ...
 Here are routines that calculates Julian data based on month, day of month
 and year.
******************************************************************************/

int is_leapyear(int year) {
  return ( (!(year%4)) && ( (year%100) || (!(year%400)) ) );
}


int day_of_year(int date, int month, int year) {
  short d[12] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  short tmp = 0;
  int a;
  if ( (date < 1) || (date > 31) || (month < 1) || (month > 12) )
    return 0;

  if (is_leapyear(year)) {
    d[1]++;
  }
  for (a=0; a < (month-1); a++) {
    tmp += d[a];
  }
  tmp += date;
  return tmp;
}


//==================================================================================
//				E N D   O F   P R O G R A M
//==================================================================================