vgrid2hvel.c

seismic software,very useful

    }
    if (!getparfloat("dtvgrid",&dtvgrid)) { 
        if(ierr==0) {
            getgval(&gh,"d1",&dtvgrid);
        } else {
            err(" dtvgrid missing ");
        }  
    }
    if (!getparint("ntvgrid",&ntvgrid)) { 
        if(ierr==0) {
            getgval(&gh,"n1",&tmp);
            ntvgrid = (int)tmp;
        } else {
            err(" ntvgrid missing ");
        }  
    }
    if (!getparint("ncdpline",&ncdpline)) ncdpline = ncdpvgrid;
    if (!getparint("cdp1",&cdp1)) cdp1 = 1;
    if (!getparint("trace1",&trace1)) trace1 = 1;
    if (!getparint("line1",&line1)) line1 = 1;
    if (!getparint("dline",&dline)) dline = 1;

    if (!getparfloat("fcdphvel",&fcdphvel)) fcdphvel = fcdpvgrid;
    if (!getparfloat("dcdphvel",&dcdphvel)) dcdphvel = dcdpvgrid;
    if (!getparint("ncdphvel",&ncdphvel)) ncdphvel = 1;
    if (!getparfloat("flinehvel",&flinehvel)) flinehvel = flinevgrid;
    if (!getparfloat("dlinehvel",&dlinehvel)) dlinehvel = dlinevgrid;
    if (!getparint("nlinehvel",&nlinehvel)) nlinehvel = 1;
    if (!getparfloat("fthvel",&fthvel)) fthvel = ftvgrid;
    if (!getparfloat("dthvel",&dthvel)) dthvel = 10. * dtvgrid;
    if (!getparint("nthvel",&nthvel)) nthvel = ntvgrid/10;

    if (!getparint("ivtype",&ivtype)) ivtype=0;
    if (!getparint("ovtype",&ovtype)) ovtype=0;
    if (!getparint("ittype",&ittype)) ittype=0;
    if (!getparint("ottype",&ottype)) ottype=0;
    if (!getparint("hcdptype",&hcdptype)) hcdptype=0;

    if(getparstring("template",&template)) {
                tmpfp = efopen(template,"r");
                if (!getparint("nvfmax",&nvfmax)) nvfmax = 4096;
                if (!getparint("ntvmax",&ntvmax)) ntvmax = 256;
        if (!getparint("nthvel",&nthvel)) nthvel= 4096; 

        /* arrays used to store all VS3D card's x,y,time and velocity */
        cdp = (int*)emalloc(nvfmax*sizeof(int));
        nps = (int*)emalloc(nvfmax*sizeof(int));
        ts = (float*) emalloc(nvfmax*ntvmax*sizeof(float));
        vs = (float*) emalloc(nvfmax*ntvmax*sizeof(float));
        ncdp = 0;
        bzero(nps,nvfmax*sizeof(int));
        hvelread(tmpfp,cdp,ts,vs,&ncdp,nps,ntvmax,nvfmax);
    }

    if(!getparfloat("nulltime",&nulltime)) nulltime = 0.0;
    twb = (float*) malloc(ncdpvgrid*nlinevgrid*sizeof(float));
    tst = (float*) malloc(ncdpvgrid*nlinevgrid*sizeof(float));
    tsb = (float*) malloc(ncdpvgrid*nlinevgrid*sizeof(float));
    for(it=0;it<ncdpvgrid*nlinevgrid;it++) {
        twb[it] = nulltime;
        tst[it] = nulltime;
        tsb[it] = nulltime;
    }
    if(getparstring("wbgrid",&wbgrid)) {
        wbfp = efopen(wbgrid,"r");
        efread(twb,sizeof(float),ncdpvgrid*nlinevgrid,wbfp);
    }
    if(getparstring("stgrid",&stgrid)) {
        stfp = efopen(stgrid,"r");
        efread(tst,sizeof(float),ncdpvgrid*nlinevgrid,stfp);
    }
    if(getparstring("sbgrid",&sbgrid)) {
        sbfp = efopen(sbgrid,"r");
        efread(tsb,sizeof(float),ncdpvgrid*nlinevgrid,sbfp);
    }
    
    ierr = fgetghdr(infp,&gh);

    nthvel0 = nthvel;
    tin = (float*)malloc(ntvgrid*sizeof(float));
    vin = (float*)malloc(ntvgrid*sizeof(float));
    tout0 = (float*)malloc(nthvel*sizeof(float));
    tout = (float*)malloc((nthvel+3)*sizeof(float));
    vout = (float*)malloc((nthvel+3)*sizeof(float));

    for(it=0;it<ntvgrid;it++) tin[it] = ftvgrid + it*dtvgrid;

    if(ncdp==1) {
        nthvel = nps[0];
        for(it=0;it<nthvel;it++) tout[it] = ts[it];
    } else if(ncdp==0) {
        for(it=0;it<nthvel;it++) tout0[it] = fthvel + it*dthvel;
    }

    if(ncdp>1) {
        for(icdp=0;icdp<ncdp;icdp++) {
            ocdp = cdp[icdp];
            if(hcdptype==0) {
                io3 = (ocdp-cdp1)/ncdpline*dline + line1;
                /* io2 = ocdp - (io3-line1)/dline*ncdpline + trace1 - 1; */
                io2 = ocdp - (io3-line1)/dline*ncdpline;
            } else {
                tmp = ocdp/10000;
                io3 = tmp;
                tmp = ocdp - io3*10000 + 0.5;
                io2 = tmp;
            }
            tmp = (io3 - flinevgrid)/dlinevgrid + 0.5;
            i3vgrid = (int) tmp;
            tmp = (io2 - fcdpvgrid)/dcdpvgrid + 0.5;
            i2vgrid = (int) tmp;
            if(i3vgrid<0) i3vgrid=0;
            if(i3vgrid>=nlinevgrid) i3vgrid=nlinevgrid-1;
            if(i2vgrid<0) i2vgrid=0;
            if(i2vgrid>=ncdpvgrid)  i2vgrid=ncdpvgrid-1;
            lpos = i2vgrid+i3vgrid*ncdpvgrid;
            lpos = lpos * ntvgrid*sizeof(float);
            fseek64(infp,lpos,0);
            fread(vin,sizeof(float),ntvgrid,infp);

            nthvel = nps[icdp];
            for(it=0;it<nthvel;it++) tout[it] = ts[it+icdp*ntvmax];
            /* insert times of wbgrid, stgrid, sbgrid */
            inserttime(tout,&nthvel,twb[i2vgrid+i3vgrid*ncdpvgrid],nulltime);
            inserttime(tout,&nthvel,tst[i2vgrid+i3vgrid*ncdpvgrid],nulltime);
            inserttime(tout,&nthvel,tsb[i2vgrid+i3vgrid*ncdpvgrid],nulltime);

            /* time/depth conversion if needed */
            vconvert(tin,vin,ntvgrid,ivtype,ittype,
                     tout,vout,nthvel,ovtype,ottype);

        /*
        for(it=0;it<ntvgrid;it++) 
            fprintf(stderr,"tin=%f vin=%f \n",tin[it],vin[it]);
        for(it=0;it<nthvel;it++) 
            fprintf(stderr,"tout=%f vout=%f \n",tout[it],vout[it]);
        */

        if( !strcmp( outfmt ,"hvel" )) {
           printhvel(ocdp,nthvel,tout,vout,outfp);

        }else if( !strcmp( outfmt ,"avf" )){

		   sl2xydb(lmk_trace1,lmk_line1,lmk_x1,lmk_y1
                  ,lmk_trace2,lmk_line2,lmk_x2,lmk_y2
			  ,lmk_trace3,lmk_line3,lmk_x3,lmk_y3
                  ,(double)io2,(double)io3,&xout,&yout);
           if( vout[1] > 0.0 ){
              printavf(ocdp,xout,yout,nthvel,tout,vout,outfp);
           }


        }else{
            fprintf( stderr ,"Unknown output format: %s\n" ,outfmt );
            exit(-1);
        }
        if(hcdptype==0) { 
fprintf(stderr,"Output at trace=%d line=%d cdp=%d i2vg=%d i3vg=%d\n",
                io2,io3,ocdp,i2vgrid+1,i3vgrid+1);
        } else {
fprintf(stderr,"Output at trace=%d line=%d cdplbl=%d i2vg=%d i3vg=%d\n",
                io2,io3,ocdp,i2vgrid+1,i3vgrid+1);
        }
        }
    } else {
        for(i3=0;i3<nlinehvel;i3++) {
            o3 = flinehvel+i3*dlinehvel;
            io3 = o3;
            tmp = (o3 - flinevgrid)/dlinevgrid + 0.5;
            i3vgrid = (int) tmp;
            if(i3vgrid<0) i3vgrid=0;
            if(i3vgrid>=nlinevgrid) i3vgrid=nlinevgrid-1;
            for(i2=0;i2<ncdphvel;i2++) {
                o2 = fcdphvel+i2*dcdphvel;
                io2 = o2;
                tmp = (o2 - fcdpvgrid)/dcdpvgrid + 0.5;
                i2vgrid = (int) tmp;
                if(i2vgrid<0) i2vgrid=0;
                if(i2vgrid>=ncdpvgrid)  i2vgrid=ncdpvgrid-1;

                if(hcdptype==0) {
                    tmp = o2 - trace1 + (o3-line1)/dline*ncdpline + cdp1 + 0.5; 
                    ocdp = tmp;
                } else {
                    ocdp = io3*10000 + io2;
                }
                lpos = i2vgrid+i3vgrid*ncdpvgrid;
                lpos = lpos * ntvgrid*sizeof(float);
                fseek64(infp,lpos,0);
                efread(vin,sizeof(float),ntvgrid,infp);

                nthvel = nthvel0;
                for(it=0;it<nthvel;it++) tout[it] = tout0[it];
                /* insert times of wbgrid, stgrid, sbgrid */
                inserttime(tout,&nthvel,
                    twb[i2vgrid+i3vgrid*ncdpvgrid],nulltime);
                inserttime(tout,&nthvel,
                    tst[i2vgrid+i3vgrid*ncdpvgrid],nulltime);
                inserttime(tout,&nthvel,
                    tsb[i2vgrid+i3vgrid*ncdpvgrid],nulltime);
                /* time/depth conversion if needed */
                vconvert(tin,vin,ntvgrid,ivtype,ittype,
                        tout,vout,nthvel,ovtype,ottype);

                if( !strcmp( outfmt ,"hvel" )) {
                   printhvel(ocdp,nthvel,tout,vout,outfp);

                }else if( !strcmp( outfmt ,"avf" )){

		   sl2xydb(lmk_trace1,lmk_line1,lmk_x1,lmk_y1
                          ,lmk_trace2,lmk_line2,lmk_x2,lmk_y2
			  ,lmk_trace3,lmk_line3,lmk_x3,lmk_y3
                          ,(double)io2,(double)io3,&xout,&yout);
                   if( vout[1] > 0.0 ){
                      printavf(ocdp,xout,yout,nthvel,tout,vout,outfp);
                   }


                }else{
                    fprintf( stderr ,"Unknown output format: %s\n" ,outfmt );
                    exit(-1);
                }
                fprintf(stderr
                ,"Output at trace=%d line=%d cdp=%d i2vg=%d i3vg=%d\n"
                ,io2,io3,ocdp,i2vgrid+1,i3vgrid+1);
            }
        }
    }
    
    free(tout0);
    free(tout);
    free(vout);
    free(tin);
    free(vin);
    free(twb);
    free(tst);
    free(tsb);
}

void inserttime(float *ts, int *np, float tgrid, float nulltime) { 

    int i, j;
    int ns;

    ns = *np;

    if(tgrid==nulltime) return;

    if(tgrid<ts[0]) {
        for(j=ns;j>0;j--) ts[j] = ts[j-1];
        ts[0] = tgrid;
        *np += 1;
    } else if(tgrid > ts[ns-1]) {
        ts[ns] = tgrid;
        *np += 1;
    } else {
        for(i=0;i<ns-1;i++) {
            if((tgrid>ts[i])&&(tgrid<ts[i+1])) {
                for(j=ns;j>i+1;j--) ts[j] = ts[j-1];
                ts[i+1] = tgrid;
                *np +=  1;
                break;
            }
        }
    }
}