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

/***************************************************************************************

                      s e g c 2 s u - v e r _ 0 _ 2 . 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-ver_0_2 segc2su-ver_0_2.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/


 University of Bergen, Dept. of Earth Science, Norway

 -------------------------------------------------------------------------------------------
 Rev.  Date          By                 Description
 -------------------------------------------------------------------------------------------
 0.2   Oct 11, 2006  O.M.            a) Printing info on all SEG-C header fields
                                     b) Fixed gain SEG-Y trace header parameter now in dB.
                                     c) Transferring info to these SEG-Y trace headers:
                                         * gain Field instrument gain type code)
                                         * duse Data use (1 = production, 2 = test)
                                     b) Treating DFS-V front panel switches different. 
                                     b) Computing size of optional header part (instead of using constant)
                                        by reading number of bytes per scan, and then dividing by 4
                                        (number of bytes per channel). The optional SEG-C header size
                                        will be (no_of_channels * 4).
                                        
 0.1   Dec 18, 2003  O.M., 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 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

								
//	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;			// field instrument gain type code: 3 = floating point
	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[512];					// DFS-V had max 60 channels + 6 AUX's = 66, so maximum header size
										// should be 24 + (66 * 4) = 288
	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;
	unsigned short serial_no, gain_mode, record_type, gain_constant, gain_word, gain_constant_dB, gain_word_dB, bit_value;
	short year, month, day_of_month, Julian_day;
	segy tr;
 
	//---------------------------------------------------------------------------------------------
	//--------   F I X E D - L E N G T H   ( 2 4   B Y T E S )  H E A D E R   B L O C K 
	//---------------------------------------------------------------------------------------------
	
	test = fread(input, 24, 1, stdin);	
	if (test < 1)
	{
		fprintf(stderr, "Error reading fixed-part (24 bytes) of SEG-C header - exit.\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. 6 bytes, from no. 5 to 10 
	//
	// 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
	//
	// N O T E :
	// Rev 0.2, new meaning of DFS-V front panel switches. All are BCD encoded.
	//		Byte #
	//		  5  = Day of month
	//		  6	 = Month
	//		  7	 = Year, two most significant digits
	//		  8  = Year, two least sign. digits
	// When cross-checking recording logs, we see that the observers sometimes forgot to set the DATE switch.
	// Regardless, we'll transfer this information to the SU file - better then nothing ..?
	
	/* From Ver 0.1
	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);
	*/
	// Remember, we index from zero ...
	day_of_month	= (get_nibble(input[4], HIGH_NIBBLE) * 10)   + (get_nibble(input[4], LOW_NIBBLE));
	month			= (get_nibble(input[5], HIGH_NIBBLE) * 10)   + (get_nibble(input[5], LOW_NIBBLE));
	year			= (get_nibble(input[6], HIGH_NIBBLE) * 1000) + (get_nibble(input[6], LOW_NIBBLE) * 100) + (get_nibble(input[7], HIGH_NIBBLE) * 10) + (get_nibble(input[7], LOW_NIBBLE) * 1);
	tr.day			= day_of_year(day_of_month, month, year);
	tr.year			= year;

	
	// ----------  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 us\n", tr.dt);


	// ------------  Get INSTRUMENT SERIAL NUMBER (for info only, no place in SEG-Y trace headers (in reel header, yes - but we're writing SU data now)
	// 6 BCD digits
	serial_no =  (get_nibble(input[13], HIGH_NIBBLE) * 100000) + (get_nibble(input[13], LOW_NIBBLE) * 10000);
	serial_no += (get_nibble(input[14], HIGH_NIBBLE) * 1000) + (get_nibble(input[14], LOW_NIBBLE) * 100);
	serial_no += (get_nibble(input[15], HIGH_NIBBLE) * 10) + (get_nibble(input[15], LOW_NIBBLE) * 1);
	fprintf(stderr, "     Instrument serial no.: %d\n", serial_no);

	
	// ------------  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
	

	// ------------  GAIN MODE, it should always be 9 to designate floating point gain control system
	gain_mode = get_nibble(input[17], HIGH_NIBBLE);
	if (gain_mode != 9)
	{
		fprintf(stderr, "Error GAIN MODE parameter in SEG-C header - exit.\n");
		exit(1);
	}
	fprintf(stderr, "     Gain mode (9): ok\n");
	tr.gain = 3;		// "3" = floating point

	
	// ------------  RECORD TYPE, Normal or test
	// In SEG-Y trace header, "duse" = Data use (1 = production, 2 = test)

	record_type = get_nibble(input[17], LOW_NIBBLE);
	switch (record_type)
	{	
		case 8:		fprintf(stderr, "     Record type (8): Normal shot\n"); tr.duse = 1; break;
		case 2:		fprintf(stderr, "     Record type (2): Test\n"); tr.duse = 2; break;
		default:	fprintf(stderr, "     Record type: %d - warning, illegal\n", record_type); break;
	}


	// -----------  Get Low-cut filter data

	// First get LC frequency
	// Two BCD digits, 00 = out, 03, 05, 08, 12, 18 or 27; 03 is actually 3.56Hz, 05 is actually 5.33Hz
	tr.lcf = (get_nibble(input[18], HIGH_NIBBLE) * 10)   + (get_nibble(input[18], LOW_NIBBLE));
	// Now get LC filter slope, in multiples of 6 dB/octave
	tr.lcs = 6 * get_nibble(input[19], HIGH_NIBBLE);
	fprintf(stderr, "     Low-cut filter frequency: %d Hz\n", tr.lcf);
	fprintf(stderr, "     Low-cut filter slope: %d dB/oct\n", tr.lcs);


	// ------------  Get Notch filter data
	// Two BCD digits, 00=out, 50 or 60 (for 50 or 60Hz)

	tr.nofilf = (get_nibble(input[22], HIGH_NIBBLE) * 10)   + (get_nibble(input[22], LOW_NIBBLE));
	if (tr.nofilf == 0)
		fprintf(stderr, "     Notch filter frequency (0): Out\n");
	else
		fprintf(stderr, "     Notch filter frequency: %d Hz\n", tr.nofilf);


	// ------------  Get Alias filter data
	// One BCD digit, 1=256Hz, 2=128Hz, 4=64Hz, 8=32Hz

	switch (get_nibble(input[23], HIGH_NIBBLE))
	{
		case 1:		tr.afilf = 256; break;
		case 2:		tr.afilf = 128; break;
		case 4:		tr.afilf =  64; break;
		case 8:		tr.afilf =  32; break;
		default:	break;
	}
	fprintf(stderr, "     Alias filter frequency: %d Hz\n", tr.afilf);


	//---------------------------------------------------------------------------------------------
	// --------   O P T I O N A L   H E A D E R   B L O C K 
	//---------------------------------------------------------------------------------------------

	char opt_header[no_of_channels * 4];	// Declare char array of proper size, 4 bytes per ch.
	test = fread(opt_header, sizeof(opt_header), 1, stdin);
	if (test < 1)
	{
		fprintf(stderr, "Error reading optional part of SEG-C header - exit.\n");
		exit(1);
	}

	
	/*----------  Get FIXED and EARLY GAIN values for each seis channel
	Both recorded as 4-bit binary words
	Interpretation:
		Gain constant:	Most sign. bit has a gain value of 2 exp 8 = 256 (in dB: 20 log 256 = 48 dB),
						least sign. bit has a gain value of 2 exp 1 = 2 (6 dB).
						The least sign. bit is recorded as a zero for the DFS-V.
		Gain word:		0000 is recorded when operating in floating point gain control.
						When operating in manual gain control, the most significant bit has a gain
						value of 2 exp 8 (48 dB) and the least significant bit has a gain value of 
						2 exp 1 (6 dB). 
						The least sign. bit is recorded as a zero for the DFS-V.
	*/

	// Put values into two arrays.
	// Also determine number of SEIS 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++)  
	{
		if (get_nibble(opt_header[i*4], HIGH_NIBBLE) == 2) 		// Channel ID code, inverse nibble reading
			no_of_seis_ch++;

		gain_constant = get_nibble(opt_header[i*4], LOW_NIBBLE);
		bit_value = (gain_constant >> 1) & 0x01;
		gain_constant_dB = 12 * bit_value;
		bit_value = (gain_constant >> 2) & 0x01;
		gain_constant_dB += 24 * bit_value;
		bit_value = (gain_constant >> 3) & 0x01;
		gain_constant_dB += 48 * bit_value;
		fixed_gain[i] = gain_constant_dB;

		gain_word = get_nibble(opt_header[1+(i*4)], LOW_NIBBLE);
		bit_value = (gain_word >> 1) & 0x01;
		gain_word_dB = 12 * bit_value;
		bit_value = (gain_word >> 2) & 0x01;
		gain_word_dB += 24 * bit_value;
		bit_value = (gain_word >> 3) & 0x01;
		gain_word_dB += 48 * bit_value;
		early_gain[i] = gain_word_dB;
	}
	
	// Print info for each seis channel
	for (int i=0; i < no_of_seis_ch; i++) 
	{
		if (early_gain[i] == 0)
			fprintf(stderr, "        Fixed, early gain seis ch %2d: %d dB, Auto (floating point gain control)\n", i+1, fixed_gain[i]);
		else
			fprintf(stderr, "        Fixed, early gain seis ch %2d: %d, %d dB\n", i+1, fixed_gain[i], early_gain[i]);
	}

	fprintf(stderr, "     No. of seis channels: %d\n", no_of_seis_ch);


	//---------------------------------------------------------------------------------------------
	// --------   D A T A   B L O C K 
	//---------------------------------------------------------------------------------------------
	
	// 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 data section - exit.\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 count, channel count: %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 traces written: %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");
			exit(1);	
		}
		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);
}



/*******************************************************************************************
 !!!!!!!!!!!   There are errors in this routine - don't use. OM, 11 Oct 2006      !!!!!!!!

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
//==================================================================================