#include <iostream>
#include <cmath>
#include <string>
using namespace std;
#define PI 3.1415926
int main()
{
float seisx[4375][1],seisy[4375][1],seisr[4375][1],seist[4375][1];
int nsamp=4375,ntraces=479;
for(int i=0;i<nsamp;i++)
{
seisx[i][1]=0.0;
seisy[i][1]=0.0;
seisr[i][1]=0.0;
seist[i][1]=0.0;
}
int sr[179][1],st[179][1];
for(i=0;i<179;i++)
{
sr[i][1]=0;
st[i][1]=0;
cout<<sr[i][1]<<" "<<st[i][1]<<" "<<i<<endl;
}
FILE *fp1,*fp2,*fp3,*fp4;
fp1=fopen("D:\\programming\\matlab\\h17609400x.segy","r");
fp2=fopen("D:\\programming\\matlab\\h17609400y.segy","r");
fp3=fopen("D:\\programming\\matlab\\h17609400t.segy","rb+");
fp4=fopen("D:\\programming\\matlab\\h17609400r.segy","rb+");
for(int seqno=345;seqno<ntraces;seqno++)
{
cout<<seqno+1<<endl;
fseek(fp1,long(3600+(seqno+1)*240+seqno*nsamp*4),0);
for(int i=0;i<nsamp;i++)
{
fread(&seisx[i][1],4,nsamp+1,fp1);
}
fseek(fp2,long(3600+(seqno+1)*240+seqno*nsamp*4),0);
fread(seisy,4,nsamp+1,fp2);
cout<<seisx[seqno][1]<<seisx[seqno][1]<<endl;
float seismaxr1=0.0,seismaxr2=0.0;
int rind1,rind2,rind;
for(i=0;i<nsamp;i++)
{
if(abs(seisy[i][1])>seismaxr1)
{
seismaxr1=abs(seisy[i][1]);
rind1=i;
}
if(abs(seisx[i][1])>seismaxr2)
{
seismaxr2=abs(seisx[i][1]);
rind2=i;
}
}
cout<<rind2<<endl;
if(seismaxr2>seismaxr1)
rind=rind2;
else
rind=rind1;
float seisrmax=0.0;
int seisindex;
if((seisx[rind][1]*seisy[rind][1])>0)
{
for(int j=1;j<=179;j++)
{
float t=(j*PI)/180;
sr[j-1][1]=pow((cos(t)*seisx[rind][1]+sin(t)*seisy[rind][1]),2.0);
st[j-1][1]=pow(((-1)*cos(t)*seisy[rind][1]+sin(t)*seisx[rind][1]),2.0);
}
for(j=1;j<=179;i++)
if(abs(sr[j-1][1])>seisrmax)
{
seisrmax=abs(sr[j-1][1]);
seisindex=j;
}
for(int k=0;k<nsamp;k++)
{
seisindex=seisindex*PI/180;
seisr[k][1]=cos(seisindex)*seisx[k][1]+sin(seisindex)*seisy[k][1];
seist[k][1]=(-1)*cos(seisindex)*seisy[k][1]+sin(seisindex)*seisx[k][1];
}
}
else
{
for(int j=1;j<=179;j++)
{
float t=(j*PI)/180;
sr[j-1][1]=pow(((-1)*cos(t)*seisx[rind][1]+sin(t)*seisy[rind][1]),2.0);
st[j-1][1]=pow((cos(t)*seisy[rind][1]+sin(t)*seisx[rind][1]),2.0);
}
float seisrmax=0.0;
int seisindex;
for(j=1;j<=179;i++)
if(abs(sr[j-1][1])>seisrmax)
{
seisrmax=abs(sr[j-1][1]);
seisindex=j;
}
for(int k=0;k<nsamp;k++)
{
seisindex=seisindex*PI/180;
seisr[k][1]=(-1)*cos(seisindex)*seisx[k][1]+sin(seisindex)*seisy[k][1];
seist[k][1]=cos(seisindex)*seisy[k][1]+sin(seisindex)*seisx[k][1];
}
}
fseek(fp3,long(3600+seqno*(nsamp+60)*4+180),0);
fwrite(&seisindex,4,1,fp3);
fseek(fp3,long(3600+seqno*(nsamp+60)*4+240),0);
for(i=0;i<nsamp;i++)
fwrite(&seist[i][1],4,1,fp3);
fseek(fp4,long(3600+seqno*(nsamp+60)*4+180),0);
fwrite(&seisindex,4,1,fp4);
fseek(fp4,long(3600+seqno*(nsamp+60)*4+240),0);
for(i=0;i<nsamp;i++)
fwrite(&seisr[i][1],4,1,fp4);
}
fclose(fp1);
fclose(fp2);
fclose(fp3);
fclose(fp4);
return 0;
}