getrecpos.c

#include "par.h"
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

/**
* read receiver positions used in green
*
*   AUTHOR:
*           Jan Thorbecke (janth@xs4all.nl)
*           The Netherlands 
**/

#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
#ifndef MAX
#define	MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
#ifndef MIN
#define	MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define SGN(x) ((x) < 0 ? -1.0 : 1.0)
#ifndef ABS
#define ABS(x) ((x) < 0 ? -(x) : (x))
#endif

void getrecpos(float *xi, float *zi, int nx, float *xrcv, float *zrcv, int verbose)
{
	int		nrx, i, j, ndeltx, np, lint, seed;
	long    idum;
	float	xprev, zprev, deltx, deltz, dxrcv, dzrcv, var, irr, maxirr;
    float rrcv, dphi, oxrcv, ozrcv;

	nrx = countparval("xrcv");
	if(!getparfloat("dxrcv",&dxrcv)) dxrcv = 15;
	if(!getparfloat("var", &var)) var=0;
	if(!getparint("lint", &lint)) lint=1;
	if(!getparint("seed", &seed)) seed=0;
    
    /* check if receiver positions on a circle are defined */
	if (getparfloat("rrcv", &rrcv)) {
		if (!getparfloat("dphi",&dphi)) dphi=2.0;
		if (!getparfloat("oxrcv",&oxrcv)) oxrcv=0.0;
		if (!getparfloat("ozrcv",&ozrcv)) ozrcv=0.0;
		
        np = 0;
		for (i=0; i<nx; i++) {
			xi[np]   = oxrcv+rrcv*cos(((i*dphi)/360.0)*(2.0*M_PI));
			zi[np++] = ozrcv+rrcv*sin(((i*dphi)/360.0)*(2.0*M_PI));
			if (verbose>4) fprintf(stderr,"Receiver Circle: xrcv[%d]=%f zrcv=%f\n", i, xi[i],zi[i]);
		}
		return;
	}


	if (var <= 0) {
		if (lint == 1) {
			xprev = xrcv[0];
			zprev = zrcv[0];
			np = 0;
			for (i = 1; i < nrx; i++) {
				deltx = xrcv[i] - xprev;
				deltz = zrcv[i] - zprev;
				ndeltx = NINT(ABS(deltx/dxrcv));
				dzrcv = deltz/ndeltx;
				for (j = 0; j < ndeltx; j++) {
					zi[np]   = zprev + j*dzrcv;
					xi[np++] = xprev + j*dxrcv;
				}
				xprev = xrcv[i];
				zprev = zrcv[i];
			}
			xi[nx-1] = xrcv[nrx-1];
			zi[nx-1] = zrcv[nrx-1];
		}
		else {
			for (i = 0; i < nrx; i++) {
				xi[i] = xrcv[i];
				zi[i] = zrcv[i];
			}
		}
	}
	else {
		xprev = xrcv[0];
		zprev = zrcv[0];
		np = 0;
		maxirr = 0;
		idum = (long) seed;
		srand48(idum);
		for (i = 1; i < nrx; i++) {
			deltx = xrcv[i] - xprev;
			deltz = zrcv[i] - zprev;
			ndeltx = NINT(ABS(deltx/dxrcv));
			dzrcv = deltz/ndeltx;
			for (j = 0; j < ndeltx; j++) {
				irr = var*((float)drand48());
				if (fabs(irr) > maxirr) maxirr = fabs(irr);
				zi[np]   = zprev + j*dzrcv;
				xi[np++] = xprev + j*dxrcv + irr;
				if (verbose==13)vmess("xrcv %d = %f (%f)",np-1,xi[np-1], irr);
			}
			xprev = xrcv[i];
			zprev = zrcv[i];
		}
		irr = var*((float)drand48());
		if (fabs(irr) > maxirr) maxirr = fabs(irr);
		xi[nx-1] = xrcv[nrx-1] + irr;
		zi[nx-1] = zrcv[nrx-1];
		if (verbose) vmess("maximum error in receiver position %f", maxirr);
		if (verbose==13) vmess("xrcv %d = %f (%f)", nx-1, xi[nx-1], irr);
	}

	if (verbose) vmess("getrecpos number of receivers = %d", np+1);

	return;
}