sustaticrrs.c

su 的源代码库

        	efclose(fpr);

		/* open and read from source statics file */
        	if((fps=efopen(sou_file,"rb"))==NULL)
                	err("cannot open stat_file=%s\n",sou_file);
        	efread(sou_statics, sizeof(float),ns,fps);
        	efclose(fps);
	}

        if (hdrs==3){
                getparstring("blvl_file",&blvl_file);
                getparstring("refr_file",&refr_file);
                getparint("nsamp",&nsamp);
                getparfloat("fx",&fx);
                getparfloat("dx",&dx);
                getparfloat("V_r",&V_r);

                blvl_arr = alloc1float(nsamp);
                if((fpb=efopen(blvl_file,"rb"))==NULL)
                        err("cannot open blvl_file=%s\n",blvl_file);
                efread(blvl_arr, sizeof(float),nsamp,fpb);

                refr_arr = alloc1float(nsamp);
                if((fpf=efopen(refr_file,"rb"))==NULL)
                        err("cannot open refr_file=%s\n",refr_file); 
                efread(refr_arr, sizeof(float),nsamp,fpf);

                g1=fopen("tcdp.b","wb");
                g2=fopen("num.b","wb");

                if (!getparstring("vfile",&vfile)) err("must specify vfile");
                if (!getparint("mx",&mx)) err("must specify mx");
                if (!getparint("mz",&mz)) err("must specify mz");
                if (!getparfloat("dzv",&dzv)) err("must specify dzv");
                vfp = fopen(vfile,"r");
                vel = alloc2float(mz,mx);
                fread(vel[0],sizeof(float),mz*mx,vfp);
                fclose(vfp);

                tcdp_stat = alloc1float(nsamp);
                num = alloc1float(nsamp);
                for(i=0;i<nsamp;i++){
                      tcdp_stat[i] = 0.0;
                      num[i] = 0.0;
                }
        }

        if (hdrs==4){
                getparint("nsamp",&nsamp);
                getparfloat("fx",&fx);
                getparfloat("dx",&dx); 
                tcdp_arr = alloc1float(nsamp);
                if((g3=efopen("tcdp.b","rb"))==NULL)
                        err("cannot open tcdp.b\n");
                efread(tcdp_arr, sizeof(float),nsamp,g3);

                num_arr = alloc1float(nsamp);
                if((g4=efopen("num.b","rb"))==NULL)
                        err("cannot open num.b\n");
                efread(num_arr, sizeof(float),nsamp,g4);

        }

	/* Loop on traces */	
	io = 0; is = 0;
	do {
		int temp = SGN(intrace.scalel)*log10(abs((int) intrace.scalel));
		float scale;
                scale = pow(10., (float)temp);
		
		/* copy and adjust header */
		memcpy( (void *) &outtrace, (const void *) &intrace, HDRBYTES);
	
		/* compute static correction if necessary */
		if(!hdrs) {
		    	tsd = scale *
			(-intrace.selev + intrace.sdel + intrace.sdepth)/v1;
			trd = tsd - intrace.sut/1000.0;
			tstat = tsd + trd +
				scale * (intrace.selev - intrace.gelev)/v0;

		/* else, read statics from headers */
		} else { 

			if (hdrs==2) {
				ir = is + io;
				if (is < ns) tsd = sou_statics[is]/1000.0;
				if (ir >= 0 && ir < nr)
					trd = rec_statics[ir]/1000.0;

				intrace.tstat =  NINT(tsd*1000.0 +  trd*1000.0);
				io ++;
				if (io > no-1) {
					io = 0; is++;
				}
			}

                        if(hdrs==3){
                              is_c = (int)((((float)intrace.sx)-fx)/dx);
                              z_blvl_s = blvl_arr[is_c];
                              z_refr_s = refr_arr[is_c];
                              
                              ir_c = (int)((((float)intrace.gx)-fx)/dx);
                              z_blvl_r = blvl_arr[ir_c];
                              z_refr_r = refr_arr[ir_c];

                              tsd = 0.0;
                              trd = 0.0;
                              is_c = (int)((((float)intrace.sx)-fx)/dx);
                              ir_c = (int)((((float)intrace.gx)-fx)/dx);
                              for(j=0;j<mz-1;j++){
                                 tsd = tsd + dzv/vel[is_c][j];
                                 trd = trd + dzv/vel[ir_c][j];
                              }
                              tsd = tsd - (z_blvl_s-z_refr_s)/V_r;
                              trd = trd - (z_blvl_r-z_refr_r)/V_r;

                              intrace.tstat = NINT(tsd*1000.0+trd*1000.0);

                              ic_c = (int)((((float)intrace.cdp)-fx)/dx);
                              tcdp_stat[ic_c] = tcdp_stat[ic_c] + tsd + trd;
                              num[ic_c] = num[ic_c] + 1.0;
                        }

                        if(hdrs==4){
                              ic_c = (int)((((float)intrace.cdp)-fx)/dx);
                              if(num_arr[ic_c]>=0.5){
                                 bulk = NINT(1000.0*tcdp_arr[ic_c]/num_arr[ic_c]);
                                 res = intrace.tstat - bulk;
                              }
                              if(num_arr[ic_c]<=0.5){
                                 bulk = 0;
                                 res = 0;
                              }
                        }


			/* if total statics not supplied, calculate it */
			if(intrace.tstat==0 && hdrs!=3 && hdrs!=4){
				outtrace.tstat = intrace.sstat+intrace.gstat;
				tstat = outtrace.tstat/1000.0;
			} 
                        /* for hdrs==5, here is where tstat is calculated */
                        if(intrace.tstat!=0 && hdrs!=3 && hdrs!=4){
				tstat = intrace.tstat/1000.0;
			}
                        if(hdrs==3){
                                tstat = 0.0;
                        }
                        if(hdrs==4){
                                tstat = res/1000.0;
                        }
		}
		
		/* Compute output times */
		for (itime=0; itime<nt; ++itime)
			t[itime] = tmin + itime*dt + sign*tstat;

		/* sinc interpolate new data */
		ints8r(nt, dt, tmin, intrace.data, 
				0.0, 0.0, nt, t, outtrace.data);
		
		/* set header field for output trace */
		if(hdrs!=1 && hdrs!=3 && hdrs!=4 && hdrs!=5){

			/* value is added to existing header values */
			/* this permits multiple static corrections */
			outtrace.sstat += (1000.0 * tsd);
			outtrace.gstat += (1000.0 * trd);
			outtrace.tstat += (1000.0 * tstat);
		}
                if(hdrs==3){
                        outtrace.sstat = 0;
                        outtrace.gstat = 0;
                        outtrace.tstat = intrace.tstat;
                }
                if(hdrs==4){
                        outtrace.sstat = 0;
                        outtrace.gstat = res;
                        outtrace.tstat = bulk;
                }
                if(hdrs==5){
                        outtrace.tstat = intrace.tstat;
                        outtrace.sstat = intrace.sstat;
                        outtrace.gstat = intrace.gstat;
                }
		
		puttr(&outtrace);
	} while (gettr(&intrace));

        if(hdrs==3){
           fwrite(tcdp_stat,sizeof(float),nsamp,g1);
           fwrite(num,sizeof(float),nsamp,g2);
           fclose(fpb);
           fclose(fpf);
           fclose(g1);
           fclose(g2);
        }
        if(hdrs==4){
           fclose(g3);
           fclose(g4);
        }

	return(CWP_Exit());
}