diff --git a/marchenko/fmute.c b/marchenko/fmute.c
index 4d63faf38fa8bb4cb489ee554bca07e52ce9a439..f5c7668a16dc382b586c2981427888c6fa256f66 100644
--- a/marchenko/fmute.c
+++ b/marchenko/fmute.c
@@ -228,6 +228,8 @@ int main (int argc, char **argv)
/* alternative find maximum at source position */
dxrcv = (hdrs_in1[nx1-1].gx - hdrs_in1[0].gx)*sclsxgx/(float)(nx1-1);
imax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
+ /* make sure that the position fits into the receiver array */
+ imax = MIN(MAX(0,imax),nx1-1);
tmax=0.0;
jmax = 0;
xmax=0.0;
diff --git a/raytime/Makefile b/raytime/Makefile
index 192a17315ed0dc71f494646357315a4ffcc2dc82..b3aa6f39b0b3a26f00e0c892e352eebbbe223768 100644
--- a/raytime/Makefile
+++ b/raytime/Makefile
@@ -47,16 +47,16 @@ SRCC = $(PRG).c \
OBJC = $(SRCC:%.c=%.o)
$(PRG): $(OBJC) raytime.h
- $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o Raytime $(OBJC) $(LIBS)
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o raytime $(OBJC) $(LIBS)
install: raytime
- cp Raytime $B
+ cp raytime $B
clean:
- rm -f core $(OBJC) $(OBJM) Raytime
+ rm -f core $(OBJC) $(OBJM) raytime
realclean:
- rm -f core $(OBJC) $(OBJM) $(PRG) $B/Raytime
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/raytime
print: Makefile $(SRC)
diff --git a/raytime3d/Makefile b/raytime3d/Makefile
index 1bed95dcdef7cd3fb664ca9017a3bad0b38788b3..be0bb6e043b09ee22c89254f9c54afceb990ad92 100644
--- a/raytime3d/Makefile
+++ b/raytime3d/Makefile
@@ -27,15 +27,15 @@ PRG = raytime3d
SRCC = $(PRG).c \
vidale3d.c \
src3d.c \
- getParameters.c \
+ getParameters3d.c \
getWaveletInfo.c \
writeSrcRecPos.c \
- readModel.c \
+ readModel3d.c \
getWaveletHeaders.c \
verbosepkg.c \
- getModelInfo.c \
+ getModelInfo3d.c \
wallclock_time.c \
- recvPar.c \
+ recvPar3d.c \
writesufile.c \
name_ext.c \
atopkge.c \
@@ -45,17 +45,17 @@ SRCC = $(PRG).c \
OBJC = $(SRCC:%.c=%.o)
-$(PRG): $(OBJC) raytime.h
- $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o Raytime $(OBJC) $(LIBS)
+$(PRG): $(OBJC) raytime3d.h
+ $(CC) $(LDFLAGS) $(CFLAGS) $(OPTC) -o raytime3d $(OBJC) $(LIBS)
-install: raytime
- cp Raytime $B
+install: raytime3d
+ cp raytime3d $B
clean:
- rm -f core $(OBJC) $(OBJM) Raytime
+ rm -f core $(OBJC) $(OBJM) raytime3d
realclean:
- rm -f core $(OBJC) $(OBJM) $(PRG) $B/Raytime
+ rm -f core $(OBJC) $(OBJM) $(PRG) $B/raytime3d
print: Makefile $(SRC)
diff --git a/raytime3d/getModelInfo.c b/raytime3d/getModelInfo.c
deleted file mode 100644
index 378a1b50ebac46e5b7b8a8bef4b5365ac15bef9d..0000000000000000000000000000000000000000
--- a/raytime3d/getModelInfo.c
+++ /dev/null
@@ -1,109 +0,0 @@
-#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 getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose)
-{
- FILE *fp;
- size_t nread, trace_sz;
- off_t bytes;
- int ret, i, one_shot, ntraces;
- float *trace, cmin;
- segy hdr;
-
- 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 );
-
- *n1 = hdr.ns;
- *d1 = hdr.d1;
- *d2 = hdr.d2;
- *f1 = hdr.f1;
- *f2 = hdr.f2;
-
- 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 dx= to set value");
- }
- *d2 = *d1;
- }
- trace_sz = sizeof(float)*(*n1)+TRCBYTES;
- ntraces = (int) (bytes/trace_sz);
- *n2 = ntraces;
-
- /* check to find out min and max values gather */
-
- one_shot = 1;
- trace = (float *)malloc(trace_sz);
- fseeko( fp, TRCBYTES, SEEK_SET );
- nread = fread( trace, sizeof(float), hdr.ns, fp );
- assert (nread == hdr.ns);
- fseeko( fp, TRCBYTES, SEEK_SET );
-
- if (hdr.trid == TRID_DEPTH) *axis = 1; /* samples are z-axis */
- else *axis = 0; /* sample direction respresents the x-axis */
-
- i=0; cmin=trace[0];
- while ( ( (cmin==0.0) && zeroch) && (i<hdr.ns) ) cmin=trace[i++];
-
- *max = cmin;
- *min = cmin;
- /* keep on reading traces until there are no more traces (nread==0) */
- while (one_shot) {
- nread = fread( trace, sizeof(float), hdr.ns, fp );
- assert (nread == hdr.ns);
- for (i=0;i<(*n1);i++) {
- *max = MAX(trace[i],*max);
- cmin = MIN(trace[i],*min);
- if (zeroch) {
- if (cmin!=0.0) *min = MIN(*min, cmin);
- }
- else {
- *min = cmin;
- }
- }
- nread = fread( &hdr, 1, TRCBYTES, fp );
- if (nread==0) break;
- }
- fclose(fp);
- free(trace);
-
- if (verbose>2) {
- vmess("For file %s", file_name);
- vmess("nz=%d nx=%d", *n1, *n2);
- vmess("dz=%f dx=%f", *d1, *d2);
- vmess("min=%f max=%f", *min, *max);
- vmess("zstart=%f xstart=%f", *f1, *f2);
- if (*axis) vmess("sample represent z-axis\n");
- else vmess("sample represent x-axis\n");
- }
- return 0;
-}
-
diff --git a/raytime3d/getParameters.c b/raytime3d/getParameters.c
deleted file mode 100644
index 89d1d34585a977f9ba6fd58eff936bf4396e4ea2..0000000000000000000000000000000000000000
--- a/raytime3d/getParameters.c
+++ /dev/null
@@ -1,305 +0,0 @@
-#include<stdlib.h>
-#include<stdio.h>
-#include<math.h>
-#include<assert.h>
-#include"par.h"
-#include"raytime.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 getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose);
-
-int recvPar(recPar *rec, float sub_x0, float sub_z0, float dx, float dz, int nx, int nz);
-
-int getParameters(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose)
-{
- int nx, nz, nsrc, ix, axis, is0;
- int idzshot, idxshot;
- int src_ix0, src_iz0, src_ix1, src_iz1;
- float cp_min, cp_max;
- float sub_x0,sub_z0;
- float srcendx, srcendz, dx, dz;
- float xsrc, zsrc, dxshot, dzshot;
- float dxrcv,dzrcv,dxspread,dzspread;
- float xmax, zmax;
- float xsrc1, xsrc2, zsrc1, zsrc2;
- float *xsrca, *zsrca;
- float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
- size_t nsamp;
- int nxsrc, nzsrc;
- int is;
- char *src_positions;
-
- if (!getparint("verbose",&verbose)) verbose=0;
-
- if (!getparstring("file_cp",&mod->file_cp)) {
- verr("parameter file_cp required!");
- }
- 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 */
-
- getModelInfo(mod->file_cp, &nz, &nx, &dz, &dx, &sub_z0, &sub_x0, &cp_min, &cp_max, &axis, 1, verbose);
- mod->cp_max = cp_max;
- mod->cp_min = cp_min;
- mod->dz = dz;
- mod->dx = dx;
- mod->nz = nz;
- mod->nx = nx;
-
- /* origin of model in real (non-grid) coordinates */
- mod->x0 = sub_x0;
- mod->z0 = sub_z0;
- xmax = sub_x0+(nx-1)*dx;
- zmax = sub_z0+(nz-1)*dz;
-
- if (verbose) {
- vmess("*******************************************");
- vmess("*************** model info ****************");
- vmess("*******************************************");
- vmess("nz = %8d nx = %8d", nz, nx);
- vmess("dz = %8.4f dx = %8.4f", dz, dx);
- vmess("zmin = %8.4f zmax = %8.4f", sub_z0, zmax);
- vmess("xmin = %8.4f xmax = %8.4f", sub_x0, xmax);
- vmess("min(cp) = %9.3f max(cp) = %9.3f", cp_min, cp_max);
- }
-
- /* define the number and type of shots to model */
- /* each shot can have multiple sources arranged in different ways */
-
- if (!getparfloat("xsrc",&xsrc)) xsrc=sub_x0+((nx-1)*dx)/2.0;
- if (!getparfloat("zsrc",&zsrc)) zsrc=sub_z0;
- if (!getparint("nxshot",&shot->nx)) shot->nx=1;
- if (!getparint("nzshot",&shot->nz)) shot->nz=1;
- if (!getparfloat("dxshot",&dxshot)) dxshot=dx;
- if (!getparfloat("dzshot",&dzshot)) dzshot=dz;
-
- shot->n = (shot->nx)*(shot->nz);
-
- 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;
- }
-
- /* calculate the shot positions */
-
- 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);
- srcendx=(shot->nx-1)*dxshot+xsrc;
- srcendz=(shot->nz-1)*dzshot+zsrc;
- 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->x = (int *)calloc(shot->nx,sizeof(int));
- shot->z = (int *)calloc(shot->nz,sizeof(int));
- 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 */
-
- nxsrc = countparval("xsrca");
- nzsrc = countparval("zsrca");
- if (nxsrc != nzsrc) {
- 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("oxsrc",&oxsrc)) oxsrc=0.0;
- if (!getparfloat("ozsrc",&ozsrc)) ozsrc=0.0;
- ncsrc = NINT(360.0/dphisrc);
- src->n = nsrc;
-
- src->x = (int *)malloc(ncsrc*sizeof(int));
- src->z = (int *)malloc(ncsrc*sizeof(int));
-
- for (ix=0; ix<ncsrc; ix++) {
- 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: xsrc[%d]=%d zsrc=%d\n", 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("xsrc1", &xsrc1)) xsrc1=sub_x0;
- if (!getparfloat("xsrc2", &xsrc2)) xsrc2=xmax;
- if (!getparfloat("zsrc1", &zsrc1)) zsrc1=sub_z0;
- if (!getparfloat("zsrc2", &zsrc2)) zsrc2=zmax;
- dxshot = xsrc2-xsrc1;
- dzshot = zsrc2-zsrc1;
- 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->x = (int *)malloc(nsrc*sizeof(int));
- src->z = (int *)malloc(nsrc*sizeof(int));
- xsrca = (float *)malloc(nsrc*sizeof(float));
- zsrca = (float *)malloc(nsrc*sizeof(float));
- getparfloat("xsrca", xsrca);
- getparfloat("zsrca", zsrca);
- for (is=0; is<nsrc; is++) {
- 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: xsrc[%d]=%f zsrc=%f\n", is, xsrca[is], zsrca[is]);
- }
- src->random = 1;
- 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->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->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("dxspread",&dxspread)) dxspread=0;
- if (!getparfloat("dzspread",&dzspread)) dzspread=0;
- if (!getparint("nt",&rec->nt)) rec->nt=1024;
-
- /* calculates the receiver coordinates */
-
- recvPar(rec, sub_x0, sub_z0, dx, dz, nx, nz);
-
- if (verbose) {
- if (rec->n) {
- 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 ", dzrcv, dxrcv);
- 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("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]);
- 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/raytime.h b/raytime3d/raytime.h
deleted file mode 100644
index 882c78b9a2eb0b945737dac10ae5c6f07aa55123..0000000000000000000000000000000000000000
--- a/raytime3d/raytime.h
+++ /dev/null
@@ -1,135 +0,0 @@
-#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;
- float *zr;
- float *xr;
- 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 dt;
- float z0;
- float x0;
- /* medium max/min values */
- float cp_min;
- float cp_max;
- int nz;
- int nx;
-} modPar;
-
-typedef struct _waveletPar { /* Wavelet Parameters */
- char *file_src; /* general source */
- int nsrcf;
- int nt;
- int ns;
- int nx;
- 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 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 nx;
- int nz;
- int *z;
- int *x;
-} 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/raytime3d.c b/raytime3d/raytime3d.c
index 9b0c88d6e59fe63c13275cc139f7352bff266501..73e57f41822da6bec2322f85b0b943e25afad73f 100644
--- a/raytime3d/raytime3d.c
+++ b/raytime3d/raytime3d.c
@@ -1,15 +1,56 @@
-#include <DELPHI_IOc.h>
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
+#include <genfft.h>
+#include"par.h"
+#include"raytime3d.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))
+
+double wallclock_time(void);
+
+void name_ext(char *filename, char *extension);
+
+void threadAffinity(void);
+
+
+int getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose);
+
+int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord r, rayPar ray, fcoord *T, float *Jr);
+
+void applyMovingAverageFilter(float *slowness, icoord size, int window, int dim, float *averageModel);
+
+int readModel3d(char *file_name, float *slowness, int n1, int n2, int n3, int nz, int nx, int ny, float h, int verbose);
+
+int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int reverse, int verbose);
+
+int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
+
+void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int xs, int ys, int zs, int NCUBE);
+
+void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube);
+
/*********************** self documentation **********************/
char *sdoc[] = {
" ",
" RAYTIME3D - modeling of one-way traveltime for operators in 3D media",
" ",
-" raytime3d file_vel= xsrc1= zsrc1= ysrc1= [optional parameters]",
+" raytime3d file_cp= xsrc1= zsrc1= ysrc1= [optional parameters]",
" ",
" Required parameters:",
" ",
-" file_vel= ................ gridded velocity file ",
+" file_cp= ................ gridded velocity file ",
+" file_src= ......... file with source signature",
+" file_rcv=recv.su .. base name for receiver files",
+" file_rcvtime= ..... receiver file in x-t",
+" h= ................ read from model file: if d1==0 then h= can be used to set it",
+" nt=1024 ........... number of time-samples in file_rcvtime",
" xsrc1= ................... x-position of the source (m)",
" ysrc1= ................... y-position of the source (m)",
" zsrc1= ................... z-position of the source (m)",
@@ -18,16 +59,17 @@ char *sdoc[] = {
" ",
" INPUT AND OUTPUT ",
" key=gy ................... input data sorting key",
-" ny=1 ..................... if 2D file number of y traces (2D model)",
-" nxmax=512 ................ maximum number of traces in input files",
-" ntmax=1024 ............... maximum number of samples/trace in input files",
+" nx=1 ..................... if 1D file number of points in x ",
+" ny=1 ..................... if 2D file number of points in y ",
" file_out= ................ output file with traveltime cube",
" file_amp= ................ output file with approximate amplitudes",
-" RAY TRACING ",
-" dT=0 ..................... put traces on one-way time grid with step dT",
-" Tmin=first shot........... minimum time of one-way time grid",
-" Tmax=last shot ........... maximum time of one-way time grid",
-" hom=1 .................... 1: draw straight rays in homogeneous layers",
+" ",
+//" RAY TRACING PARAMETERS:",
+//" dT=0 ..................... put traces on one-way time grid with step dT",
+//" Tmin=first shot........... minimum time of one-way time grid",
+//" Tmax=last shot ........... maximum time of one-way time grid",
+//" hom=1 .................... 1: draw straight rays in homogeneous layers",
+//" ",
" SOURCE POSITIONS ",
" xsrc2=xsrc1 .............. x-position of last source",
" dxsrc=0 .................. step in source x-direction",
@@ -36,11 +78,11 @@ char *sdoc[] = {
" zsrc2=zsrc1 .............. z-position of last source",
" dzsrc=0 .................. step in source z-direction",
" RECEIVER POSITIONS ",
-" xrcv=0,(nx-1)*dx ......... x-position's of receivers (array)",
-" yrcv=0,(ny-1)*dy ......... y-position's of receivers (array)",
+" xrcv=-(nx-1/2)*h,(nx-1/2)*h .. x-position's of receivers (array)",
+" yrcv=-(ny-1)/2*h,(ny-1/2)*h .. y-position's of receivers (array)",
" zrcv=0,0 ................. z-position's of receivers (array)",
-" dxrcv=dx ................. step in receiver x-direction",
-" dyrcv=dy ................. step in receiver y-direction",
+" dxrcv=h ................. step in receiver x-direction",
+" dyrcv=h ................. step in receiver y-direction",
" dzrcv=0 .................. step in receiver z-direction",
" dxspr=0 .................. step of receiver spread in x-direction",
" dyspr=0 .................. step of receiver spread in y-direction",
@@ -65,11 +107,13 @@ NULL};
#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
#define s0(x,y,z) slow0[nxy*(z) + nx*(y) + (x)]
-void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int xs, int ys, int zs, int NCUBE);
-void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube);
-
void main(int argc, char *argv[])
{
+ modPar mod;
+ recPar rec;
+ srcPar src;
+ shotPar shot;
+ rayPar ray;
int
nx, /* x-dimension of mesh (LEFT-TO-RIGHT) */
ny, /* y-dimension of mesh (FRONT-TO-BACK) */
@@ -80,13 +124,11 @@ void main(int argc, char *argv[])
h, /* spatial mesh interval (units consistant with vel) */
*slow0, *time0;
-/* to read the delphi velocity file */
- int32 type, dom1, dom2;
- int error, n1, n2, ret, size, nkeys, verbose;
- int ntmax, nxmax;
- float d1, d2, f1, f2, *tmpdata, c, scl, ox, oz, oy;
- char *file_vel, *file_out, *key;
- segyhdr *hdrs;
+/* to read the velocity file */
+ int error, n1, n2, n3, ret, size, nkeys, verbose;
+ float d1, d2, d3, f1, f2, f3, *tmpdata, c, scl, ox, oz, oy;
+ char *file_cp, *file_out;
+ segy *hdrs;
/*---------------------------------------------------------------------------*
* Read input parameters and query for any that are needed but missing.
@@ -97,122 +139,57 @@ void main(int argc, char *argv[])
if (!getparint("verbose",&verbose)) verbose=0;
if (verbose) {
- samess("Hole, J.A., and B.C. Zelt, 1995. \"3-D finite-difference");
- samess("reflection traveltimes\". Geophys. J. Int., 121, 427-434");
- }
- if(!getparstring("file_out",&file_out)) saerr("file_out not given");
- if(!getparstring("file_vel", &file_vel)) {
- sawarn("file_vel not defined, assuming homogeneous model");
- if(!getparfloat("c",&c)) saerr("c not defined");
- if(!getparint("nx",&nx)) saerr("nx not defined");
- if(!getparint("ny",&ny)) saerr("ny not defined");
- if(!getparint("nz",&nz)) saerr("nz not defined");
- if(!getparfloat("h",&h)) saerr("h not defined");
+ vmess("Hole, J.A., and B.C. Zelt, 1995. \"3-D finite-difference");
+ vmess("reflection traveltimes\". Geophys. J. Int., 121, 427-434");
}
- if(!getparint("ny",&ny)) saerr("ny not defined");
+ if(!getparstring("file_out",&file_out)) verr("file_out not given");
+
+ getParameters3d(&mod, &rec, &src, &shot, &ray, verbose);
+
/*---------------------------------------------------------------------------*
* Open velocity file
*---------------------------------------------------------------------------*/
- if (file_vel != NULL) {
- error = open_file(file_vel, GUESS_TYPE, DELPHI_READ);
- if (error < 0 ) saerr("error in opening file %s", file_vel);
- error = get_dims(file_vel, &n1, &n2, &type);
- if (ret >= 0) {
- if (!getparint("ntmax", &ntmax)) ntmax = n1;
- if (!getparint("nxmax", &nxmax)) nxmax = n2;
- if (verbose>=2 && (ntmax!=n1 || nxmax!=n2))
- samess("dimensions overruled: %d x %d",ntmax,nxmax);
+ if (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 (!getparint("ntmax", &ntmax)) ntmax = 1024;
- if (!getparint("nxmax", &nxmax)) nxmax = 512;
- if (verbose>=2) samess("dimensions used: %d x %d",ntmax,nxmax);
+ 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;
}
- size = ntmax * nxmax;
- tmpdata = alloc1float(size);
- hdrs = (segyhdr *) malloc(nxmax*sizeof(segyhdr));
-
- if (!getparstring("key", &key)) {
- ret = get_sort(file_vel);
- if (ret < 0) key = "gy";
- else key = getkey(ret);
+ else { /* Full 3D model */
+ nz = n1; nx = n2; nz = n3;
+ oz = f1; ox = f2; oy = f3;
+ dz = d1; dx = d1; dy = d1;
}
- if (verbose) samess("input sorting key is %s",key);
- set_sukey(key);
-
- ret = read_data(file_vel,tmpdata,size,&n1,&n2,&f1,&f2,&d1,&d2,
- &type,hdrs);
- if (ret < 0) saerr("error in reading data from file %s", file_vel);
- if (hdrs[0].scalco < 0) scl = 1.0/fabs(hdrs[0].scalco);
- else if (hdrs[0].scalco == 0) scl = 1.0;
- else scl = hdrs[0].scalco;
- get_axis(&dom1, &dom2);
- if (d1 != d2)
- saerr("d1 != d2; this is not allowed in the calculation");
+
h = d1;
- if (dom1 == SA_AXIS_Z) {
- nx = n2; nz = n1;
- ox = hdrs[0].gx*scl; oy = hdrs[0].gy*scl; oz = f1;
- }
- else {
- nx = n1; nz = n2;
- ox = f1; oy = hdrs[0].gy*scl; oz = f1;
- }
+ slow0 = (float *)malloc(nz*nx*ny*sizeof(float));
+ if (slow0 == NULL) verr("Out of memory for slow0 array!");
- slow0 = (float *)malloc(ny*n1*n2*sizeof(float));
- if (slow0 == NULL) saerr("Out of memory for slow0 array!");
- nxy = nx * ny;
- if (verbose) samess("h = %.2f nx = %d nz = %d ny = %d", h, nx, nz, ny);
-
- yy = 0;
- while (ret >= 0) {
- if (verbose==2) disp_info(file_vel,n1,n2,f1,f2,d1,d2,type);
-
- if (dom1 == SA_AXIS_Z) {
- if (n2 != nx || n1 != nz) saerr("dimensions changed");
- for (i = 0; i < n2; i++) {
- for (j = 0; j < n1; j++)
- slow0[j*nxy+yy*nx+i] = h/tmpdata[i*n1+j];
- }
- }
- else {
- if (n1 != nx || n2 != nz) saerr("dimensions changed");
- for (i = 0; i < n2; i++) {
- for (j = 0; j < n1; j++)
- slow0[i*nxy+yy*nx+j] = h/tmpdata[i*n1+j];
- }
- }
+ readModel3d(file_cp, slow0, n1, n2, n3, nz, nx, ny, h, verbose);
- yy += 1;
- ret = read_data(file_vel, tmpdata, size, &n1, &n2, &f1, &f2,
- &d1, &d2, &type, hdrs);
- }
- ret = close_file(file_vel);
- if (ret < 0) sawarn("err %d on closing input file",ret);
- free1float(tmpdata);
- if (yy == 1) {
- if(!getparint("ny",&ny)) samess("2D model defined");
- else {
- slow0 = (float *)realloc(slow0, ny*nx*nz*sizeof(float));
- if (slow0 == NULL) saerr("Out of memory for slow0 array!");
-
- samess("3D model defined from 2D model");
- for (zz = 0; zz < nz; zz++) {
- for (yy = 1; yy < ny; yy++) {
- for (xx = 0; xx < nx; xx++)
- slow0[zz*nxy+yy*nx+xx] = slow0[zz*nxy+xx];
- }
- }
- }
- }
+ if (verbose) vmess("h = %.2f nx = %d nz = %d ny = %d", h, nx, nz, ny);
}
else {
nxy = nx * ny;
+ 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");
+ oz = 0; ox = ((nx-1)/2)*d1; oy = ((ny-1)/2)*d1;
+
slow0 = (float *)malloc(nx*nz*ny*sizeof(float));
- if (slow0 == NULL) saerr("Out of memory for slow0 array!");
+ if (slow0 == NULL) verr("Out of memory for slow0 array!");
scl = h/c;
ox = 0; oy = 0; oz = 0;
for (zz = 0; zz < nz; zz++) {
@@ -228,7 +205,7 @@ void main(int argc, char *argv[])
/* ALLOCATE MAIN GRID FOR TIMES */
time0 = (float *) malloc(sizeof(float)*nxyz);
- if(time0 == NULL) saerr("error in allocation of array time0");
+ if(time0 == NULL) verr("error in allocation of array time0");
/*---------------------------------------------------------------------------*
* Input the source locations.
@@ -237,7 +214,7 @@ void main(int argc, char *argv[])
*---------------------------------------------------------------------------*/
src3d(time0, slow0, nz, nx, ny, h, ox, oy, oz, &xs, &ys, &zs, &cube);
- if (verbose) samess("source positions xs = %d ys = %d zs = %d", xs,ys,zs);
+ if (verbose) vmess("source positions xs = %d ys = %d zs = %d", xs,ys,zs);
/* for (zz = 0; zz < nz; zz++) {
for (yy = 0; yy < ny; yy++) {
@@ -255,7 +232,6 @@ void main(int argc, char *argv[])
* Check and set parameters
*---------------------------------------------------------------------------*/
-
/*---------------------------------------------------------------------------*
* Compute traveltimes.
*---------------------------------------------------------------------------*/
@@ -274,17 +250,17 @@ void main(int argc, char *argv[])
}
}
*/
- ret = open_file(file_out, GUESS_TYPE, DELPHI_CREATE);
- if (ret < 0 ) saerr("error in creating output file %s", file_out);
+// ret = open_file(file_out, GUESS_TYPE, DELPHI_CREATE);
+// if (ret < 0 ) verr("error in creating output file %s", file_out);
- hdrs = (segyhdr *) malloc(ny*sizeof(segyhdr));
- tmpdata = alloc1float(nxy);
+ hdrs = (segy *) malloc(ny*sizeof(segy));
+ tmpdata = (float *)malloc(nxy*sizeof(float));
f1 = ox;
f2 = oy;
d1 = h;
d2 = h;
- gen_hdrs(hdrs,nx,ny,f1,f2,d1,d2,TRID_ZX);
+// gen_hdrs(hdrs,nx,ny,f1,f2,d1,d2,TRID_ZX);
for (i = 0; i < ny; i++) {
hdrs[i].scalco = -1000;
hdrs[i].scalel = -1000;
@@ -300,18 +276,17 @@ void main(int argc, char *argv[])
}
}
- if (ret < 0 ) sawarn("error on writing keys.");
- ret = set_axis(dom1, dom2);
- if (ret < 0 ) saerr("error on writing axis.");
+/*
ret = write_data(file_out,tmpdata,nx,ny,f1,f2,d1,d2,type,hdrs);
- if (ret < 0 ) saerr("error on writing output file.");
+ if (ret < 0 ) verr("error on writing output file.");
ret = close_file(file_out);
- if (ret < 0) saerr("err %d on closing output file",ret);
+ if (ret < 0) verr("err %d on closing output file",ret);
+*/
free(time0);
free(slow0);
free(hdrs);
- free1float(tmpdata);
+ free(tmpdata);
exit(0);
diff --git a/raytime3d/readModel.c b/raytime3d/readModel.c
deleted file mode 100644
index 27f1da715a50dbc36eac4827efe68434288ee9e7..0000000000000000000000000000000000000000
--- a/raytime3d/readModel.c
+++ /dev/null
@@ -1,80 +0,0 @@
-#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 "raytime.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 readModel(modPar mod, float *velocity, float *slowness, int nw)
-{
- FILE *fpcp;
- size_t nread;
- int i, tracesToDo, j;
- int nz, nx;
- segy hdr;
-
-
- /* grid size and start positions for the components */
- nz = mod.nz;
- nx = mod.nx;
-
-/* open files and read first header */
-
- fpcp = fopen( mod.file_cp, "r" );
- assert( fpcp != NULL);
- nread = fread(&hdr, 1, TRCBYTES, fpcp);
- assert(nread == TRCBYTES);
-
-/* read all traces */
-
- tracesToDo = mod.nx;
- i = 0;
- while (tracesToDo) {
- nread = fread(&velocity[i*nz], sizeof(float), hdr.ns, fpcp);
- assert (nread == hdr.ns);
- for (j=0;j<nz;j++) {
- if (velocity[i*nz+j]!=0.0) {
- slowness[(i+nw)*nz+j+nw] = 1.0/velocity[i*nz+j];
- }
- }
- for (j=0;j<nw;j++) slowness[(i+nw)*nz+j] = slowness[(i+nw)*nz+nw];
- for (j=nz+nw;j<nz+2*nw;j++) slowness[(i+nw)*nz+j] = slowness[(i+nw)*nz+nz+nw-1];
-
- nread = fread(&hdr, 1, TRCBYTES, fpcp);
- if (nread==0) break;
- i++;
- }
- fclose(fpcp);
-
- for (i=0;i<nw;i++) {
- for (j=0;j<nz+2*nw;j++) {
- slowness[(i)*nz+j] = slowness[(nw)*nz+j];
- slowness[(nx+nw+i)*nz+j] = slowness[(nx+nw-1)*nz+j];
- }
- }
-
- return 0;
-}
-
-
diff --git a/raytime3d/recvPar.c b/raytime3d/recvPar.c
deleted file mode 100644
index 4ad9eae9be2f3c7bcbeb1ee061dc97d5a4ff9d77..0000000000000000000000000000000000000000
--- a/raytime3d/recvPar.c
+++ /dev/null
@@ -1,519 +0,0 @@
-#include <stdio.h>
-#include <assert.h>
-#include <math.h>
-
-#include "raytime.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;
-}
diff --git a/raytime3d/src3d.c b/raytime3d/src3d.c
index 1dadfa5746b00a8bb5504f16de26188c06da5c97..545c1974d14ab60f4bd28b259e68d2d4b4ee0e21 100644
--- a/raytime3d/src3d.c
+++ b/raytime3d/src3d.c
@@ -1,5 +1,8 @@
-#include <cwp.h>
-#include <su.h>
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
#include <fcntl.h>
#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
@@ -7,10 +10,10 @@
#define SQR(x) ((x) * (x))
#define DIST(x,y,z,x1,y1,z1) sqrt(SQR(x-(x1))+SQR(y-(y1)) + SQR(z-(z1)))
-/* definitions from saerrpkge.c */
-extern void saerr(char *fmt, ...);
-extern void sawarn(char *fmt, ...);
-extern void samess(char *fmt, ...);
+/* definitions from verbose.c */
+extern void verr(char *fmt, ...);
+extern void vwarn(char *fmt, ...);
+extern void vmess(char *fmt, ...);
void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox, float oy, float oz, int *pxs, int *pys, int *pzs, int *cube)
{
@@ -48,9 +51,9 @@ void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox
if(!getparint("srctype",&srctype)) srctype=1;
if(srctype==1) {
- if(!getparfloat("xsrc1",&xsrc1)) saerr("xsrc1 not given");
- if(!getparfloat("ysrc1",&ysrc1)) saerr("ysrc1 not given");
- if(!getparfloat("zsrc1",&zsrc1)) saerr("zsrc1 not given");
+ if(!getparfloat("xsrc1",&xsrc1)) verr("xsrc1 not given");
+ if(!getparfloat("ysrc1",&ysrc1)) verr("ysrc1 not given");
+ if(!getparfloat("zsrc1",&zsrc1)) verr("zsrc1 not given");
fxs = (xsrc1-ox)/h;
fys = (ysrc1-oy)/h;
fzs = (zsrc1-oz)/h;
@@ -58,30 +61,30 @@ void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox
ys = (int)(fys + 0.5);
zs = (int)(fzs + 0.5);
if(xs<2 || ys<2 || zs<2 || xs>nx-3 || ys>ny-3 || zs>nz-3){
- sawarn("Source near an edge, beware of traveltime errors");
- sawarn("for raypaths that travel parallel to edge ");
- sawarn("while wavefronts are strongly curved, (JV, 8/17/88)\n");
+ vwarn("Source near an edge, beware of traveltime errors");
+ vwarn("for raypaths that travel parallel to edge ");
+ vwarn("while wavefronts are strongly curved, (JV, 8/17/88)\n");
}
*pxs = xs; *pys = ys, *pzs = zs, *cube = NCUBE;
}
else if (srctype==2) {
- if (!getparint("srcwall",&srcwall)) saerr("srcwall not given");
- if (!getparstring("wallfile",&wallfile)) saerr("wallfile not given");
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
+ if (!getparstring("wallfile",&wallfile)) verr("wallfile not given");
if((wfint=open(wallfile,O_RDONLY,0664))<=1) {
fprintf(stderr,"cannot open %s\n",wallfile);
exit(-1);
}
}
else if (srctype==3) {
- if (!getparint("srcwall",&srcwall)) saerr("srcwall not given");
- if (!getparstring("oldtfile",&oldtfile)) saerr("oldtfile not given");
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
+ if (!getparstring("oldtfile",&oldtfile)) verr("oldtfile not given");
if((ofint=open(oldtfile,O_RDONLY,0664))<=1) {
fprintf(stderr,"cannot open %s\n",oldtfile);
exit(-1);
}
}
else {
- saerr("ERROR: incorrect value of srctype");
+ verr("ERROR: incorrect value of srctype");
}
nxy = nx * ny;
@@ -127,8 +130,8 @@ void src3d(float *time0, float *slow0, int nz, int nx, int ny, float h, float ox
/* can't handle asin(1.) ! */
if (fabs(rz1-rzc)>=rho) rho=1.0000001*fabs(rz1-rzc);
theta2 = asin((rz1-rzc)/rho);
- if (rxyc<0) theta1=PI-theta1;
- if (rxyc<rxy1) theta2=PI-theta2;
+ if (rxyc<0) theta1=M_PI-theta1;
+ if (rxyc<rxy1) theta2=M_PI-theta2;
t0(xx,yy,zz) = log(tan(theta2/2)/tan(theta1/2)) / dv;
}
}
diff --git a/raytime3d/vidale3d.c b/raytime3d/vidale3d.c
index 06b4b388274def5752a95ab1f0bbb1b881b47c11..ff8a74c1c4f75f4179ce168c514b6ee1949ab837 100644
--- a/raytime3d/vidale3d.c
+++ b/raytime3d/vidale3d.c
@@ -1,5 +1,9 @@
-#include <cwp.h>
-#include <su.h>
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include<string.h>
+#include"par.h"
#define SQR2 1.414213562
#define SQR3 1.732050808
@@ -7,10 +11,10 @@
#define t0(x,y,z) time0[nxy*(z) + nx*(y) + (x)]
#define s0(x,y,z) slow0[nxy*(z) + nx*(y) + (x)]
-/* definitions from saerrpkge.c */
-extern void saerr(char *fmt, ...);
-extern void sawarn(char *fmt, ...);
-extern void samess(char *fmt, ...);
+/* definitions from verbose.c */
+extern void verr(char *fmt, ...);
+extern void vwarn(char *fmt, ...);
+extern void vmess(char *fmt, ...);
struct sorted
{ float time; int i1, i2;};
@@ -82,7 +86,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
/* SET MAXIMUM RADIUS TO COMPUTE */
if (maxoff > 0.) {
maxrad = maxoff/h + 1;
- sawarn("WARNING: Computing only to max radius = %d",maxrad);
+ vwarn("WARNING: Computing only to max radius = %d",maxrad);
}
else maxrad = 99999999;
@@ -106,7 +110,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
}
wall = (float *) malloc(4*nwall);
if(sort == NULL || wall == NULL)
- saerr("error in allocation of arrays sort and wall");
+ verr("error in allocation of arrays sort and wall");
if(!getparint("srctype",&srctype)) srctype=1;
if(srctype==1) {
@@ -126,7 +130,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
else{ z2 = nz; dz2 = 0;}
}
else {
- if (!getparint("srcwall",&srcwall)) saerr("srcwall not given");
+ if (!getparint("srcwall",&srcwall)) verr("srcwall not given");
/* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE */
radius = 1;
if (srcwall == 1) x2=1;
@@ -1245,7 +1249,7 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
/* UPDATE RADIUS */
radius++;
- if(radius%10 == 0 && verbose) samess("Completed radius = %d",radius);
+ if(radius%10 == 0 && verbose) vmess("Completed radius = %d",radius);
if(radius == maxrad) rad0 = 0;
} /* END BIG LOOP */
@@ -1259,42 +1263,42 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
else {
head=0;
if (headw[1]>0) {
- if(verbose) samess("Head waves found on left: %d",headw[1]);
+ if(verbose) vmess("Head waves found on left: %d",headw[1]);
if (headw[1]>head) {
head = headw[1];
srcwall = 1;
}
}
if (headw[2]>0) {
- if(verbose) samess("Head waves found on right: %d",headw[2]);
+ if(verbose) vmess("Head waves found on right: %d",headw[2]);
if (headw[2]>head) {
head = headw[2];
srcwall = 2;
}
}
if (headw[3]>0) {
- if(verbose) samess("Head waves found on front: %d",headw[3]);
+ if(verbose) vmess("Head waves found on front: %d",headw[3]);
if (headw[3]>head) {
head = headw[3];
srcwall = 3;
}
}
if (headw[4]>0) {
- if(verbose) samess("Head waves found on back: %d",headw[4]);
+ if(verbose) vmess("Head waves found on back: %d",headw[4]);
if (headw[4]>head) {
head = headw[4];
srcwall = 4;
}
}
if (headw[5]>0) {
- if(verbose) samess("Head waves found on top: %d",headw[5]);
+ if(verbose) vmess("Head waves found on top: %d",headw[5]);
if (headw[5]>head) {
head = headw[5];
srcwall = 5;
}
}
if (headw[6]>0) {
- if(verbose) samess("Head waves found on bottom: %d",headw[6]);
+ if(verbose) vmess("Head waves found on bottom: %d",headw[6]);
if (headw[6]>head) {
head = headw[6];
srcwall = 6;
@@ -1302,32 +1306,32 @@ void vidale3d(float *slow0, float *time0, int nz, int nx, int ny, float h, int x
}
if (headpref>0 && headw[headpref]>0) {
if(verbose)
- samess("Preference to restart on wall opposite source");
+ vmess("Preference to restart on wall opposite source");
srcwall = headpref;
}
/* SET LOCATIONS OF SIDES OF THE CUBE SO THAT CUBE IS A FACE */
dx1=1; dx2=1; dy1=1; dy2=1; dz1=1; dz2=1; rad0=1;
radius = 1;
if (srcwall == 1) { x2=1;
- samess("RESTART at left side of model"); }
+ vmess("RESTART at left side of model"); }
else { x2=nx; dx2=0; }
if (srcwall == 2) { x1=nx-2;
- samess("RESTART at right side of model"); }
+ vmess("RESTART at right side of model"); }
else { x1= -1; dx1=0; }
if (srcwall == 3) { y2=1;
- samess("RESTART at front side of model"); }
+ vmess("RESTART at front side of model"); }
else { y2=ny; dy2=0; }
if (srcwall == 4) { y1=ny-2;
- samess("RESTART at back side of model"); }
+ vmess("RESTART at back side of model"); }
else { y1= -1; dy1=0; }
if (srcwall == 5) { z2=1;
- samess("RESTART at top side of model"); }
+ vmess("RESTART at top side of model"); }
else { z2=nz; dz2=0; }
if (srcwall == 6) { z1=nz-2;
- samess("RESTART at bottom side of model"); }
+ vmess("RESTART at bottom side of model"); }
else { z1= -1; dz1=0; }
if (reverse == 0)
- sawarn("RESTART CANCELLED by choice of input parameter `reverse`");
+ vwarn("RESTART CANCELLED by choice of input parameter `reverse`");
}
reverse--;
diff --git a/raytime3d/writeSrcRecPos.c b/raytime3d/writeSrcRecPos.c
index faf297811e469314c8de0265f78b87543b27af9b..d9d03da16dae329b7628e3420224bda8d9a6cce4 100644
--- a/raytime3d/writeSrcRecPos.c
+++ b/raytime3d/writeSrcRecPos.c
@@ -3,7 +3,7 @@
#include<math.h>
#include<assert.h>
#include"par.h"
-#include"raytime.h"
+#include"raytime3d.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
diff --git a/raytime3d/writesufile.c b/raytime3d/writesufile.c
index 6eac57d300aaa653e621b99fe2731c9ed3ddb60a..7ebcbbd198858193b67a2b9c2fbb5670feb8b5d8 100644
--- a/raytime3d/writesufile.c
+++ b/raytime3d/writesufile.c
@@ -3,7 +3,7 @@
#include <assert.h>
#include <string.h>
#include "par.h"
-#include "raytime.h"
+#include "raytime3d.h"
#include "SUsegy.h"
#include "segy.h"