the one is not in the list of SConstruct

user/pyang/Mlsprtm2d.c

@@ -1,5 +1,4 @@
-/* least-squares prestack RTM in 2-D
- */
+/* least-squares prestack RTM in 2-D */
 /*
   Copyright (C) 
   - 2014  Xi'an Jiaotong University, UT Austin (Pengliang Yang)
@@ -26,145 +25,145 @@
 
 int main(int argc, char* argv[])
 {   
-    bool verb, fromBoundary;
-    int nb, nz, nx, nt, ns, ng, niter, csd, sxbeg, szbeg, jsx, jsz, gxbeg, gzbeg, jgx, jgz;
-    float dz, dx, dt, fm, o1, o2, amp;
-    float **v0, *mod, *dat; 
-    int testadj;
+	bool verb, fromBoundary;
+	int nb, nz, nx, nt, ns, ng, niter, csd, sxbeg, szbeg, jsx, jsz, gxbeg, gzbeg, jgx, jgz;
+	float dz, dx, dt, fm, o1, o2, amp;
+	float **v0, *mod, *dat; 
+	int testadj;
 
-    sf_file shots, imglsm, imgrtm, velo;/* I/O files */
+	sf_file shots, imglsm, imgrtm, velo;/* I/O files */
 
-    /* initialize Madagascar */
-    sf_init(argc,argv);
+	/* initialize Madagascar */
+	sf_init(argc,argv);
 
-    shots = sf_input ("in"); /* shot records, data 	*/
-    velo = sf_input ("vel"); /* velocity */
-    imglsm = sf_output("out"); /* output LSRTM imglsme, model */
-    imgrtm = sf_output("imgrtm"); /* output RTM imglsme */
+	shots = sf_input ("in"); /* shot records, data 	*/
+	velo = sf_input ("vel"); /* velocity */
+	imglsm = sf_output("out"); /* output LSRTM imglsme, model */
+	imgrtm = sf_output("imgrtm"); /* output RTM imglsme */
 
-    if (!sf_histint(velo,"n1",&nz)) sf_error("n1");
-    /* 1st dimension size */
-    if (!sf_histint(velo,"n2",&nx)) sf_error("n2");
-    /* 2nd dimension size */
-    if (!sf_histfloat(velo,"d1",&dz)) sf_error("d1");
-    /* d1 */
-    if (!sf_histfloat(velo,"d2",&dx)) sf_error("d2");
-    /* d2 */
-    if (!sf_histfloat(velo,"o1",&o1)) sf_error("o1");
-    /* o1 */
-    if (!sf_histfloat(velo,"o2",&o2)) sf_error("o2");
-    /* o2 */
-    if (!sf_getbool("verb",&verb)) verb=true;
-    /* verbosity */
-    if (!sf_getint("niter",&niter)) niter=10;
-    /* totol number of least-squares iteration*/
-    if (!sf_getint("nb",&nb)) nb=20;
-    /* number (thickness) of ABC on each side */
-    if (!sf_getbool("fromBoundary",&fromBoundary)) fromBoundary=true;
-    /* if fromBoundary=true, reconstruct source wavefield from stored boundary */
-    if (!sf_getint("testadj",&testadj)) testadj=0;
-    /* if testadj = 1 then program only testadj without calculating */
+	if (!sf_histint(velo,"n1",&nz)) sf_error("n1");
+	/* 1st dimension size */
+	if (!sf_histint(velo,"n2",&nx)) sf_error("n2");
+	/* 2nd dimension size */
+	if (!sf_histfloat(velo,"d1",&dz)) sf_error("d1");
+	/* d1 */
+	if (!sf_histfloat(velo,"d2",&dx)) sf_error("d2");
+	/* d2 */
+	if (!sf_histfloat(velo,"o1",&o1)) sf_error("o1");
+	/* o1 */
+	if (!sf_histfloat(velo,"o2",&o2)) sf_error("o2");
+	/* o2 */
+	if (!sf_getbool("verb",&verb)) verb=true;
+	/* verbosity */
+	if (!sf_getint("niter",&niter)) niter=10;
+	/* totol number of least-squares iteration*/
+	if (!sf_getint("nb",&nb)) nb=20;
+	/* number (thickness) of ABC on each side */
+	if (!sf_getbool("fromBoundary",&fromBoundary)) fromBoundary=true;
+	/* if fromBoundary=true, reconstruct source wavefield from stored boundary */
+	if (!sf_getint("testadj",&testadj)) testadj=0;
+	/* if testadj = 1 then program only testadj without calculating */
 
-    if (!sf_histint(shots,"n1",&nt)) sf_error("no nt");
-    /* total modeling time steps */
-    if (!sf_histint(shots,"n2",&ng)) sf_error("no ng");
-    /* total receivers in each shot */
-    if (!sf_histint(shots,"n3",&ns)) sf_error("no ns");
-    /* number of shots */
-    if (!sf_histfloat(shots,"d1",&dt)) sf_error("no dt");
-    /* time sampling interval */
-    if (!sf_histfloat(shots,"amp",&amp)) sf_error("no amp");
-    /* maximum amplitude of ricker */
-    if (!sf_histfloat(shots,"fm",&fm)) sf_error("no fm");
-    /* dominant freq of ricker */
-    if (!sf_histint(shots,"sxbeg",&sxbeg)) sf_error("no sxbeg");
-    /* x-begining index of sources, starting from 0 */
-    if (!sf_histint(shots,"szbeg",&szbeg)) sf_error("no szbeg");
-    /* x-begining index of sources, starting from 0 */
-    if (!sf_histint(shots,"gxbeg",&gxbeg)) sf_error("no gxbeg");
-    /* x-begining index of receivers, starting from 0 */
-    if (!sf_histint(shots,"gzbeg",&gzbeg)) sf_error("no gzbeg");
-    /* x-begining index of receivers, starting from 0 */
-    if (!sf_histint(shots,"jsx",&jsx)) sf_error("no jsx");
-    /* source x-axis  jump interval  */
-    if (!sf_histint(shots,"jsz",&jsz)) sf_error("no jsz");
-    /* source z-axis jump interval  */
-    if (!sf_histint(shots,"jgx",&jgx)) sf_error("no jgx");
-    /* receiver x-axis jump interval  */
-    if (!sf_histint(shots,"jgz",&jgz)) sf_error("no jgz");
-    /* receiver z-axis jump interval  */
-    if (!sf_histint(shots,"csdgather",&csd)) sf_error("csdgather or not required");
-    /* default, common shot-gather; if n, record at every point*/
-
-    sf_putint(imglsm,"n1",nz);
-    sf_putint(imglsm,"n2",nx);
-    sf_putint(imglsm,"n3",1);
-    sf_putfloat(imglsm,"d1",dz);
-    sf_putfloat(imglsm,"d2",dx);
-    sf_putfloat(imglsm,"o1",o1);
-    sf_putfloat(imglsm,"o2",o2);
-    sf_putstring(imglsm,"label1","Depth");
-    sf_putstring(imglsm,"label2","Distance");
-
-    sf_putint(imgrtm,"n1",nz);
-    sf_putint(imgrtm,"n2",nx);
-    sf_putint(imgrtm,"n3",1);
-    sf_putfloat(imgrtm,"d1",dz);
-    sf_putfloat(imgrtm,"d2",dx);
-    sf_putfloat(imgrtm,"o1",o1);
-    sf_putfloat(imgrtm,"o2",o2);
-    sf_putstring(imgrtm,"label1","Depth");
-    sf_putstring(imgrtm,"label2","Distance");
-
-    /* In rtm, vv is the velocity model [modl], which is input parameter; 
-       mod is the imglsme/reflectivity [imglsm]; dat is seismogram [data]! */
-    v0=sf_floatalloc2(nz,nx);
-    mod=sf_floatalloc(nz*nx);
-    dat=sf_floatalloc(nt*ng*ns);
-
-    /* initialize velocity, model and data */
-    sf_floatread(v0[0], nz*nx, velo);
-    memset(mod, 0, nz*nx*sizeof(float));
-    sf_floatread(dat, nt*ng*ns, shots);
-    prtm2d_init(verb, csd, fromBoundary, dz, dx, dt, amp, fm, nz, nx, nb, nt, ns, ng, 
-		sxbeg, szbeg, jsx, jsz, gxbeg, gzbeg, jgx, jgz, v0, mod, dat);
-
-    // run adjoint test first.
-    sf_warning("adjoint test \n");
-    if (testadj){
-	prtm2d_adjtest();
-	sf_warning("exiting after testadj\n");
-	mod=NULL;
-	sf_floatwrite(mod,nz*nx,imgrtm);
-	sf_floatwrite(mod,nz*nx,imglsm);
+	if (!sf_histint(shots,"n1",&nt)) sf_error("no nt");
+	/* total modeling time steps */
+	if (!sf_histint(shots,"n2",&ng)) sf_error("no ng");
+	/* total receivers in each shot */
+	if (!sf_histint(shots,"n3",&ns)) sf_error("no ns");
+	/* number of shots */
+	if (!sf_histfloat(shots,"d1",&dt)) sf_error("no dt");
+	/* time sampling interval */
+	if (!sf_histfloat(shots,"amp",&amp)) sf_error("no amp");
+	/* maximum amplitude of ricker */
+	if (!sf_histfloat(shots,"fm",&fm)) sf_error("no fm");
+	/* dominant freq of ricker */
+	if (!sf_histint(shots,"sxbeg",&sxbeg)) sf_error("no sxbeg");
+	/* x-begining index of sources, starting from 0 */
+	if (!sf_histint(shots,"szbeg",&szbeg)) sf_error("no szbeg");
+	/* x-begining index of sources, starting from 0 */
+	if (!sf_histint(shots,"gxbeg",&gxbeg)) sf_error("no gxbeg");
+	/* x-begining index of receivers, starting from 0 */
+	if (!sf_histint(shots,"gzbeg",&gzbeg)) sf_error("no gzbeg");
+	/* x-begining index of receivers, starting from 0 */
+	if (!sf_histint(shots,"jsx",&jsx)) sf_error("no jsx");
+	/* source x-axis  jump interval  */
+	if (!sf_histint(shots,"jsz",&jsz)) sf_error("no jsz");
+	/* source z-axis jump interval  */
+	if (!sf_histint(shots,"jgx",&jgx)) sf_error("no jgx");
+	/* receiver x-axis jump interval  */
+	if (!sf_histint(shots,"jgz",&jgz)) sf_error("no jgz");
+	/* receiver z-axis jump interval  */
+	if (!sf_histint(shots,"csdgather",&csd)) sf_error("csdgather or not required");
+	/* default, common shot-gather; if n, record at every point*/
+
+	sf_putint(imglsm,"n1",nz);
+	sf_putint(imglsm,"n2",nx);
+	sf_putint(imglsm,"n3",1);
+	sf_putfloat(imglsm,"d1",dz);
+	sf_putfloat(imglsm,"d2",dx);
+	sf_putfloat(imglsm,"o1",o1);
+	sf_putfloat(imglsm,"o2",o2);
+	sf_putstring(imglsm,"label1","Depth");
+	sf_putstring(imglsm,"label2","Distance");
+
+	sf_putint(imgrtm,"n1",nz);
+	sf_putint(imgrtm,"n2",nx);
+	sf_putint(imgrtm,"n3",1);
+	sf_putfloat(imgrtm,"d1",dz);
+	sf_putfloat(imgrtm,"d2",dx);
+	sf_putfloat(imgrtm,"o1",o1);
+	sf_putfloat(imgrtm,"o2",o2);
+	sf_putstring(imgrtm,"label1","Depth");
+	sf_putstring(imgrtm,"label2","Distance");
+
+	/* In rtm, vv is the velocity model [modl], which is input parameter; 
+	   mod is the imglsme/reflectivity [imglsm]; dat is seismogram [data]! */
+	v0=sf_floatalloc2(nz,nx);
+	mod=sf_floatalloc(nz*nx);
+	dat=sf_floatalloc(nt*ng*ns);
+
+	/* initialize velocity, model and data */
+	sf_floatread(v0[0], nz*nx, velo);
+	memset(mod, 0, nz*nx*sizeof(float));
+	sf_floatread(dat, nt*ng*ns, shots);
+	prtm2d_init(verb, csd, fromBoundary, dz, dx, dt, amp, fm, nz, nx, nb, nt, ns, ng, 
+				sxbeg, szbeg, jsx, jsz, gxbeg, gzbeg, jgx, jgz, v0, mod, dat);
+
+	// run adjoint test first.
+	sf_warning("adjoint test \n");
+	if (testadj){
+		prtm2d_adjtest();
+		sf_warning("exiting after testadj\n");
+		mod=NULL;
+		sf_floatwrite(mod,nz*nx,imgrtm);
+		sf_floatwrite(mod,nz*nx,imglsm);
+		prtm2d_close();
+		free(*v0); free(v0);
+		free(mod);
+		free(dat); 	  
+		exit (0);
+	}
+
+	sf_warning("start migration\n");
+	/* original RTM is simply applying adjoint of prtm2d_lop once!*/
+	prtm2d_lop(true, false, nz*nx, nt*ng*ns, mod, dat); 
+	sf_floatwrite(mod, nz*nx, imgrtm);/* output RTM image */
+	sf_warning("migration ends");
+
+	if (niter>0){
+		memset(mod, 0, nz*nx*sizeof(float));
+
+		sf_warning("inside LS loop niter= %d\n",niter);
+		/* least squares migration */
+		sf_solver(prtm2d_lop, sf_cgstep, nz*nx, nt*ng*ns, mod, dat, niter, "verb", verb, "end");
+		sf_floatwrite(mod, nz*nx, imglsm);  /* output inverted image */
+		sf_cgstep_close();
+	}
+
 	prtm2d_close();
 	free(*v0); free(v0);
 	free(mod);
-	free(dat); 	  
-	exit (0);
-    }
-
-    sf_warning("start migration\n");
-    /* original RTM is simply applying adjoint of prtm2d_lop once!*/
-    prtm2d_lop(true, false, nz*nx, nt*ng*ns, mod, dat); 
-    sf_floatwrite(mod, nz*nx, imgrtm);/* output RTM image */
-    sf_warning("migration ends");
-
-    if (niter>0){
-	memset(mod, 0, nz*nx*sizeof(float));
-
-	sf_warning("inside LS loop niter= %d\n",niter);
-	/* least squares migration */
-	sf_solver(prtm2d_lop, sf_cgstep, nz*nx, nt*ng*ns, mod, dat, niter, "verb", verb, "end");
-	sf_floatwrite(mod, nz*nx, imglsm);  /* output inverted image */
-	sf_cgstep_close();
-    }
-
-    prtm2d_close();
-    free(*v0); free(v0);
-    free(mod);
-    free(dat); 
+	free(dat); 
 
-    exit(0);
+	exit(0);
 }
 

user/pyang/SConstruct

@@ -6,7 +6,7 @@ import bldutil
 ######################################################################
 # 1. C main programs
 progs = '''
-cohn checkptdemo fwi2d fbrec2d istseislet lsrtm2d  
+cohn checkptdemo fwi2d fbrec2d istseislet lsrtm2d lsprtm2d 
 mythresh mshots mcaseislet modeling2d mysnr
 pocsseislet rtm2d rtmadcig rtmodcig rtmva2d
 Testfd2d Testfd3d Testeb Testelastic2d Testaniso Testspml