//#############################################################a
//##s
//##s 2D Acoustic VTI Medium RTM & Pick Up ADCIG
//##s
//##s----------------------------------------------------------
//##s PML , P&SV , VTI , Acoustic , RTM , -SV , Adcig,
//##s Migration , Both sides receive , SU or dat(v,e,d) ,
//##s Adcig( poynting , s-r , offset-domain ),
//##s
//##s----------------------------------------------------------
//##s Rong Tao
//##s 2016
//############################################################a
#include<stdio.h>
#include<malloc.h>
#include<math.h>
#include<stdlib.h>
#include "mpi.h"
#include "/home/Toa/hc/alloc.c"
#include "/home/Toa/hc/alloc.h"
#include "/home/Toa/hc/complex.c"
#include "/home/Toa/hc/complex.h"
#include "/home/Toa/hc/bhdr.h"
#include "/home/Toa/hc/hdr.h"
#include "/home/Toa/hc/fft.h"
#include "/home/Toa/hc/fft.c"
#include "/home/Toa/hc/mute_direct.h"
#include "/home/Toa/hc/cjbsegy.h"
/********* SEG-Y header *********/
typedef cjbsegy segy;
/********** SU & SEG-Y **********/
#define SU_NKEYS 80 /* Number of key header words */
#define HDRBYTES 240 /* Bytes in the trace header */
#define EBCBYTES 3200 /* Bytes in the card image EBCDIC block */
#define BNYBYTES 400 /* Bytes in the binary coded block */
#define SEGY_HDRBYTES 240 /* Bytes in the tape trace header */
#define SEGY_NKEYS 71 /* Number of mandated header fields */
#define BHED_NKEYS 27 /* Number of mandated binary fields */
#define pi 3.1415926535898
/********* SEG-Y header *********/
Y_3200 y3200;
bhed head_400;
cjbsegy tr, vtr;
/********* SEG-Y header *********/
/********** SU function *********/
void swap_short_2(short *tni2);
void swap_u_short_2(unsigned short *tni2);
void swap_int_4(int *tni4);
void swap_u_int_4(unsigned int *tni4);
void swap_long_4(long *tni4);
void swap_u_long_4(unsigned long *tni4);
void swap_float_4(float *tnf4);
void swap_double_8(double *tndd8);
void swaphval(segy *tr, int index);
/********* SU function *********/
/* MAIN */
main(int argc,char *argv[])
{
/*************************************** function ********************************************/
void model_vti_get_boundry(int nx,int nz,int vnx,int vnz,int nt,int npd,float dx,float dz,
float vdx,float vdz,float favg,float tmax,float dt,float dtout,float pfac,
float **vp,float **epsilu,float **deta,float **rho,
int ns_sxd,int ds_sxd,int fs_sxd,int zs_sxd,int is,
float **p_cal_x,float **p_cal_z,
float **p_top_x,float **p_bottom_x,float **p_left_x,float **p_right_x,
float **p_top_z,float **p_bottom_z,float **p_left_z,float **p_right_z,int _Circle_,
int mm,int wtype,int hsx,int myid,float *mu_v,int flag_snap,int seismic);
void mute_directwave(int flag_mu,int nx,int nt,float dt,float favg,
float dx,float dz,int fs_sxd,int ds_sxd,int zs_sxd,int is,
float mu_v,float **p_cal,int tt);
void RTM_corr_adcig(int nx,int nz,int vnx,int vnz,int nt,int npd,float dx,float dz,
float vdx,float vdz,float favg,float tmax,float dt,float dtout,
float pfac,float **vp,float **epsilu,float **deta,float **rho,char FN5[],
int ns_sxd,int ds_sxd,int fs_sxd,int ds_initial,int fs_initial,int zs_sxd,int is,int flag_cdp,
float **p_top_x,float **p_bottom_x,float **p_left_x,float **p_right_x,////////////////////////////////
float **p_top_z,float **p_bottom_z,float **p_left_z,float **p_right_z,float **p_obs_x,float **p_obs_z,
int mm,int wtype,int hsx,int myid,float **mig_is,float **mig_ns0,
float ***adcig_is,float ***adcig_ns0,float ***Ixhz_is,float ***Ixhz_ns0,int nh,
int flag_snap,int seismic,int flag_adcig);
void smooth1float(float *v,int r,int n);
void smooth2float(int nx,int rx,int nz,int rz,float **v);
void pad_vv(int nx,int nz,int npd,float **ee);
void read_v_e_d_r(char FN1[],char FN2[],char FN3[],int nx,int nz,float **vv,float **epsilu,float **deta,
float **rho0,int npd,int seismic);
/*************************** parameter statement ***********************/
int i,j,k,l,m,ih,is,nx,nz,nt,vnx,vnz,i_start,i_end,mm,wtype,hsx,ia,nxs,flag_cdp,flag_adcig,nh;
int ns_sxd,ds_sxd,fs_sxd,zs_sxd,fs,ds,npd;
float dx,dz,vdx,vdz,tmax,dt,dtout,pfac,favg;
int myid,numprocs,_Circle_,flag_mu,flag_snap,seismic,nangle,dangle,fangle;
float mu_v;
/***************** wave float **************/
float **p_cal_x,**p_cal_z,**p_top_x,**p_bottom_x,**p_left_x,**p_right_x;
float **p_top_z,**p_bottom_z,**p_left_z,**p_right_z;
/***************** image float *************/
float **mig_is,**mig_ns,**mig_ns0,***adcig_is,***adcig_ns,***adcig_ns0,**adcig2d;
float ***Ixhz_ns0,***Ixhz_ns,***Ixhz_is;
float **vp,**rho,**deta,**epsilu;
float **vps,**rhos,**detas,**epsilus;
//a###########################################################################
//#### input the parameter and confirm ####
//a###########################################################################
/************************ dat document **********************/
/* Input velocity */char FN1[250]={"vel_1801_862.dat"};
/* Input epsilu */char FN2[250]={"eps_1801_862.dat"};
/* Input deta */char FN3[250]={"del_1801_862.dat"};
/* Cal data */char FN4[250]={"shot_cal.dat"};
/* Obs data */char FN5[250]={"shot_obs.dat"};
/* Migration */char FN6[250]={"mig_ns.dat"};
/* IS tempor migration */char FN7[250]={"mig_is_tempor.dat"};
/* Adcig initial */char FN8[250]={"adcig.dat"};
/* Adcig smooth */char FN9[250]={"adcig_smooth.dat"};
/* Migration adcig stack */char FN10[250]={"mig_stack_adcig.dat"};
/* Half offset I(x,h,z) */char FN11[250]={"Image_x_h_z.dat"};
/*************************** type ***************************/
/* Wavelet type (1-3) */ wtype=1;
/* 1-mute , 0->don't mute */ flag_mu=1;
/* 1-snap , 0->nosnap */ flag_snap=0;
/* 1.dat, 2.su (v,e,d) */ seismic=1;
/*************************** cdp ****************************/
/* Activate "nxs"(=1) */ flag_cdp=1; /* 1-activate ;0-invaliable */
/* CDP of each shot */ nxs=501; /* Must be odd number */
/*************************** cdp ****************************/
/* choice adcig type */ flag_adcig=1; /* 1-poynting */
/* 2-source-receivers domain */
/* 3-half-offset-domain */
/*************************** wave ***************************/
hsx=1;mm=4;npd=20;_Circle_=15;
/******************** observation system ********************/
favg=20; pfac=1000.0;
nx=1801; dx=10.0;
nz=862; dz=5.0;
tmax=6.5;
dt=0.5;//ms
nt=11001;
ns_sxd=250;
fs_sxd=255;
ds_sxd=5;
zs_sxd=1;
nangle=90;
fangle=0;
dangle=1;
/*************************v****************************/
vdz=dz;vdx=dx;vnx=nx;vnz=nz;dtout=dt;
/************************FILE**************************/
FILE *fp4,*fp6,*fp7,*fp8,*fp9,*fp10,*fp11;
fp4=fopen(FN4,"wb"); /* Cal data */
fp6=fopen(FN6,"wb"); /* Migration */
fp7=fopen(FN7,"wb"); /* IS tempor migration */
fp8=fopen(FN8,"wb"); /* Adcig initial */
fp9=fopen(FN9,"wb"); /* Adcig smooth */
fp10=fopen(FN10,"wb"); /* Migration adcig stack */
if(flag_adcig==3)
fp11=fopen(FN11,"wb"); /* Image half-offset */
/************************** read_file **************
