Commit a6156e24 authored by niklas.thiel's avatar niklas.thiel
Browse files

added moving streamer geometry caluclation

parent a19a1cf8
...@@ -322,7 +322,9 @@ The locations of the receivers may either be specified in a separate file REC\_F ...@@ -322,7 +322,9 @@ The locations of the receivers may either be specified in a separate file REC\_F
\end{verbatim}}} \end{verbatim}}}
These receiver coordinates in REC\_FILE are shifted by REFREC[1], REFREC[2] into the x- and y-direction, respectively. This allows for completely moving the receiver spread without modifying REC\_FILE. This may be useful for the simulation of moving profiles in reflection seismics. These receiver coordinates in REC\_FILE are shifted by REFREC[1], REFREC[2] into the x- and y-direction, respectively. This allows for completely moving the receiver spread without modifying REC\_FILE. This may be useful for the simulation of moving profiles in reflection seismics.
If READREC=0 the receiver locations must be specified in the parameter file. In this case, it is assumed that the receivers are located along a straight line. The first receiver position is defined by (XREC1, YREC1), and the last receiver position by (XREC1, YREC1). The spacing between receivers is NGEOPH grid points. If READREC=0 the receiver locations must be specified in the parameter file. In this case, it is assumed that the receivers are located along a straight line. The first receiver position is defined by (XREC1, YREC1), and the last receiver position by (XREC1, YREC1). The spacing between receivers is NGEOPH grid points.
READREC=2 is made for a moving streamer aquisition geometry. The receiver positions for shot 1 must be specified in the parameter file (as described for READREC=0). For all other shots the receiver positions are moved according to the movement of the following source. This implementation allows a memory saving usage of moving geometries. The receivers must be in the same depth and all receivers must be located in the model for all shots.
Receivers are always located on full grid indices, i.e. a receiver that is located between two grid points will be shifted by the FD program to the closest next grid point. It is not yet possible Receivers are always located on full grid indices, i.e. a receiver that is located between two grid points will be shifted by the FD program to the closest next grid point. It is not yet possible
to output seismograms for arbitrary receiver locations since this would require a certain wavefield interpolation. to output seismograms for arbitrary receiver locations since this would require a certain wavefield interpolation.
......
...@@ -243,10 +243,14 @@ int main(int argc, char **argv){ ...@@ -243,10 +243,14 @@ int main(int argc, char **argv){
/* In the following, NX and NY denote size of the local grid ! */ /* In the following, NX and NY denote size of the local grid ! */
NX = IENDX; NX = IENDX;
NY = IENDY; NY = IENDY;
/* Reading source positions from SOURCE_FILE */
srcpos=sources(&nsrc);
nsrc_glob=nsrc;
ishot=0;
if (SEISMO){ if (SEISMO){
recpos=receiver(FP, &ntr); recpos=receiver(&ntr, srcpos, ishot);
recswitch = ivector(1,ntr); recswitch = ivector(1,ntr);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch); recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr; ntr_glob=ntr;
...@@ -711,105 +715,14 @@ int main(int argc, char **argv){ ...@@ -711,105 +715,14 @@ int main(int argc, char **argv){
} }
if (ntr>0){ if (ntr>0){
switch (SEISMO){ alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
case 1 : /* particle velocities only */ &sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
switch (WAVETYPE) { &sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
case 1:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
sectionvz=matrix(1,ntr,1,ns);
break;
}
break;
case 2 : /* pressure only */
sectionp=matrix(1,ntr,1,ns);
sectionpnp1=matrix(1,ntr,1,ns);
sectionpn=matrix(1,ntr,1,ns);
break;
case 3 : /* curl and div only */
sectioncurl=matrix(1,ntr,1,ns);
sectiondiv=matrix(1,ntr,1,ns);
break;
case 4 : /* everything */
switch (WAVETYPE) {
case 1:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
sectionvz=matrix(1,ntr,1,ns);
break;
}
sectioncurl=matrix(1,ntr,1,ns);
sectiondiv=matrix(1,ntr,1,ns);
sectionp=matrix(1,ntr,1,ns);
break;
case 5 : /* everything except curl and div*/
switch (WAVETYPE) {
case 1:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
sectionvx=matrix(1,ntr,1,ns);
sectionvy=matrix(1,ntr,1,ns);
sectionvz=matrix(1,ntr,1,ns);
break;
}
sectionp=matrix(1,ntr,1,ns);
break;
}
} }
/* Memory for seismic data */ /* Memory for seismic data */
sectionread=matrix(1,ntr_glob,1,ns); sectionread=matrix(1,ntr_glob,1,ns);
sectionpdata=matrix(1,ntr,1,ns);
sectionpdiff=matrix(1,ntr,1,ns);
sectionpdiffold=matrix(1,ntr,1,ns);
switch (WAVETYPE) {
case 1:
sectionvxdata=matrix(1,ntr,1,ns);
sectionvxdiff=matrix(1,ntr,1,ns);
sectionvxdiffold=matrix(1,ntr,1,ns);
sectionvydata=matrix(1,ntr,1,ns);
sectionvydiff=matrix(1,ntr,1,ns);
sectionvydiffold=matrix(1,ntr,1,ns);
break;
case 2:
sectionvzdata=matrix(1,ntr,1,ns);
sectionvzdiff=matrix(1,ntr,1,ns);
sectionvzdiffold=matrix(1,ntr,1,ns);
break;
case 3:
sectionvxdata=matrix(1,ntr,1,ns);
sectionvxdiff=matrix(1,ntr,1,ns);
sectionvxdiffold=matrix(1,ntr,1,ns);
sectionvydata=matrix(1,ntr,1,ns);
sectionvydiff=matrix(1,ntr,1,ns);
sectionvydiffold=matrix(1,ntr,1,ns);
sectionvzdata=matrix(1,ntr,1,ns);
sectionvzdiff=matrix(1,ntr,1,ns);
sectionvzdiffold=matrix(1,ntr,1,ns);
break;
}
/* Memory for inversion for source time function */ /* Memory for inversion for source time function */
if((INV_STF==1)||(TIME_FILT==1) || (TIME_FILT==2)){ if((INV_STF==1)||(TIME_FILT==1) || (TIME_FILT==2)){
sectionp_conv=matrix(1,ntr_glob,1,NT); sectionp_conv=matrix(1,ntr_glob,1,NT);
...@@ -857,10 +770,6 @@ int main(int argc, char **argv){ ...@@ -857,10 +770,6 @@ int main(int argc, char **argv){
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
/* Reading source positions from SOURCE_FILE */
srcpos=sources(&nsrc);
nsrc_glob=nsrc;
if(FORWARD_ONLY==0&&USE_WORKFLOW){ if(FORWARD_ONLY==0&&USE_WORKFLOW){
read_workflow(FILE_WORKFLOW,&workflow, &workflow_lines,workflow_header); read_workflow(FILE_WORKFLOW,&workflow, &workflow_lines,workflow_header);
} }
...@@ -1176,6 +1085,28 @@ int main(int argc, char **argv){ ...@@ -1176,6 +1085,28 @@ int main(int argc, char **argv){
/*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/
for (ishot=1;ishot<=nshots;ishot+=SHOTINC){ for (ishot=1;ishot<=nshots;ishot+=SHOTINC){
if (SEISMO && READREC==2){
if (ntr>0) {
dealloc_sections(ntr,ns,recpos_loc,sectionvx,sectionvy,sectionvz,sectionp,sectionpnp1,sectionpn,sectioncurl,sectiondiv,
sectionpdata,sectionpdiff,sectionpdiffold,sectionvxdata,sectionvxdiff,sectionvxdiffold,sectionvydata,
sectionvydiff,sectionvydiffold,sectionvzdata,sectionvzdiff,sectionvzdiffold);
}
free_imatrix(recpos,1,3,1,ntr_glob);
recpos=receiver(&ntr, srcpos, ishot);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr;
ntr=ntr_loc;
if (ntr>0){
alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
&sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
&sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
}
if (ntr) group_id = 1;
else group_id = 0;
MPI_Comm_split(MPI_COMM_WORLD, group_id, MYID, &MPI_COMM_NTR);
MPI_Comm_rank(MPI_COMM_NTR, &myid_ntr);
}
SOURCE_SHAPE = SOURCE_SHAPE_OLD; SOURCE_SHAPE = SOURCE_SHAPE_OLD;
if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_SH=SOURCE_SHAPE_OLD_SH; if(WAVETYPE==2 || WAVETYPE==3) SOURCE_SHAPE_SH=SOURCE_SHAPE_OLD_SH;
...@@ -3325,7 +3256,26 @@ int main(int argc, char **argv){ ...@@ -3325,7 +3256,26 @@ int main(int argc, char **argv){
fprintf(FP,"\n=================================================================================================\n"); fprintf(FP,"\n=================================================================================================\n");
fprintf(FP,"\n ***** Starting simulation (test-forward model) no. %d for shot %d of %d (rel. step length %.5f) \n",itest,ishot,nshots,eps_scale); fprintf(FP,"\n ***** Starting simulation (test-forward model) no. %d for shot %d of %d (rel. step length %.5f) \n",itest,ishot,nshots,eps_scale);
fprintf(FP,"\n=================================================================================================\n\n"); fprintf(FP,"\n=================================================================================================\n\n");
if (SEISMO && READREC==2){
if (ntr>0) {
dealloc_sections(ntr,ns,recpos_loc,sectionvx,sectionvy,sectionvz,sectionp,sectionpnp1,sectionpn,sectioncurl,sectiondiv,
sectionpdata,sectionpdiff,sectionpdiffold,sectionvxdata,sectionvxdiff,sectionvxdiffold,sectionvydata,
sectionvydiff,sectionvydiffold,sectionvzdata,sectionvzdiff,sectionvzdiffold);
}
free_imatrix(recpos,1,3,1,ntr_glob);
recpos=receiver(&ntr, srcpos, ishot);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr;
ntr=ntr_loc;
if (ntr>0){
alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
&sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
&sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
}
}
for (nt=1;nt<=8;nt++) srcpos1[nt][1]=srcpos[nt][ishot]; for (nt=1;nt<=8;nt++) srcpos1[nt][1]=srcpos[nt][ishot];
/*-----------------------------------*/ /*-----------------------------------*/
...@@ -4308,53 +4258,9 @@ int main(int argc, char **argv){ ...@@ -4308,53 +4258,9 @@ int main(int argc, char **argv){
break; break;
} }
if ((ntr>0) && (SEISMO)){ if ((ntr>0) && (SEISMO)){
free_imatrix(recpos_loc,1,3,1,ntr); dealloc_sections(ntr,ns,recpos_loc,sectionvx,sectionvy,sectionvz,sectionp,sectionpnp1,sectionpn,sectioncurl,sectiondiv,
switch (SEISMO){ sectionpdata,sectionpdiff,sectionpdiffold,sectionvxdata,sectionvxdiff,sectionvxdiffold,sectionvydata,
case 1 : /* particle velocities only */ sectionvydiff,sectionvydiffold,sectionvzdata,sectionvzdiff,sectionvzdiffold);
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvx,1,ntr,1,ns);
free_matrix(sectionvy,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvz,1,ntr,1,ns);
}
break;
case 2 : /* pressure only */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionp,1,ntr,1,ns);
free_matrix(sectionpn,1,ntr,1,ns);
free_matrix(sectionpnp1,1,ntr,1,ns);
}
break;
case 3 : /* curl and div only */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectioncurl,1,ntr,1,ns);
free_matrix(sectiondiv,1,ntr,1,ns);
}
break;
case 4 : /* everything */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvx,1,ntr,1,ns);
free_matrix(sectionvy,1,ntr,1,ns);
free_matrix(sectionp,1,ntr,1,ns);
free_matrix(sectioncurl,1,ntr,1,ns);
free_matrix(sectiondiv,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvz,1,ntr,1,ns);
}
break;
case 5 : /* everything except curl and div */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvx,1,ntr,1,ns);
free_matrix(sectionvy,1,ntr,1,ns);
free_matrix(sectionp,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvz,1,ntr,1,ns);
}
break;
}
} }
...@@ -4372,23 +4278,7 @@ int main(int argc, char **argv){ ...@@ -4372,23 +4278,7 @@ int main(int argc, char **argv){
free_matrix(We_SH,-nd+1,NY+nd,-nd+1,NX+nd); free_matrix(We_SH,-nd+1,NY+nd,-nd+1,NX+nd);
} }
} }
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvxdata,1,ntr,1,ns);
free_matrix(sectionvxdiff,1,ntr,1,ns);
free_matrix(sectionvydata,1,ntr,1,ns);
free_matrix(sectionvydiff,1,ntr,1,ns);
free_matrix(sectionvydiffold,1,ntr,1,ns);
free_matrix(sectionvxdiffold,1,ntr,1,ns);
free_matrix(sectionpdata,1,ntr,1,ns);
free_matrix(sectionpdiff,1,ntr,1,ns);
free_matrix(sectionpdiffold,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvzdata,1,ntr,1,ns);
free_matrix(sectionvzdiff,1,ntr,1,ns);
free_matrix(sectionvzdiffold,1,ntr,1,ns);
}
free_matrix(sectionread,1,ntr_glob,1,ns); free_matrix(sectionread,1,ntr_glob,1,ns);
......
...@@ -70,6 +70,7 @@ IFOS2D= \ ...@@ -70,6 +70,7 @@ IFOS2D= \
IFOS2D.c \ IFOS2D.c \
stf.c \ stf.c \
window_cos.c \ window_cos.c \
alloc_sections.c \
calc_mat_change_test.c \ calc_mat_change_test.c \
calc_res.c \ calc_res.c \
calc_misfit.c \ calc_misfit.c \
...@@ -205,4 +206,4 @@ clean: ...@@ -205,4 +206,4 @@ clean:
install: clean all install: clean all
-include $(IFOS_OBJ:.o=.d) -include $(IFOS_OBJ:.o=.d)
\ No newline at end of file
...@@ -277,7 +277,7 @@ void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus, ...@@ -277,7 +277,7 @@ void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus,
therefore we determine the minimum/maximum position in y-direction by the ZREC1 variable and vice versa. therefore we determine the minimum/maximum position in y-direction by the ZREC1 variable and vice versa.
this has to be considered for the receiver line coordinates specified in both the input file and separate source/receiver files*/ this has to be considered for the receiver line coordinates specified in both the input file and separate source/receiver files*/
if (READREC==0) { if (READREC==0 || READREC==2) {
if (XREC1>XREC2) { if (XREC1>XREC2) {
srec_maxx=XREC1; srec_maxx=XREC1;
srec_minx=XREC2; srec_minx=XREC2;
......
...@@ -47,6 +47,10 @@ void spat_filt(float ** waveconv, int iter, int sws); ...@@ -47,6 +47,10 @@ void spat_filt(float ** waveconv, int iter, int sws);
float norm(float ** waveconv, int iter, int sws); float norm(float ** waveconv, int iter, int sws);
void alloc_sections(int ntr,int ns,float ***sectionvx,float ***sectionvy,float ***sectionvz,float ***sectionp,float ***sectionpnp1,float ***sectionpn,float ***sectioncurl,float ***sectiondiv,
float ***sectionpdata,float ***sectionpdiff,float ***sectionpdiffold,float ***sectionvxdata,float ***sectionvxdiff,float ***sectionvxdiffold,float ***sectionvydata,
float ***sectionvydiff,float ***sectionvydiffold,float ***sectionvzdata,float ***sectionvzdiff,float ***sectionvzdiffold);
void av_mat(float ** pi, float ** u, void av_mat(float ** pi, float ** u,
float ** ppijm, float ** puip, float ** pujm); float ** ppijm, float ** puip, float ** pujm);
...@@ -95,6 +99,10 @@ void count_killed_traces(int ntr, int swstestshot, int ntr_glob, int **recpos_lo ...@@ -95,6 +99,10 @@ void count_killed_traces(int ntr, int swstestshot, int ntr_glob, int **recpos_lo
void create_trkill_table(int ** killtable, int ntr_glob, int **recpos, int nsrc_glob, float **srcpos, int ishot, float kill_offset_lower, float kill_offset_upper); void create_trkill_table(int ** killtable, int ntr_glob, int **recpos, int nsrc_glob, float **srcpos, int ishot, float kill_offset_lower, float kill_offset_upper);
void dealloc_sections(int ntr,int ns,int **recpos_loc,float **sectionvx,float **sectionvy,float **sectionvz,float **sectionp,float **sectionpnp1,float **sectionpn,float **sectioncurl,float **sectiondiv,
float **sectionpdata,float **sectionpdiff,float **sectionpdiffold,float **sectionvxdata,float **sectionvxdiff,float **sectionvxdiffold,float **sectionvydata,
float **sectionvydiff,float **sectionvydiffold,float **sectionvzdata,float **sectionvzdiff,float **sectionvzdiffold);
float exchange_L2(float L2, int sw, int bcast_l2); float exchange_L2(float L2, int sw, int bcast_l2);
void exchange_rsg(float ** vx, float ** vy, float ** vz, void exchange_rsg(float ** vx, float ** vy, float ** vz,
...@@ -226,7 +234,7 @@ void readmod_elastic(float ** rho, float ** pi, float ** u); ...@@ -226,7 +234,7 @@ void readmod_elastic(float ** rho, float ** pi, float ** u);
void readmod_elastic_es(float ** rho, float ** pi, float ** u, float ** matmod, int is); void readmod_elastic_es(float ** rho, float ** pi, float ** u, float ** matmod, int is);
int **receiver(FILE *fp, int *ntr); int **receiver(int* ntr, float** srcpos, int shotno);
void save_checkpoint(int nx1, int nx2, int ny1, int ny2, void save_checkpoint(int nx1, int nx2, int ny1, int ny2,
float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy); float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy);
......
...@@ -354,7 +354,7 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -354,7 +354,7 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("READREC",number_readobjects,&READREC,varname_list, value_list)) if (get_int_from_objectlist("READREC",number_readobjects,&READREC,varname_list, value_list))
declare_error("Variable READREC could not be retrieved from the json input file!"); declare_error("Variable READREC could not be retrieved from the json input file!");
else { else {
if (READREC==0) { if (READREC==0 || READREC==2) {
if (get_float_from_objectlist("XREC1",number_readobjects,&XREC1,varname_list, value_list)) if (get_float_from_objectlist("XREC1",number_readobjects,&XREC1,varname_list, value_list))
declare_error("Variable XREC1 could not be retrieved from the json input file!"); declare_error("Variable XREC1 could not be retrieved from the json input file!");
if (get_float_from_objectlist("XREC2",number_readobjects,&XREC2,varname_list, value_list)) if (get_float_from_objectlist("XREC2",number_readobjects,&XREC2,varname_list, value_list))
...@@ -366,7 +366,7 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -366,7 +366,7 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("NGEOPH",number_readobjects,&NGEOPH,varname_list, value_list)) if (get_int_from_objectlist("NGEOPH",number_readobjects,&NGEOPH,varname_list, value_list))
declare_error("Variable NGEOPH could not be retrieved from the json input file!"); declare_error("Variable NGEOPH could not be retrieved from the json input file!");
} }
else { if (READREC>0) {
if (get_string_from_objectlist("REC_FILE",number_readobjects,REC_FILE,varname_list, value_list)) if (get_string_from_objectlist("REC_FILE",number_readobjects,REC_FILE,varname_list, value_list))
declare_error("Variable REC_FILE could not be retrieved from the json input file!"); declare_error("Variable REC_FILE could not be retrieved from the json input file!");
} }
...@@ -1052,7 +1052,7 @@ void read_par_json(FILE *fp, char *fileinp){ ...@@ -1052,7 +1052,7 @@ void read_par_json(FILE *fp, char *fileinp){
}*/ }*/
/* receiver file */ /* receiver file */
if (READREC) if (READREC==1)
{ {
if (access(REC_FILE,0) != 0) if (access(REC_FILE,0) != 0)
{ {
......
...@@ -24,42 +24,44 @@ ...@@ -24,42 +24,44 @@
#include "fd.h" #include "fd.h"
int **receiver(FILE *fp, int *ntr){ int **receiver(int *ntr, float **srcpos, int shotno){
/* declaration of extern variables */ /* declaration of extern variables */
extern char REC_FILE[STRING_SIZE]; extern char REC_FILE[STRING_SIZE];
extern float DH, REFREC[4], REC_ARRAY_DEPTH, REC_ARRAY_DIST, FW; extern float DH, REFREC[4], REC_ARRAY_DEPTH, REC_ARRAY_DIST, FW;
extern int READREC, NGEOPH, DRX, REC_ARRAY, NX; extern int READREC, NGEOPH, DRX, REC_ARRAY, NXG, NYG;
extern int MYID; extern int MYID, VERBOSE,TRKILL;
extern FILE *FP;
int **recpos1, **recpos = NULL, nxrec=0, nyrec=0, nzrec=0; int **recpos1, **recpos = NULL, nxrec=0, nyrec=0, nzrec=0;
int itr=1, itr1=0, itr2=0, recflag=0, c, ifw, n, i, j; int itr=1, itr1=0, itr2=0, recflag=0, c, ifw, n, i, j, ntr1=0;
int nxrec1, nxrec2, nyrec1, nyrec2, nzrec1=1, nzrec2=1; int nxrec1, nxrec2, nxrec3, nyrec1, nyrec2, nzrec1=1, nzrec2=1, offset, shotdist;
extern float XREC1, YREC1, XREC2, YREC2; extern float XREC1, YREC1, XREC2, YREC2;
float xrec, yrec; float xrec, yrec;
FILE *fpr; FILE *fpr,*f;
char filename[STRING_SIZE];
if (MYID==0) if (MYID==0)
{ {
if (READREC){ /* read receiver positions from file */ if (READREC==1){ /* read receiver positions from file */
fprintf(fp,"\n Reading receiver positions from file: \n\t%s\n",REC_FILE); fprintf(FP,"\n Reading receiver positions from file: \n\t%s\n",REC_FILE);
fpr=fopen(REC_FILE,"r"); fpr=fopen(REC_FILE,"r");
if (fpr==NULL) declare_error(" Receiver file could not be opened !"); if (fpr==NULL) declare_error(" Receiver file could not be opened !");
*ntr=0; *ntr=0;
while ((c=fgetc(fpr)) != EOF) while ((c=fgetc(fpr)) != EOF)
if (c=='\n') ++(*ntr); if (c=='\n') ++(ntr1);
rewind(fpr); rewind(fpr);
recpos1=imatrix(1,3,1,*ntr); recpos1=imatrix(1,3,1,ntr1);
for (itr=1;itr<=*ntr;itr++){ for (itr=1;itr<=ntr1;itr++){
fscanf(fpr,"%f%f\n",&xrec, &yrec); fscanf(fpr,"%f%f\n",&xrec, &yrec);
recpos1[1][itr]=iround((xrec+REFREC[1])/DH); recpos1[1][itr]=iround((xrec+REFREC[1])/DH);
recpos1[2][itr]=iround((yrec+REFREC[2])/DH); recpos1[2][itr]=iround((yrec+REFREC[2])/DH);
recpos1[3][itr]=iround((0.0+REFREC[3])/DH); recpos1[3][itr]=iround((0.0+REFREC[3])/DH);
} }
fclose(fpr); fclose(fpr);
fprintf(fp," Message from function receiver (written by PE %d):\n",MYID);/***/ fprintf(FP," Message from function receiver (written by PE %d):\n",MYID);/***/
fprintf(fp," Number of receiver positions found: %i\n",*ntr); fprintf(FP," Number of receiver positions found: %i\n",ntr1);
/* recpos1[1][itr] X in m /* recpos1[1][itr] X in m
* recpos1[2][itr] Y in m * recpos1[2][itr] Y in m
...@@ -68,15 +70,15 @@ int **receiver(FILE *fp, int *ntr){ ...@@ -68,15 +70,15 @@ int **receiver(FILE *fp, int *ntr){
/* check if more than one receiver is located /* check if more than one receiver is located
at the same gridpoint */ at the same gridpoint */
for (itr=1;itr<=(*ntr-1);itr++) for (itr=1;itr<=(ntr1-1);itr++)
for (itr1=itr+1;itr1<=*ntr;itr1++) for (itr1=itr+1;itr1<=ntr1;itr1++)
if ((recpos1[1][itr]==recpos1[1][itr1]) if ((recpos1[1][itr]==recpos1[1][itr1])
&& (recpos1[2][itr]==recpos1[2][itr1]) && (recpos1[2][itr]==recpos1[2][itr1])
&& (recpos1[3][itr]==recpos1[3][itr1])) && (recpos1[3][itr]==recpos1[3][itr1]))
recpos1[1][itr1]=-(++recflag); recpos1[1][itr1]=-(++recflag);
recpos=imatrix(1,3,1,*ntr-recflag); recpos=imatrix(1,3,1,ntr1-recflag);
for (itr=1;itr<=*ntr;itr++) for (itr=1;itr<=ntr1;itr++)
if (recpos1[1][itr]>0){ if (recpos1[1][itr]>0){
recpos[1][++itr2]=recpos1[1][itr]; recpos[1][++itr2]=recpos1[1][itr];
recpos[2][itr2]=recpos1[2][itr]; recpos[2][itr2]=recpos1[2][itr];
...@@ -85,36 +87,72 @@ int **receiver(FILE *fp, int *ntr){ ...@@ -85,36 +87,72 @@ int **receiver(FILE *fp, int *ntr){
*ntr=itr2; *ntr=itr2;
if (recflag>0){ if (recflag>0){
fprintf(fp,"\n\n"); fprintf(FP,"\n\n");
fprintf(fp," Warning:\n"); fprintf(FP," Warning:\n");
fprintf(fp," Several receivers located at the same gridpoint !\n"); fprintf(FP," Several receivers located at the same gridpoint !\n");
fprintf(fp," Number of receivers reduced to %i\n", *ntr); fprintf(FP," Number of receivers reduced to %i\n", *ntr);
fprintf(fp,"\n\n"); fprintf(FP,"\n\n");
} }
free_imatrix(recpos1,1,3,1,*ntr); free_imatrix(recpos1,1,3,1,ntr1);
} }
else if (REC_ARRAY>0){ else if (READREC==2){
ifw=iround(FW/DH); /* frame width in gridpoints */ fprintf(FP,"\n Moving streamer acquisition geometry choosen. \n");
*ntr=(1+(NX-2*ifw)/DRX)*REC_ARRAY; nxrec1=iround(XREC1/DH); /* (nxrec1,nyrec1) and (nxrec2,nyrec2) are */
recpos=imatrix(1,3,1,*ntr); nyrec1=iround(YREC1/DH); /* the positions of the first and last receiver*/
itr=0; nxrec2=iround(XREC2/DH); /* in gridpoints */
for (n=0;n<=REC_ARRAY-1;n++){ nyrec2=iround(YREC2/DH);
j=iround((REC_ARRAY_DEPTH+REC_ARRAY_DIST*(float)n)/DH); offset=iround((XREC1-srcpos[1][1])/DH);