kmmerge.c

seismic software,very useful

                tmp = (cdp2 - cdp1)/(ncdppl-1.) + .5;
                cdppds = tmp;
        } else {
                cdppdds = nlines;
                cdppds = nlines;
        }
        if(cdppdds<1) cdppdds = 1;
        if(cdppds<1) cdppds = 1;
        if(nlines>1) {
                tmp = (cdp3 - cdp1)/(nlines-1.) + .5;
                cdppddl = tmp;
                tmp = (cdp3 - cdp1)/(nlines-1.) + .5;
                cdppdl = tmp;
        } else {
                cdppddl = ncdppl;
                cdppdl = ncdppl;
        }
        if(cdppddl<1) cdppddl = 1;
        if(cdppdl<1) cdppdl = 1;
	if(cdppdds!=1 && cdppddl!=1 ) {
                warn("cdppdds=%d cdppddl=%d ",cdppdds,cdppddl);
                err("check cdp1, cdp2, cdp3 ");
        }
        if(cdppdds==1) {
                cdpnum = 0;
        } else {
                cdpnum = 1;
        }
	if ( getparfloat("send",&send)
          && getparfloat("lend",&lend)
          && getparint("ntrend",&ntrend) ) {
                if(ntrend<1) err("ntrend must be greater than 0 ");
                nsout = ntrend;
                nlout = 1;
        } else {
                ntrend = 0;
                nsout = ns;
                nlout = nl;
        }
	ntmp = nl*ndl;
        nl = ntmp;
        ntmp = ns*nds;
        ns = ntmp;

        if(ntrend==0) {
                nsout = ns;
                nlout = nl;
        } else {
                nsout = ntrend;
                nlout = 1;
	}

	s = (float*) emalloc(nsout*nlout*sizeof(float));
        l = (float*) emalloc(nsout*nlout*sizeof(float));

	if(ntrend==0) {
                for(il=0;il<nlout;il++) {
                        for(is=0;is<nsout;is++) {
                                s[is+il*nsout] = sstart+is/nds*ds+is%nds*dds;
                        }
                }
                for(il=0;il<nlout;il++) {
                        for(is=0;is<nsout;is++) {
                                l[is+il*nsout] = lstart+il/ndl*dl+il%ndl*ddl;
                        }
                }
        } else {
                for(is=0;is<ntrend;is++) {
                        s[is] = sstart+is*(send-sstart)/(ntrend-1.);
                        l[is] = lstart+is*(lend-lstart)/(ntrend-1.);
                }
        }




	if(!getparint("flexbin",&flexbin)) flexbin = 0;

	
	if(kzorkt==1) {
		nzo = ntau;	
		nz = ntau;
		fzo = tau0;
		dzo = dtau;
		zpow = -tpow;
	}

	mig = (float*) emalloc(nzo*sizeof(float));
	migs = (float*) emalloc(nzo*nsout*nofo*sizeof(float));
	fold = (float*) emalloc(nz*sizeof(float));
	folds = (float*) emalloc(nz*nsout*nofo*sizeof(float));
	zgain = (float*) emalloc(nzo*sizeof(float));

	fprintf(stderr," nofo=%d nsout=%d nlout=%d \n",nofo,nsout,nlout);
	fprintf(stderr," start output ... \n");

	for(iz=0;iz<nzo;iz++) {
		tmp = fzo+iz*dzo;
		if(zpow!=0.0) {
			if(tmp!=0.0) {
				zgain[iz] = pow(tmp,zpow);
			} else {
				zgain[iz] = 0.0;
			}
		} else {
			zgain[iz] = 1.0;
		}
	}

	for(il=0;il<nlout;il++) {
		bzero(migs,nzo*nsout*nofo*sizeof(float));
		bzero(folds,nz*nsout*nofo*sizeof(float));
		for(io=0;io<nofo;io++) {
			itmp = io*nsout*nlout + il*nsout;
			for(i=0;i<ni;i++) {
				lpos =  itmp*nzo*sizeof(float);
				efseek(difp[i],lpos,0);
				for(is=0;is<nsout;is++) {
					efread(mig,sizeof(float),nzo,difp[i]);
					itmp1 = (is*nofo+io)*nzo;
					for(iz=0;iz<nzo;iz++) {
						migs[itmp1+iz] += mig[iz];
					}
				}
				lpos =  itmp*nz*sizeof(float);
				efseek(dffp[i],lpos,0);
				for(is=0;is<nsout;is++) {
					efread(fold,sizeof(float),nz,dffp[i]);
					itmp1 = (is*nofo+io)*nz;
					for(iz=0;iz<nz;iz++) {
						folds[itmp1+iz] += fold[iz];
					}
				}
			}
		}
		for(is=0;is<nsout;is++) {
			for(io=0;io<nofo;io++) {
				for(i=0;i<ni;i++) {
					efread(&tr,sizeof(char),240,dhfp[i]);
					if(i==0) {
						bcopy(&tr,&tra,240);
					} else {
						if(tr.ungpow<tra.ungpow) 
							bcopy(&tr,&tra,240);
						if(tr.trid==1 && tra.trid!=1) 
							bcopy(&tr,&tra,240);
					}
				}
				itmp = (is*nofo+io)*nzo;
				itmp1 = (is*nofo+io)*nz;
				if(flexbin==0) {
					if(kzorkt==0) {
                                        	for(iz=0;iz<nz;iz++)
                                                	fold[iz] =
                                                	sqrt(fold0+
							folds[itmp1+iz]);
                                        	for(iz=0;iz<nzo;iz++){
                                                	zi=(fzo+iz*dzo-fz)/dzt;
                                                	izi = zi;
                                                	if(izi>nz-2) izi = nz-2;
                                                	scale = (1.0-zi+izi)*
                                                        	fold[izi]+
                                                        	(zi-izi)*
                                                        	fold[izi+1];
                                                	tra.data[iz] = 
								migs[itmp+iz]/
                                                        	scale;
                                        	}
					} else {
						for(it=0;it<ntau;it++) {
							scale = folds[itmp+it];
							if(scale>1.) {
                                                	tra.data[it] = 
								migs[itmp+it]/
                                                        	scale;
							} else {
                                                	tra.data[it] = 
								migs[itmp+it];
							}
						}
					}
                                } else {
                                        for(iz=0;iz<nzo;iz++)
                                                tra.data[iz] = migs[itmp+iz];
                                }

				for(iz=0;iz<nzo;iz++)
					tra.data[iz] *= zgain[iz];

				/* update header */
				tra.offset = ofo[io];
                                tra.tracl = is + 1;
                                tra.tracr = il + 1;
                                tra.cdpt = io + 1;
				if(kzorkt==0) {
                                	tra.dz = dzo;
                                	tra.fz = fzo;
					tra.ns = nzo;
					tra.mute = 0;
					tra.f1 = fzo;
					tra.d1 = dzo;
				} else {
					tra.delrt = tau0;
                                	tra.ns = ntau;
                                	tra.dt = dtau * 1000000.;
					if(tra.trid!=1) {
                                        tra.d1 = tra.dt;
                                        tmp = (tau0 + (ntau-1)*dtau)*1000;
                                        tra.mute = tmp;
					}
				}
                                ofs = ofo[io];
                                if(tra.trid==1) {
                                        mx = (tra.gx+tra.sx)/2;
                                        my = (tra.gy+tra.sy)/2;
                                        tmp = tra.gx-tra.sx;
                                        dd = tra.gy-tra.sy;
                                        tmp = tmp*tmp + dd*dd;
                                        dd = sqrt(tmp);
                                        if(dd>0.) {
                                                tra.sx = mx+ofs*(tra.sx-mx)/dd;
                                                tra.gx = mx+ofs*(tra.gx-mx)/dd;
                                                tra.sy = my+ofs*(tra.sy-my)/dd;
                                                tra.gy = my+ofs*(tra.gy-my)/dd;
                                        } else {
						tra.sx = mx-ofs/2;
                                                tra.gx = mx+ofs/2;
                                                tra.sy = my;
                                                tra.gy = my;
                                        }
                                }


				if(ntrend>0) {
                                        tmp = (l[il*nsout+is] - l1)/ddl + 0.5;
                                        iiy = tmp;
                                        tmp = (s[il*nsout+is] - s1)/dds + 0.5;
                                        iix = tmp;
                                        if(cdpnum==0) {
                                                tra.cdp = iiy*ncdppl+iix+cdp1;
                                        } else {
                                                tra.cdp = iix*nlines+iiy+cdp1;
                                        }
                                }

                                /* output */
                                fputtr(outfp,&tra);
			}
		}
	}
	
	fprintf(stderr," merge output done \n");

	free(migs);
	free(mig);
	free(folds);
	free(fold);
	
	return 0;
}