/*
	seismic.c

	cl /AL seismic.c getfile

	Processes seg-y data in mainframe IBM float format

	R. Ambroziak USGS/OMG

*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "font.h"
#include "nomouse.h"

#define MSK_MANT_IBM 0xffffff
#define MSK_EXP_IBM  0x7f000000
#define MSK_SIGN_IBM 0x80000000
#define MSK_NORM_IEEE 0x800000

#define MAX_BLOCK 20000
#define NUM_SEG 38
#define WIDTH 640
#define HEIGHT 480
#define D_SCALE 5
#define W_SCALE 6
#define MAX_SHOTS 5000
#define NUMLINE 8
#define MENU_C 5
#define MENUXC 550

int Mouse=0;
float ZeroVal=0.0;
unsigned char InBuff[BUFSIZ],RmsBuff[BUFSIZ],ValBuff[BUFSIZ],HdrBuff[BUFSIZ];
int White=255;
int BlkSize;
float Gain=W_SCALE,Vf,Hf;
int T1,P1;
char *menu_c[MENU_C]= 
	{"intensity","hue","saturation",	"reset","exit color"}; 
int Map_Xmin,Map_Ymin;
float Map_Scale;
int Shot_X[MAX_SHOTS],Shot_Y[MAX_SHOTS];
int Num_Shots=0;
unsigned char S_Buff[MAX_BLOCK],Buff[12][640],ImBuff[20][640];
unsigned char Inbuff1[BUFSIZ],Inbuff2[BUFSIZ],Inbuff3[BUFSIZ],Inbuff4[BUFSIZ];
unsigned char P_Out[WIDTH][4],P_In[HEIGHT][8];
float Data[(MAX_BLOCK-240)/4];

int Map=0,Color=0;
int Cxc,Cyc,Cxs,Cys;
int Pxc,Pyc,Pxs,Pys;
float H1=0.0,H2=155.0,H3=0.0;
float S1=1.0,S2=1.0,S3=1.0;
float I1=0.5,I2=0.5,I3=1.0;
double decode_line();
FILE *Fpin,*FpRms,*FpVal,*FpPlot,*FpHdr,*FpOut;
int FileNum;
int LabelTrace=0;
unsigned char Pbuff[4][NUMLINE][800];
int CutOff[NUMLINE];
int PixBrite[12][WIDTH];
unsigned char ValA[WIDTH],ValB[WIDTH];

unsigned char getpt();
void seismic();

main(argc,argv)

int argc;
char *argv[];
{
	int i,j,k;
	int vscale=10;
	char ans;
	char filename[100],name[100];
	int numpts,val,blksize,plotpts,numfile;
	double rms,dnum=0.0,srms,max_rms,min_rms;
	double dval,ddat;
	long offskip,numseg=NUM_SEG,offset,line;

	if(argc<2)
	{
		numfile=get_file_name("*.dat",name,0);
		for(k=0;k<strlen(name);k++)
			if(name[k]=='.')
				name[k]='\0';
	}
	else
	{
		strcpy(name,argv[1]);
	}
	strcpy(filename,name);
	strcat(filename,".dat");
	Fpin=fopen(filename,"rb");
	if(!Fpin)
	{
		printf("Could not open '%s'.\n\n",filename);
	}
	setbuf(Fpin,InBuff);
	fread((char *)S_Buff,sizeof(char),3600,Fpin);

	if(argc<3)
	{
		printf("Give blocksize in bytes\n");
		scanf("%d",&blksize);
	}
	else
	{
		sscanf(argv[2],"%d",&blksize);
	}
	BlkSize=blksize;
	numpts=(blksize-240)/4;

	FileNum-=1;
	if(argc<3)
	{
		do
		{
			printf("   r -- read data and decode 1 data line\n");
			printf("   m -- compute rms for all lines and put in ____.rms\n");
			printf("\n   x -- exit\n\n");
			ans=getch();
			if(ans=='r')
			{
				printf("Give number of line to read.\n\n");
				scanf("%ld",&line);
				offset=(line-1)*blksize+3600;
				fseek(Fpin,offset,SEEK_SET);
				fread((char *)S_Buff,sizeof(char),240,Fpin);
				fread((char *)S_Buff,sizeof(char),blksize-240,Fpin);
				printf("blksize = %d  numpts = %d\n",blksize,numpts);
				rms=decode_line(numpts,S_Buff,Data);

				for(i=0;i<numpts;i++)
				{
					if(Data[i]<0.0)
						printf("%4d %20.5f\n",i+1,Data[i]);
					else
						printf("%4d                     %0.5f\n",i+1,Data[i]);
				}

				printf("rms = %lf\n",rms);
			}
			if(ans=='m')
			{
				strcpy(filename,name);
				strcat(filename,".rms");
				FpRms=fopen(filename,"wb");
				setbuf(FpRms,RmsBuff);
				strcpy(filename,name);
				strcat(filename,".val");
				FpVal=fopen(filename,"wb");
				setbuf(FpVal,ValBuff);
				strcpy(filename,name);
				strcat(filename,".hdr");
				FpHdr=fopen(filename,"wb");
				setbuf(FpHdr,HdrBuff);
				set_zero_val(numpts,blksize);
				rewind(Fpin);
				fread((char *)S_Buff,sizeof(char),3600,Fpin);
				fwrite((char *)S_Buff,sizeof(char),3600,FpHdr);
				max_rms=-1;
				min_rms=999999999999.0;
				srms=0.0;
				i=0;
				while((fread((char *)S_Buff,sizeof(char),240,Fpin))>0)
				{
					fwrite((char *)S_Buff,sizeof(char),240,FpHdr);
					fread((char *)S_Buff,sizeof(char),blksize-240,Fpin);
					rms=decode_line(numpts,S_Buff,Data);
					if(rms>0.0)
					{
						dnum+=1.0;
						srms+=rms;
					}
					i+=1;
					if(rms<min_rms)min_rms=rms;
					if(rms>max_rms)max_rms=rms;
					printf("%4d -- rms = %lf\n",i,rms);
					fwrite((char *)&rms,sizeof(double),1,FpRms);
					fwrite((char *)Data,sizeof(float),numpts,FpVal);
				}
				if(dnum>0.0)
					srms/=dnum;
				printf("mean rms = %lf\n",srms);
				printf("max = %lf   min = %lf\n",max_rms,min_rms);
				fclose(FpRms);
				fclose(FpVal);
				fclose(Fpin);
				fclose(FpHdr);
			}
		}while(ans!='x');
	}
	else
	{
		FpOut=fopen("seismic.out","at");
		fprintf(FpOut,"%s ",name);
		strcpy(filename,name);
		strcat(filename,".rms");
		FpRms=fopen(filename,"wb");
		setbuf(FpRms,RmsBuff);
		strcpy(filename,name);
		strcat(filename,".val");
		FpVal=fopen(filename,"wb");
		setbuf(FpVal,ValBuff);
		strcpy(filename,name);
		strcat(filename,".hdr");
		FpHdr=fopen(filename,"wb");
		setbuf(FpHdr,HdrBuff);
		set_zero_val(numpts,blksize);
		rewind(Fpin);
		fread((char *)S_Buff,sizeof(char),3600,Fpin);
		fwrite((char *)S_Buff,sizeof(char),3600,FpHdr);
		max_rms=-1;
		min_rms=999999999999.0;
		srms=0.0;
		i=0;
		while((fread((char *)S_Buff,sizeof(char),240,Fpin))>0)
		{
			fwrite((char *)S_Buff,sizeof(char),240,FpHdr);
			fread((char *)S_Buff,sizeof(char),blksize-240,Fpin);
			rms=decode_line(numpts,S_Buff,Data);
			i+=1;
			if(rms<min_rms)min_rms=rms;
			if(rms>max_rms)max_rms=rms;
			srms+=rms;
			printf("%4d -- rms = %lf\n",i,rms);
			fwrite((char *)&rms,sizeof(double),1,FpRms);
			fwrite((char *)Data,sizeof(float),numpts,FpVal);
		}
		dnum=i;
		srms/=dnum;
		printf("mean rms = %lf\n",srms);
		fprintf(FpOut,"number of traces = %d  ",i);
		fprintf(FpOut,"mean rms = %lf\n",srms);
		printf("max = %lf   min = %lf\n",max_rms,min_rms);
		fclose(FpRms);
		fclose(FpVal);
		fclose(Fpin);
		fclose(FpHdr);
		fclose(FpOut);
	}
	fclose(Fpin);
	fclose(FpHdr);
}

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

int set_zero_val(numpts,blksize)

int numpts,blksize;

{
	int i,j,k;
	int num=10;
	float rms=0.0,fnum=num;
	long offset=blksize;

	offset*=20;

	rewind(Fpin);
	fread((char *)S_Buff,sizeof(char),3600,Fpin);
	fseek(Fpin,offset,SEEK_CUR);	/* skip first 20 traces */
	for(i=0;i<num;i++)
	{
		fread((char *)S_Buff,sizeof(char),240,Fpin);
		fread((char *)S_Buff,sizeof(char),blksize-240,Fpin);
		rms+=decode_line(numpts,S_Buff,Data);
	}
	rms/=fnum;
	printf("\nrms bar = %f   ",rms);
	ZeroVal=rms/100000.0;
	printf("ZeroVal = %30.15f\n",ZeroVal);
}



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

double decode_line(numpts,buff,data)

int numpts;
float *buff;
float *data;

{
	int i,j,k;
	float fdata,fval;
	double rms=0.0,pts=numpts,ddata,used=0.0;
	double dval,pval,base=16.0;
	int val1,val2,val3,val4;
	double frac,sign,power;

	if(numpts<=0)return(-1.0);

	flip(buff,numpts);
	ibm_iee(buff,data,numpts);

	for(i=0;i<numpts;i++)
	{
		if(data[i]<-ZeroVal||data[i]>ZeroVal)
		{
			ddata=data[i];
			rms+=ddata*ddata;
			used+=1.0;
		}
	}
	if(used>0.0)
		rms/=used;
	if(rms>0.0)
		rms=sqrt(rms);
	else
		rms=0.0;
	printf("Used %5.0f  ",used);
	return(rms);
}


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

int ibm_iee(ibm_fpn,ieee_fpn,numpts)

long *ibm_fpn,*ieee_fpn;
int numpts;

{
	int i;
	long exponent;
	long mantissa;
	long sign;

	for(i=0;i<numpts;i++)
	{
		if ((ibm_fpn[i] & MSK_MANT_IBM) == 0)
	  	{
			ieee_fpn[i] = 0;
		}
		else
		{
			mantissa = (ibm_fpn[i] & MSK_MANT_IBM); /* get mantissa */
			exponent = ((((ibm_fpn[i] & MSK_EXP_IBM) >> 24) - 64) << 2 ); 
				/* get exponent */
			sign     = (ibm_fpn[i] & MSK_SIGN_IBM); /* get sign */
			while ((mantissa & MSK_NORM_IEEE) == 0) 
			{
				mantissa <<= 1;  /* normalize */
				exponent-- ;
			}
			mantissa = mantissa & 0x7fffff;    /* shift understood one out */
			exponent += 126;
			if (( exponent < 0 ) || (exponent > 255))
			{
				printf ("IBM floating point exponent out of range \n");
				sign=mantissa=exponent = 0;
			}
			exponent <<= 23;
			ieee_fpn[i] =  sign | exponent | mantissa;
		}
	}
}

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

int flip(buffer,numpts)

unsigned char *buffer;

{
	int i,j;
	unsigned char val;

	for(i=0;i<numpts;i++)
	{
		j=i*4;
		val=buffer[j+0];
		buffer[j+0]=buffer[j+3];
		buffer[j+3]=val;
		val=buffer[j+1];
		buffer[j+1]=buffer[j+2];
		buffer[j+2]=val;
	}
}