vmerge.c

seismic software,very useful

    }

    gh.n4 = gh.scale;
    gh.n5 = gh.scale;
    gh.d4 = sbin * gh.scale;
    gh.d5 = vbin * gh.scale;
    gh.o4 = shigh * gh.scale;;
    gh.o5 = vclip[0] * gh.scale;
    gh.dcdp2 = dv * gh.scale;
    gh.dline3 = ddv * gh.scale;
    gh.oline3 = 0.;
    gh.ocdp2 = 0.;
    gh.gtype=5 * gh.scale;

    write(1, &gh, sizeof(gh));

    return( 0 );
}

/*--------------------------------------------------------------------*\
  Read a pair of data vectors from the two input files.
\*--------------------------------------------------------------------*/

static int readvec(FILE *fp, char *type, float *x, int n)
{
    int      nread;

    if( !strcmp(type, "segy") ){
        fseeko(fp, (long long)240, SEEK_CUR);
        if( fread(x, 4, n,fp) != n ){
            return (0);

        }else{
            return (n);

        }

    }else if( !strcmp(type, "vgrid") ){

        if( (nread = fread(x,4,n,fp)) == n ){
            return (n);

        }else{
            return (nread - 25);

        }

    }else if( !strcmp(type, "float") || !strcmp(type, "floats") ){
        return (fread(x,4,n,fp));

    }else if( !strcmp(type, "byte") || !strcmp(type, "bytes") ){
        return (fread(x,1,n,fp));

    }

    return (0);
}

/*--------------------------------------------------------------------*\
   Determine the total size of the input dataset.
\*--------------------------------------------------------------------*/

static void findsize(FILE *fp, char *type, long *size, int *ns)
{
    segy     t1;

    *ns = N;

    if( !strcmp(type, "segy") ){

        fseeko(fp, (long long)3600, SEEK_SET);
        fread(&t1,1,240,fp);

        fseeko(fp, (long long)0, SEEK_END);

        *ns = ftello(fp) / (t1.ns * 4 + 240);
        *size = t1.ns * *ns;
        *ns = t1.ns;

    }else if( !strcmp(type, "vgrid") ){
        fseeko(fp, (long long)0, SEEK_END);
        *size = (ftello(fp) - 100) / 4;

    }else if( !strcmp(type, "float") || !strcmp(type, "floats") ){
        fseeko(fp, (long long)0, SEEK_END);
        *size = ftello(fp) / 4;

    }else if( !strcmp(type, "byte") || !strcmp(type, "bytes") ){
        fseeko(fp, (long long)0, SEEK_END);
        *size = ftello(fp);

    }
}

/*--------------------------------------------------------------------*\
\*--------------------------------------------------------------------*/

static void findclip(FILE *fp, char *type, long size, float *low, float *high,
                  float perc, float nan)
{
#define NSORT 10000000

    int      nsort;
    int      dsort;
    int      i;
    int      n;
    float   *buf;
    unsigned char cbuf[N];
    segy     tr;

    if( !strcmp(type, "byte") || !strcmp(type, "bytes") ){
        return;

    }

    nsort = size < NSORT ? size : NSORT;
    nsort = nsort > (int) (.05 * size) ? nsort : (int) (.05 * size);
    dsort = ceil((float) size / (float) nsort);
    buf   = (float *) malloc((nsort + N) * 4);

    if( !strcmp(type, "segy") ){
        fseeko(fp, (long long)3600, SEEK_SET);
        fread(&tr,1,240,fp);
        nsort = 0;

        while( fread(tr.data,1,tr.ns * 4 + 240,fp) > 0 ){

            for( i = 0; i < tr.ns; i += dsort ){
                if( tr.data[i] != nan ){
                    buf[nsort++] = tr.data[i];
                }
            }
        }

    }else if( !strcmp(type, "vgrid") ){
        fseeko(fp, (long long)0, SEEK_SET);
        nsort = 0;

        while( (n = fread(tr.data,1,N * 4,fp) / 4) > 0 ){

            for( i = 0; i < n; i += dsort, nsort++ ){
                buf[nsort] = tr.data[i];
            }
        }
        nsort -= 25;

    }else if( !strcmp(type, "float") || !strcmp(type, "floats") ){
        fseeko(fp, (long long)0, SEEK_SET);
        nsort = 0;
        while( (n = fread(tr.data,1, N * 4,fp) / 4) > 0 ){

            for( i = 0; i < n; i += dsort, nsort++ ){
                buf[nsort] = tr.data[i];
            }
        }

    }else if( !strcmp(type, "float") || !strcmp(type, "floats") ){
        fseeko(fp, (long long)0, SEEK_SET);
        nsort = 0;
        while( (n = fread(cbuf,1,N,fp)) > 0 ){

            for( i = 0; i < n; i += dsort, nsort++ ){
                buf[nsort] = cbuf[i];
            }
        }
    }

    perc = perc < 100. ? perc : 100.;
    *low = cent(buf, nsort, 100. - perc);
    *high = cent(buf, nsort, perc);
    free(buf);
}

/*--------------------------------------------------------------------*\
  percentile subroutine based on Canales, SEP-10
  this routine changes data order, so sort a copy

  p - percentile <0.,99.999999999999>
  x - data
  n - vector length
 \*--------------------------------------------------------------------*/

static float cent(float *x, int n, float p ){

    int      q;
    register float *i, *j, ak;
    float   *low, *hi, buf, *k;

    p = p < 99.999 ? p : 99.999;
    p = p > 0.0 ? p : 0.0;
    q = (p * n) / 100.;

    if( q == n)
        q--;

    for( low = x, hi = x + n - 1, k = x + q; low < hi; ){
        ak = *k;
        i = low;
        j = hi;
        do{
            while( *i < ak)
                i++;
            while( *j > ak)
                j--;
            if( i <= j ){
                buf = *i;
                *i++ = *j;
                *j-- = buf;
            }
        }while( i <= j);
        if( j < k)
            low = i;
        if( k < i)
            hi = j;
    }
    return (*k);
}

/*--------------------------------------------------------------------*\
\*--------------------------------------------------------------------*/

static int findbin(float *bin, int nbin, float low, float high)
{
    float    interval;
    int      nbin1;

    interval = (high - low) / nbin;

    interval = pow(10., rint(log10(interval)));

    /* fprintf (stderr,"interval=%g\n",interval); */

    if( (nbin1 = findbin2(bin, nbin, low, high, 0.2 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 0.25 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 0.3 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 0.4 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 0.5 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 0.6 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 0.75 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 1.0 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 2.0 * interval)) > 0 ){
        return (nbin1);
    }

    if( (nbin1 = findbin2(bin, nbin, low, high, 5.0 * interval)) > 0 ){
        return (nbin1);
    }

    return (0);
}

/*--------------------------------------------------------------------*\
\*--------------------------------------------------------------------*/

static int findbin2(float *bin, int nbin, float low, float high, float interval)
{
    int      ilow, ihigh, i;

    ilow = ceil(low / interval);
    ihigh = floor(high / interval);

#if 1

    fprintf (stderr ,"interval=%g ilow=%d ihigh=%d nbin=%d\n"
                ,interval,ilow,ihigh,ihigh-ilow+2); 

#endif

    if( ihigh - ilow + 2 > nbin)
        return (0);

    for( i = 0; ilow <= ihigh; i++, ilow++ ){
        bin[i] = ilow * interval;
    }
    return (i);
}

/*--------------------------------------------------------------------*\
\*--------------------------------------------------------------------*/

static int findbin3(float *bin, int nbin, float low, float high, float range)
{
    int      i;
    float    ratio, length, interval;

    ratio    = pow(range, 1. / (nbin - 3));
    interval = 1.0;
    length   = 0.0;

    for( i = 0; i < nbin - 2; i++ ){
        interval *= ratio;
        length   += interval;
    }

    interval = (high - low) / length;

    fprintf( stderr 
            ,"range=%g ratio=%g length=%g low=%g high=%g interval=%g-%g " 
            ,range ,ratio ,length ,low ,high ,interval ,interval * range);

    bin[0] = low;

    for (i = 1; i < nbin - 1; i++ ){
        interval *= ratio;
        bin[i] = bin[i - 1] + interval;

    }

    fprintf(stderr, "bin1=%g bin$=%g\n", bin[0], bin[nbin - 2]);

    return (nbin - 1);
}