-
Jan Thorbecke authored
removing unused variables in all code, bug fixes in marchenko to handle different acquisition geometries, re-organizing calls within marchenko algorithm in preperation for software release
Jan Thorbecke authoredremoving unused variables in all code, bug fixes in marchenko to handle different acquisition geometries, re-organizing calls within marchenko algorithm in preperation for software release
recvPar.c 17.58 KiB
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include "fdelmodc.h"
#include "par.h"
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
/**
* Calculates the receiver positions based on the input parameters
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
*
* Ammendments:
* Max Holicki changing the allocation receiver array (2-2016)
* The Netherlands
**/
void name_ext(char *filename, char *extension);
int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx, int nz)
{
float *xrcv1, *xrcv2, *zrcv1, *zrcv2;
int i, ix, ir, verbose;
float dxrcv, dzrcv, *dxr, *dzr;
float rrcv, dphi, oxrcv, ozrcv, arcv;
double circ, h, a, b, e, s, xr, zr, dr, srun, phase;
float xrange, zrange, sub_x1, sub_z1;
int Nx1, Nx2, Nz1, Nz2, Ndx, Ndz, iarray, nrec, nh;
int nxrcv, nzrcv, ncrcv, nrcv, ntrcv, *nlrcv;
float *xrcva, *zrcva;
char* rcv_txt;
FILE *fp;
if (!getparint("verbose", &verbose)) verbose = 0;
/* Calculate Model Dimensions */
sub_x1=sub_x0+(nx-1)*dx;
sub_z1=sub_z0+(nz-1)*dz;
/* Compute how many receivers are defined and then allocate the receiver arrays */
/* Receivers on a Circle */
if (getparfloat("rrcv",&rrcv)) {
if (!getparfloat("dphi",&dphi)) dphi=2.0;
ncrcv=NINT(360.0/dphi);
if (verbose) vmess("Total number of receivers on a circle: %d",ncrcv);
}
else {
ncrcv=0;
}
/* Receivers from a File */
ntrcv=0;
if (!getparstring("rcv_txt",&rcv_txt)) rcv_txt=NULL;
if (rcv_txt!=NULL) {
/* Open text file */
fp=fopen(rcv_txt,"r");
assert(fp!=NULL);
/* Get number of lines */
while (!feof(fp)) if (fgetc(fp)=='\n') ntrcv++;
fseek(fp,-1,SEEK_CUR);
if (fgetc(fp)!='\n') ntrcv++; /* Checks if last line terminated by /n */
if (verbose) vmess("Number of receivers in rcv_txt file: %d",ntrcv);
rewind(fp); }
/* Receiver Array */
nxrcv=countparval("xrcva");
nzrcv=countparval("zrcva");
if (nxrcv!=nzrcv) verr("Number of receivers in array xrcva (%d), zrcva(%d) are not equal",nxrcv,nzrcv);
if (verbose&&nxrcv) vmess("Total number of array receivers: %d",nxrcv);
/* Linear Receiver Arrays */
Nx1 = countparval("xrcv1");
Nx2 = countparval("xrcv2");
Nz1 = countparval("zrcv1");
Nz2 = countparval("zrcv2");
if (Nx1!=Nx2) verr("Number of receivers starting points in 'xrcv1' (%d) and number of endpoint in 'xrcv2' (%d) are not equal",Nx1,Nx2);
if (Nz1!=Nz2) verr("Number of receivers starting points in 'zrcv1' (%d) and number of endpoint in 'zrcv2' (%d) are not equal",Nz1,Nz2);
if (Nx1!=Nz2) verr("Number of receivers starting points in 'xrcv1' (%d) and number of endpoint in 'zrcv2' (%d) are not equal",Nx1,Nz2);
rec->max_nrec=ncrcv+ntrcv+nxrcv;
/* no receivers are defined use default linear array of receivers on top of model */
if (!rec->max_nrec && Nx1==0) Nx1=1; // Default is to use top of model to record data
if (Nx1) {
/* Allocate Start & End Points of Linear Arrays */
xrcv1=(float *)malloc(Nx1*sizeof(float));
xrcv2=(float *)malloc(Nx1*sizeof(float));
zrcv1=(float *)malloc(Nx1*sizeof(float));
zrcv2=(float *)malloc(Nx1*sizeof(float));
if (!getparfloat("xrcv1",xrcv1)) xrcv1[0]=sub_x0;
if (!getparfloat("xrcv2",xrcv2)) xrcv2[0]=sub_x1;
if (!getparfloat("zrcv1",zrcv1)) zrcv1[0]=sub_z0;
if (!getparfloat("zrcv2",zrcv2)) zrcv2[0]=zrcv1[0];
/* check if receiver arrays fit into model */
for (iarray=0; iarray<Nx1; iarray++) {
xrcv1[iarray] = MAX(sub_x0, xrcv1[iarray]);
xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
xrcv2[iarray] = MAX(sub_x0, xrcv2[iarray]);
xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
zrcv1[iarray] = MAX(sub_z0, zrcv1[iarray]);
zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
zrcv2[iarray] = MAX(sub_z0, zrcv2[iarray]);
zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
}
/* Crop to Fit Model */
/* Max's addtion still have to check if it has the same fucntionality */
for (iarray=0;iarray<Nx1;iarray++) {
if (xrcv1[iarray]<sub_x0) {
if (xrcv2[iarray]<sub_x0) {
verr("Linear array %d outside model bounds",iarray);
}
else {
vwarn("Cropping element %d of 'xrcv1' (%f) to model bounds (%f)",iarray,xrcv1[iarray],sub_x0);
xrcv1[iarray]=sub_x0;
}
}
else if (xrcv1[iarray] > sub_x1) {
verr("Linear array %d outside model bounds",iarray);
}
if ( (xrcv2[iarray] < xrcv1[iarray]) ) {
verr("Ill defined linear array %d, 'xrcv1' (%f) greater than 'xrcv2' (%f)",iarray,xrcv1[iarray],xrcv2[iarray]);
}
else if (xrcv2[iarray]>sub_x1) {
vwarn("Cropping element %d of 'xrcv2' (%f) to model bounds (%f)",iarray,xrcv2[iarray],sub_x1);
xrcv2[iarray]=sub_x1;
}
if (zrcv1[iarray] < sub_z0) { if (zrcv2[iarray] < sub_z0) {
verr("Linear array %d outside model bounds",iarray);
}
else {
vwarn("Cropping element %d of 'zrcv1' (%f) to model bounds (%f)",iarray,zrcv1[iarray],sub_z0);
zrcv1[iarray]=sub_z0;
}
}
else if (zrcv1[iarray] > sub_z1) {
verr("Linear array %d outside model bounds",iarray);
}
if ( (zrcv2[iarray] < zrcv1[iarray]) ) {
verr("Ill defined linear array %d, 'zrcv1' (%f) greater than 'zrcv2' (%f)",iarray,zrcv1[iarray],zrcv2[iarray]);
}
else if (zrcv2[iarray]>sub_z1) {
vwarn("Cropping element %d of 'xrcv2' (%f) to model bounds (%f)",iarray,zrcv2[iarray],sub_z1);
zrcv2[iarray]=sub_z1;
}
}
/* Get Sampling Rates */
Ndx = countparval("dxrcv");
Ndz = countparval("dzrcv");
dxr = (float *)malloc(Nx1*sizeof(float));
dzr = (float *)malloc(Nx1*sizeof(float));
if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
for (i=1; i<Nx1; i++) {
dxr[i] = dxr[0];
dzr[i] = dzr[0];
}
Ndx=1;
Ndz=1;
}
else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
for (i=1; i<Nx1; i++) {
dxr[i] = dxr[0];
dzr[i] = dzr[0];
}
Ndz=1;
Ndx=1;
}
else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
if (Ndx!=Ndz) {
verr("Number of 'dxrcv' (%d) is not equal to number of 'dzrcv' (%d) or 1",Ndx,Ndz);
}
if (Ndx!=Nx1 && Ndx!=1) {
verr("Number of 'dxrcv' (%d) is not equal to number of starting points in 'xrcv1' (%d) or 1",Ndx,Nx1);
}
}
/* check consistency of receiver steps */
for (iarray=0; iarray<Ndx; iarray++) {
if (dxr[iarray]<0) {
dxr[i]=dx;
vwarn("'dxrcv' element %d (%f) is less than zero, changing it to %f'",iarray,dxr[iarray],dx);
}
}
for (iarray=0;iarray<Ndz;iarray++) {
if (dzr[iarray]<0) {
dzr[iarray]=dz;
vwarn("'dzrcv' element %d (%f) is less than zero, changing it to %f'",iarray,dzr[iarray],dz);
}
}
for (iarray=0;iarray<Ndx;iarray++){
if (dxr[iarray]==0 && dzr[iarray]==0) {
xrcv2[iarray]=xrcv1[iarray];
dxr[iarray]=1.; vwarn("'dxrcv' element %d & 'dzrcv' element 1 are both 0.",iarray+1);
vmess("Placing 1 receiver at (%d,%d)",xrcv1[iarray],zrcv1[iarray]);
}
}
for (iarray=0;iarray<Ndx;iarray++){
if (xrcv1[iarray]==xrcv2[iarray] && dxr[iarray]!=0) {
dxr[iarray]=0.;
vwarn("Linear array %d: 'xrcv1'='xrcv2' and 'dxrcv' is not 0. Setting 'dxrcv'=0",iarray+1);
}
}
for (iarray=0;iarray<Ndx;iarray++){
if (zrcv1[iarray]==zrcv2[iarray] && dzr[iarray]!=0.){
dzr[iarray]=0.;
vwarn("Linear array %d: 'zrcv1'='zrcv2' and 'dzrcv' is not 0. Setting 'dzrcv'=0",iarray+1);
}
}
/* Calculate Number of Receivers */
nrcv = 0;
nlrcv=(int *)malloc(Nx1*sizeof(int));
for (iarray=0; iarray<Nx1; iarray++) {
xrange = (xrcv2[iarray]-xrcv1[iarray]);
zrange = (zrcv2[iarray]-zrcv1[iarray]);
if (dxr[iarray] != 0.0) {
nlrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
}
else {
if (dzr[iarray] == 0) {
verr("For receiver array %d: receiver distance dzrcv is not given", iarray);
}
nlrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
}
nrcv+=nlrcv[iarray];
}
/* Calculate Number of Receivers */
/*
nlrcv=(int *)malloc(Nx1*sizeof(int));
if (!isnan(*xrcv1)) *nlrcv=MIN(NINT((*xrcv2-*xrcv1)/(*dxr)),NINT((*zrcv2-*zrcv1)/(*dzr)))+1;
else *nlrcv=0;
nrcv=*nlrcv;
if (verbose>4 && nlrcv[iarray]!=0) vmess("Linear receiver array 1 has final bounds: (X: %f -> %f,Z: %f ->
%f)",xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*(*dxr),zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*(*dzr));
if (Ndx>1) {
for (iarray=1;iarray<Nx1;iarray++) {
if (!isnan(xrcv1[iarray])) {
nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/dxr[iarray]),NINT((zrcv2[iarray]-zrcv1[iarray])/dzr[iarray]))+1;
}
else {
nlrcv[iarray]=0;
}
nrcv+=nlrcv[iarray];
if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %d has final bounds: (X: %f -> %f,Z: %f ->
%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]*dxr[iarray],zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]*dzr[iarray]);
}
}
else {
for (iarray=1;iarray<Nx1;iarray++) {
if (!isnan(xrcv1[iarray])) nlrcv[iarray]=MIN(NINT((xrcv2[iarray]-xrcv1[iarray])/(*dxr)),NINT((zrcv2[iarray]-zrcv1[iarray])/(*dzr)))+1;
else nlrcv[iarray]=0;
nrcv+=nlrcv[iarray];
if (verbose>4&&nlrcv[iarray]!=0) vmess("Linear receiver array %d has final bounds: (X: %f -> %f,Z: %f ->
%f)",iarray,xrcv1[iarray],xrcv1[iarray]+nlrcv[iarray]**dxr,zrcv1[iarray],zrcv1[iarray]+nlrcv[iarray]**dzr);
}
}
*/
if (verbose) vmess("Total number of linear array receivers: %d",nrcv);
if (!nrcv) {
free(xrcv1);
free(xrcv2); free(zrcv1);
free(zrcv2);
free(dxr);
free(dzr);
free(nlrcv);
}
rec->max_nrec+=nrcv;
}
else {
nrcv=0;
}
/* allocate the receiver arrays */
/* Total Number of Receivers */
if (verbose) vmess("Total number of receivers: %d",rec->max_nrec);
/* Allocate Arrays */
rec->x = (int *)calloc(rec->max_nrec,sizeof(int));
rec->z = (int *)calloc(rec->max_nrec,sizeof(int));
rec->xr = (float *)calloc(rec->max_nrec,sizeof(float));
rec->zr = (float *)calloc(rec->max_nrec,sizeof(float));
/* read in the receiver postions */
nrec=0;
/* Receivers on a Circle */
if (ncrcv) {
if (!getparfloat("oxrcv",&oxrcv)) oxrcv=0.0;
if (!getparfloat("ozrcv",&ozrcv)) ozrcv=0.0;
if (!getparfloat("arcv",&arcv)) {
arcv=rrcv;
for (ix=0; ix<ncrcv; ix++) {
rec->xr[ix] = oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI));
rec->zr[ix] = ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI));
rec->x[ix] = NINT(rec->xr[ix]/dx);
rec->z[ix] = NINT(rec->zr[ix]/dz);
//rec->x[ix] = NINT((oxrcv-sub_x0+rrcv*cos(((ix*dphi)/360.0)*(2.0*M_PI)))/dx);
//rec->z[ix] = NINT((ozrcv-sub_z0+arcv*sin(((ix*dphi)/360.0)*(2.0*M_PI)))/dz);
if (verbose>4) fprintf(stderr,"Receiver Circle: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
}
}
else { /* an ellipse */
/* simple numerical solution to find equidistant points on an ellipse */
nh = (ncrcv)*1000; /* should be fine enough for most configurations */
h = 2.0*M_PI/nh;
a = MAX(rrcv, arcv);
b = MIN(rrcv, arcv);
e = sqrt(a*a-b*b)/a;
//fprintf(stderr,"a=%f b=%f e=%f\n", a, b, e);
circ = 0.0;
for (ir=0; ir<nh; ir++) {
s = sin(ir*h);
circ += sqrt(1.0-e*e*s*s);
}
circ = a*h*circ;
//fprintf(stderr,"circ = %f circle=%f\n", circ, 2.0*M_PI*rrcv);
/* define distance between receivers on ellipse */
dr = circ/ncrcv;
ix = 0;
srun = 0.0;
if (arcv >= rrcv) phase=0.0;
else phase=0.5*M_PI;
for (ir=0; ir<nh; ir++) {
s = sin(ir*h);
srun += sqrt(1.0-e*e*s*s);
if (a*h*srun >= ix*dr ) {
xr = rrcv*cos(ir*h+phase);
zr = arcv*sin(ir*h+phase);
rec->xr[ix] = oxrcv-sub_x0+xr; rec->zr[ix] = ozrcv-sub_z0+zr;
rec->x[ix] = NINT(rec->xr[ix]/dx);
rec->z[ix] = NINT(rec->zr[ix]/dz);
if (verbose>4) fprintf(stderr,"Receiver Ellipse: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
ix++;
}
if (ix == ncrcv) break;
}
}
/* check if receivers fit into the model otherwise clip to edges */
for (ix=0; ix<ncrcv; ix++) {
rec->x[ix] = MIN(nx-1, MAX(rec->x[ix], 0));
rec->z[ix] = MIN(nz-1, MAX(rec->z[ix], 0));
}
nrec += ncrcv;
}
/* Receiver Text File */
if (ntrcv) {
/* Allocate arrays */
xrcva = (float *)malloc(nrcv*sizeof(float));
zrcva = (float *)malloc(nrcv*sizeof(float));
/* Read in receiver coordinates */
for (i=0;i<nrcv;i++) {
if (fscanf(fp,"%e %e\n",&xrcva[i],&zrcva[i])!=2) vmess("Receiver Text File: Can not parse coordinates on line %d.",i);
}
/* Close file */
fclose(fp);
/* Process coordinates */
for (ix=0; ix<nrcv; ix++) {
rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
if (verbose>4) vmess("Receiver Text Array: xrcv[%d]=%f zrcv=%f", ix, rec->xr[nrec+ix]+sub_x0, rec->zr[nrec+ix]+sub_z0);
}
free(xrcva);
free(zrcva);
nrec += ntrcv;
}
/* Receiver Array */
if (nxrcv != 0) {
/* receiver array is defined */
xrcva = (float *)malloc(nxrcv*sizeof(float));
zrcva = (float *)malloc(nxrcv*sizeof(float));
getparfloat("xrcva", xrcva);
getparfloat("zrcva", zrcva);
for (ix=0; ix<nxrcv; ix++) {
rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
rec->x[nrec+ix] = NINT((xrcva[ix]-sub_x0)/dx);
rec->z[nrec+ix] = NINT((zrcva[ix]-sub_z0)/dz);
if (verbose>4) fprintf(stderr,"Receiver Array: xrcv[%d]=%f zrcv=%f\n", ix, rec->xr[nrec+ix]+sub_x0, rec->zr[nrec+ix]+sub_z0);
}
free(xrcva);
free(zrcva);
nrec += nxrcv;
}
/* Linear Receiver Arrays */
if (nrcv!=0) {
xrcv1 = (float *)malloc(Nx1*sizeof(float));
xrcv2 = (float *)malloc(Nx1*sizeof(float));
zrcv1 = (float *)malloc(Nx1*sizeof(float));
zrcv2 = (float *)malloc(Nx1*sizeof(float));
if(!getparfloat("xrcv1", xrcv1)) xrcv1[0]=sub_x0; if(!getparfloat("xrcv2", xrcv2)) xrcv2[0]=(nx-1)*dx+sub_x0;
if(!getparfloat("zrcv1", zrcv1)) zrcv1[0]=sub_z0;
if(!getparfloat("zrcv2", zrcv2)) zrcv2[0]=zrcv1[0];
Ndx = countparval("dxrcv");
Ndz = countparval("dzrcv");
dxr = (float *)malloc(Nx1*sizeof(float));
dzr = (float *)malloc(Nx1*sizeof(float));
if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
if ( (Ndx<=1) && (Ndz==0) ){ /* default values are set */
for (i=1; i<Nx1; i++) {
dxr[i] = dxr[0];
dzr[i] = dzr[0];
}
Ndx=1;
}
else if ( (Ndz==1) && (Ndx==0) ){ /* default values are set */
for (i=1; i<Nx1; i++) {
dxr[i] = dxr[0];
dzr[i] = dzr[0];
}
Ndz=1;
}
else { /* make sure that each array has dzrcv or dxrcv defined for each line or receivers */
if (Ndx>1) assert(Ndx==Nx1);
if (Ndz>1) assert(Ndz==Nx1);
}
/*
if ( (Ndx!=0) && (Ndz!=0) ) {
vwarn("Both dzrcv and dxrcv are set: dxrcv value is used");
Ndz=0;
for (i=0; i<Nx1; i++) dzr[i] = 0.0;
}
*/
/* check if receiver arrays fit into model */
for (iarray=0; iarray<Nx1; iarray++) {
xrcv1[iarray] = MAX(sub_x0, xrcv1[iarray]);
xrcv1[iarray] = MIN(sub_x0+nx*dx,xrcv1[iarray]);
xrcv2[iarray] = MAX(sub_x0, xrcv2[iarray]);
xrcv2[iarray] = MIN(sub_x0+nx*dx,xrcv2[iarray]);
zrcv1[iarray] = MAX(sub_z0, zrcv1[iarray]);
zrcv1[iarray] = MIN(sub_z0+nz*dz,zrcv1[iarray]);
zrcv2[iarray] = MAX(sub_z0, zrcv2[iarray]);
zrcv2[iarray] = MIN(sub_z0+nz*dz,zrcv2[iarray]);
}
/* calculate receiver array and store into rec->x,z */
for (iarray=0; iarray<Nx1; iarray++) {
xrange = (xrcv2[iarray]-xrcv1[iarray]);
zrange = (zrcv2[iarray]-zrcv1[iarray]);
if (dxr[iarray] != 0.0) {
nrcv = nlrcv[iarray];
dxrcv=dxr[iarray];
dzrcv = zrange/(nrcv-1);
if (dzrcv != dzr[iarray]) {
vwarn("For receiver array %d: calculated dzrcv=%f given=%f", iarray, dzrcv, dzr[iarray]);
vwarn("The calculated receiver distance %f is used", dzrcv);
}
}
else {
if (dzr[iarray] == 0) {
verr("For receiver array %d: receiver distance dzrcv is not given", iarray);
}
nrcv = nlrcv[iarray];
dxrcv = xrange/(nrcv-1); dzrcv = dzr[iarray];
if (dxrcv != dxr[iarray]) {
vwarn("For receiver array %d: calculated dxrcv=%f given=%f", iarray, dxrcv, dxr[iarray]);
vwarn("The calculated receiver distance %f is used", dxrcv);
}
}
// calculate coordinates
for (ir=0; ir<nrcv; ir++) {
rec->xr[nrec]=xrcv1[iarray]-sub_x0+ir*dxrcv;
rec->zr[nrec]=zrcv1[iarray]-sub_z0+ir*dzrcv;
rec->x[nrec]=NINT((rec->xr[nrec])/dx);
rec->z[nrec]=NINT((rec->zr[nrec])/dz);
nrec++;
}
}
free(xrcv1);
free(xrcv2);
free(zrcv1);
free(zrcv2);
free(dxr);
free(dzr);
free(nlrcv);
}
rec->n=rec->max_nrec;
return 0;
}