kirmig.c

seismic software,very useful

		}
	}
	if (!getparint("nxt", &nxt)) {
		if(ierr==0) {
			getgval(&gh,"n2",&fhread);
			fhread = fhread + 0.5;
			nxt = (int) fhread;
		} else {
			nxt = nx;
		}
	} else {
		if(ierr==0) {
			getgval(&gh,"n2",&fhread);
			fhread = fhread + 0.5;
			ihread = (int) fhread;
			if(ihread!=nxt) err(" check nxt & timefile ");	
		}
	}
	if (!getparint("nst", &nst)) {
		if(ierr==0) {
			getgval(&gh,"n3",&fhread);
			fhread = fhread + 0.5;
			nst = (int) fhread;
		} else {
			nst = nx;
		}
	} else {
		if(ierr==0) {
			getgval(&gh,"n3",&fhread);
			fhread = fhread + 0.5;
			ihread = (int)fhread;
			if(ihread!=nst) err(" check nst & timefile ");	
		}
	}

	if (smint>x0) warn("=====> smint > x0 <=====\n");
        if (zmint>z0) warn("=====> zmint > z0 <=====\n");
        if (xmint>x0) warn("=====> xmint > x0 <=====\n");
        if ((smint+(nxt-1)*dxt)<(x0+(nx-1)*dx))
         warn ("=====> (smint+(nxt-1)*dxt) < (x0+(nx-1)*dx) <=====\n");
        if ((xmint+(nxt-1)*dxt)<(x0+(nx-1)*dx))
         warn ("=====> (xmint+(nxt-1)*dxt) < (x0+(nx-1)*dx) <=====\n");
        if ((zmint+(nzt-1)*dzt)<(z0+(nz-1)*dz))
         warn ("=====> (zmint+(nzt-1)*dzt) < (z0+(nz-1)*dz) <=====\n");

	if (!getparint("nm", &nm)) nm = 1;   
	if (!getparint("i2p5", &i2p5)) i2p5 = 1;   
	if (!getparint("intype", &intype)) intype = 0 ;   
	if (!getparfloat("dcdp", &dcdp)) dcdp = 0. ;   
	if (!getparint("sort", &sort)) sort = 0 ;   
	if (!getparint("amptype", &amptype)) amptype = 0 ;   
	if (!getparint("mlimit", &mlimit)) mlimit = 48 ;   
	if (!getparint("fcdpout", &fcdpout)) fcdpout = 1 ;   
	if (!getparint("dcdpout", &dcdpout)) dcdpout = 1 ;   
	if (!getparint("lagc", &lagc)) lagc = 0 ;
        if (!getparfloat("tpow", &tpow)) {
                if(lagc==0) {
                        tpow = 0. ;
                } else {
                        tpow = -2. ;
                }
        }
        if (!getparfloat("aper", &aper)) aper=nx*dx/4.;
        if (!getparfloat("apan", &apan)) apan=60.;
	if (!getparfloat("apanl", &apanl)) apanl=-apan;
        if (!getparfloat("apanr", &apanr)) apanr=apan;

	if (apanl<-89.9) apanl = -89.9;
        if (apanl>89.9) apanl = 89.9;
        if (apanl>-0.1 && apanl <0.1) {
		if(apanl>=0) {
			apanl = 0.1;
		} else {
			apanl = -0.1;
		}
	}
	if (apanr<-89.9) apanr = -89.9;
        if (apanr>89.9) apanr = 89.9;
        if (apanr>-0.1 && apanr <0.1) {
		if(apanr>=0) {
			apanr = 0.1;
		} else {
			apanr = -0.1;
		}
	}

        if (!getparfloat("v0", &v0)) v0=1500.;
        if (!getparfloat("dxtol", &dxtol)) dxtol=1.;
	if(dxtol<0.) dxtol = - dxtol;


/* update binary header */
        sgybin.hns = nz ;
        sgybin.ntrpr = nx ;
        sgybin.hdt = (short) ( dz < 1. ? dz*1000. : dz );
        sgybin.fold = nofo;
        sgybin.tsort = 3;
	
	
        fputhdr(outfp,&sgycrd,&sgybin) ;

	/* set length of migration trace to lt to intepolate input traces */
	if (!getparint("lt", &lt)) lt = 4*nt;   
	if (lt<nt) lt=nt;
        lt = lt/nt*nt;
	
/* number of offsets per subprocess to migrate */
	iofpm = (nofo+nm-1)/nm;

  	intpx = 1;
        if(x0==xmint && dx==dxt && nx<=nxt) intpx=0;
        intpz = 1;
        if(z0==zmint && dz==dzt && nz==nzt) intpz=0;
        intps = 1;
        if(x0==smint && dx==dst && nx<=nst) intps=0;
	
	fprintf(stderr,"x0=%f z0=%f nx=%d nz=%d \n",x0,z0,nx,nz);
        fprintf(stderr,"dx=%f dz=%f \n",dx,dz);
        fprintf(stderr,"timefile=%s \n",timefile);
        fprintf(stderr,"ampfile=%s \n",ampfile);
        fprintf(stderr,"nt=%d dt=%f tmin=%f \n",nt,dt,tmin);
        fprintf(stderr,"ofomin=%f dofo=%f nofo=%d \n",ofomin,dofo,nofo);
        fprintf(stderr,"xmint=%f zmint=%f smint=%f \n",xmint,zmint,smint);
        fprintf(stderr,"dxt=%f dzt=%f dst=%f \n",dxt,dzt,dst);
        fprintf(stderr,"nxt=%d nzt=%d nst=%d \n",nxt,nzt,nst);
        fprintf(stderr,"lt=%d i2p5=%d \n",lt,i2p5);
	fprintf(stderr,"intps=%d intpx=%d intpz=%d\n",intps,intpx,intpz);


	tmp = nx * dx / dxdip;
	nxdip = (int)tmp;
	/* compute memory requirements */
	lmem = (2*nzt*nxt + iofpm*(nx*nz+1) + 4*nx*nz + nx + 2*nz
                + nzt + 2*lt + nt)*sizeof(float)
              +(nx + nz) * sizeof(int) + idip*2*nxdip*nz*sizeof(float);
	lmem = lmem + nzt*nxt*nst*sizeof(short) +
		iamp * (nzt*nxt*nst*sizeof(short)+nzt*(nxt+1)*sizeof(float));

        fprintf(stderr," \n");
        fprintf(stderr," KIRMIG needs memory size=%d \n",lmem);
        fprintf(stderr," \n");
	
        llimit = mlimit * 1024 * 1024;
        if ( lmem > llimit ) {
	   fprintf(stderr,"travel time table disk i/o used \n");
	   if(iamp==1) fprintf(stderr,"amplitude table disk i/o used \n");
	   lmem = lmem - nzt*nxt*nst*sizeof(short) -
		iamp * (nzt*nxt*nst*sizeof(short)+nzt*(nxt+1)*sizeof(float));
           if ( lmem > llimit ) {
           	err("Need at least memory size mlimit=%d (Bytes)\n",lmem);
	   }
        }



	mtrace = (int *)malloc(nm*sizeof(int));
	lseekini = (int *)malloc(nm*sizeof(int));
	dataini = (char**) malloc(nm*sizeof(char *));
	dataouti = (char**) malloc(nm*sizeof(char *));

       	for(im=0;im<nm;im++) {
              dataouti[im] = (char*) malloc(strlen(dataout)+10);
              sprintf(dataouti[im],"%s%d",dataout,im+1) ;
        }

/* test to see if output offsets were created in the previous run */
        for(im=0;im<nm;im++) {
              nofoi = iofpm;
              kofo = im * iofpm + 1;
              if (kofo+nofoi-1 > nofo ) nofoi = nofo - kofo + 1;
              isize=0;
              if((infp = fopen(dataouti[im],"r"))!=NULL) {
                 fseek(infp,0L,SEEK_END);
                 isize= (int) ftell(infp);
                 fclose(infp);
              }
              if( isize!=(nz*sizeof(float)+HDRBYTES)*nx*nofoi ) {
                        isize = 0;
                        break;
              }
        }
        if(isize!=0) goto no_migration;


/* sort data into common offsets for each machine*/
	if ( nm > 1 && intype == 0 ) {
	   FILE ** fdataini  ;

	   fdataini = (FILE**) malloc(sizeof(FILE *)*nm);
	   for(kofo=1;kofo<=nm;kofo++) {
	      dataini[kofo-1] = (char*) malloc(strlen(datains)+10);
	      sprintf(dataini[kofo-1],"%s%d",datains,kofo) ;
	      if(sort==0 || sort==2) 
			fdataini[kofo-1] = efopen(dataini[kofo-1],"w") ;
           }

	   if ( sort == 0 || sort==2 ) {
	      infp = efopen(datain,"r");
	      fseek(infp,EBCBYTES+BNYBYTES,SEEK_SET);
	      while( fgettr(infp,&tra) ) {
	         float offset = tra.offset;
                 /* offset = tra.gx - tra.sx; */
                 iofo =  (int) ( ( fabs(offset) - ofomin ) /dofo  + 1.5 );
	         im = (iofo+iofpm-1)/iofpm;
	
	         if(iofo>0 && iofo<=nofo) fputtr(fdataini[im-1],&tra);
	      }
	      fclose(infp);
	   }
	   for(im=0; im<nm ; im++) {
	      if(sort==0 || sort==2) fclose(fdataini[im]) ;
              lseekini[im] = 0 ;
           }
	} 
	else if(nm>1 && intype == 1 ) {

	/* find out the average number of traces per machine to mig */
	   int ntrpm, trlen,jm;
	   infp = efopen(datain,"r");
	   fseek(infp,0L,SEEK_END);
	   trlen = ntmax*sizeof(float)+HDRBYTES; 
	   ntrpm = ((int)ftell(infp)-(EBCBYTES+BNYBYTES) )/trlen/nm; 
           ntrpm -= 1 ; /* to be "in average" before the offset to point to */

	/* find out where to seek */
	   lseekini[0] = EBCBYTES+BNYBYTES;
	   for( im=0; im<nm; im++ ) {
               int before = (im==0) ? 1 : 0 ;
	       fseek(infp,lseekini[im],SEEK_SET);
	       while( fgettr(infp,&tra) ) {
	          float offset = tra.offset;
		  if (offset<0.) 
			err("Negative offsets not allowed for intype=1 \n");
		  /* offset = tra.gx - tra.sx; */
                  iofo = (int)( ( fabs(offset) - ofomin ) /dofo + .5 ) ;
	          jm = iofo/iofpm;
                  if( jm < im ) {
                     /* smaller offset: increment the seek position */
                     before = 1 ;
                     lseekini[im] += trlen ;
                  } else if( jm == im && before ) {
                     /* initialize next machine seek position */
                     if( im < nm-1 ) lseekini[im+1] = lseekini[im]+trlen*ntrpm;
                     break ;
                  } else {
                     /* bigger or equal offset: I have to go back of 2 traces */
                     lseekini[im] -= trlen ;
	             fseek(infp,-2*trlen,SEEK_CUR);
                  } 
               }
               dataini[im] = datain ;
	   }
	   fclose(infp);
	} 
	else {
           dataini[0] = datain ;
           lseekini[0] = EBCBYTES+BNYBYTES ;
        }

        kofo = 0;

        dispatcher_attributes(argv[0],0) ;
	init_dispatch();
	
	for(im=0; im<nm ; im++) {
	      kofo = im * iofpm + 1;
              ofomini = ofomin + (kofo-1)*dofo;
              mtrace[im] = 0 ;
	      nofoi = iofpm;
	      if (kofo+nofoi-1 > nofo ) nofoi = nofo - kofo + 1; 
/* temporary dataset */
	      dataouti[im] = (char*) malloc(strlen(dataout)+10);
              sprintf(dataouti[im],"%s%d",dataout,im+1) ;
	      lseekout=0;
	      isize=0;
	      if((infp = fopen(dataouti[im],"r"))!=NULL) {
	         fseek(infp,0L,SEEK_END);
		 isize= (int) ftell(infp);
		 fclose(infp);
	      }

	      if ( nofoi<= 0 ) continue;	

	      if( isize!=(nz*sizeof(float)+HDRBYTES)*nx*nofoi ) {
                 
	         fprintf(stderr,"start dispatch im=%d \n",im+1);

	         dispatch("kirmigs",dataini[im],dataouti[im],
                                    lseekini[im],lseekout,
                                    timefile,ampfile,iamp,
                                    nt,nx,nz,lt,
                                    tmin,x0,z0,
                                    dt,dx,dz,
                                    smint,xmint,zmint,
                                    dst,dxt,dzt,
                                    nst,nxt,nzt,
                                    ofimin,ofimax,
                                    ofomini,dofo,nofoi,
                                    mtrace+im,i2p5,
				    intps,intpx,intpz,intype,dcdp,
				    amptype,mlimit,lagc,tpow,aper,apanl,apanr,
				    v0,fcdpout,dcdpout,dipfile,dxdip,dxtol);
	      }

	}
	wait_dispatched();

	no_migration:

   	fprintf(stderr,"Starting merge of %d offset datasets \n",nm);

	ntrace = 0;
        for(im=0; im<nm; im++ ) {
	    if(mtrace[im]==0) {
		warn(" zero trace migrated at im=%d \n",mtrace[im]);
	    } else if(mtrace[im]<0) {
		err(" error occured at im=%d",mtrace[im]);
	    }
	    ntrace = ntrace + mtrace[im];
        }
/* remove temporary input common-offset data sets */

       if( nm > 1 && intype == 0 && sort<=1 ) {
  	   for(kofo=0;kofo<nm;kofo++) {
            	unlink(dataini[kofo]); 
            /*  eremove(dataini[kofo]); */
           }
        }

	/* if( ntrace == 0 ) return 1 ; */

/* copy data from temporary datasets to output */
	fseek(outfp,EBCBYTES+BNYBYTES,SEEK_SET);

	tracl = 1;
	for(im=0;im<nm;im++) {
	      infp = efopen(dataouti[im],"r");
	      ltrace = nz*sizeof(float) + HDRBYTES;
	      nofoi = iofpm;
	      kofo = im * iofpm + 1;
	      if (kofo+nofoi-1 > nofo ) nofoi = nofo - kofo + 1; 
	      for(iofo=0;iofo<nofoi;iofo++) {
	         for(ix=0;ix<nx;ix++) {
		    fread((char *)&tra,sizeof(char),ltrace,infp);
		    tra.cdp = fcdpout + ix * dcdpout;
		    tra.cdpt = iofo + 1;
		    tra.mute = 0;
		    tra.muts = 0;
		    tra.mutb = 0;
		    tra.dz = dz;
		    tra.fz = z0;
		    tra.tracl = tracl;
		    fwrite((char *)&tra,sizeof(char),ltrace,outfp);
		    tracl = tracl + 1;
	         }
	      }
	      fclose(infp);
	      eremove(dataouti[im]);
	}
	      fclose(outfp);

   fprintf(stderr,"KIRMIG COMPLETED for total number of traces =%d \n",ntrace);
	
        return 0;
}