kmmerge.c

seismic software,very useful

        if(ddl==0.) { ddl = 1.0; dl = ddl; }
	if(ncdppl>1) {
                tmp = (cdp2 - cdp1)/(ncdppl-1.) + .5;
                cdppdds = tmp;
                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;
	}

	if ( getparstring("tracekey",&tracekey) && getparstring("linekey",&linekey) )
		ikey=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 = 1;
		fzo = tau0;
		dzo = dtau;
		zpow = tpowaf;
	}

	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);
	if(kzorkt==0) fprintf(stderr," nzo=%d dzo=%g fzo=%g \n",nzo,dzo,fzo);
	if(kzorkt==1) fprintf(stderr," ntau=%d dtau=%g tau0=%g \n",
		ntau,dtau*1000,tau0*1000);
	fprintf(stderr," start output ... \n");

	bh.hns = nzo;
	itmp = dzo*1000;
	if(itmp>32000) itmp = itmp/1000;
	bh.hdt = itmp;
	if(kzorkt==1) bh.hdt = dtau * 1000000;
	bh.fold = nofo;
	if(nofo>1) {
		bh.tsort = 2;
	} else {
		bh.tsort = 4;
	}

	if(noinput==1) idhdrs(&ch,&bh,nzo);
	fputhdr(outfp,&ch,&bh);


	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 && tr.trid==1 && tra.trid==1) 
							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 {
						scale = folds[itmp1];
						if(scale>1.) {
						    for(it=0;it<ntau;it++)
                               	tra.data[it] = migs[itmp+it]/scale;
						} else {
						    for(it=0;it<ntau;it++)
                                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];
				if(ikey!=1) {
                	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;
					if(dzo>32) { 
						tra.dt = dzo;
					} else {
						tra.dt = dzo*1000.+0.5;
					}
					if(ihz==0) tra.delrt = fzo;
				} else {
					if(ihz==0) tra.delrt = tau0*1000;
                   	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];
                ms = s[il*nsout+is];
		ml = l[il*nsout+is];
		sl2xy(s1,l1,x1,y1,s2,l2,x2,y2,s3,l3,x3,y3,ms,ml,&mx,&my);

                if(tra.trid==1) {
		    if(tra.scalco>1) {
			tra.gx = tra.gx * tra.scalco;
			tra.gy = tra.gy * tra.scalco;
			tra.sx = tra.sx * tra.scalco;
			tra.sy = tra.sy * tra.scalco;
		    } else if(tra.scalco<0) {
			tra.gx = tra.gx / (-tra.scalco);
			tra.gy = tra.gy / (-tra.scalco);
			tra.sx = tra.sx / (-tra.scalco);
			tra.sy = tra.sy / (-tra.scalco);
		    }
                    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;
                    }
                } else {
		       tra.sx = mx-ofs/2;
                       tra.gx = mx+ofs/2;
                       tra.sy = my;
                       tra.gy = my;
		}
		tra.scalco = 0;

/*
				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;
}