PDA

View Full Version : azimuth search



khawar_geo
11-24-2016, 02:38 PM
Dear friends.

Good morning.
I hope you are fine. I need help in the azimuth search in the code. The code is only working for single azimuth. Please if you can help me in resolving the issue. Humble request.

smartkhawar@gmail.com
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\\\\\\\\\\\\\\\\\\\\\\\\
#include"main.h"

#define mymin(a,b) (a < b ? a : b)
#define mymax(a,b) (a > b ? a : b)
#define mysqr(a) (a*a)

#define PI 3.1415927

using namespace std;

int main (int argc, char** argv) {

///MPI begin
MPI_Init(&argc,&argv);
int mpi_rank=0;
MPI_Comm_rank(MPI_COMM_WORLD,&mpi_rank);
int mpi_size;
MPI_Comm_size(MPI_COMM_WORLD,&mpi_size);

if(argc==1){
PipeToMore(doc);
return 0;
}

///declarationen of variables///////////////

float Gx, Gy;
float Sx, Sy,s;
double starttime = 0.0;
float velmin, velmax, *vel, dv, v;
double endtime = 0.0;
float dt, dx, apert;
float x, xm,tt, t0;
int ns, nv, s_count, rcoher, time, flag;
int idown, iup;
string str;
string Inputfile, predictionfile;
string Outputfile;
float Value, Val_up, Val_down, Value_temp, Val_temp;
float *sembl_up, *sembl_down;
string Output_name[mpi_size];
float apery, aperx;

///Getting parameter /////////////////////////////////////////////////////////////////////
Parser* parser = new Parser(argc, argv);
int *intdummy;
float *fdummy;


Inputfile = parser->parse("input");
if (!Inputfile.length()) {
if (mpi_rank==0)
cerr << "ERROR: No In-filename found." << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}

if (!parser->parse("coherR",intdummy)) {
if (mpi_rank==0) cerr << "ERROR: No coherence band width found. Use coherR=" << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}
rcoher=*intdummy;

if (!parser->parse("mode",intdummy)) {
if (mpi_rank==0) cerr << "ERROR: No mode found. Use mode= 1 for 3D od 0 for 2D" << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}
int moda=*intdummy;


if (!parser->parse("aperx",fdummy)) {
if (mpi_rank==0) cerr << "ERROR: No X-aperture found. Use aperx =" << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}
aperx=apert=*fdummy;

if (moda ==1){
if (!parser->parse("apery",fdummy)) {
if (mpi_rank==0) cerr << "ERROR: No Y-aperture found. Use apery =" << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}

apery=*fdummy;
}

if (!parser->parse("cmpdx",fdummy)) {
if (mpi_rank==0) cerr << "ERROR: No cmp distance found. Use cmpdx =" << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}
dx=*fdummy;


Outputfile = parser->parse("output");
if (!Outputfile.length()) {
if (mpi_rank==0)
cerr << "ERROR: No output-filename found." << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}

float *vtarget;
if(!parser->parse("vmin",vtarget))
{
if (mpi_rank==0)
cerr << "ERROR: No velocity found." << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}
velmin=*vtarget;
if(!parser->parse("vmax",vtarget))
{
if (mpi_rank==0)
cerr << "ERROR: No velocity found." << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}
velmax=*vtarget;

float vel_temp;
if( velmax < velmin)
{
vel_temp=velmax; velmax=velmin; velmin=vel_temp;
}

if(!parser->parse("dv",vtarget))
{
if (mpi_rank==0)
cerr << "ERROR: No velocity found." << endl;
MPI_Abort(MPI_COMM_WORLD,1);
}
dv=*vtarget;

delete parser;
parser=0;


/// SU data read in...
SUdata *Input_data, *Output;
Input_data = new SUdata(Inputfile.c_str());
Input_data->readFile();
SUdata *Output_data[mpi_size];

if ( moda==0 ) apert=aperx;

///allocate all arrays and variable's assignment
vector<Trace>::iterator Input_it, Input_temp, it_left, it_right;
dt = float(Input_data->dt);
ns = int(Input_data->nt);
short fac1=short(Input_data->scalco);
nv= static_cast<int>((velmax -velmin)/dv);
vel = new float[nv+1];
float dtnew=dt/1000000.;
int aper = static_cast<int>(apert/dx);
int width = rcoher*2 ;
sembl_up = new float[width+1];
sembl_down = new float[width+1];
float fac, ym, y, t_aper;

if ( ( fac1 ) || ( fac1 != 0 ) ) fac=1./(1.*abs(fac1));
else fac=1.;
///create string array for output name
for (int i=mpi_rank; i<mpi_size;i++)
{
Output_name[i]+=Outputfile.c_str();
str=ItoA(i);
Output_name[i]+=str;
}


if (mpi_rank==0) {
cout << endl
<< "INFO: Velocity model building" << endl
<< "INFO: Copyright (C) Stefan Duemmong & Sergius Dell - 2009" << endl
<< "INFO: sergius.dell@zmaw.de" << endl
<< endl;
cout << endl<<"Starting program on "<<mpi_size<<" nodes."<<endl<<endl;
cout << endl
<< "INFO: Using: " << endl
<< "INFO: kind of processing: time-migration velocity analysis" <<endl
<< "INFO: infile: " << Inputfile << endl
<< "INFO: outfile: " << Outputfile << endl
<< "INFO: cmp distance: " <<dx<<endl
<< "INFO: scaling: " <<fac<<endl
<< "INFO: radius of coherence band: " <<rcoher<<endl
<< "INFO: X aperture: " <<2*aperx<<endl;
if (moda==1) cout<< "INFO: Y aperture: " <<2*apery<<endl
<< "INFO: mimimal velocity: " <<velmin<<endl
<< "INFO: maximal velocity: " <<velmax<<endl
<< "INFO: search step: " <<dv<<endl
<< endl;
}
///allocate output data
for (int i=mpi_rank; i<mpi_size;i++)
{
Output_data[i]=new SUdata(Output_name[i].c_str());
}


for (int i=0;i<=nv;i++) {
vel[i]=velmin + i*int(dv);
}

///Time start
starttime=MPI_Wtime();


if (moda==0 ) {
///Loop ovel all velocities

for (int j=0; j<=nv;j++) {
v=vel[j];
if (mpi_rank==0)
cout << "INFO: 2D semblance analysis: velocity " << v<<" is done"<<endl;
/// Loop over all traces: //////////////////////////////////////////////////
for ( int i=mpi_rank; i<Input_data->_traces.size(); i+=mpi_size) {

Input_it = Input_data->_traces.begin() + i;

/// Output preparing: ///////////////////////////////////////////////////////////
Gx = Input_it->gx;
Sx = Input_it->sx;
xm=0.5*(Gx+Sx)*fac;
Output_data[mpi_rank]->_traces.push_back( Trace(ns));
Output_data[mpi_rank]->_traces.back().gx = int(Gx);
Output_data[mpi_rank]->_traces.back().sx = int(Sx);
Output_data[mpi_rank]->_traces.back().dt = int(dt);
Output_data[mpi_rank]->_traces.back().offset = int(v);
Output_data[mpi_rank]->_traces.back().cdp = int(Input_it->cdp);
Output_data[mpi_rank]->_traces.back().scalco = short(fac1);
( i < aper )? it_left=Input_data->_traces.begin() +1 :it_left=Input_it - aper;
( (Input_data->_traces.size() - i) > aper )? it_right=Input_it + aper:it_right=Input_data->_traces.end()-1;

/// Loop over all samples: /////////////////////////////////////////////////
for (int ii=0; ii<ns; ii++)
{
t0=ii*dtnew;
Value=0.;
s_count=0;
for (int is=0; is < width;is++){
sembl_up[is]=0.;sembl_down[is]=0.;
}
for ( Input_temp=it_left;Input_temp != it_right;Input_temp++) {
x=0.5*(Input_temp->gx + Input_temp->sx)*fac ; /// xm difference
tt=sqrt(t0*t0+4.*(xm-x)*(xm-x)/(v*v)); /// calculate tt hyperbola over all times
idown=static_cast<int>(tt/dtnew); ///calculate the current index in arrays
if (Input_temp > Input_data->_traces.end()) break;
( (idown -rcoher ) > 0)? idown=idown-rcoher: idown=idown;
( (idown + rcoher) > ns)? idown=ns - rcoher: idown=idown;
s_count++; /// trace counter
for (int is=0; is<width;is++){ /// semblance within the coherence band
Value_temp=float(Input_temp->data[idown+is]);
sembl_up[is]+=Value_temp;
sembl_down[is]+=Value_temp*Value_temp;
}
s=1.*s_count;
}
Val_up= Val_down=0.0;
for (int is=0; is<width;is++){
Val_up+=sembl_up[is]*sembl_up[is];
Val_down+=sembl_down[is];
}
Value=Val_up/Val_down/s;
Output_data[mpi_rank]->_traces.back().data[ii]=Value;
}
}
}

}
if (moda==1 ) {
///Loop ovel all velocities

for (int j=0; j<=nv;j++) {
v=vel[j];
if (mpi_rank==0)
cout << "INFO: 3D semblance analysis: velocity " << v<<" is done"<<endl;
/// Loop over all traces: //////////////////////////////////////////////////
for ( int i=mpi_rank; i<Input_data->_traces.size(); i+=mpi_size) {
Input_it = Input_data->_traces.begin() + i;

/// Output preparing: ///////////////////////////////////////////////////////////
Gx = Input_it->gx;
Sx = Input_it->sx;
Gy = Input_it->gy;
Sy = Input_it->sy;
xm=0.5*(Gx+Sx)*fac;
ym=0.5*(Gy+Sy)*fac;
Output_data[mpi_rank]->_traces.push_back( Trace(ns));
Output_data[mpi_rank]->_traces.back().gx = int(Gx);
Output_data[mpi_rank]->_traces.back().sx = int(Sx);
Output_data[mpi_rank]->_traces.back().gy = int(Gy);
Output_data[mpi_rank]->_traces.back().sy = int(Sy);
Output_data[mpi_rank]->_traces.back().dt = int(dt);
Output_data[mpi_rank]->_traces.back().offset = int(v);
Output_data[mpi_rank]->_traces.back().cdp = int(Input_it->cdp);
Output_data[mpi_rank]->_traces.back().scalco = short(fac1);
it_left=Input_data->_traces.begin() +1;
it_right=Input_data->_traces.end()-1;

/// Loop over all samples: /////////////////////////////////////////////////
for (int ii=0; ii<ns; ii++)
{
t0=ii*dtnew;
Value=0.;
s_count=0;
for (int is=0; is < width;is++){
sembl_up[is]=0.;sembl_down[is]=0.;
}
for ( Input_temp=it_left;Input_temp != it_right;Input_temp++) {

y=0.5*(Input_temp->gy + Input_temp->sy)*fac; /// ym difference
x=0.5*(Input_temp->gx + Input_temp->sx)*fac; /// xm difference
if ( (xm-x)*(xm-x)/(aperx*aperx) + (ym-y)*(ym-y)/(apery*apery) > 1.) continue; /// aperture
tt=sqrt(t0*t0+4.*(xm-x)*(xm-x)/(v*v) +4.*(ym-y)*(ym-y)/(v*v)); /// calculate tt hyperbola over all times
idown=static_cast<int>(tt/dtnew); ///calculate the current index in arrays
if (Input_temp > Input_data->_traces.end()) break;
( (idown -rcoher ) > 0)? idown=idown-rcoher: idown=idown;
( (idown + rcoher) > ns)? idown=ns - rcoher: idown=idown;
s_count++; /// trace counter
for (int is=0; is<width;is++){ /// semblance within the coherence band
Value_temp=float(Input_temp->data[idown+is]);
sembl_up[is]+=Value_temp;
sembl_down[is]+=Value_temp*Value_temp;
}
s=1.*s_count;
}
Val_up= Val_down=0.0;
for (int is=0; is<width;is++){
Val_up+=sembl_up[is]*sembl_up[is];
Val_down+=sembl_down[is];
}
Value=Val_up/Val_down/s;
if ( isinf( Value) || isnan(Value ) ) Value=0.0;
Output_data[mpi_rank]->_traces.back().data[ii]=Value;
}
}
}
}

/// Write Data////////////////////////////////////////////////////////////////////////////////

// MPI_Barrier(MPI_COMM_WORLD);
// Output_data[mpi_rank]->write();

int loopn;
Output=new SUdata(Outputfile.c_str());

for (int i=mpi_rank; i<mpi_size;i++) {
Input_it = Output_data[mpi_rank]->_traces.begin();
loopn= Output_data[mpi_rank]->_traces.size();

for (int i=1; i<loopn; i++)
{

Output->_traces.push_back( Trace(ns));
Output->_traces.back().gx = int(Input_it->sx);
Output->_traces.back().sx = int(Input_it->gx);
Output->_traces.back().gy = int(Input_it->gy);
Output->_traces.back().sy = int(Input_it->sy);
Output->_traces.back().dt = int(Input_it->dt);
Output->_traces.back().offset = int(Input_it->offset);
Output->_traces.back().cdp = int(Input_it->cdp);
Output->_traces.back().scalco = short(Input_it->scalco);
for (int ii=0; ii<ns; ii++) Output->_traces.back().data[ii]=Input_it->data[ii];
Input_it++;
}
}

MPI_Barrier(MPI_COMM_WORLD);
Output->write();

///Time stop
if (mpi_rank==0) {
endtime=MPI_Wtime();
time=static_cast<int>(endtime-starttime);
sec_to_time(time);
}
///MPI end
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
}
/// end main
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ \\\\\\\\\\\\\\\\\\\\\\\\\