diff --git a/raytime3d/getModelInfo3d.c b/raytime3d/getModelInfo3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..000200c5cf3bf2d917cfa892b07618aacc6d0d96
--- /dev/null
+++ b/raytime3d/getModelInfo3d.c
@@ -0,0 +1,125 @@
+#define _FILE_OFFSET_BITS 64
+#define _LARGEFILE_SOURCE
+#define _LARGEFILE64_SOURCE
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <math.h>
+#include "par.h"
+#include "segy.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* reads gridded model file to compute minimum and maximum values and sampling intervals
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int getModelInfo3d(char *file_name, int *n1, int *n2, int *n3, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, int *axis, int verbose)
+{
+ FILE *fp;
+ size_t nread, trace_sz;
+ off_t bytes, pos;
+ int ret, i, one_shot, ntraces, model, i2, i3;
+ float *trace, scl;
+ segy hdr, lasthdr;
+
+ if (file_name == NULL) return;
+
+ fp = fopen( file_name, "r" );
+ assert( fp != NULL);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ ret = fseeko( fp, 0, SEEK_END );
+ if (ret<0) perror("fseeko");
+ bytes = ftello( fp );
+ rewind(fp);
+ pos = bytes-hdr.ns*sizeof(float)-TRCBYTES;
+ ret = fseeko( fp, pos, SEEK_SET );
+ if (ret<0) perror("fseeko");
+ nread = fread( &lasthdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+
+ if (hdr.trid == TRID_DEPTH) *axis = 1; /* samples are z-axis */
+ else *axis = 0; /* sample direction respresents the x-axis */
+
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+
+ *n1 = hdr.ns;
+ *d1 = hdr.d1;
+ *d2 = hdr.d2;
+ *f1 = hdr.f1;
+ *f2 = hdr.gx*scl;
+ *f3 = hdr.gy*scl;
+
+ if ( NINT(100.0*((*d1)/(*d2)))!=100 ) {
+ verr("dx and dz are different in the model !");
+ }
+ if ( NINT(1000.0*(*d1))==0 ) {
+ if(!getparfloat("dx",d1)) {
+ verr("dx is equal to zero use parameter h= to set value");
+ }
+ *d2 = *d1;
+ }
+ trace_sz = sizeof(float)*(*n1)+TRCBYTES;
+ ntraces = (int) (bytes/trace_sz);
+
+ if (ntraces == 1) { /* 1D medium */
+ model = 1;
+ *n2 = 1;
+ *n3 = 1;
+ *d2 = *d1;
+ *d3 = *d1;
+ }
+ else { /* find out if this is a 2D or 3D model */
+ if (hdr.gy == lasthdr.gy) { /* assume 2D model */
+ *n3 = 1;
+ *n2 = ntraces;
+ *d3 = *d1;
+ }
+ else { /* 3D model */
+ /* find the number of traces in the x-direction */
+ rewind(fp);
+ one_shot = 1;
+ i3=0;
+ while (one_shot) {
+ i2=0;
+ lasthdr.gy = hdr.gy;
+ while (hdr.gy == lasthdr.gy) { /* number of samples in x */
+ pos = i2*trace_sz;
+ ret = fseeko( fp, pos, SEEK_SET );
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread==0) break;
+ i2++;
+ }
+ fprintf(stderr,"3D model gy=%d %d traces in x = %d\n", lasthdr.gy, i3, i2-1);
+ if (nread==0) break;
+ i3++;
+ }
+ *n3=i3;
+ *n2=ntraces/i3;
+ }
+ }
+
+ if (verbose>2) {
+ vmess("For file %s", file_name);
+ vmess("nz=%d nx=%d ny=%d ", *n1, *n2, *n3);
+ vmess("dz=%f dx=%f *dy=%f", *d1, *d2, *d3);
+ vmess("zstart=%f xstart=%f ystart=%f", *f1, *f2, *f3);
+ if (*axis) vmess("sample represent z-axis\n");
+ else vmess("sample represent x-axis\n");
+ }
+
+
+ return 0;
+}
+
diff --git a/raytime3d/getParameters3d.c b/raytime3d/getParameters3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..10522518777dc64deee8ade8ca0e0f0009936949
--- /dev/null
+++ b/raytime3d/getParameters3d.c
@@ -0,0 +1,347 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"par.h"
+#include"raytime3d.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+*
+* The routine getParameters reads in all parameters to set up a FD modeling.
+* Model and source parameters are used to calculate stability and dispersion relations
+* Source and receiver positions are calculated and checked if they fit into the model.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int getModelInfo3d(char *file_name, int *n1, int *n2, int *n3, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, int *axis, int verbose);
+
+int recvPar3d(recPar *rec, float sub_x0, float sub_z0, float sub_y0, float dx, float dz, float dy, int nx, int nz, int ny);
+
+int getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose)
+{
+ int nx, nz, ny, nsrc, ix, iy, axis, is0;
+ int n1, n2, n3;
+ int idzshot, idxshot, idyshot;
+ int src_ix0, src_iz0, src_iy0, src_ix1, src_iz1, src_iy1;
+ float cp_min, cp_max;
+ float sub_x0,sub_z0,sub_y0;
+ float srcendx, srcendz, srcendy, dy, dx, dz;
+ float xsrc, zsrc, ysrc, dxshot, dzshot, dyshot;
+ float dxrcv,dzrcv,dyrcv,dxspread,dzspread,dyspread;
+ float xmax, zmax, ymax;
+ float xsrc1, xsrc2, zsrc1, zsrc2, ysrc1, ysrc2;
+ float *xsrca, *zsrca, *ysrca;
+ float rsrc, oxsrc, ozsrc, oysrc, dphisrc, ncsrc;
+ size_t nsamp;
+ int nxsrc, nzsrc, nysrc;
+ int is;
+ char *src_positions, *file_cp=NULL;
+
+ if (!getparint("verbose",&verbose)) verbose=0;
+
+ if (!getparstring("file_cp", &file_cp)) {
+ vwarn("file_cp not defined, assuming homogeneous model");
+ }
+
+ if (!getparstring("file_rcv",&rec->file_rcv)) rec->file_rcv="recv.su";
+ if (!getparint("src_at_rcv",&src->src_at_rcv)) src->src_at_rcv=1;
+
+ /* read model parameters, which are used to set up source and receivers and check stability */
+
+ getModelInfo3d(mod->file_cp, &n1, &n2, &n3, &dz, &dx, &dy, &sub_z0, &sub_x0, &sub_y0, &axis, verbose);
+ mod->dz = dz;
+ mod->dx = dx;
+ mod->dy = dx;
+ mod->nz = nz;
+ mod->nx = nx;
+ mod->ny = nx;
+
+ /* origin of model in real (non-grid) coordinates */
+ mod->x0 = sub_x0;
+ mod->z0 = sub_z0;
+ mod->y0 = sub_y0;
+ xmax = sub_x0+(nx-1)*dx;
+ zmax = sub_z0+(nz-1)*dz;
+ ymax = sub_y0+(ny-1)*dy;
+
+ if (verbose) {
+ vmess("*******************************************");
+ vmess("*************** model info ****************");
+ vmess("*******************************************");
+ vmess("nz = %8d nx = %8d ny = %8d", nz, nx, ny);
+ vmess("dz = %8.4f dx = %8.4f dy = %8.4f", dz, dx, dy);
+ vmess("zmin = %8.4f zmax = %8.4f", sub_z0, zmax);
+ vmess("xmin = %8.4f xmax = %8.4f", sub_x0, xmax);
+ vmess("ymin = %8.4f ymax = %8.4f", sub_y0, ymax);
+ }
+
+ /* define the number and type of shots to model */
+ /* each shot can have multiple sources arranged in different ways */
+
+ if (!getparfloat("ysrc",&ysrc)) ysrc=sub_y0+((ny-1)*dy)/2.0;
+ if (!getparfloat("xsrc",&xsrc)) xsrc=sub_x0+((nx-1)*dx)/2.0;
+ if (!getparfloat("zsrc",&zsrc)) zsrc=sub_z0;
+ if (!getparint("nyshot",&shot->ny)) shot->ny=1;
+ if (!getparint("nxshot",&shot->nx)) shot->nx=1;
+ if (!getparint("nzshot",&shot->nz)) shot->nz=1;
+ if (!getparfloat("dyshot",&dyshot)) dyshot=dy;
+ if (!getparfloat("dxshot",&dxshot)) dxshot=dx;
+ if (!getparfloat("dzshot",&dzshot)) dzshot=dz;
+
+ shot->n = (shot->nx)*(shot->nz)*(shot->ny);
+
+ if (shot->nx>1) {
+ idxshot=MAX(0,NINT(dxshot/dx));
+ }
+ else {
+ idxshot=0.0;
+ }
+ if (shot->nz>1) {
+ idzshot=MAX(0,NINT(dzshot/dz));
+ }
+ else {
+ idzshot=0.0;
+ }
+ if (shot->ny>1) {
+ idyshot=MAX(0,NINT(dyshot/dy));
+ }
+ else {
+ idyshot=0.0;
+ }
+
+ /* calculate the shot positions */
+
+ src_iy0=MAX(0,NINT((ysrc-sub_y0)/dy));
+ src_iy0=MIN(src_iy0,ny);
+ src_ix0=MAX(0,NINT((xsrc-sub_x0)/dx));
+ src_ix0=MIN(src_ix0,nx);
+ src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
+ src_iz0=MIN(src_iz0,nz);
+ srcendy=(shot->ny-1)*dyshot+ysrc;
+ srcendx=(shot->nx-1)*dxshot+xsrc;
+ srcendz=(shot->nz-1)*dzshot+zsrc;
+ src_iy1=MAX(0,NINT((srcendy-sub_y0)/dy));
+ src_iy1=MIN(src_iy1,ny);
+ src_ix1=MAX(0,NINT((srcendx-sub_x0)/dx));
+ src_ix1=MIN(src_ix1,nx);
+ src_iz1=MAX(0,NINT((srcendz-sub_z0)/dz));
+ src_iz1=MIN(src_iz1,nz);
+
+ shot->y = (int *)calloc(shot->ny,sizeof(int));
+ shot->x = (int *)calloc(shot->nx,sizeof(int));
+ shot->z = (int *)calloc(shot->nz,sizeof(int));
+ for (is=0; is<shot->ny; is++) {
+ shot->y[is] = src_iy0+is*idyshot;
+ if (shot->y[is] > ny-1) shot->ny = is-1;
+ }
+ for (is=0; is<shot->nx; is++) {
+ shot->x[is] = src_ix0+is*idxshot;
+ if (shot->x[is] > nx-1) shot->nx = is-1;
+ }
+ for (is=0; is<shot->nz; is++) {
+ shot->z[is] = src_iz0+is*idzshot;
+ if (shot->z[is] > nz-1) shot->nz = is-1;
+ }
+
+ /* check if source array is defined */
+
+ nysrc = countparval("xyrca");
+ nxsrc = countparval("xsrca");
+ nzsrc = countparval("zsrca");
+ if (nxsrc != nzsrc || nxsrc != nysrc) {
+ verr("Number of sources in array xsrca (%d), zsrca(%d) are not equal",nxsrc, nzsrc);
+ }
+
+ /* check if sources on a circle are defined */
+
+ if (getparfloat("rsrc", &rsrc)) {
+ if (!getparfloat("dphisrc",&dphisrc)) dphisrc=2.0;
+ if (!getparfloat("oysrc",&oysrc)) oysrc=0.0;
+ if (!getparfloat("oxsrc",&oxsrc)) oxsrc=0.0;
+ if (!getparfloat("ozsrc",&ozsrc)) ozsrc=0.0;
+ ncsrc = NINT(360.0/dphisrc);
+ src->n = nsrc;
+
+ src->y = (int *)malloc(ncsrc*sizeof(int));
+ src->x = (int *)malloc(ncsrc*sizeof(int));
+ src->z = (int *)malloc(ncsrc*sizeof(int));
+
+ for (ix=0; ix<ncsrc; ix++) {
+ src->y[ix] = NINT((oysrc-sub_y0+rsrc*cos(((iy*dphisrc)/360.0)*(2.0*M_PI)))/dy);
+ src->x[ix] = NINT((oxsrc-sub_x0+rsrc*cos(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dx);
+ src->z[ix] = NINT((ozsrc-sub_z0+rsrc*sin(((ix*dphisrc)/360.0)*(2.0*M_PI)))/dz);
+ if (verbose>4) fprintf(stderr,"Source on Circle: ysrc[%d]=%d xsrc[%d]=%d zsrc=%d\n", iy, src->y[iy], ix, src->x[ix], src->z[ix]);
+ }
+
+ }
+
+ /* TO DO propagate src_positions parameter and structure through code */
+
+ if (!getparstring("src_positions",&src_positions)) src_positions="single";
+ src->random=0;
+ src->plane=0;
+ src->array=0;
+ src->single=0;
+ if (strstr(src_positions, "single")) src->single=1;
+ else if (strstr(src_positions, "array")) src->array=1;
+ else if (strstr(src_positions, "random")) src->random=1;
+ else if (strstr(src_positions, "plane")) src->plane=1;
+ else src->single=1;
+
+ /* to maintain functionality of older parameters usage */
+ if (!getparint("src_random",&src->random)) src->random=0;
+ if (!getparint("plane_wave",&src->plane)) src->plane=0;
+
+ if (src->random) {
+ src->plane=0;
+ src->array=0;
+ src->single=0;
+ }
+ if (src->plane) {
+ src->random=0;
+ src->array=0;
+ src->single=0;
+ }
+
+
+ /* number of sources per shot modeling */
+
+ if (!getparint("src_window",&src->window)) src->window=0;
+ if (!getparint("distribution",&src->distribution)) src->distribution=0;
+ if (!getparfloat("amplitude", &src->amplitude)) src->amplitude=0.0;
+ if (src->random && nxsrc==0) {
+ if (!getparint("nsrc",&nsrc)) nsrc=1;
+ if (!getparfloat("ysrc1", &ysrc1)) ysrc1=sub_y0;
+ if (!getparfloat("ysrc2", &ysrc2)) ysrc2=ymax;
+ if (!getparfloat("xsrc1", &xsrc1)) xsrc1=sub_x0;
+ if (!getparfloat("xsrc2", &xsrc2)) xsrc2=xmax;
+ if (!getparfloat("zsrc1", &zsrc1)) zsrc1=sub_z0;
+ if (!getparfloat("zsrc2", &zsrc2)) zsrc2=zmax;
+ dyshot = ysrc2-ysrc1;
+ dxshot = xsrc2-xsrc1;
+ dzshot = zsrc2-zsrc1;
+ src->y = (int *)malloc(nsrc*sizeof(int));
+ src->x = (int *)malloc(nsrc*sizeof(int));
+ src->z = (int *)malloc(nsrc*sizeof(int));
+ nsamp = 0;
+
+ }
+ else if (nxsrc != 0) {
+ /* source array is defined */
+ nsrc=nxsrc;
+ src->y = (int *)malloc(nsrc*sizeof(int));
+ src->x = (int *)malloc(nsrc*sizeof(int));
+ src->z = (int *)malloc(nsrc*sizeof(int));
+ ysrca = (float *)malloc(nsrc*sizeof(float));
+ xsrca = (float *)malloc(nsrc*sizeof(float));
+ zsrca = (float *)malloc(nsrc*sizeof(float));
+ getparfloat("ysrca", ysrca);
+ getparfloat("xsrca", xsrca);
+ getparfloat("zsrca", zsrca);
+ for (is=0; is<nsrc; is++) {
+ src->y[is] = NINT((ysrca[is]-sub_y0)/dy);
+ src->x[is] = NINT((xsrca[is]-sub_x0)/dx);
+ src->z[is] = NINT((zsrca[is]-sub_z0)/dz);
+ if (verbose>3) fprintf(stderr,"Source Array: ysrc[%d]=%f xsrc=%f zsrc=%f\n", is, ysrca[is], xsrca[is], zsrca[is]);
+ }
+ src->random = 1;
+ free(ysrca);
+ free(xsrca);
+ free(zsrca);
+ }
+ else {
+ if (src->plane) { if (!getparint("nsrc",&nsrc)) nsrc=1;}
+ else nsrc=1;
+
+ if (nsrc > nx) {
+ vwarn("Number of sources used in plane wave is larger than ");
+ vwarn("number of gridpoints in X. Plane wave will be clipped to the edges of the model");
+ nsrc = mod->nx;
+ }
+
+ /* for a source defined on mutliple gridpoint calculate p delay factor */
+
+ src->y = (int *)malloc(nsrc*sizeof(int));
+ src->x = (int *)malloc(nsrc*sizeof(int));
+ src->z = (int *)malloc(nsrc*sizeof(int));
+ is0 = -1*floor((nsrc-1)/2);
+ for (is=0; is<nsrc; is++) {
+ src->y[is] = is0 + is;
+ src->x[is] = is0 + is;
+ src->z[is] = 0;
+ }
+
+ }
+
+ src->n=nsrc;
+
+ if (verbose) {
+ if (src->n>1) {
+ vmess("*******************************************");
+ vmess("*********** source array info *************");
+ vmess("*******************************************");
+ vmess("Areal source array is defined with %d sources.",nsrc);
+ vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
+ }
+ if (src->random) {
+ vmess("Sources are placed at random locations in domain: ");
+ vmess(" x[%.2f : %.2f] z[%.2f : %.2f] ", xsrc1, xsrc2, zsrc1, zsrc2);
+ }
+ }
+
+ /* define receivers */
+
+ if (!getparint("sinkdepth",&rec->sinkdepth)) rec->sinkdepth=0;
+ if (!getparint("sinkdepth_src",&src->sinkdepth)) src->sinkdepth=0;
+ if (!getparint("sinkvel",&rec->sinkvel)) rec->sinkvel=0;
+ if (!getparint("max_nrec",&rec->max_nrec)) rec->max_nrec=15000;
+ if (!getparfloat("dyspread",&dyspread)) dyspread=0;
+ if (!getparfloat("dxspread",&dxspread)) dxspread=0;
+ if (!getparfloat("dzspread",&dzspread)) dzspread=0;
+ if (!getparint("nt",&rec->nt)) rec->nt=1024;
+
+ /* calculates the receiver coordinates */
+
+ recvPar3d(rec, sub_x0, sub_z0, sub_y0, dx, dz, dy, nx, nz, ny);
+
+ if (verbose) {
+ if (rec->n) {
+ dyrcv = rec->yr[MIN(1,rec->n-1)]-rec->yr[0];
+ dxrcv = rec->xr[MIN(1,rec->n-1)]-rec->xr[0];
+ dzrcv = rec->zr[MIN(1,rec->n-1)]-rec->zr[0];
+ vmess("*******************************************");
+ vmess("************* receiver info ***************");
+ vmess("*******************************************");
+ vmess("ntrcv = %d nrcv = %d ", rec->nt, rec->n);
+ vmess("dzrcv = %f dxrcv = %f dyrcv = %f ", dzrcv, dxrcv, dyrcv);
+ vmess("Receiver array at coordinates: ");
+ vmess("zmin = %f zmax = %f ", rec->zr[0]+sub_z0, rec->zr[rec->n-1]+sub_z0);
+ vmess("xmin = %f xmax = %f ", rec->xr[0]+sub_x0, rec->xr[rec->n-1]+sub_x0);
+ vmess("ymin = %f ymax = %f ", rec->yr[0]+sub_y0, rec->yr[rec->n-1]+sub_y0);
+ vmess("which are gridpoints: ");
+ vmess("izmin = %d izmax = %d ", rec->z[0], rec->z[rec->n-1]);
+ vmess("ixmin = %d ixmax = %d ", rec->x[0], rec->x[rec->n-1]);
+ vmess("iymin = %d iymax = %d ", rec->y[0], rec->y[rec->n-1]);
+ fprintf(stderr,"\n");
+ }
+ else {
+ vmess("*************** no receivers **************");
+ }
+ }
+
+ /* Ray tracing parameters */
+ if (!getparint("smoothwindow",&ray->smoothwindow)) ray->smoothwindow=0;
+ if (!getparint("useT2",&ray->useT2)) ray->useT2=0;
+ if (!getparint("geomspread",&ray->geomspread)) ray->geomspread=1;
+ if (!getparint("nraystep",&ray->nray)) ray->nray=5;
+
+ return 0;
+}
+
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
index 73e57f41822da6bec2322f85b0b943e25afad73f..b1735d85dc72268229f9e5d3e9c0a93329e7b8c2 100644
--- a/raytime3d/raytime3d.c
+++ b/raytime3d/raytime3d.c
@@ -146,46 +146,44 @@ void main(int argc, char *argv[])
getParameters3d(&mod, &rec, &src, &shot, &ray, verbose);
-
/*---------------------------------------------------------------------------*
* Open velocity file
*---------------------------------------------------------------------------*/
- if (file_cp != NULL) {
+ if (mod.file_cp != NULL) {
if (n2==1) { /* 1D model */
if(!getparint("nx",&nx)) verr("for 1D medium nx not defined");
if(!getparint("ny",&nx)) verr("for 1D medium ny not defined");
nz = n1;
oz = f1; ox = ((nx-1)/2)*d1; oy = ((ny-1)/2)*d1;
- dz = d1; dx = d1; dy = d1;
}
else if (n3==1) { /* 2D model */
if(!getparint("ny",&nx)) verr("for 2D medium ny not defined");
nz = n1; nx = n2;
oz = f1; ox = f2; oy = ((ny-1)/2)*d1;
- dz = d1; dx = d1; dy = d1;
}
else { /* Full 3D model */
nz = n1; nx = n2; nz = n3;
oz = f1; ox = f2; oy = f3;
- dz = d1; dx = d1; dy = d1;
}
- h = d1;
+ h = mod.dx;
slow0 = (float *)malloc(nz*nx*ny*sizeof(float));
if (slow0 == NULL) verr("Out of memory for slow0 array!");
- readModel3d(file_cp, slow0, n1, n2, n3, nz, nx, ny, h, verbose);
+ readModel3d(mod.file_cp, slow0, n1, n2, n3, nz, nx, ny, h, verbose);
if (verbose) vmess("h = %.2f nx = %d nz = %d ny = %d", h, nx, nz, ny);
}
else {
- nxy = nx * ny;
+ if(!getparfloat("c",&c)) verr("c not defined");
+ if(!getparfloat("h",&h)) verr("h not defined");
if(!getparint("nx",&nx)) verr("for homogenoeus medium nx not defined");
if(!getparint("ny",&nx)) verr("for homogenoeus medium ny not defined");
if(!getparint("nz",&nx)) verr("for homogenoeus medium nz not defined");
+ nxy = nx * ny;
oz = 0; ox = ((nx-1)/2)*d1; oy = ((ny-1)/2)*d1;
slow0 = (float *)malloc(nx*nz*ny*sizeof(float));
@@ -264,7 +262,6 @@ void main(int argc, char *argv[])
for (i = 0; i < ny; i++) {
hdrs[i].scalco = -1000;
hdrs[i].scalel = -1000;
-/* hdrs[i].offset = xi[0]*dx + is*ispr*dx - xsrc;*/
hdrs[i].sx = (int)(ox+xs*h)*1000;
hdrs[i].sy = (int)(oy+ys*h)*1000;
hdrs[i].gy = (int)(oy+i*d2)*1000;
diff --git a/raytime3d/raytime3d.h b/raytime3d/raytime3d.h
new file mode 100644
index 0000000000000000000000000000000000000000..9900d8b13315c9b5ec7662a053d8eb25973e7915
--- /dev/null
+++ b/raytime3d/raytime3d.h
@@ -0,0 +1,144 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<limits.h>
+#include<float.h>
+#include<math.h>
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+typedef struct _dcomplexStruct { /* complex number */
+ double r,i;
+} dcomplex;
+#endif/* complex */
+
+
+typedef struct _icoord { /* 3D coordinate integer */
+ int z;
+ int x;
+ int y;
+} icoord;
+
+typedef struct _fcoord { /* 3D coordinate float */
+ float z;
+ float x;
+ float y;
+} fcoord;
+
+struct s_ecount {
+ int corner,corner_min,side;
+};
+
+typedef struct _receiverPar { /* Receiver Parameters */
+ char *file_rcv;
+ int n;
+ int nt;
+ int max_nrec;
+ int *z;
+ int *x;
+ int *y;
+ float *zr;
+ float *xr;
+ float *yr;
+ int scale;
+ int sinkdepth;
+ int sinkvel;
+ float cp;
+ float rho;
+} recPar;
+
+typedef struct _modelPar { /* Model Parameters */
+ int sh;
+ char *file_cp;
+ float dz;
+ float dx;
+ float dy;
+ float dt;
+ float z0;
+ float x0;
+ float y0;
+ /* medium max/min values */
+ float cp_min;
+ float cp_max;
+ int nz;
+ int nx;
+ int ny;
+} modPar;
+
+typedef struct _waveletPar { /* Wavelet Parameters */
+ char *file_src; /* general source */
+ int nsrcf;
+ int nt;
+ int ns;
+ int nx;
+ int ny;
+ float dt;
+ float ds;
+ float fmax;
+ int random;
+ int seed;
+ int nst;
+ size_t *nsamp;
+} wavPar;
+
+typedef struct _sourcePar { /* Source Array Parameters */
+ int n;
+ int type;
+ int orient;
+ int *z;
+ int *x;
+ int *y;
+ int single;
+ int plane;
+ int circle;
+ int array;
+ int random;
+ int multiwav;
+ float angle;
+ float velo;
+ float amplitude;
+ int distribution;
+ int window;
+ int injectionrate;
+ int sinkdepth;
+ int src_at_rcv; /* Indicates that wavefield should be injected at receivers */
+} srcPar;
+
+typedef struct _shotPar { /* Shot Parameters */
+ int n;
+ int ny;
+ int nx;
+ int nz;
+ int *z;
+ int *x;
+ int *y;
+} shotPar;
+
+typedef struct _raypar { /* ray-tracing parameters */
+ int smoothwindow;
+ int useT2;
+ int geomspread;
+ int nray;
+} rayPar;
+
+#ifndef TRUE
+# define TRUE 1
+#endif
+
+#ifndef FALSE
+# define FALSE 0
+#endif
+
+#define equal(x,y) !strcmp(x,y)
+#define min2(a,b) (((a) < (b)) ? (a) : (b))
+#define max2(a,b) (((a) > (b)) ? (a) : (b))
+
+#define Infinity FLT_MAX
+
+#if __STDC_VERSION__ >= 199901L
+ /* "restrict" is a keyword */
+#else
+#define restrict
+#endif
+
diff --git a/raytime3d/readModel3d.c b/raytime3d/readModel3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..e8f4f9a98148717742210409df1f2d6358982f77
--- /dev/null
+++ b/raytime3d/readModel3d.c
@@ -0,0 +1,87 @@
+#define _FILE_OFFSET_BITS 64
+#define _LARGEFILE_SOURCE
+#define _LARGEFILE64_SOURCE
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <math.h>
+#include "segy.h"
+#include "par.h"
+#include "raytime3d.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* Reads gridded model files and compute from them medium parameters used in the FD kernels.
+* The files read in contain the P (and S) wave velocity and density.
+* The medium parameters calculated are lambda, mu, lambda+2mu, and 1/ro.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int readModel3d(char *file_name, float *slowness, int n1, int n2, int n3, int nz, int nx, int ny, float h, int verbose)
+{
+ FILE *fpcp;
+ size_t nread;
+ int i, j, k, tracesToDo;
+ int nxy, nxz;
+ float *tmp;
+ segy hdr;
+
+ nxy = nx * ny;
+ tmp = (float *)malloc(n1*sizeof(float));
+
+/* open files and read first header */
+
+ fpcp = fopen( file_name, "r" );
+ assert( fpcp != NULL);
+ nread = fread(&hdr, 1, TRCBYTES, fpcp);
+ assert(nread == TRCBYTES);
+ assert(hdr.ns == n1);
+
+ if (n2==1) { /* 1D model */
+ nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
+ for (j = 0; j < nz; j++) {
+ for (k = 0; k < ny; k++) {
+ for (i = 0; i < nx; i++) {
+ slowness[j*nxy+k*nx+i] = h/tmp[j];
+ }
+ }
+ }
+ }
+ else if (n3==1) { /* 2D model */
+ for (i = 0; i < nx; i++) {
+ nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
+ for (j = 0; j < nz; j++) {
+ for (k = 0; k < ny; k++) {
+ slowness[j*nxy+k*nx+i] = h/tmp[j];
+ }
+ }
+ nread = fread(&hdr, 1, TRCBYTES, fpcp);
+ }
+ }
+ else { /* Full 3D model */
+ /* read all traces */
+ for (k = 0; k < ny; k++) {
+ for (i = 0; i < nx; i++) {
+ nread = fread(&tmp[0], sizeof(float), hdr.ns, fpcp);
+ for (j = 0; j < nz; j++) {
+ slowness[j*nxy+k*nx+i] = h/tmp[j];
+ }
+ nread = fread(&hdr, 1, TRCBYTES, fpcp);
+ }
+ }
+ }
+
+ fclose(fpcp);
+
+ return 0;
+}
+
+
diff --git a/raytime3d/recvPar3d.c b/raytime3d/recvPar3d.c
new file mode 100644
index 0000000000000000000000000000000000000000..65dd258e5c95c0086fbdbf247ecfd084b22647c8
--- /dev/null
+++ b/raytime3d/recvPar3d.c
@@ -0,0 +1,616 @@
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+
+#include "raytime3d.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 recvPar3d(recPar *rec, float sub_x0, float sub_z0, float sub_y0, float dx, float dz, float dy, int nx, int nz, int ny)
+{
+ float *yrcv1, *yrcv2, *xrcv1, *xrcv2, *zrcv1, *zrcv2;
+ int i, ix, iy, ir, verbose;
+ float dxrcv, dzrcv, dyrcv, *dxr, *dzr, *dyr;
+ float rrcv, dphi, oxrcv, ozrcv, oyrcv, arcv;
+ double circ, h, a, b, e, s, xr, zr, yr, dr, srun, phase;
+ float xrange, zrange, yrange, sub_x1, sub_z1, sub_y1;
+ int Nx1, Nx2, Nz1, Nz2, Ny1, Ny2, Ndx, Ndz, Ndy, iarray, nrec, nh;
+ int nxrcv, nzrcv, nyrcv, ncrcv, nrcv, ntrcv, *nlrcv;
+ float *xrcva, *zrcva, *yrcva;
+ 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;
+ sub_y1=sub_y0+(ny-1)*dy;
+
+/* 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 */
+ nyrcv=countparval("yrcva");
+ 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 (nxrcv!=nyrcv) verr("Number of receivers in array xrcva (%d), yrcva(%d) are not equal",nxrcv,nyrcv);
+ if (verbose&&nxrcv) vmess("Total number of array receivers: %d",nxrcv);
+
+ /* Linear Receiver Arrays */
+ Ny1 = countparval("yrcv1");
+ Ny2 = countparval("yrcv2");
+ Nx1 = countparval("xrcv1");
+ Nx2 = countparval("xrcv2");
+ Nz1 = countparval("zrcv1");
+ Nz2 = countparval("zrcv2");
+ if (Ny1!=Ny2) verr("Number of receivers starting points in 'yrcv1' (%d) and number of endpoint in 'yrcv2' (%d) are not equal",Ny1,Ny2);
+ 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);
+ if (Nx1!=Ny2) verr("Number of receivers starting points in 'xrcv1' (%d) and number of endpoint in 'yrcv2' (%d) are not equal",Nx1,Ny2);
+
+ 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 */
+ yrcv1=(float *)malloc(Nx1*sizeof(float));
+ yrcv2=(float *)malloc(Nx1*sizeof(float));
+ 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("yrcv1",yrcv1)) yrcv1[0]=sub_y0;
+ if (!getparfloat("yrcv2",yrcv2)) yrcv2[0]=sub_y1;
+ 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++) {
+ yrcv1[iarray] = MAX(sub_y0, yrcv1[iarray]);
+ yrcv1[iarray] = MIN(sub_y0+ny*dy,yrcv1[iarray]);
+ yrcv2[iarray] = MAX(sub_y0, yrcv2[iarray]);
+ yrcv2[iarray] = MIN(sub_y0+ny*dy,yrcv2[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 (yrcv1[iarray]<sub_y0) {
+ if (yrcv2[iarray]<sub_y0) {
+ verr("Linear array %d outside model bounds",iarray);
+ }
+ else {
+ vwarn("Cropping element %d of 'yrcv1' (%f) to model bounds (%f)",iarray,yrcv1[iarray],sub_y0);
+ yrcv1[iarray]=sub_y0;
+ }
+ }
+ else if (yrcv1[iarray] > sub_y1) {
+ verr("Linear array %d outside model bounds",iarray);
+ }
+ if ( (yrcv2[iarray] < yrcv1[iarray]) ) {
+ verr("Ill defined linear array %d, 'yrcv1' (%f) greater than 'yrcv2' (%f)",iarray,yrcv1[iarray],yrcv2[iarray]);
+ }
+ else if (yrcv2[iarray]>sub_y1) {
+ vwarn("Cropping element %d of 'yrcv2' (%f) to model bounds (%f)",iarray,yrcv2[iarray],sub_y1);
+ yrcv2[iarray]=sub_y1;
+ }
+
+ 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 */
+ Ndy = countparval("dyrcv");
+ Ndx = countparval("dxrcv");
+ Ndz = countparval("dzrcv");
+
+ dyr = (float *)malloc(Nx1*sizeof(float));
+ dxr = (float *)malloc(Nx1*sizeof(float));
+ dzr = (float *)malloc(Nx1*sizeof(float));
+ if(!getparfloat("dyrcv", dyr)) dyr[0]=dy;
+ if(!getparfloat("dxrcv", dxr)) dxr[0]=dx;
+ if(!getparfloat("dzrcv", dzr)) dzr[0]=0.0;
+/* TODO logic extended to 3D */
+ 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;
+ }
+/* end TODO */
+ 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);
+ }
+ if (Ndy!=Ndz) {
+ verr("Number of 'dyrcv' (%d) is not equal to number of 'dzrcv' (%d) or 1",Ndy,Ndz);
+ }
+ }
+
+ /* check consistency of receiver steps */
+ for (iarray=0; iarray<Ndy; iarray++) {
+ if (dyr[iarray]<0) {
+ dyr[i]=dy;
+ vwarn("'dyrcv' element %d (%f) is less than zero, changing it to %f'",iarray,dyr[iarray],dy);
+ }
+ }
+ 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++) {
+ yrange = (yrcv2[iarray]-yrcv1[iarray]);
+ xrange = (xrcv2[iarray]-xrcv1[iarray]);
+ zrange = (zrcv2[iarray]-zrcv1[iarray]);
+ if (dyr[iarray] != 0.0) {
+ nlrcv[iarray] = NINT(fabs(yrange/dyr[iarray]))+1;
+ }
+ 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(yrcv1);
+ free(yrcv2);
+ free(xrcv1);
+ free(xrcv2);
+ free(zrcv1);
+ free(zrcv2);
+ free(dyr);
+ 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->y = (int *)calloc(rec->max_nrec,sizeof(int));
+ rec->x = (int *)calloc(rec->max_nrec,sizeof(int));
+ rec->z = (int *)calloc(rec->max_nrec,sizeof(int));
+ rec->yr = (float *)calloc(rec->max_nrec,sizeof(float));
+ 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("oyrcv",&oxrcv)) oyrcv=0.0;
+ 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->yr[ix] = oyrcv-sub_y0+rrcv*cos(((iy*dphi)/360.0)*(2.0*M_PI));
+ 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->y[ix] = NINT(rec->yr[ix]/dy);
+ 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: yrcv[%d]=%f xrcv=%f zrcv=%f\n", ix, rec->yr[ix]+sub_y0, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
+ }
+ }
+ else { /* an ellipse */
+/*TODO fix for 3D */
+ /* 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 ) {
+ yr = rrcv*cos(ir*h+phase);
+ xr = rrcv*cos(ir*h+phase);
+ zr = arcv*sin(ir*h+phase);
+ rec->yr[ix] = oyrcv-sub_y0+xr;
+ rec->xr[ix] = oxrcv-sub_x0+xr;
+ rec->zr[ix] = ozrcv-sub_z0+zr;
+ rec->y[ix] = NINT(rec->yr[ix]/dy);
+ 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->y[ix] = MIN(ny-1, MAX(rec->y[ix], 0));
+ 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 */
+ yrcva = (float *)malloc(nrcv*sizeof(float));
+ 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 %e\n",&yrcva[i],&xrcva[i],&zrcva[i])!=3) 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->yr[nrec+ix] = yrcva[ix]-sub_y0;
+ rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
+ rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
+ rec->y[nrec+ix] = NINT((yrcva[ix]-sub_y0)/dy);
+ 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 Text Array: yrcv[%d]=%f xrcv=%f zrcv=%f\n", ix, rec->yr[ix]+sub_y0, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
+ }
+ free(yrcva);
+ free(xrcva);
+ free(zrcva);
+ nrec += ntrcv;
+ }
+
+ /* Receiver Array */
+ if (nxrcv != 0) {
+ /* receiver array is defined */
+ yrcva = (float *)malloc(nxrcv*sizeof(float));
+ xrcva = (float *)malloc(nxrcv*sizeof(float));
+ zrcva = (float *)malloc(nxrcv*sizeof(float));
+ getparfloat("yrcva", yrcva);
+ getparfloat("xrcva", xrcva);
+ getparfloat("zrcva", zrcva);
+ for (ix=0; ix<nxrcv; ix++) {
+ rec->yr[nrec+ix] = yrcva[ix]-sub_y0;
+ rec->xr[nrec+ix] = xrcva[ix]-sub_x0;
+ rec->zr[nrec+ix] = zrcva[ix]-sub_z0;
+ rec->y[nrec+ix] = NINT((yrcva[ix]-sub_y0)/dy);
+ 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: yrcv[%d]=%f xrcv=%f zrcv=%f\n", ix, rec->yr[ix]+sub_y0, rec->xr[ix]+sub_x0, rec->zr[ix]+sub_z0);
+ }
+ free(yrcva);
+ free(xrcva);
+ free(zrcva);
+ nrec += nxrcv;
+ }
+
+ /* Linear Receiver Arrays */
+ if (nrcv!=0) {
+ yrcv1 = (float *)malloc(Nx1*sizeof(float));
+ yrcv2 = (float *)malloc(Nx1*sizeof(float));
+ 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("yrcv1", yrcv1)) yrcv1[0]=sub_y0;
+ if(!getparfloat("yrcv2", yrcv2)) yrcv2[0]=(ny-1)*dy+sub_y0;
+ 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];
+
+ Ndy = countparval("dyrcv");
+ Ndx = countparval("dxrcv");
+ Ndz = countparval("dzrcv");
+
+ dyr = (float *)malloc(Nx1*sizeof(float));
+ dxr = (float *)malloc(Nx1*sizeof(float));
+ dzr = (float *)malloc(Nx1*sizeof(float));
+ if(!getparfloat("dyrcv", dyr)) dyr[0]=dy;
+ 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++) {
+ dyr[i] = dyr[0];
+ 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++) {
+ dyr[i] = dyr[0];
+ 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 (Ndy>1) assert(Ndy==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++) {
+ yrcv1[iarray] = MAX(sub_y0, yrcv1[iarray]);
+ yrcv1[iarray] = MIN(sub_y0+ny*dy,yrcv1[iarray]);
+ yrcv2[iarray] = MAX(sub_y0, yrcv2[iarray]);
+ yrcv2[iarray] = MIN(sub_y0+ny*dy,yrcv2[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++) {
+ yrange = (yrcv2[iarray]-yrcv1[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->yr[nrec]=yrcv1[iarray]-sub_y0+ir*dyrcv;
+ rec->xr[nrec]=xrcv1[iarray]-sub_x0+ir*dxrcv;
+ rec->zr[nrec]=zrcv1[iarray]-sub_z0+ir*dzrcv;
+
+ rec->y[nrec]=NINT((rec->yr[nrec])/dy);
+ rec->x[nrec]=NINT((rec->xr[nrec])/dx);
+ rec->z[nrec]=NINT((rec->zr[nrec])/dz);
+ nrec++;
+ }
+ }
+ free(yrcv1);
+ free(yrcv2);
+ free(xrcv1);
+ free(xrcv2);
+ free(zrcv1);
+ free(zrcv2);
+ free(dyr);
+ free(dxr);
+ free(dzr);
+ free(nlrcv);
+ }
+
+ rec->n=rec->max_nrec;
+ return 0;
+}