2dimaging.mpi.c

主从模式粗粒级并行算法C程序：这是我以前研究生期间编写的叠前地震成像C源码

	step = 1;

	/* if x values increasing */
	if (ax[nx-1]>ax[0]) {

		/* find indices such that ax[lower] <= x < ax[upper] */
		while (lower>0 && ax[lower]>x) {
			upper = lower;
			lower -= step;
			step += step;
		}
		if (lower<0) lower = 0;
		while (upper<nx && ax[upper]<=x) {
			lower = upper;
			upper += step;
			step += step;
		}
		if (upper>nx) upper = nx;

		/* find index via bisection */
		while ((middle=(lower+upper)>>1)!=lower) {
			if (x>=ax[middle])
				lower = middle;
			else
				upper = middle;
		}

	/* else, if not increasing */
	} else {

		/* find indices such that ax[lower] >= x > ax[upper] */
		while (lower>0 && ax[lower]<x) {
			upper = lower;
			lower -= step;
			step += step;
		}
		if (lower<0) lower = 0;
		while (upper<nx && ax[upper]>=x) {
			lower = upper;
			upper += step;
			step += step;
		}
		if (upper>nx) upper = nx;

		/* find index via bisection */
		while ((middle=(lower+upper)>>1)!=lower) {
			if (x<=ax[middle])
				lower = middle;
			else
				upper = middle;
		}
	}

	/* return lower index */
	*index = lower;
}
/*****************************************************************************
INTLIN - evaluate y(x) via linear interpolation of y(x[0]), y(x[1]), ...

intlin		evaluate y(x) via linear interpolation of y(x[0]), y(x[1]), ...

******************************************************************************
Function Prototype:
void intlin (int nin, float xin[], float yin[], float yinl, float yinr,
	int nout, float xout[], float yout[]);

******************************************************************************
Input:
nin		length of xin and yin arrays
xin		array[nin] of monotonically increasing or decreasing x values
yin		array[nin] of input y(x) values
yinl		value used to extraplate y(x) to left of input yin values
yinr		value used to extraplate y(x) to right of input yin values
nout		length of xout and yout arrays
xout		array[nout] of x values at which to evaluate y(x)

Output:
yout		array[nout] of linearly interpolated y(x) values

******************************************************************************
Notes:
xin values must be monotonically increasing or decreasing.

Extrapolation of the function y(x) for xout values outside the range
spanned by the xin values in performed as follows:

	For monotonically increasing xin values,
		yout=yinl if xout<xin[0], and yout=yinr if xout>xin[nin-1].

	For monotonically decreasing xin values, 
		yout=yinl if xout>xin[0], and yout=yinr if xout<xin[nin-1].

If nin==1, then the monotonically increasing case is used.

******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/02/89
*****************************************************************************/
/**************** end self doc ********************************/

void intlin (int nin, float xin[], float yin[], float yinl, float yinr, 
	int nout, float xout[], float yout[])
/*****************************************************************************
evaluate y(x) via linear interpolation of y(x[0]), y(x[1]), ...
******************************************************************************
Input:
nin		length of xin and yin arrays
xin		array[nin] of monotonically increasing or decreasing x values
yin		array[nin] of input y(x) values
yinl		value used to extraplate y(x) to left of input yin values
yinr		value used to extraplate y(x) to right of input yin values
nout		length of xout and yout arrays
xout		array[nout] of x values at which to evaluate y(x)

Output:
yout		array[nout] of linearly interpolated y(x) values
******************************************************************************
Notes:
xin values must be monotonically increasing or decreasing.

Extrapolation of the function y(x) for xout values outside the range
spanned by the xin values in performed as follows:

	For monotonically increasing xin values,
		yout=yinl if xout<xin[0], and yout=yinr if xout>xin[nin-1].

	For monotonically decreasing xin values, 
		yout=yinl if xout>xin[0], and yout=yinr if xout<xin[nin-1].

If nin==1, then the monotonically increasing case is used.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/02/89
*****************************************************************************/
{
	static int idx;
	int jout;
	float x;

	/* if input x values are monotonically increasing, then */
	if (xin[0]<=xin[nin-1]) {
		for (jout=0; jout<nout; jout++) {
			x = xout[jout];
			if (x<xin[0])
				yout[jout] = yinl;
			else if (x>xin[nin-1])
				yout[jout] = yinr;
			else if (x==xin[nin-1] || nin==1)
				yout[jout] = yin[nin-1];
			else {
				xindex(nin,xin,x,&idx);
				yout[jout] = yin[idx]+(x-xin[idx])
					*(yin[idx+1]-yin[idx])
					/(xin[idx+1]-xin[idx]);
			}
		}
	
	/* else, if input x values are monotonically decreasing, then */
	} else {
		for (jout=0; jout<nout; jout++) {
			x = xout[jout];
			if (x>xin[0])
				yout[jout] = yinl;
			else if (x<xin[nin-1])
				yout[jout] = yinr;
			else if (x==xin[nin-1] || nin==1)
				yout[jout] = yin[nin-1];
			else {
				xindex(nin,xin,x,&idx);
				yout[jout] = yin[idx]+(x-xin[idx])
					*(yin[idx+1]-yin[idx])
					/(xin[idx+1]-xin[idx]);
			}
		}
	}
}

void keyin_par(char *iwfile, char *ivfile, char *owfile1, char *owfile2, 
     char *owfile3, int *seismic, int *operator, int *nxv, int *nzv,
     int *ncdp0, float *dmx, float *dhx,  float *depth, float *dz, 
     float *f1, float *f2, float *f3, float *f4, int *iflag_fft,
     float *dxv, float *dzv, float *xv0, int *ncdpstep_cig, 
     float *gamamin, float *gamamax, int *ngama)
{
        FILE *fp;

        fp=fopen("2DcmpPreSDM_ODCIG2ADCIG_Gama.PAR","w");  

        printf("***************************************************\n");
        printf("*             Parameters Inputing                 *\n");
        printf("***************************************************\n");
        printf("Keyin shot gather filename (*su): iwfile\n");
        scanf("%s",iwfile);
        fprintf(fp,"%s\n",iwfile);
        printf("Keyin velocity file name(pc no su head): ivfile\n");
        scanf("%s",ivfile);
        fprintf(fp,"%s\n",ivfile);
        printf("Keyin filename of image by 't=0,h=0' before output CIGs: owfile1\n");
        scanf("%s",owfile1);
        fprintf(fp,"%s\n",owfile1);
        printf("Keyin filename of offset domain image-point gather file name: owfile2\n");
        scanf("%s",owfile2);
        fprintf(fp,"%s\n",owfile2);
        printf("Keyin filename of angle domain image-point gather file name: owfile3\n");
        scanf("%s",owfile3);
        fprintf(fp,"%s\n",owfile3);
 
        printf("Keyin input/ouput seismic or velocity data type: seismic\n");
        printf("SEG-Y type: seismic=1\n");
        printf("SU type:    seismic=2\n");
        scanf("%d",seismic);
        fprintf(fp,"%d\n",*seismic);

        printf("Keyin imaging operator type: operator\n");
        printf("Split-step Fourier:           1\n");
        printf("Quasi-linear Fourier Born:    2\n");
        scanf("%d",operator);
        fprintf(fp,"%d\n",*operator);

        printf("Keyin velocity field lateral/vertical gridpoint number: nxv nzv\n");
        scanf("%d %d",nxv, nzv);
        fprintf(fp,"%d %d\n",*nxv, *nzv);

        printf("Keyin wavefield first cdp number: ncdp0\n");
        printf("Default value: ncdp0=1\n");
        scanf("%d",ncdp0);
        fprintf(fp,"%d\n",*ncdp0);

        printf("Keyin midpoint interval: dmx (m)\n");
        scanf("%f",dmx);
        fprintf(fp,"%f\n",*dmx);

        printf("Keyin half-offset interval: dhx (m)\n");
        scanf("%f",dhx);
        fprintf(fp,"%f\n",*dhx);

        printf("Keyin maximum extrapolating depth: depth \n");
        scanf("%f",depth);
        fprintf(fp,"%f\n",*depth);

        printf("Keyin depth extrapolating interval: dz\n");
        scanf("%f",dz);
        fprintf(fp,"%f\n",*dz);

        printf("Keyin imaging frequency range: f1-f2--f3-f4 (Hz)\n");
        scanf("%f %f %f %f",f1, f2, f3, f4);
        fprintf(fp,"%f %f %f %f\n",*f1, *f2, *f3, *f4);

        printf("If FFTed Wavefield Temp File Existing,  Keyin: 1\n");
        printf("If FFTed Wavefield Temp File not Exist, Keyin: 0\n");
        scanf("%d",iflag_fft);
        fprintf(fp,"%d\n",*iflag_fft);

        printf("Keyin velocity lateral/vertical grid interval: dxv dzv (m)\n");
        printf("Requred:  dxv = dx, dzv=dz \n");
        scanf("%f %f",dxv, dzv);
        fprintf(fp,"%f %f\n",*dxv, *dzv);

        printf("Keyin velocity lateral starting coordinate: xv0 (m)\n");
        scanf("%f",xv0);
        fprintf(fp,"%f\n",*xv0);

	printf("Keyin CIGs cdp number interval: ncdpstep_cig\n");
        scanf("%d",ncdpstep_cig);
        fprintf(fp,"%d\n",*ncdpstep_cig);

	printf("Keyin Minimum & maximum Incidence Angle: gamamin, gamamax\n");
	printf("Example:  (0,80), (-80, 0), (-80, +80)\n");
	scanf("%f %f", gamamin, gamamax);
	fprintf(fp,"%f %f\n",*gamamin, *gamamax);

        printf("Keyin Incidence Angle Numer: ngama\n");
        printf("Suggested: 81 \n");
        scanf("%d",ngama);
        fprintf(fp,"%d\n",*ngama);

        fclose(fp);
} 
void read_par(char *iwfile, char *ivfile, char *owfile1, char *owfile2, 
     char *owfile3, int *seismic, int *operator, int *nxv, int *nzv,
     int *ncdp0, float *dmx, float *dhx,  float *depth, float *dz, 
     float *f1, float *f2, float *f3, float *f4, int *iflag_fft,
     float *dxv, float *dzv, float *xv0, int *ncdpstep_cig, 
     float *gamamin, float *gamamax, int *ngama)
{
        FILE *fp;

        fp=fopen("2DcmpPreSDM_ODCIG2ADCIG_Gama.PAR","r");  

        fscanf(fp,"%s\n",iwfile);
        fscanf(fp,"%s\n",ivfile);
        fscanf(fp,"%s\n",owfile1);
        fscanf(fp,"%s\n",owfile2);
        fscanf(fp,"%s\n",owfile3);
        fscanf(fp,"%d\n",seismic);
        fscanf(fp,"%d\n",operator);
        fscanf(fp,"%d %d\n",nxv, nzv);
        fscanf(fp,"%d\n",ncdp0);
        fscanf(fp,"%f\n",dmx);
        fscanf(fp,"%f\n",dhx);
        fscanf(fp,"%f\n",depth);  
        fscanf(fp,"%f\n",dz);   
        fscanf(fp,"%f %f %f %f\n",f1, f2, f3, f4);
        fscanf(fp,"%d\n",iflag_fft);
        fscanf(fp,"%f %f\n",dxv, dzv);
        fscanf(fp,"%f\n",xv0);
	fscanf(fp,"%d\n",ncdpstep_cig);
        fscanf(fp,"%f %f\n",gamamin, gamamax);
        fscanf(fp,"%d",ngama);
        fclose(fp);
}