csmodel.c

seismic software,very useful

	fprintf(paramfp,
"geometry-file_os=%g                       :receiver geometry\n",os);  
	if(plot==-1) {
	fprintf(paramfp,
"                                    :second plot descriptor (sgq)\n");
	fprintf(paramfp,
"l                                   :job descriptor (rlt)\n");
	} else {
	fprintf(paramfp,
"sg                                  :second plot descriptor (sgq)\n");
	fprintf(paramfp,
"rl                                  :job descriptor (rlt)\n");
	}
	fprintf(paramfp,
"csm_os=%g                                 :output filename(s)\n",os);
	fprintf(paramfp,
"%-5.1f  %-5.1f                        :range of takeoff angles\n",amin,amax);
	fprintf(paramfp,
"%-5.1f                               :increment in takeoff angle\n",da);
	for(ih=0;ih<nhs;ih=ih+3) {
		imax = 3;
		if(ih+imax >nhs) imax = nhs - ih;
		for(i=0;i<imax;i++) {
			jh = ih + i;
			tmp = 0.;
			j1 = 1;
			j2 = npicks[jh];
			nn = npicks[jh];
			if(ns==1) {
				tmp1 = sx[0]+r1;
				if(r1>0.) tmp1 = sx[0]; 
				tmp2 = sx[0]+r4;
				if(r4<0.) tmp2 = sx[0];
				if(tmp1>tmp2) {
					tmp = tmp2;
					tmp1 = tmp2;
					tmp2 = tmp;
				}
				if(tmp1<xpicks[jh*maxpks]) {
					j1 = 1;
				} else if(tmp1>xpicks[jh*maxpks+nn-1]) {
					j1 = nn;
				} else {
					bisear_(&nn,&one,
						xpicks+jh*maxpks,&tmp1,&j1);
				}
				if(tmp2<xpicks[jh*maxpks]) {
					j2 = 1;
				} else if(tmp2>xpicks[jh*maxpks+nn-1]) {
					j2 = nn;
				} else {
					bisear_(&nn,&one,
						xpicks+jh*maxpks,&tmp2,&j2);
				}
			}
		 	for(j=j1-1;j<j2;j++) {
				tmp = tmp + veltops[jh*maxpks+j] 
				    + velbots[jh*maxpks+j];
			}	
			tmp = tmp /(j2-j1+1);
			fprintf(paramfp,"%-7.1f ",0.5*tmp); 
		}
		if(ih+imax!=nhs) {
			fprintf(paramfp,"\n");
		} else {
			fprintf(paramfp,"                       :velocities\n");
		}
	}
	fprintf(paramfp,
"%1s                                   :direct wave? (y or n)\n",direct);
	if(nhead==0) {
		fprintf(paramfp,
"                                    :headwave interface numbers (1,2,...)\n");
	} else {
		for(ih=0;ih<nhead;ih++) {
			fprintf(paramfp,"%d ",headhn[ih]);
		}
		fprintf(paramfp,
"                                    :headwave interface numbers (1,2,...)\n");
	}
	fprintf(paramfp,
"n                                   :all primaries? (y or n)\n");
	for(ih=0;ih<nh;ih++) {
		fprintf(paramfp,
"%-2d                                  :reflection from interface %d\n",
			hn[ih],hn[ih]);
	}
	for(i=0;i<nm;i++) fprintf(paramfp, "%s",mn[i]);
	
	fclose(paramfp);

	/* output param1_2 file */
	paramfp = fopen(param1_2,"w");
	fprintf(paramfp,
"model-file_os=%g                          :model file\n",os);
	fprintf(paramfp,
"%-2d                                  :#interfaces in model\n",nhs-1);
	fprintf(paramfp,
"plotcolors_os=%g                          :model colors file\n",os);
	fprintf(paramfp,
"                                    :first plot descriptor (mwq)\n");
	fprintf(paramfp,
"don't care                          :well coordinates\n");
	fprintf(paramfp,
"s                                   :shooting mode (sd)\n");
	fprintf(paramfp,
"geometry-file_os=%g                       :receiver geometry\n",os);  
	fprintf(paramfp,
"                                    :second plot descriptor (sgq)\n");
	fprintf(paramfp,
"t                                   :job descriptor (rlt)\n");
	fprintf(paramfp,
"csm_os=%g                                 :output filename(s)\n",os);
	fprintf(paramfp,
"%-5.1f  %-5.1f                        :range of takeoff angles\n",amin,amax);
	fprintf(paramfp,
"%-5.1f                               :increment in takeoff angle\n",da);
	for(ih=0;ih<nhs;ih=ih+3) {
		imax = 3;
		if(ih+imax >nhs) imax = nhs - ih;
		for(i=0;i<imax;i++) {
			jh = ih + i;
                        tmp = 0.;
                        j1 = 1;
                        j2 = npicks[jh];
                        nn = npicks[jh];
                        if(ns==1) {
                                tmp1 = sx[0]+r1;
                                if(r1>0.) tmp1 = sx[0];
                                tmp2 = sx[0]+r4;
                                if(r4<0.) tmp2 = sx[0];
                                if(tmp1>tmp2) {
                                        tmp = tmp2;
                                        tmp1 = tmp2;
                                        tmp2 = tmp;
                                }
                                if(tmp1<xpicks[jh*maxpks]) {
                                        j1 = 1;
                                } else if(tmp1>xpicks[jh*maxpks+nn-1]) {
                                        j1 = nn;
                                } else {
                                        bisear_(&nn,&one,
                                                xpicks+jh*maxpks,&tmp1,&j1);
                                }
                                if(tmp2<xpicks[jh*maxpks]) {
                                        j2 = 1;
                                } else if(tmp2>xpicks[jh*maxpks+nn-1]) {
                                        j2 = nn;
                                } else {
                                        bisear_(&nn,&one,
                                                xpicks+jh*maxpks,&tmp2,&j2);
                                }
                        }
                        for(j=j1-1;j<j2;j++) {
                                tmp = tmp + veltops[jh*maxpks+j]
                                    + velbots[jh*maxpks+j];
                        }
                        tmp = tmp /(j2-j1+1);
                        fprintf(paramfp,"%-7.1f ",0.5*tmp);
			/*
			jh = ih + i;
			tmp = 0.;
		 	for(j=0;j<npicks[jh];j++) {
				tmp = tmp + veltops[jh*maxpks+j] 
				    + velbots[jh*maxpks+j];
			}	
			tmp = tmp /(npicks[jh]*2);
			fprintf(paramfp,"%-7.1f ",tmp); 
			*/
		}
		if(ih+imax!=nhs) {
			fprintf(paramfp,"\n");
		} else {
			fprintf(paramfp,"                       :velocities\n");
		}
	}
	fprintf(paramfp,
"%1s                                   :direct wave? (y or n)\n",direct);
	if(nhead==0) {
		fprintf(paramfp,
"                                    :headwave interface numbers (1,2,...)\n");
	} else {
		for(ih=0;ih<nhead;ih++) {
			fprintf(paramfp,"%d ",headhn[ih]);
		}
		fprintf(paramfp,
"                                    :headwave interface numbers (1,2,...)\n");
	}
	fprintf(paramfp,
"n                                   :primaries? (y or n)\n");
	for(ih=0;ih<nh;ih++) {
		fprintf(paramfp,
"%-2d                                  :reflection from interface %d\n",
			hn[ih],hn[ih]);
	}

	for(i=0;i<nm;i++) fprintf(paramfp, "%s",mn[i]);
		
	fclose(paramfp);

	/* output param2 file */
	tmp = (r2-r1)/dr + (r4-r3)/dr + 1 + 1 + 0.5;
	nr  = tmp;
	
	paramfp = fopen(param2,"w");
	fprintf(paramfp,
"s                                :job option (s,r)\n");
	fprintf(paramfp,
"1  %d                             :first, last shot for sort\n",ns);
	fprintf(paramfp,
"1  %d                           :first, last trace OR first last receiver\n",
		nr);
	fprintf(paramfp,
"%-4.1f %-4.1f %-4.1f %-4.1f              :frequency spectrum of wavelet\n",
				f1,f2,f3,f4);
	fprintf(paramfp,
"%-4.3f                            :wavelet length (secs)\n",wl*0.001);
	fprintf(paramfp,
"%-4.3f                            :sample rate (secs)\n",dt*0.001);
	fprintf(paramfp,
"%-6.3f                            :record length (secs)\n",tmax*0.001);
	fprintf(paramfp,
"csm_os=%gshot                          :input filename\n",os);
	fprintf(paramfp,
"csm_os=%gtraces                        :output filename\n",os);
	fclose(paramfp);

	/* output plotcolors file */
	sprintf(fname,"plotcolors_os=%g\0",os);
	paramfp = fopen(fname,"w");

	fprintf(paramfp, "0          receivers\n");
	fprintf(paramfp, "4          sources\n");
	fprintf(paramfp, "6          well color\n");
	fprintf(paramfp, "%d          caustic rays\n",ccolor);
	fprintf(paramfp, "%d          rays\n",rcolor);
	fprintf(paramfp, "%d          interfaces\n",icolor);
	fprintf(paramfp,"\n");
	fprintf(paramfp,"\n");
	fprintf(paramfp,"\n");
	fprintf(paramfp,"key:  (CWP's xgraph colors)\n");
	fprintf(paramfp,"0      black\n");
	fprintf(paramfp,"1      white\n");
	fprintf(paramfp,"2      red\n");
	fprintf(paramfp,"3      green\n");
	fprintf(paramfp,"4      dark blue\n");
	fprintf(paramfp,"5      light blue\n");
	fprintf(paramfp,"6      violet\n");
	fprintf(paramfp,"7      yellow\n");
	fclose(paramfp);	


	/* plot rays */
	bzero(cmd,1024);
	sprintf(cmd,"cp %s %s",param1_1,param1);
	system(cmd);
	if(comp==0 || comp==1) {
		bzero(cmd,1024);
		sprintf(cmd,
	"cshot1 param1=%s | cshotplot >csmplot_os=%g outpar=csmpar_os=%g",
			param1,os,os);
		system(cmd);
	}else if(comp==-1) {
		bzero(cmd,1024);
		sprintf(cmd,"cshot1 param1=%s >/dev/null",param1);
		system(cmd);
	}

	if(plot==0) {
		bzero(cmd,1024);
		if(rfsave==0) {
			sprintf(cmd,
"( xgraph <csmplot_os=%g par=csmpar_os=%g style=seismic title='Csmodel raypaths' label1=Depth label2=Distance x1beg=%g x1end=%g x2beg=%g x2end=%g ; /bin/rm -f csmpar_os=%g csmplot_os=%g ) &",
			os,os,zmini,zmaxi,xmini,xmaxi,os,os); 
		} else {
			sprintf(cmd,
"xgraph <csmplot_os=%g par=csmpar_os=%g style=seismic title='Csmodel raypaths' label1=Depth label2=Distance x1beg=%g x1end=%g x2beg=%g x2end=%g &",
			os,os,zmini,zmaxi,xmini,xmaxi); 
		}
	} else if(plot==1) {
		sprintf(cmd,"psgraph <csmplot_os=%g par=csmpar_os=%g title=Csmodel label1=Depth label2=Distance x1beg=%g x1end=%g x2beg=%g x2end=%g | lpr &",
			os,os,zmini,zmaxi,xmini,xmaxi); 
	}

	if(comp!=2 && comp!=-1 && plot!=-1) {
		system(cmd);
		if(plot==1 && rfsave==0) {
			bzero(cmd,1024);
			sprintf(cmd,"/bin/rm -f csmpar_os=%g csmplot_os=%g");
			system(cmd);
		}
	}

	if(comp==0 || comp==-1) goto notrace;

	bzero(cmd,1024);
	sprintf(cmd,"cp %s %s",param1_2,param1);
	system(cmd);
	bzero(cmd,1024);
	sprintf(cmd,"cshot1 param1=%s >/dev/null",param1);
	system(cmd);
	bzero(cmd,1024);
	sprintf(cmd,"cshot2 param2=%s",param2);
	system(cmd);
	nt = 1 + ( tmax + dt / 2. ) / dt;
	j = dt*1000;
	bzero(cmd,1024);
	sprintf(cmd,"suaddhead <csm_os=%gtraces ftn=1 ns=%d | sushw key=dt a=%d >csmtraces.segy_os=%g", os,nt,j,os);
	system(cmd);

	bzero(fname,80);
	sprintf(fname,"csm_os=%gtraces\0",os);
	if(rfsave==0) unlink(fname); 

	bzero(fname,80);
	sprintf(fname,"csmtraces.segy_os=%g\0",os);
	datafp = fopen(fname,"r");
	i = 0;
	dcdp = dr * 0.5;
	if (ds!=0. && ds<dr) dcdp = ds * 0.5; 

	fgethdr(datafp,&ch,&bh);
	nt = nt - 1;
	bh.hns = nt;
	bh.ntrpr = nr;
	bh.tsort = 1;
	bh.hdt = dt * 1000;
	bh.format = 1;
	fputhdr(outfp,&ch,&bh);

	while( fgettr(datafp,&tr) ) {
		i = i + 1;
		js = (i-1)/nr;
		tr.tracl = i;
		tr.sx = (int) sx[js];
		tmp = r1 + (i-js*nr-1)*dr + 0.0001; 
		if(tmp<=r2) {
			tr.gx = (int)tmp + tr.sx;
			j = i;
		} else {
			tmp = r3 + (i-j-1)*dr;
			tr.gx = (int) tmp + tr.sx;
		}
		tr.ep = js + 1;
		tr.fldr = js + 1;
		tr.tracf = i-js*nr;
		tr.offset = tr.gx - tr.sx;
		tmp = tr.gx + tr.sx;
		tmp = tmp * 0.5;
		tmp = tmp/dcdp + 1.5; 
		js  = tmp;
		tr.cdp = js;
		tr.scalco = 1;
		tr.trid = 1;
		tr.delrt = (int) dt;
		tr.ns = nt;
		for(it=0;it<nt;it++) tr.data[it] = tr.data[it+1]; 
		fputtr(outfp,&tr);
	}
	
	if(plot==0) {
		bzero(cmd,1024);
		if(rfsave==0) {
			sprintf(cmd,
	"( <csmtraces.segy_os=%g suxwigb title='Modeled gather' label1=Time label2=Trace xcur=2 grid1=solid perc=99 ; /bin/rm -f csmtraces.segy_os=%g ) & ",os,os,os); 
		} else {
			sprintf(cmd,
	"<csmtraces.segy_os=%g suxwigb title='Modeled gather' label1=Time label2=Trace xcur=2 grid1=solid perc=99 & ",os); 
		}
		system(cmd);
	} else {
		bzero(cmd,1024);
		sprintf(cmd,"<csmtraces.segy_os=%g supswigb title='Modeled gather' label1=Time label2=Trace xcur=2 grid1=solid perc=99 | lpr ",os);
		system(cmd);
		if(rfsave==0) {
			bzero(cmd,1024);
			sprintf(cmd,"/bin/rm -f csmtraces.segy_os=%g ",os);
			system(cmd);
		}
	}


	notrace:

	if(rfsave==0) {
		bzero(fname,80);
		sprintf(fname,"model-file_os=%g\0",os);
		unlink(fname);
		bzero(fname,80);
		sprintf(fname,"geometry-file_os=%g\0",os); 
		unlink(fname);
		unlink(param1_1);
		unlink(param1_2);
		unlink(param1);
		unlink(param2);
		bzero(fname,80);
		sprintf(fname,"plotcolors_os=%g\0",os); 
		unlink(fname);
		if(comp!=0) {
			bzero(fname,80);
			sprintf(fname,"csm_os=%gshot\0",os); 
			unlink(fname);
		}
		bzero(fname,80);
		sprintf(fname,"csm_os=%gdata\0",os); 
		unlink(fname);
		bzero(fname,80);
		sprintf(fname,"csm_os=%glisting\0",os); 
		unlink(fname);

		system("sleep 5");
		if(comp==1 || comp==2) {
			bzero(fname,80);
			sprintf(fname,"csmtraces.segy_os=%g\0",os); 
			unlink(fname);
		}
	}

	return(0);
}