From 6bc8fb75c68243a046e629f03084b2f3214597d3 Mon Sep 17 00:00:00 2001
From: Rui Gong <rg49534@BEGW-SF-A09289.beg.utexas.edu>
Date: Tue, 21 May 2024 20:15:26 -0500
Subject: [PATCH] reproducible version

---
 book/tccs/shapestacks/Fig/schematic.pdf  | Bin 0 -> 42655 bytes
 book/tccs/shapestacks/SConstruct         |   8 +
 book/tccs/shapestacks/elf/SConstruct     | 422 ++++++++++++++++++
 book/tccs/shapestacks/paper.bib          | 309 +++++++++++++
 book/tccs/shapestacks/paper.tex          | 529 +++++++++++++++++++++++
 book/tccs/shapestacks/synseis/SConstruct | 476 ++++++++++++++++++++
 6 files changed, 1744 insertions(+)
 create mode 100644 book/tccs/shapestacks/Fig/schematic.pdf
 create mode 100644 book/tccs/shapestacks/SConstruct
 create mode 100644 book/tccs/shapestacks/elf/SConstruct
 create mode 100644 book/tccs/shapestacks/paper.bib
 create mode 100644 book/tccs/shapestacks/paper.tex
 create mode 100644 book/tccs/shapestacks/synseis/SConstruct

diff --git a/book/tccs/shapestacks/Fig/schematic.pdf b/book/tccs/shapestacks/Fig/schematic.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..a8040967c4d9761cc71ad1e0f5994e64beecf7b8
GIT binary patch
literal 42655
zcmeEvby$|!x;G`AO6LpG{nFjtE#2MHozmT*lpsihba!`%NFyy0(g?^0=s0_rIp=Km
ze6wet>pK6)<NMV8yVu?8Uh5%|6B3~T(lWx5bgmt&ZRZ_kj&^p#G6CoTHu~nUTwDM;
zDI;qW2U7s^%})vdIuSEV2P3<izm|FqMnXmgHikw39v)bG2RkD@D_9rMp_GA0+#VF+
z(_^$c;s<5Pdp#<Zdec6j1D1Sn18{Ei3>5IvNa7kgAFn^xsYO|<ylte0rH4{*7B&(Q
zF`f&qjq$W-W~}Kxz>s8aW;^^m<6-kvPy6d`9HaNOp?33!!;f)Kul6Pqv|XCYJdZAn
z7HlV4)a)m!=96pIdu)u?6dDi>eG=XDL<_yyyFY5}POva2H8azatTf^Xnxl|h0jGhp
ziODMh%Pn*S1~~JVD)t+q@yB)f63Q7GVD_`Bc|A4j7dF$UWAc|<Usm28^YUJA4@I1x
z4bLH#pXl{5X4pSn&}?F*^S)k4S@hV{a4JEvO<V&KK6CZSfV7PCP>uB9F|UrZEM!T|
zwy9l8XB=Z8@FeBn!yEgYa<tMD#3O7mG?JUjD`V<%jNNf*{^9HF_LFv{Q=5De^|?ZW
zwE2c$ZCfv3yl=(f`*+<$)JpbZXRa~NW?l>kBy1Sz4*C-3ePZ8!u9sFNHteSG`oYj<
zem?OQ1<5L9pT;3}skDq*7RTmncGgT4#9mK#i)171n`N+9ZQZVOSal}vv6=_w&^n6V
z^ajG^1(;%Zw!jZ*SC^V1Ql+kOWe%pW2MqeHpE-^;T5eshOyeV+E#pGSf;HOlQ@<CL
zafiILO?<J8r?`}*1>%)Zw6A^iOqii47Hz*skVj|)a-WXXY)I%KqW1!)A~fp*QgNCc
zy`{&^=JqlJ6YeHn4+X?fyqlX{56|aYN_9MI&-XviUu=&MENUUFN^?RQRgIn{S{9%&
zJ5i6-#m2dgqx4r!S(1ixuVs@;_M3KfuLqxDg5|#=!7&HSodf!4IyB&80>ob$#(nH>
zy>8>qpHtM+By4|Cv-y;VZ7m16@tDFQutuIzUM~AI*BIlt*BCeBE!jaopX28kwmODa
z)x-gHOa-_n0WVysv0Ci>a)dt;>_+gj%c8%+4lWt>c+HG7ysmi>j9Q*}8X!#}i5@8D
z5WP<kV&m(tnGrX@VC;&YxN#URY&T)Q=qbk+(NiB8z(D=LS{mk4@FYu1=_7PR!cobM
zxB|bh_vS4;=-c&STH#Aj1WD)>lop+=cE0vQ!pgILHLHvo<Es^`j377ymicIngM+z@
zGWz&~<N=S3A)9+qq#pO-HqK2Ey#h-^C#TP58fJ>^4Z{sugq{9a3q8g>WJ;l`@+$Yk
zSmqKWKNn-8J4_Bc(;Jb<x>W|DOH`}36=X9AX9Ld|M*8{+FcsiSvYGv5rT`CIGTPd?
zh}vv?uFG7F5#BQO;?Sc!nIITkMn}&p#U*?nrrKu^wPGXNEuE9=-DvB`w)}YlMk-|@
zYW@9I9s8A`mwW2wkQRsA(CK<7;uTZiqbE-#xLPu~awF6&B9GcdcdLd5s}WA8t0z7i
z&dz07(OlcjI&q|1K5r?r!e}(>Um;8lzv$&t+^=@`tJv7sEhk&kGdp47{~R1dO06M*
zl}Y-75(V7#@x&0jv0jZQR+*~4NpLbW43cnBsUNa3!vj2GHe$9-nZRajkqz8c=@p@E
zWJ4Rc@<2ee2nofM;wz4>Jb%oRlQ-hZ`B`s>Ra{<5W%rgb6h;QvcKP%>#czB}hNA$P
z9Z0q~DZ=QVrq%)x0d;8CSy6r7v`Au~9Yta0OAi@)EjPc}d$j7i)boK{dy^;*tpMWT
zu*jMr-K=<E+N!N~&0EDN2M(`F-!G~&VpV2w2+^%90GRa3nI7VP>UuXtXIYg#F9&*E
z;XP<)Uq_#48`%!hI%-fZlBAYo6J9R&L(MDSjcxm~t0(qfI^wOWH}S<T`gQdw47FA;
zSYXem(a7kU>lx7M-|&1W<GO-{Cy_($Ow@CWw1a2{FR3g;vR&sM^?SV^DXp{tr!~xx
zO+w(Y*4hT82U6%0JDFP$Kh`hy(aYch75g)u@dwKGG%Q_ZjuZ-LAcKjcuCIcg1)g7C
z*f@L+XbXKul~|Ur0+m-{Le-No1uoB&vq)t}Khb=QPby<&Yh2H|VJBCR;c(sthublA
zN@S@#c=kfi|Aapvo>&z+wdi>|p;%?F_;FoJMh2ZRb!<=y>=qQ;U}b_|MUa30IToV5
zzddwN{+PEV5><Q$6kuA%1?LN&OuTcJ!CBb{FIO;8R9*OnAUUsP&S{E;#L`|aG;Prj
z{RenT{#g6)J}80F+-E3sKETNZG(icX5Jn}25@RK(IzGjGG==;sB4ULNb--(G9Sle}
z4^`88#^wX+?K1~ql@L4q<|!X6M;1SG`VK}V2}VsSsi_OfLJ&=J8QjpTHKMOzj|x&e
zr5+D{7}z=NRI#S(702Pn4T9l+iZFS?^k$7*CpV^nt-rh=NOD)?0;v2Mh^x*q?<l3T
zDXCv2@BU>31-xEzC+ls!bnxnyOc~3^C~^2cqw;!!Qt!$0ult`yC!ri)_3@ZwCg2=?
zJ!vs>!1{<xRTIn^(d0p`<fy>~vZ~IdbU5-}zj*+}@6nU>$8t@yZ@T87ykA+i;i2|{
z!A;GyGH$6v?^WrE1x#*^tK&>zSy+h)5gd0=R}sdkn?X&@M5{k0c(iO`1+ZS0HcYQ$
z;j5aArLJltRsf%_eT@*p0*jM0<D>ks)g|2mT&-WxxR*}SO6;gD5#Kl;;iH+EIwtUK
zppch-#Ih!if^>qwo2kGx*FN;w4XTQm)`xF*UK~P^D4eCILa#-I1P-!Bp9-&yE8vQ(
zX_*s5_KQ|=6Tj>FNI_{*$hvBlz0F2RFy)AqfRpe<RQ|k7q=xjR)`g3=QhNldWB+Mq
z+HR@%X}#>b`EfoRcZXLwu-siA<J8fOI4+^19tLwK@u_uE`#ww<LEW0;9`v-%xUeCZ
zV5p>q=)tTd|M+A*L?wvxC5l7cqzG*}>|rFaS9o+f5A41mv{(wl@FWS!6t^6Y=q?*m
zZe(X>3EByRnTe*NMjoeog2@)8%)B_|5ruF)U^zgF2qtNVa@mv^MzMk<n`V5EX1^8X
zA;tEVtvXvgiaK^4Q??!utppOqmgKR~r-j3YaANfFcO>QSn^#o=^=9SWL#mtb0$|dK
z7BXL|$6U{W0P&Xb;}9CVA&J$JNU-C~@ZX<ImRRnw*D<RXwvrf7yLZ2S#QCWE+zqY-
zvF7|7Gx6hN{Ogq<A8jrvO^lLSkWv^WX<TnrVmtn%){I>DrsqcEYV=pZ<y0;kP2<b$
zogQ?@p~dGMOK(27mcYN5iNyExvVOfA-mGEWx*Xti&M&HGcZn0lG)^x^%x1#@*rdK#
zrIq~T7z0uz{KnABiD!gE62E*5B=KFuIN2^XXzlKr$cD-?l~w%3$-|JOLy?Ll@_L|E
zS_EzEVRfRKNgp0g6G*9Ci{^y?i+Pqh{AM4zG%<qm%XH`3d7UhTT*;n2HlC>k7hC>L
z1Zs!w(KqV#sHWU~Ppo92Wiav_*z<Kf1bHT1qET!8mBmJfRd**~p!4uU;PH_`2(5iC
zNR&wYWf^saugq0x3TTbuX_8Tj6`o*ix`I0wF#&zDMwVbxjLdAr-ZY)sI+S8Bmup(s
z5frOX1&RpH6~81LPk17wxpeyRfys)rT{I~VH~j8_tGRvi1)ZaQi5aYswc+<L^ycU`
zdc6%<Z~ro|urb}9y*vJg$X3bK)(Ai+r)P5aubq*#1Ay^cl&fH5Z{uiZU}O(q`SydL
zjkUwgd3(TZWPB3_OB)%Q=?U1l05s@teqdr{12C`xwP0@|;qULcz3-O@TF%bKK+(tn
zpmFmuArSzbl97u8fKJ@%=EVZv{|bKpD-O^C&<WaD+Sn=D>KPaTZUb;Zdmw=I_O-A#
zhYSFQZzS;W0O*8W97GizZUXVU>qKu(GP8Z7PxR&<X7;<|+Y3c*$iKTrNeKvGxC`ZP
zXa&%z%IcdN8GJh>ZUtllFn-5-^B*99@fO<+SaB;Rz|Fbu$3VdMcZgdt0)Py+Kf~T!
z!f?ZeA14?9z#sSCu;%Wf+vl(XZZ7-fIrIRAA2$Hm08DppxP26m9l-P}9gG0RAMd_-
z67!FTF)#v{e?7wlVE%Cr3*hE?-{}If0GNJ#Co6#I$NfNN0OJn=7;gCW<1%Kz&0W8c
z$pB#d0g@fS@a@D6^L`!cJEr|I*x%D7#=8-_8S6V4`-RVRijMjYcfZ_?K7j6tnc+=t
z#l#Gt6EiY1F?G1%74U{pH^OFTYjc}o!8+JE-iYTvh{X2?D@yCxTL4(T(=1@-U@vE6
zCun13Yh!&Y$ZUU@s@)3ck5p}E6chmj1f)IlRMrZSHJDGY)dv(D5+0LY_wPAqZEFQ_
z2m=BAy&~T+^OnIt<{yfDC!F`04*WK1{}j`2B=i@i-yZ%1&%aC1zru3^UCalafgm8D
z)PIQcmi{}OKVaVJ-94B$f_o>Ne=*GOy7m7b<{kaNg?X!o_fUSbUjL*5-H7{l1-h{h
zKQor!_55D}`49TD^9dA%{D07#C2(*&zCRlETcEetnOS~u9(Pvt9y@^FUDH1`=r@-5
zU*yMkgZOvy;|}R>8FJ@Y?lFYvX9VP5R(94uVPOAzWxu2U78CGSWxsvdJ(xFs_$QS8
z#(8|#gP*f$tbfv`O~QbJfJ*$4bl=hcTbOt8z&)5Z&ip4}{;Q^&^)It&cl6(b`3IwU
z<K6#xm^Ts8FQfUBVFBx(GHXE4yRhI`cE<+t?+y$0!NH#%{gGVVaq^Dg%)npl=3P8_
zU$B10O#JJC*i8m-H?V)waQ=q}Vhq3ds(%Q??&9`)OuymdKQHG0vTbJj>umEcspGG<
z`7R^55A$b&eBe#y@vA!jRFuN@C#>_%=nr+aLO}Yv)%kBqE&krb=AV6gGbH3w_dl|A
zOt*4&D|<|gKVrnYr1+jZGX8ti?k{tqcl6(4`c>`j(*1idng1=Z0^6VPq5pm4cSrv{
zn72H=5A$b6^Ebxjr%DKHf68<D-h>Ju9hv8f1gimZe^Uja^<687W%nN`4)d>e{6~;|
zS0=j8>Yp;?|7!Hf{+IbH=3flNuZnY50lNqDrm*&JaLep}!dw0KigQQ*Zx!dRws#LE
z)6WFiHzM+jS7v1VnbHCKpD=?vqd)W}V+r!VR2*krG~quI>pOtApqW^I34QO%n)mte
zQ+D!Sw;=3)!czVNC;mnY^3Qtc|DnY5_tZ`Q9mM-4`GVSv{h=)03fY}3GBW-mjCW=A
zd&0>0@3B39!kPW|w&#xiTTH*2*t-dXdoY=Q#`gS)lmJNoC!ERu<ph7%9LG=V>VKIN
z+zHw@z59_9+|63u6G-6Cgf{<r#R*9NC%pbY$m4%<$?1-jcP8N1lGEMX;60Y##1a2|
znDu8^{-+|4A7z36WR~Bt@*d0YreN=}oZ)Ar{ASkpN62?`_|MP$0qOr#r1OL0R<OMP
zWg({27ZmN;A354Ph_|qr*naV|EVth79w%6S#?RiHC;D}k8J14S$jQvWNYqZx^*0;>
ze4CH{C(E3_X$^jxsHWIXkL-ivmWgYpB8n^&4W>;84_HG__NB~}Ge5uf?!+acK_KN~
z-xb79-kbfD^B%!(p8yQRRw{+>V4(Ri1Td#t_Ur+GOft2f4rJ}kR6XYgx=*Y%)|LC)
z5w(YMB^(4{UL<`K052mNw??paEAsd_uk<}Wm)TV30hj<HkPObwI7??En`pUE_XmJ9
zNWo#~o>1_K1wdA|>jPjsC2Crh@<&AraDPe=HP4R-qQjUtlF9gR6;8l)S{N@pPrVW?
zJrvyCvJ>c4CR>hvuz)G6$-0g`uDNAtI-#}sat0S7At$gfQ2|DRm}*XXSMo-dLlX_(
zp~uVT##FZZa`f2m0cKG6UIwMzNH7X0WV~*{hlkbd#YSLr$9#=mcwcz{PTbXIOkTFW
zec2s-Ruynk3yx>DVhjg+EAU=>2TPX}p3MwedY)c2VD&?;J!?XZIXq5~bi<f3NkaRY
zUQDZ=ZDE=VXzGq{L~vTy$8FM@{A@vkzxq)5jR_EpJemx1pO*>0a!qJE?4R9EKA~QH
zW@RxoW%?dz|7dyLcNvNN+D<wxfC=|-??}zvM|L|7yy69E>uJf%ss~jx(yCB->xBy;
zKRofv7dCk9%Lxj^!!@Xf?J?J^!AE4<7!Z$8+iyklw_P}Q0{KmDe`G6n^Aq<3m-S~<
z_`cQrYbAFlwJ))N12qsP%@6j3gk{<w%qWQb*`()nUOmn>;xP2rxhnHwY+O<1GpQZo
zYwPTnmjzJ;+mxo!=2aTX@W#eAIF*ipMS}%i9t6pl7Ifp&9jaB;AEJDHlQUmO_$>0$
ztJqgx&Op6ikIZP`F!Z0|*jI!yICkodvWnjEI`CBVw;XlHXtoDVZR(tNw;k!i?|NuK
zsqM8TF$Wv1@AqM_kmA|L@25O?6BmD>9A13Jm8L{r&9T;p8An&Ke#zdh`c`}21^SI6
zR*h<hbruJXm*9N+{DGuwdw^T=WYv1Q_t)SCJlEiQjonn5!GNRT5CS=s?uZ7HxwQ~s
zE$zL&GK^+nL70HFV*j+1&v4URQZL=(3`mAqTUeQj224do;vTBMQ?Er5KQEW0K#xE|
z%^saBo_fb6FaYHa)1zZNbOsl`UBn=+$b^)Vjget+E)68_xoNVVYcVdHp4AJ4lo*;f
z!r1>X$z|qBs?<SmkxIIsg<x<7=#68nj3#<uQZ{Zb`SnHS+ivo8esH*ic<s(Q8@^W|
z>N>f7xHJurPd>lzF~=`1ml)54IR=|#dC_Sw<)ix+f3a(*^fCR!3s<Q!e_Tk$)WW>l
z5e>9VXd`2bbL6L^HK3fy%fmKmZmCRE>RZM1O+Z9N6e~ysMVjms<5NHKl_aAU<M{D{
zEtyEV7?0tGYWZsMT#fm2c9c@v8H&Xm^T~@s3(9oS7LPb7W*q;e?uwo%sus~)Wn7y$
z8@rE*i;u^4&-NdVwF_H|RB#N42{mgu;kv@VJbd1zLpW6w^WsFO{E78v_q7=bVpS(y
zj$%m!>l7wawb=0^!py)Jtl=Zlm*VH_^u;sUbHa@d{WttJwmC(OGOFlQ&&94<IOwk!
zr!p?Eb~erv;yNW6Cmc7w=*bfsIh3jD-D~a>q&~B{lBctSe>1Ur*2A`P8nJ;|XVRQK
zd~`-^$LSC@$3xaou}@(t{fw)rSB40#M##y@#D-xX_L|eI^P`lY`XRQ9OiUWABGnK~
z8d;3Lek^)10t-^iyrg)*vyF^87t0_Li*-Av^|2?|GN}RXBUj?`R=i(A-*Fd8<#eu}
zGiWnO10c%!&nT*kAF{X5$SD<^ofXYxEv&06rqgW8d70#5T-l#cw`vNH%^{2klydpR
zp$&c_a5-Ur8usj}iIrpNV-F9?SiKH7WQciHL!{6bb}J3x$g{<|loD@=H!i%CD!b2(
z>*qD(N$a4x(k~DSa~bEjh1@7YpcgAutj-BZT04+mOZ7(Y9;6z*V+Rl)oI4wa`y<Lp
zuc0`5J?hK#-;BhH_{{EJoY|aeRy8I(Cz_dTKEcMR3yzXo<D?e8rW~Lg<V-C)HWZyx
z;Q~K3^!77&y2{5ixlBb3^G6k!lGBx(yy9&rIr{4~*eX#}=j7;V8JkRSY?=1{CnLw7
zY`4iKDkQ@{A0|HB-OCpg%)m>Qw%4#lr(%Z7eE&3bIFjno@K{PEyUaQ25s|$VqeFbx
z28x}%hxvO9RhC!qSlKo(dnwRoyP25aFdMi+`@wAhbT6_P{jIB4plGn6p9y;)_v=|x
zH8WW*jSgoz&d%HsTQ+caAAU(N@a|+aTF`tooe?T;_w-@9j%1g4NNr?`;=n3;$_tXA
zc2`Bxo@O;&*aX%@hP~rXzdfu$+`h7XPMG40c5g;NNB^XJXaA%L8<pbe*dCz|VuHS-
z7sd#%7kmDyIf$vEvdrAMLG^u7hsf9&tt^GN2PW|XZ}YzQ^(97wbz|}Goxq)Ozcu(k
zMn?Ky!k4?=iu+#Rr}F3DR$~8JCt#qP$$eAnHvFU5jP;i4J9rGgfMdDMhwg!6`I*Gv
z-n8i72kwsQ--7!V2;V~or2iQYelH0BYvJw0>d3a!qjuF<y=#+zq`geBBvJ&2;C~<Z
zT!ux6F+^XA7`S}l{g5+pYkU5A$jK{lds(lhjTaJY*O$^zApJWC&lGZPhb1<B@&$O!
zpL18VjWH}qZ+qD7KrJ`YL7&74aegMysh%Bjdm8$Le+Y5$^kRPf>w)SLtHIFci4CvQ
zFH1Vf5tNT>=Y48Q(|H%hk0}pCstMVIx5SD!^0tZ_f@9fUuxsH#<-GK+IDkxKL|B19
zfGFTR6c>FI$bCRepPZ)WleuYJGYnn$QOg9II)?8(60<!Xc~WZ@P8T_#*M9SXVXF*`
z?ieY|Bx0)}mYC!i)>n{~IfVz$Js~RtcjpN^R1jvjaoM61yO*cs<O2GNvm_v2(2*8D
zq{MolSFlh$$U)O6u)vqVCAI1SJi#;myl=5bhY(79`FmadHqN($&j|bx(XszQD#HFV
zsmOnfxZhFzTXDaexVb0o?Efnb)g9IM!2N-tV*g)hsP3r#Ex5bBwtFMT`cvS3+fnz|
znyDKXLw4h0tpCu(+*18paCgn$_rS6Ij9A`pclvA5?WC$D6|kWO!0_{-nslL{NKV?J
zqs2sjxxOwyTBtvjcgFu%^0sLv&+Sa*V=|+miKz621%U{`Q3lXj4N_uTrbb85EW2p|
z!OGIi7_oG#623k0+C_>QnRcYp(@?)mTPxV2<=GBiX=!JESo=`%bPa-~QMNsD20~=p
zs|7V`3v#a7kf;L0;=w_AN@H6YQnmH$_tkkHfROg>FXV^?k*5j?T)-TiDfE$~bv>ou
zLq6)o$@nnJulH))bf>fAL5XGMss1d4WVR10FXW2eYY+Z~xnS@aWla+QdYlkAbsRUN
z$RzTLJS?83!7;5tTuV1IlD$|ut)d3ze&=Qj->x=`LR@RK`4l4xF1Vs)nkXed9*5S%
zxARe-c9Y7lbf@+yS!wodS0WtKrKpz1PxgtB@LhC2&kzkHN>d{npExn3)ZrE&&tvs4
z;B(52V>x#S;eVO_TCTw87Hoxd5j*(FBh+%fEw|p<z2UOl($~S&ws%4}8^XG8lQ_;5
zZ`N8HD`U+nP#n-*_*GaN97Fv%mbGm{xlVvWU---?(h*b2i0;dD0!SCi36F`c>gOpt
zRW3qHkP4JY2kV@wE1Scsl&7Wu5zm-cgw)oJaa@@j?0uo1c=+E&FxQ7XzzkCPn10}h
z_fafj$siK8LqBuvHGNLj<8crU6U<U6TAdzzhG+S&$tJZy%di}f%!c(-M22c2296Y}
zS|plx<u9}%9)5Iaujr-IJ46_8=rEH#jVg8EENO~d#;9+Se2U1qmV_V)sx<SJC59;S
zi}jZBJmY+kgCAb8n;N3URG|?T>%4ZYY=-F-f_dFS&u794Hta@W3;!<mxTuMj*DO9W
zfKbL1e1Q>9EFY2g2k7O|r&mQ9%sQSl5l<VqImZhcG8HMd=*ZX@3;c()axXG>8u^*A
zlfF`5Y`C7z7dS9}u(pla+Gu6gfs)gF>EMn!FxzO<a;d)fu2_I+ZsYh})$tI+xSh^(
z^iK0GX_|~|^Bf;RTQCg0UPrt}oBeLmmhdt~k%<AXSeu4q-O=!M3B!S!tFIQP$!j}i
zv6X{Qc21hc;st|<`LL^Kf^)Z=6SoR19y5vid<KrNXk%gTH!{q@mE)*Hn~&E{EE=9a
z_1l7*-z`zVf>J%-Qn85XwFu59pSK(=nS{g~lBhDHjzd<G57c-~PEfcjAUwZ{GC31Y
z+##4b6Ob~W)+p0DJcf6jxF^qT$@G$1BkSej$DIRrpB?C9tccg#mLiicxqK5Mxg~Tt
z^W`{d3%iu1P_#8d8HvSSZ)Kp5fymQAFruACzi;*{Hj3<{2hYn$MKeKs=$j6yp6{)U
zT&PWL>+gHeB#=-%#{Quz%UC%kt5O(`vFZ)g99eC7??bNkvQTB=1^8XTZHZz5Qq0{i
zwQvu5q`FL>Fgo^yHDv8vG3`Ry&_cEZ_?74jVU!%>ITmLqiad}RgF;Y2gg~72ZNF4F
z*oPU`?N~f`s|(ulK<&|Ie|nxFvXs{4r=+>BkNJW^VK3$U!P(&5*kMPLvOlDV!vrqo
zq_4B%L$dGCb16L&?A0J_(?fie1^t*chL<Xv-VFLI_BAQnCd|VpWV}$PTO_FQqr|cC
zOND;qq~`%Vy*yPjBW<-t?_0>lL`bHy6UJ%)Z%Dn612b|$=~0)?yw7Sf$0fJ}wx1Dc
zUXp`7cSP-OK?~F*abuDD64L8Uz+;5+#UCrXscE=yiO2rkHOQi5*x__B=jep)7+v_`
zsr6At_m@KNOYo}K)t`PJblutYJM+)T!2C;&@om!Oo(BNZ|BQe8Hy6^>%<K!<kcB?b
zFTrJWLa{h`JnCT$LUFsk*4gQ8*fS7Qtjc%Os5(zA`ApAfeq714Gjd@tN6(g=xa-8p
z0%^?bU9*b`A2K*EZ1-6jes(VOxebvX^L*HCJxDeCl!AwY;~6s)M!ALmbWM;dYi26Y
zA~rMtJi+&^p9xO|u;@{!r8&Q<wl_u}bSDd|3)PJKQ+(axQO=mNPkO3vtD!BQ1KGm3
zXR4Q#Ej4t*t{PNXv(PQpHK<h_p$CfOi59ewIc2)W1JpIm1`EP-gxND!EFvEDV|-y5
z!fj)okh<2kY}Ll}m-HPLL~|Q$0+-@!tfSbq!At`>^KaA!GS{&Tb%NmpKN)kPDZ&yt
zC6@){9jw;}_vp4eMTa<lGBg!xJ%DC-wlEXyu+FE`V;<LVxv9q*vuh<>Z_0<;W}aQx
z7+y|&w9bGgPe7z+>~#8}?sH~*%9T@W-<t>X{c~vxtDY*N%~4xJ!qjw+Hu)Gk+-1Wi
zpeZ34IiDA?HS6gX)(mt7A&H}^;NTuy())DB&3X|Z;7vWcDG(9OI%PrHnK)3sL7U>+
zl^HJZ<!A?W<8H+<WolPV7md!;6>(1p5yW`TUg%6dj?6hdhp)3{28|{&@w{-uK18B$
zS?y)20GUi=kLtt0^o=sa$3vt7Qj1T?gjN``?W2*zDu7Ku+0fVErc}mEe%qxzM;HD*
zE<<oO3>sA}YiQP!S5^m9^B$a>gSPd{;=Mf=&OO0f)vOT;Dc+F92O&zmeLRU>p4<)t
zc1DAM)m?ei&t7wH&feukB$DMZ&htK#df8ebIZU#WJNsN1*cI+qqIhultmw4EK%iG}
zltNTJqU{4?+=u4pl&eW03WMFYuKk8yvyt2ltM8bNYp>SMDGs;YCQA}LVnmKV4MwJT
zXfIqDQ^fjBEuM(SJbI&-&@~A?lbZJ6xYs68S)0Zj0%7j}OIm%glbG7s#2KsfixKyi
zrpuHoSuDC4^h~T|ZppCXCn}?h^32*gA%?K>Mcp_SzG-I1W&>OZWwklBAM}IP3bE#2
zyXeC5tqdRW!$K*5q}KC=%F2PhmP~7<B@(KwpA4JG(dIGmu6Z^K>j9H0c1nIg8&^}N
zH4MBKDsg-V^YjT67NNE4K!X%(6E=g5v4s=On=6npw$9~=xZZ3+h^hH8`Cw>Np~!1~
zzfBE%@3;Pwk6>OJJ=y{wjGqyeCO;fI>Lx;SA2cS)$3+tlhH=`JvOZAfFT)`4G-P1v
z*C%`PAU7yk+M=TzPV=3TNmCV#T)*z{<+;$#<Ep}r=kNC+FGxkMi4pOcP=$mEz-XU|
zYCo!dsW8^*7`(GcCj;qxOv!_7Z(5qs>ZZ7wnIwv@M{7Lll?}A`f@;5ySWdc1qsx&<
zp$N^A+q>xt72)v~jr41u!()HMt}m><*kbDg#1N{_Wi|#FArXWba&@!EH{<d{@Ly6w
zl{mLYyrm?yJcYjaw9333)U4<Wpu`+h9yB$&%vw-Ii<2uAFzq)L?TaeIn2fjCv$d4e
z11b~?wJP_jH$BSg)Ltg^hu1h86V%&V(m->4#g>=|KFL(#V8*QS!B^W>xvcW6(^S*V
zCscI$(`(6USah%Hf&=Jv+Rspo0m2cl2s(335s9U|;KN|?69gYnQ|C-(ljEKSl{t~8
zzfQgYergmz^oB%7Z#}`q7*H;&!ook5ZG&9$%Q`6bK?AMTm7BXpXL%jHcdcr46kV5`
z`^Xl5@awXV8fpLXr&FLw%$ud$kYdBq_?F7Sx8B!7EiE7?+|>)e_X*!D{Wqt;^rP@}
z*Ux#+C$RlYCVsyW?61uVhSc9Dz7#_yi4zu<*R<_;F!pBBPEfOtxZ(P`kkr%!AC7jz
zAtT)_KZ%XsxH|#|g>|y>x{YAv1%HojuGRDC{IS%tNJKWHgwN+ogfewaV%Cp)-+U&t
z5qxx3zi=6!#ACR%yIN-g!V@v9jn5Hnv`6erdmxdE(I*Y}23K=!X1Qy_aB7n<$+36b
zyH6M)n6pqSAXyoC;G-t#Be?QP!E@ZliH*&%&g7rw)COEY0-CB4VuBzPZM^#wG0K_h
zMeA%Hw=uby)OCT3xAcyXewe^rSQQ0_<pd?JB$Nab+kXe5pwXKb1!_pK_qaC}a}}09
zAHp9wo*uF<)1fnSlD+rJW*ytV);R2`G%x8WD9G60G8%ZAapW#JRRp685xjoeekG~E
zaAELBzmi1w79}G^F7&1<Gaq{l;Mxx`B2l)Xz|*?w<BCq(IhXa!9}{cCIuhe_9nuP1
z9j6k1DIq(44XaQD>-{CDWEuLny9mCL&T>*uwi*~A){w#L4tKbIHtHx<mD<eVTnI)b
zV*#lw#KBZzz&o`SvhK>xLtXo<l*uy7%)qYA^i_;W+GQF{<sGKLFJ9VY>Ja`+kNV&#
z@zE1#Ump*rQ*>c_H#N!1j{3wt-LW!6i-N{39M*=JoHknr$ds4xJ3E_kJ%fFT*Vy@R
zOPP6sv+xk4joao~h-!n#(-;$#z06&H@Xdlc?xq5lnNheVh9)c?B__vEC^h41-gHXs
zN~!?sl(*X4lLm~I$->!OjPt@T`pJ*5ReW`f*}i^#6{|SyOes=gko&G)k67tp@s-Sc
zKl8e=L0RZSlD!Mt$r6n(r|Buy>{t6|(GqjL(^d8*qSeQ{M~rT-@@se39h&wy6}uzM
zD|lxOyGs=}-@DSaIaEx28mR$l?JaHeBpY?_+M<c0omcH4?NgoiPTK`)WHrxp**~5r
zTv#US!wx{6V47E>Dtd#ty9>1%hwP=;`ns)_n<Kxc(bn2f#M48n%cOz)QQg}Y!uiv#
z58g{7O=`?fT5KyHE~Y&nKBJP4in1no7fnA8x55+w#tJ2o+FkoLRW6;&Ia_4cs8)@j
zs5IU^Ooo}jw#wQT{bMthJdgLtMEDc+xKgq#yY+RBK?d53eWR`Su4Fkl&gdpF<3VY2
zREkoU^A{a;tMS~+q0{ScKCOAF8sx4Jkdw~F#5o&CuOv8OQn@0|)UlFw6r7$bjSQ8g
zI~%xZs(f}5&<Zwg*~oIPPfd#ull1~R9)9_dADOMG;u`dJGtH~hOYe<g-ciODwpp?Z
zO|*XGYoH;;Dp(vy&nW>E9k0#;>-g-flI6wuvR8|VqZov3a&?hkP@e8gw!fIowwFA8
z$c6dUu~x{$91R!EtC3I1a;XTgkGtS!HRK`j$YelYY?i6mzmZkKWrC|qJ?PNjowORG
zJd&`U-_N!7jT6P&!=$gknwxx81`@-}n7Yh3f9m5E{v{^DBS4I|o3cYJEej@^Q!)l3
zIg&A4QHt)J#w_Y1mpZ)QcSp|7JK7K|9hKB#z<gO7D+NpCHk_Q1@aQeh=Y0Hl`4aXd
zxy}i4v!e)|8oWZ)*e&Z4QqVnJ$5LfIq;ejof{M|d-OTL_lxs`+u#F(gwa-N37$7+#
zC5c(6W8!@ns4K=}C3GbzHCKn88AxTek*rHqL&iNqua)c-Bqn}<>Hk51tI@9VlC1W{
zpdqzogxc^!O@XED(AY%wN9_%GnUon6GE20jIVALvUoVbZOknzx9EKSZ0g1LAiLG4q
zHZ=0#pQ*b83|?4WKJhDm$sSHSnkvDle%MNbuJVR-+MLM;l5K$WtcKzG)k4PCYa}n%
zdZXWajyr$w&2g~(SUq>Q;NiaG_?eRae@qFN?Uw30c)xnZKd5rD{Y;|tyI#$|)-CSL
z4_N!jVG8T*7841W6`|D$mG51Ll2GC~M#xG{@jA%%CZXk#g*LuJ365ie^EA6W3O#&d
z4W`pX72C6s^VSSID3euYgS2J4^L&E~No#bQ?(nI7Pe@otv5k{aZ($4Fywg?c`B^$C
zVX0|eiDt{X8~IdD$K@l#q0M(7O)MZ4mqigQ)U9FjrX667)Bkrk-<UIW@WL=Lu(8I8
z<GN|a1!u++?SQA=^(wWoY1qIU(Y52!eugW(Txy`S8R>%gYV>VXl3X_1?56_H>B8#s
z!-Q#5SANzr61@zZHC}=w?+0F>FONTY{ef=&Qs`x3|1nbvr2m8UK|uPSS|Idm)7#x*
z;=fpY)dvy_#E^4bU??dmiLAn2zuw9A%p#H-J-pUo4)J`4TxE-XDrMvLC{HXjJp45|
z$<~#3IB=g(Y78ee+}CIV%bD5IsH){WJsDHQzO?V<1GKiIsARfSfK%)~y9)ROFC!0Q
zrkh&F#`{WNy$GMY430o6+OtFOj{ag2Zx`&N`PGxq#cn`##qv9pK;?Q3ucu0{Az#(3
z<PiocUqF4Ghbs4fyf_Q1<j=u%sD;#FJ{;if$X|r>x#~exPEN$0$6Ltq1%J8#Q{So-
znM&vd{xFIS9<eZtfHD=6Zuq=ttbIeU5-Ne3_1GSpEaToia^{e3`YEFR2FwT}i=k|1
z-3KBem#hZG2vMyulqWP0&E;Zyl>%S+++*2L-2sb(?9?_=LM+0@eUj^KSjxlV(>yLf
zjzW^5rZ%^%SVSv+JOyz$q)hZgHxJ0;<WvUHAQ-A3N)8_c-XXb>R}+vy$(RzU_CsUR
zil00YZ6L}wK?^__FNxvkvnB9QV9-1q`TMKZ&0sk!#GUE7#@f(TI0~r|XnG0C1Knv~
z!c?Q#Q3Z!(TL2J4@=r)tE@EJOafZ_mJ8JVmNX1#^z&G+_;n=bIavLNMO~+V$UKttG
zMur)lj*YF&2pBeDXdZB>N2sgQq~=@M;9$Hr<0%ql?AjcWQH#O#=T|JIdS*sL)@ayu
zxHMZu*Jn77T}qVR(;)t&^nyvs%%O28VxuO(V5<62R>LE=zM~~(i?)YNIXeUum#c($
zM1|7M&)0)MS5sHJZFRWvlEV_EdHVzhj~Q1Ohsz&<Jw`upm@2yPaA146zHigzn)8kz
z-mr<W0sX0*4z1##o(tWFvL{I}a=oKzxUnCaD=ji7xviW%2hLqp;+W`D>g~j}u)HHQ
z5@;?3hGorSJT9EFL&er6(#|O;#t-#X`@$Zwt0OZsJj-@M<J_7r9y@rEDE6*E-!sIa
z1I~{)5``g~Bk8eZX)@>0?hbEn*K4g~2XPgR;SW+4>F6@d2h7}&W_6M>YYi>fBeUFF
z^GVFIbXjRw(halqZb{mBMXGI0VyAI3ZyT;pI2d_eX0YK-BGUR%&6jeiokZu;<!5h?
zJPq!9#o~3D*#};HfrNqvIJRkeKU0$SW=)2^qv0#pJHd*V8eYytbHUXFV_&@hU}6<U
zZtLfcKqe})GxZk~#lT}hyW|t><A7+1Av!K^57&gwbDHNIo{e=Im)IN$lUItO*__;5
zVMm*%*p5WP68eZ>u>(Eb4=%wGa>qVu$A1N#MpoPN@-h)d1<|F<*^2hAJuh}Uk0VyA
z!HD|sfxnz|!cDKeh9$lM=vg*Gz_Y%U46Ivoryr%F!1S+|fLqjCESk9>vI7Op`KmkQ
znmA1VG|3@~h@1XvHZ+6iio~m$k+EJST1H5w4vLXeU587K;WjdS-`5(%0`Ky$UQm6^
zpg_v>N=oK&{jlL)@1&}~Lu(d4+b*iy4@O@@CO$*{0jF$C``i%<>tU?Y`iIS@YIx~O
z@J)k{M~$CcEZREB>yKT4c6IMPl=$=%WQaFJADr+LjMkH-bfHxlVd6b`Hw`wNhwbB@
zDA-l;&XNMENB9Mrr{2`;^DQ0(t(J<fM;5x^vK&9pse~V%T@`K4@k>jwD=xEB1bdEX
zVvNB7@dS>gH%m16Sasas)tv};&}*}i50Hqo4Q$&zLxLAk+CQ`iyPPxbJHrWR52NgR
z4ztf@dfAn=(3LlUikDy7hIdx1<xADnACp}y()-BtODK<@LLKW|Uto}lH2VJD8GN(x
z-y*eN*Z$ruJ-ZjFvHy%SxVPx-FZKm-TDQ#?foAWV2-<xJ!66lW=|e1g0qrc)xox^H
zYpMJ{UtgneJ|tU6({_2aVQzM?zIAv~r1T}DWgg_i<4F(!RigMlh7K@(p?JIEu{WQz
z6}xotiO?)ozYLgqZcIXRkv^Ze*r0ZIg%=?Ly6tk~@gQC}IH{qPb)@1M?%7;evDKNr
z0K%xk#Zq>#xzBUwh%X_iD+jc~3Rm*TI6r^CV;5JnVSp1VRJ7GqsF$~KriWYF$d41)
zX};0Dlc+r3n7Z2A<v8M_Dq-F+Ns;(kWsk!vepM`(BsY8?7KM~Sv>q30h_F{b21OJb
z7i+!OOEL(ZDw3Mz0H;pPPEre$t_?}DF!p)JIvf_sc0S3z>{?=|^~wB$C||Q4PQnl~
zk%nT$^6AesDezA?DKcFq<0=-fFCoy`NkD!t5O>3TJM{n1d%^Mt!2!$9h{JzO;C@H-
zZv%Ih+m`!#vd;1|S;@b-$9-1erpKKR#|;#8O&w1d*%w{zC;;@8kIzzgrq8&-!1Z--
zolhlG^8RWqlB@3U`(K;f#ZvfS`+dP7fCZ$2mtfHoMp|naK6VSHi!bcp`1&fzt~8i?
zZ<^d`8!`^Pn1?6;3zZw(+MtJssFYiP5obW7!ALqKMF5eH4{N8~Zl8mS7&+M$Sac+j
zyfZFYrSDo%^yu3wKg=d_^WD+Hq4Fj_el(w=h10;t&W%#&^xWvVcr?p2H%;ya)->FD
zk8YaWk6^wxxjRwpE-X9a@9tZpt{j}GVsWA9Yxs6O&QcB*e~y>Tg}R^~5Es$99<%S#
zgu>=RE`@K0GLFb~*$>n)r?_o$FCaO0?ZIY!qqpk@SZL(XCO~e$?@8slTmp*|4T@Vd
z#7J)3*k6t|TwUzKx*T3>JEsdzF>`JIs$K)uJXA6e;v(L(jhH?U#H9;fYPdA*g(5OJ
z_ai&{vi;Bs<oN)3o%@4jFl2u{S)mVIJn3YwF`VnuFdGgj!!{CH-T`+sdSN(EDzhiv
z>@9C~3_*1e$bIdK<$yC|Tj$-Za6-xeYM&0;rGhLrxP=`S^YKJ0s>hsu7VyKA>0)2x
zc-BGnO6-LV#d-kvwGHx??34|Z9BRk!wen6NZ-tbJ?bnW|KM-8@pAp>OW}G*DD1Wh;
zp6^p==ex-`Ihm08Qra<}brLsUQx(7`!Xa`NR(-sb93ee^n_$Ju?HUCkkoZ-g)=Mw=
zocK5*sA*v1+1VTJ_=}t;iQ7+|x?#?%EeZ|#Ip1nk1e?6DRH!cP;##9od3)hzdy!hu
zYtX+P<+OA9?o-s^*{1jII}OYtoR)Q-3YJr?=C1{}nP%5*(y2N3e~0szvd<@8W}F<4
zEwA@*h8(ZE#vkPxIlEUqt=wDa*Ba-noG<Y9(qV*8R91XKQCjkRrN~*$HE|?VKjkbp
zDNmGlQEE{C6Irv1=bJ?D2h`pYOW^nqw@rU9zjtzZC%^2!R!nXenBFtIY(LXo^dI9@
z?x?;+0{j8)+w$3a=z#P;1@AZWBscr{{MG9De+8lcs7lnotswMGEWc^zkBZRUPImX?
zm-T;T%y&oiZ?*GoGqwBRergKy*G)bC9mTGHXl#68BHa0~>V>c{an9H4mXH25yPlJ2
zLG>1nr6>9L%?Y9Gf=(y$6JKY7#-OfD!Zi|o$LlAy+{vvfry8;%Xod!tnbQ*arw<6C
zLLRiQrPM}A!A}*!J>5wvt5~T@U~`nKEl0aZ^ewimMEX!M46L2x4to{U?*IJk^~cfv
zLEKmRkUfU6Mm*S7?v(fj^rV9@!VD04O<#PW<eZa;cC$c>iYaDvL0qpCac9Vj;<1O)
zrkUFkBBr!0adgd5_t$%LJ#Z*>qFGr!_^6;g3CY|Q`QW!6sGYJx(f`RT(3#HO0fP~0
zTU$Nfxi5)UWgfeWbCJ1ib!mYzNK)9@SPk|R%)({hae_3#jBek@T~VjgfoWHI3zvlF
zq%L)zzi`OrqS>0cyTk<HyXZZ&i7^_0tVS_*a=vYB?Ej^)arAp*<9WHl<D14tjGM+r
zOYfV;#_sxh>r<z5gZB_y2TRd;VInYGv18(-huBVYjSM^kv+!|z;ZLr@8X;*bSs3S*
zC=FU*17fijZW|lPBfi`<Hj>{oHnv4RPbt1_Yy?#xziVt<zH4mU6bYXpk`z>(JCZr@
zS=b<^dUFAx%sRD6ioV2Csjv}WY9S|V4nbhAgo%zoBN{-+j28nr%H82&WOvTt6hkT{
z&Yg7t>(i0xEM|Z{z~rV=D@CV#QMkH|U$wACo+ZP{6?XZMWuZ)p^~IFK^_fJC!Jgyc
zS%CGP!ix^4bECF#NJ<*-=VB<Lo;635dnwi`Gi{M7P1Hf@$L4C?(lt1{6J+Y{{_56y
zD`(MaxCjd`7j&NvudY;Xw#@sm27i#s4JzAcp7RkI`$(ftxObtJXWslp9ec&e!YN1D
znZH3{T*t}w=A#ME8>CM82VDKlCuy$~X@?Jo8<e+<I_H;&CPsWMniu#!2aoT_Q`;Qo
zUL9C?^47qP*=D@c-k4$?5F7BcN<}V2qpR|WEQkR^Y(vzv6fCh=a`%77dT#V=o*`BH
zB+g^D4Z`1;NrgzGx%gBkucv4sokM)WL)m%}N*d@O29_f~;Ai^L%Os_QgGC??0+Zm-
za)L?_A^l09JE2ntY^*f4^Oxt;yfXBjS+i-%lq(>Tyz}O?4_y(q<-k4~dpf;hXf;(i
z$Ss!=kQ2tT&U)!V-2jK?7$+)c!59mV<ZF&$?BsO70ZJw6e5uhadR!(S^~N$F)T^*z
zU=P%ZNXcaK<w+UVval-{6`p&1tCK>OgYA>B!epYkZ=H>sf7IF7@V&Ef;d^Ielj;&A
zc$a$srAD@AXPn?r0g<3t7>BXB*^^pXTfG<(muGrtg`tgDb3s!H5&=^5bH`6HU?WpP
z$4ZPp#$q@N4qTYcAuMpgd%-=>>I$aS+~)kKUg(HrE>6KQs=LMqS`8lX_0{U~qZe&~
zEo+5sVCmmF8yALpQsn&dZu%ONLgjo2%xytz{7FSYIoK(?r6lON;@na<&RBZQf|j$0
za0M1%VG!+;p@lJBLoU^@0>lMdqlB2-_#caDrd{IB<%yG%XNeP2v^{wsFLeNi_uIb4
zT2ZR(if?_5q6%FsKl&QA0ye(&HEugrMtD=f*GVw_wy&}3wy%->wy$yHRhQ?dDO^(~
z)k(OO%O;$wyt1z!39dmVzgF-3KJsL}HS2fQpON9m7F}$AP|#reU+LWLsQ#^UyW67w
z7w7hEH|JmX@uxG;yWz*LTg%_=;s0M=%+D_;#lg=HFxIoQzxxc=_m2+%>HlCg4v_x4
z2bH<`$^jk!?Pt8KZ$BOIeWkyG=*_nX{->{d1t>4t%qpSsgxDzLpz{J$MT0`SN|hBg
z>Ts}&yCKEnUhPPrzuKCnBwpoU4avuWn%$}{kahQEg}emidk>Bz#zlw#sqg;&q_t99
zF^IB2T+Ds56N!*?=YjHagOfg)VyEWumjU)g_6^Pr?O7pi9kZ%bEal+D!jiNR;)J+7
z^!Rk?(D8~%N$|<SEg71afuV=AZD^;%+sk;N&;vU$yI*Ba9=EuK&5eZ+$+lXf;dx@6
z!OIU$T2xD-wm@B+nbknf9pdg!ut0znf+cyrc<3_LL&hyb$O{uWO1X_z4BZIH)i`ZT
zSEHu`j@cWxtZa=_Ph=cZQ5>shH?m9@^|S_fjYRw0@#vG5zv{({S9J!pNFfVUCrZsB
zyzi=xaQ7w4bV}>Q^}HXPMSU6f0CORbKgy5wbXeFB)zK+NU4UG8;67%X2Pb`D?g(_p
z#Vj8?CfQEzea!aG(d1!84u;R3%30;K!Dh_sUPZGdy`?q@&w>=4)Tb)z&*ldz)Kdn*
z9!D9r<PEL*K>6@i7B}8}|7)CsmGLF|f#eel>ymu_nUYjBdAKWv6i)4^UJ9#*g!XuX
z0mOQzEUW=^iSDs2Qtx78HxrbZ=?4wz7&Fls^W&c!ATeF1c^f(C)**3cwI&txtFv20
zH-&^A%b?Tfzs;B6%#llHi0^C?Vc{ZS=zSK@F`x$|k{y3&Iq%9IwbJt}{0&AE%}e=L
zoIRzsEq&#xyo&EZFA4N3KGEZtJjy9!oqqlaTC$C!4e<f2ofQWf<%B2s@nNH^qY^rY
zjNHK`pC@^+{P6od>@P!6UjtLiB&m)b%Z>t>G<F^VpdOhzvX1lQSy=FXLgea)Y&aX<
zX|}L%xO{pTF&}2t!%}A5UEO{*U0P0aX{7+hsp}S+nAN}U)DGkF7>`>Mbb<t7+h$2^
zbVgyOJp0whbKlX&q;0&tJcHwP%cBKfUw+NEey#HwV{(zx{23N!F0m)2C)E|XCz&VR
zl~kyMT0PMfu~tHJ4DAE<hP|c?l4<u+9lb~Mp<1SIEi#zN6lMC|S25P&FB>o4aN(4p
z*6HBK<*hdDe}Xv~g1^FR!gPSGFqf`!r}+$l$%I4#nFjf;TdC_=p#!?wY(`wJ0!wfA
z^wHxyC(Rc(<zN?%^iVRfki9V3&Ylm9U^VHDXWD0^Wx3$TU*jIpb8swKQA?5ZGJZ(o
zc&Z|N1ncm8bV$n!auM;Rc%jo2OF!C^_oW4;6fhqBwG3KsMT5|Wn9xSn$XT3go7%&&
ziJSe2A~*RLE}Dg%=Hab&HH3n$+c~#YOf_CV^4vi_Zo!@j+?VlKXM7n3bzI-{VIx~?
z6*5MVqipd4Z%gvSyhF{Hi_fX@2$>Kk0W0=(2`e_I4$(FGn9(7x8pu<Z-y8LBrsmG<
z(f_(d-7ifUH=jWf7ZeoGvo|ul`Ao>IrBVQB-9F^*hMTIv@5&CqUp7Ji6Q(OI(OF6H
zJ$lec6Km7xLm~jl5w;r)V4PD74MILd+(!;9@ENI-&;t=v6bcnXrB+bUAR+@{R5*V>
z{dX`2$a39n@<PIE^WhdSd&SI$3#nH<YK!kTYTms!+4Fo?4<dKz2hr}J3hKj@rHZ<{
z-KP2Q(Lm=Z7$!0(qAF<XvavCQ=$R1M<Yg0gO3G`o?uwILS|7uPS?U%!Qv0cE?-(Lh
zL1IV{SmC!UTwbIYYhdw;6<t^i*AV?j2^0*0KnC<In3qLFE1jj5Cwc2s4SmT;q^nuj
zYQtEKOzGhbato;>%h}-&S(?ctOD5SKP)1Ek;;A5Or{oh8XD_5c=uuO21~*8vEJGNC
zQL%PI4!}U5<2%_K?GVGx*-g_yTyjlTLm@)M^>b0)C7Z2n_8NcrT%mZb54`zcm-KS5
zi`2+bC43}3hEJ(YqA)0lnpKKv%uX@x&W_sDy6Jh#`$c>&#VhV!w+rRZq@PC8`NPxI
z7A-#&(KVrioDi!&$E#>gPzI^Ag5+CUE|T2?fE+vksa$w*h#zRF9ozW<o(wMR(<e)N
z&<qP&s!i653_qLxmO%E?l!4s~70z4?Uf;*!XpBa}Yq@#3l0t)n@Q6txhF7P^*DI~;
zW<z}))GsnxF+K-uD$8%i!#oQyZ3y7<RgI=Kgnj1ycG(+~g$xSIE6oDc)KLe-Gk{&L
zRvfC=h8sliF*va`uN#=>wgxf~wzs*_v^b{JK8nHZ>&6!XIL%P&Whs6l#q|^TC__7v
zm(*C<i_a|Ez#c|JIb^~z`O?yZQR^b$`Y_spF91NE`G_HdCxX+8ftm3^yy54W27%=V
zPJ=r6-vnxtzHXxsGJY(!keU!rS`lpF;C)2-U<2DE=|LBK&7l4Ek<niA7eu4yAf*UQ
zA>*h+y9)Blv7v%@35G-?$MT(zlD@{I^e@N~l;b{n?0{Vlt>ynE6KWK|1bGGO`3xl)
zY(fuV8l>9ieHkY%v}pUll8_793n+}%*=31ZbRNW-c8q15eGuvZo%T3rzFuhXF%ea0
zoX<i6(W=CtD1;U<L^)vj!l}`SIXHC@Vlhuw!6AdP^pY{ryRlY+)WVnbRrIFxCFtKK
z8ziepVXU(lBQyG*`@PbCTb`>*R>bQB!-DcWf~HeNH%eEdEU6rPTJIAlE@DqB)(2c$
zuUbq^25aULh#I7K;SE099SF;$pIsU$7f|>7klO6`^*soiVZCv^@fV{Ax?9MiVbmck
z{E!F}0mPQXO32vI8_-wJP}-T~eP7Eu5CtN;c53QVlt)hUH;Utu%EY0H2NGi_5L2X7
z`YMT=6KRohlQt2?^>G-A+F`2kY0H)nX8?3akYah|;>k%-WhHQh1>VcBOVshX^0^AW
zQb2fZJ*irdT%V>c+a%;E;YlYN$Cz}Ow9;1`R~TcFu$|;iA{+lM$vH`fHiJxu&|Bhb
z?qw~Qd>dbzqHlJYN|wlSWjl9(xNz#r1pQ$BRQ(Wr{}0Sv@+6J1zOgkAW%K#A%En1v
z5n>BWymB3{t(30hm@%56#}j4ExtNIC>N+CWY&wE5A~DifBVWT^!(OZCSAcA%*F9rY
zz<L!P7v2yab;z_{H?lG^%Tj8lececra87AW>5)K`V4MJ&5LF<hZ1WcNZCZg<f$Su`
zIjOmx`Frz=x3q;3BY7jv=_6@6OxLTIK_AQNL)>b=m>+k1<o+0kor>LzjfJg@oy645
zjFoDg+LyYPy1~>_X{rua1);vCzEP!3pOCUUsANo15yqZXT(Xl#GscrMU&vi-@VY|P
zBwMCp_>opDtp+Q9jb@Favw*B5ohqG9sr+6Rk9<aPMox=pi-vc9>nLDGH;n&z+v}{*
ztQNg=z3gR!#|1HkW0cwI*__>M2GvhPKl7|;kd}p{iN}gZ^hqmDPcA6w){j(<z+0@F
zk}?%CH}t$Ln^RLO5}C?m@8Ih7$aZ~xP6fvurW%%pK98<Ly+f@?vrSE)MxpjV)vEGi
z*_V!)&QuF61J9m7^A&@N;oSN4s^WZ_MW>b>PtH`WvU&6dnTC-Q#1o<m(hDFIcaTd^
z6XNN#gA;kfGtaju+ndX^!=J-fpvc_G%Ku+&UmX`!+qEmApdg^qC5Tc|lQ5%%fPfOx
zNJtFbAf18;0@5uh(%m5-EnU(eEuA7D9p4^(p7(ih-t+s;`Qr@CkJ%S@?RBp^X0LUx
ztN%&*RJ=+zYd1l6GX6F4P4WQ_3aj21GaqTJLO69f{WTpdKMW|2Ma-(M9(GC%2Q1os
zHv2p=!Z&*I@raE`f+$4HuO&P^Y>}!5U&ALV|5?bh{PNH}?eE6l)A{f7`|=0y6W7?)
zWV>?hqfRcK*d0!64{t6W4j)NhLR`wZx`dZ<N%WEuE<Zlx>dnixmm2G_>Z?5HxcW5j
zJ|7^vD0D1r=f!&!eodY9(3*mi9hSqi=J-+3;iZH9B@S(oN5Y@09tXAdwTk%#kV%l)
zi?)e1i42NKi8zZyiAY2#v#S)(iTbPh)2H!~(NU8+-?g*K-0Pc=@Sa`!hBv#x34I<}
zy!O?qbUw4ACofpDUggETThw&4Q`gUYi5oKKCG&PPE9~bl4Tn{EAB8S9vn>YvP~X=(
z4!IQNS>)OKDC)Diug;UaC!|kiFG>rRPTLk~ZoQhz3=a-ZW`$<<W_B5Qmh6}G3o)YJ
z7O<&9`IxJ@8ez#G;mUTS!3{;?q=&sBk}TFL`c>}&PbGgo_P-m;_*w<O*tlrYKuqs_
zyKLu;g&n{>$h?ELZM*I0C0rkFZ);FvclzeI9=>5hnCmr^_k674sl9rA<;>5T)gUu^
zv&W;G9iklpdel?qZzl(26J-l(v<!>CFEk5P2+7V=b9Ac_nl_n|nyydi7v~y}MKNgf
z{?zXH_V$}J<9O)ml#1)KgNe>v&-Lom-P9xdX}c@-;}gr3+se_}DjK6-@}|{u=R3EO
zw#uuXt<@-@)Lv*283pMEl|V9zs>Sn-#Rk(?QdVA31;s2>uQ->l$A4`Z8?K!)bL+<n
zz;43jzePgx$*pxy#w1tJtO+$flhfF<nzT|`Yv88o_Gta_kB}ZQFR|2?h9KL)2Lr=A
z89d?K4FRFeg|jbOEKa6>3Of^LFUvraEA}rTZX<tuYuTs!o^XwK)bXur)-v<S``xlH
zs&7?442=)@mZfkTIZHZEExY%}({|Hx2*jSb?D+3S#}v04Nc4nslg7m&-`l_35^DDP
z;pa>4egByae`U#BUwmFc-ekwCsC<nd6;_#J)knvQQP~E0Nc9J<@vhbDZ%zV+ZdD~j
zYOcLBS}y%=R=sRld=fdE8(4GMb@5RBh--4U(se%YES?(M59e@$v&Q*MW1Mf<b@HSa
zH<L(1pwkI)VtAsnHk-tZCw(s=OCa}TaDP}ovwCMOaagMduptJ$`6$eF7IBa--FUet
zvxiNzM07<gSFFeXz<<K2aHnar!#}Ej)~lBFv}`Y;H~ald`=e?D?Y7UGCcUfC0^*Yy
zHORxEU9&~qrMZ54qwi_yZzr>B1a@DYd^pnCX$e_glirt(cpo5O?{f8c%Yn&Ee_?9R
zbNlGO>}&r+W&`2=$CmVqkPp3gJ-4y_wU7PR#`bToPBf$6?{EER;z8^G6O|y^oMa2|
zxB+Hg2_s!=fIIMad`<O#Jsbd^<B9rRs4&U~WolsnFb)1H{Vxw^B|R%^BMWnSE-)wT
zUj_VP0zCh8;TKcjKQwg;0Q=nUg3<f+-~SX5mljl!V-vJ8LYYb_{3aKCq-U*TWn^h%
zVTA|&o2g2~9AF1GG6&ish0>9;FhiODgNPbgS=$I3qO8#5dR&sIzb+s^-rt;fz#Zd3
ze%S%g-v|EwfdcU$(eJ;n(SPAUK=;4j|Hl9I8Sh`8&*jC$xqo#rCXPmcTJV0=^egu6
ze^2v&js8Eh>3=En|5Hg2fTroZwqdF;S|R@F`roSX6lDfHRKidaQ<Q=Aug*s25&A`)
z$-xcy`)~l~nE_TvK%{_8mRlX*e*_fRNJr4zz*G;AIhUZd4w}yp2{=0gPpW=ja6kaf
z=MqL){>6300|dV1lDE+_Q=$j_qX5Qw6k6}_07c<a{+kvKh#mhc$zyup=^~e)oxxw!
zi$E-OKtCORU;S#wFSbnJNthnMV)>i)@?X53zmdPk%wJ@W;J+F9eqV!uTY&>14*%E3
zX%kM4E(#)LXPK4d7WRE(>uVCs6Y)=}9^+!)zJBr5c<`k!?xwzBq6}BB6Zi>TNn{bS
z)CF^by06{Iczu84?Kry@uBDb|FBUH9a|<sesRB(@(2Zug4W^SC+XsoA+tpKawP!EO
z5999Jt_;|%#Hntq>q$iMga)WSVKCbWAdqNvko$gH+rC@9u#;M9Ox`DaWi@KtQU^v`
z>(jb(K;WmOK_`{Icw_XeZ(sz(6_ZAHJ;?Clofv&DpW;~a=XMgX`1|-jf<^omcNf#M
zAE!@@SW37Z?Qr)^#!{J^`nMBNG4x50m9zV$;c7g5UZiWJQ_GvpSab@G*(2_JIJSg;
z?}K-nWSl~Kx{1#m6Y{m>V?P(-DVsC$Js0?x%2q1lgf_{`wbxh`9tak$(c_8kNos|n
zkzjxDoCr>Pbjz-~)cify^6e_MUPsw}X*x9?MG5Q&*sKY&M=63fICZIlwDs4d%V7`O
zmnhbI8$xTLoP~8N+{#-?E;X;J>-?}DU>#vWx*{PV*%$M%N-o-6+{F&WCSlMDS-s4T
zV~QJx#Ul4+Ze!DFpa@IuUVEDMr>CsZ-B=74e0l2*+*#d6-KE^uFEj-c$2Gab-A~-N
zFL=RND`rnoOU_m;6QTnZTC%S~OnsOP1}WRRH9eDl%o+85{MMc^l94g;rGe$>?TcQA
zhZ?ZJ4pu{v`S`+J4{W_BayN0JJp;M0^pCG2Den;RF))S_PT#CFwuG_fK_AEzxYph!
z4*U4o@RlP3?%g4jJ!yE;))fJ<7jrZ3w-XErzc1T-#4C$ZZ<<PvcR#t{TN|`_%PnYV
z6DqV=m+HO_G*ZXi`hwScK}5hNr(0CqF>O0;_(gr}GVERKZmiB@TCIp$j<v=*56;V$
zr`!eYXRlBv*duMyukJTVix6%S(~2}*p1`@$`JI^;dGg50N#I75!u16U8<z)H$VlrX
z8<Td$RbY{Co~Y0wh51mdsmZ!E%$)*DJ2-qn;RG`lz`ZQB#Tp(+ZCy;exPm2d0WTvg
z@u?h}7IgZCJE407R`vxhfs6ZW({~B54Bay?;0sK?(r-9c7SoDLXZNk`)=bqkJZ3rb
zyKxp#tNl6uhXC~t)aU3{fd<Xx)*YE>(gRP(J#Eq#{997T(($7&5*|nUefP%U1A0rY
zj-<}Cj>)}YZ#yT@F#(}a6ct$K&G+dwKlmH-G8o?opI0@J24{(s1!R1UC2@Un9V_p0
znW>0eVx8;W>zgc=i+qWf@Ww02!XLg8VS9L$<X!m%;A3>LAn;L+0FR~9?M!?3f^=`V
zSEt&-9FMTXU%z3U7G$q#B1Guu)>^x~n#+vp1uGJU_;PINOyav$>2z`fa%|xy>K&>Q
z%B!|H@K)7{9C%o%agHq)^;N|ZNbYKpZee7DYEZ!tZBD6IDUPxUp-Gp1u7*;IS4IA;
zxC3pH_?;DY;-@XH@jA`gKVMkU!B$>(9ThOG;qv47P_PmRP3B}+8mIXb-?x@XUS<;|
zW#^1-B}$qWlx!5dPY?_i#6^kX%C-qAEMDTg+=)YTArAY|1-z7EP@oi!-Ni+0ajar2
zlM9O%-uA{@FX5#()g}2o2nd|Vj=GrULn7c4>@_{h&mnA^Eohrm($*ChZ7QyRH;i0X
z?&i<}BC3g0e|JK2toGpRQ?MTRn7_ss0=f}<G%H0rrSH(N_A0PWLt2I5h*qz?;J%V1
zR^)4kyQLJ%GjnBEUf!-au=c^fKR%rJuJ}&y;%FmH(tvqS`9^Vq^f-!eOJ`4`XQ~<%
zZR>D~9k*!bJZUpEkBp-$b1&w)uJ5z!K+tQx_qZqU&YLAwUDMUg9I{^Xp7o8Z30~%R
zlEE=A)bZtF`k{QFe&~&#iXRspO2`H(#^hIOy5(<1P;s`e?Oqck_mG8<$J5SnQ~flc
zA{`>}QVEMRnAFaLOM`<z<nOjS`HNfPP)~W7mA|oPHK7Vsv^mK9v_m+?g(5vI=lsb?
zsx&;bO<&?&iwQCPa!-kZDUwp_TW7=ay}sGhdm3;Zn(RKJ&eSZV#z{8iDabFLTb}2j
z1<aTCHY=-T`Oa%al(lZ921!m|mwbM%nZa7k^^f<bDD)N7sr7Hvt48IuL=4*z>*U$)
zrTIBbS2_#d7Js5D8zV=>Go;8{^UXk}x;`iW%ISW}@nal$56Z*HEXvZuV#?dq53h&#
z%S4>6QckLrPORbE)f6Sgo=80U@YURi{V2(*a#})%P3S`wsH!>m4%ZR+oEF%J-k`pj
zeXiA{&A;8uuYEp|cdkg}QpG!pcJQ2)w)Qn-^wo;p$hG#jyQcK5cbzkvf{mW+o+>On
zCe*cZrw<zt7rPdH=Si-RD4mJItex;(;wud70po9cVpAy@M0>^vLXKaCRB`%_>y_LU
zan@<Q;%U5Qz}Y_1^KgOr?eSDw^NV`=RT|ou2iZlvnkwJ(O)m|s+UFaJ_o^*ARIghL
zGc;yPygEc$$37bWlHq+EKzjJx&8ytCM0R7U$iR24_w=n@jG~JA7{7!`&Jsm#K;90I
zg6Dl^Ne=ePoZKqTQboj1gRe37?y-@WNrBo=)YL>6oIMP)AD1ZJ{Th=f1^UX+*qlC8
zre^p&gJzrg+f(xZuj--nBca*Q0Y)8yS=x08mG~Om^oaB}s%y*KXNU>H>vnt}O)aBa
zJ*0ivaYwpObHWFcU6-9l`cG8ujc}eJ@Ho{bi{j+&Ek1kX#^F$>TJ~I0gar5UcLO>f
zg%wp{FNE5&<qBxlPM;CfWUNcya+OVJI-s((Se30Y1Y352R4eeNN2tvcnN9IaTuf=2
zo3&#ez03t4;>~GbK#(YoFplKN;1mNTMQirUgU}CeRG=n%)9j4gzLD2)MjAYF=C#eV
zX>g=u2~i^)<fzfDmsy(R9fvx?C5(^jBmM4GD67k`bUb{M7kRm`Bzjlo9(~siNoqS*
z!wK`&jcJ`5dtE^;>{}w80!>2X8@Ng7M!2!k4-5sWJG+wyp01>E-V{s3bMjH@wJ^jk
zD~+BW7wNwtIeqpfRkYF)>PSMpkZeOzRGT<h;vK3iK1Zx7yqZGHQX04xE}{pEAnUI-
zmD5FbCG~G{W%@AnN|kt)iW=Y96WIvjeJU+ez7zGBONBLHfq^E##r&|Rbqq$q_)<b`
z8CE5rPi|zGN*YSkFZY#cu#MeHc<#YYofuCt5kx92dF=vG(yr*CPildp9Mbqe_)!pz
zM|kW?hHqK#E$%p}Y5xA4*11+$aUtsR#G%{)!Kxb#J2&ICx%&@VXG6Bb#T{)JoBKub
zi?QW2zkjWfTk5*RL6D?9$DSt8mu(YwF#3L9JFH7Vt533AbTZ`<l{kE9U7qDYdKa!4
zP$23ga*K{qjBsiGnb6pdfm2Xwp$oT{dL!<SrsKUnxuV3T4_CGg`FBn(7`u_{m&A_k
zgO`c6the^WMW>>Du%R~-ahFkt#Kk%(8bb{mTFj%kpW=EJh_?1pr&@bjE#6UsR#ulp
zE29#U^Hpf5+tLb6REX7@xC=B3bJtiNbn`FheTw0lTTEffRB?K8%NrzY9^z)CYEh|n
z@G-yayUf(oVND>(qoID-NTce`Y59(`e@<**uiBOLWmkP|rG<x^5V~ws`RyG2y`q+p
z{x!-WFQvfZy;0_o&>Gz->$Q^SGs}o$pF6?txtgwS^BmuPd2I%9=I_SFD4`-JL&X{K
z__h()sZ{W*e}ROP^wY4Qwi`<obtO$D6FVi#$o`JB?YCVnxV!tcr<>lB<CFGB85TNI
zku`xed9BXjKFH4-uh^8YI!(D3T-dI@OjaKD6hv(k)oeAY%zZT`x{Itp-wq$<$UN4Z
z06ILqIEKq{5+5)8etq0*=KzmCE@)Q%dfsfnCra<aBS&(9w%V}HP(>?oU%5P(U{)Aa
zs$QSCkzb00Tf}6@WK>P)%ZEO#g+*fBq1aY33G_XYQ605W0-*wsE=aZhp^UGhxxx?|
zYF-mD=nt2MOGAgFd7|xOr%DDI?-TWLI^7_1SVkG7l**-+_-)KK&M7=cX?XMXq&fEK
z5=zWV7Y7{M4K)=@m9r!t7}~+(qp-jHWOI3N_BJSzsrVlEYqHG`y?7Cq4pd|qw**&0
z_E4tz$qQ;E)M{NltlTSgl$jI;DFPD0bBQ}t@u5$l#m6is!hw2iIo<+{QeLE=?gqgJ
zL#3U3G=gJB*Gv4K-E4YkQGsjc0(T}`6@s&AO+CIr^E&?xDd`isFI+~q&!ksg(OINS
z6*_Ja*w7K^WG~%ux~TzfwFwK&zFXAjQZwE>t3v+bfYLwjogTC;{Gol(eWO^bz*>=x
z<##j^-^F{9_gePx`*aUJ8QnMWz#V_Onu#omREja4=NTcmTWG=}#{YWIEJA_)=z4yA
z76jkqtGUD8tsHE<K3z*WEo7J$Ib8X6*==5TTcJ0XyS#(WE@qFc4x}n<gC?}BKsomv
z<hNC)Y+@)KL(4?X6Jr_n`n7U!sTs0mX8U(;n<b@>R@RIeUd;=&s4r)ZZr;{O`w>$A
zDWJB=V3Q-NhpQyMQEo}zlN<!qk#f$II(^WU*w@Mza^`MJ^ck-zVqv<@i7fe;hC?!(
z@Jm5WwxhG|r;<yvd{X@k(jvU<kueN@hIs4A%$-ceAY5?EeP`_1>Hw;PomNO@p#}d=
z!By2SiMjL2ZX12=CRrUNeT6#WqCxj=slYy2D%XFxeempBH~s=46@orEY`tg8kuREZ
zW2a;wd@kj)=RP88#b(58XoSK~1>APiApF@fo#yvJ?4Ol1o6>Oz<&+ZwE!eTAY`$|n
zm0N}X2I)Liux)e-PTT64l$@Wqan#-D6^;*~^^TL`wa)}pXZGHVb{dr?e$U<X{^S%S
zD9ydpg=2=hzWu?y9Fm&ZkX`I#k;?y4L|#qLwP~7n$UJ~!FN&v`useBIne8zf3*TsM
zN<3qb9!`B2n^>uJ(FphvUzetAX8NHl$>3vT1i|)BX1_5TTc2yqmucSkh@PK4y-(u8
zJ>Fc2pzc}{uQ2D1`O)1%DVAPUz3zznb>nI_ES=vlcJci%C|fCtN=7Vo121U;kMCfp
zm{ovD7TK*!{W2yNSas7n_wZ)Q$_O-0LX1iUSYtPc;hsVn*m{oL`-f#e3=eO0Dnuas
zR$bIHlnhqC+$D@eHfJQw`XxA7uaG<Zq&P#k$&YzBNlK~4KJ9;Yc)zc&WTTZ*JbR~1
zarn^H$*59bKxW2VI<ah-?SVYpVr3rrdi$&5Gljx+LIVL*@maW`XJ4`-j-C&GmhmBD
z7o=D0ECO+LWwJtJCuJxI96*xs;%;ddt?$)3Vr`bxo<{%9Vc^`t*nq4@&~3;>c=BPj
zndVLQx~c8a`M5VP14ENvK4Fj_X=9Y{S7{VcCBpeR{v5@yVAaWbh!}cIL>E`i(0Ovo
zcyFWUTi-+Z4`Xk>WVmEZEfN#E*efR^A1~$h8z&oabv^k`&p~HkqQPZF_`LH;&JAQP
zxiu$g{$68f)o0hpo3;|ueLa~>P0zRk`c(SPLMJG?@?&Q58+u37MwO;4N)-pAQj-pZ
zK766xI%MUm`J^PiIDPePx`e3tj;D+4rIvZsSkA@ibP7>EF4MBQdjy2rfy0A`zOpoL
ze;ypVE}SW`&>*?QDSLyO4!u7F;qBVPzM8r2=q-0t8kmG-Wfjg@#;uKg)#mDWohjn%
zsONUcuSH~Pkz+U$Cn5SZT^`My%btC<1owoW!aucYdG{TKMO{tXBgdP$F@!`|@D6pB
zmHl}6eZTCm?2ToLZQxY);aKHC$MJj_V}$a0(i`HywTcJobZ*jDs7$nzr*FkImmZBa
z8?94~A$0W3cr;!(@F)%2SvGg9ws?Im+1y(B`m&UQm_d@x+J}p&fTFNv+rTK|LWP&{
zg>)>&EnL{8AByfjKXys!<=?EyOnkT_;SoSmcF~B-E}R9|_tbm#_@rN~nUc!F@a#R^
z_L=hQ%TxT?XD#DrZjhX2nVT<eX4JQ7^HtF|J^s0cl2@+KMEdcX)_!ZyZV)4Pk!E<`
zu?XymV!00M*=Y9#p}l)^5kx+Y&j)z|hGU;j$M#cj*trwGC;rw+rbtkXiZ7W!3FcIE
zj?oT>q=;)g2~Du)bCDL8Imyb~%f9t9`_?@URu0Pho`Z!K>`BQ2Sg6w8-MFdP{Pb-1
zXB%`cR~@P4FVO1Y+P_<9@osio9(b}tRBASQB=<hO$64hlAjXWIE@SJRa)j)2|0TgS
zJKx(I4SNZzG@F~8it}*#nq{l|1m8@fCd{wDH|icy)|j<%pc((+^4yKCsq0|skn%h8
zvl;gb4Fz}O)D<SHA}_yI(ZLz7w43T|YdK6SgsxS)PN;p=uleEnJ+tSiW~q33MUMF3
zv0XIbLkY^9w0QhyecJsHJ$a$Di3A$rKo`@U$JbS}mgM~XOp}7bb&IjaE(LA%#wR}%
z%PdJRF<1)Tv74_Qw!wn4!hv{_6tRLXd8w+oiO*B;2s623>pbi;+v!%<<UVC)QMa#5
zYpl{0k}_j`yYTeWwfH-1)cc~VqKhbR)@D<wnTQ`@H}l0gr1JycM<Cd2oeqVzg;vM+
zMSSt)SNorB#YT`0sl`)=<94bIG_t1eKbcSYN!MK>fo<X2>vGy~8|<$=Kyj;hbk_Jb
z`>Po2=Ujy0aoh$qPUR?0L1hVb>+3i=*Ww4^v&Q^#yw|SdHD$lJ^M>>;V-rf8#BIp^
z;ini_j+Unc!G3Y4cvk~)P|lJvw|Tv<*wpib><t5@T&IQT!R)-OS>|`0&Ab9t#avcd
z`T5_c193HUiS|x5#vRA)$zAx}stY#zEq%U&?F&ofx5Z;76hF#nH1hn&B|C78ZVJ}^
zCg<-vig<9-WQSkL-}ha?(@zh%O$k4z$8@}tDF=Tb)Dili$Dt@W+=K)#Fn8(qt#u3a
zi%4AqgSjSeSjuPF60xCHSr{EkoygkVsejgv2`1@USzL{XWFm(?E3{<E!aqFb6af>_
zVda0DvX6+*IGs2;=$IR8&fEF~nrSr!m)#ua>!eU)axR^U)PA1D$H=ddS?GO-yz)8c
z>;30`vWha%SJNa>No_*(Q29CiO4`;o&sP_{4>Y=(b;81a@C`Lp49~u0>vr4Ds1swY
z6|*T*m=^wcJ!P-zmEjjM@|FOcQi^nZg_3w^44gAzSeNNaLPiy~BL1^NRYwA!9+<G7
zz;Mjwjot*qE<;Ib@K>1XsL9v-HF&Kta*bQl2!Aw-7T;`}ro4HA%I(>B%}1hseo2A%
z`ePQZ;Dq$qZSglTnfKg2<I_U@r#uii#SxK|*{US-J`XtR*Pj-D?7!R0VPxKxkoo3d
zHeUYuH9Ac>JKr0@9|)Oe%SIySh`*E-*p@}*H<_mwg~W7~W|MGfkZr@CFE~4J^B0;L
z2T>herf@EJR81e*G8x*JFj-lghX5Bpq^5axiMBl3JGhRQnv<Sa>+WD-W$Fi+r%7Bg
zo|JT8ZHL>EnqVS|+s6!GtVs9wrl%i#m=}*xV-@M;lXfk_IS)=_YvcEgzsYk6aI;yr
z1Qvx>;40A_ZXSBozRch_9Z`F~k>P4LHA_wXra|sao!suUSt8}$gDPKxZ=s^idK^rh
zYWJ%lVLx)`MVdMsYChlnQl(sr>yo3=R3t&xc-$DmRe6o*ir6F%ndZl4qUk!^h7f+8
zHU`;by*;IwJ=vqY*PIus4OUMxBQqL0KSX{IE?_3uV?2D$Tf-%6zAC{pYQMI>aaO0k
zoYmS=q;|_M*8B!t&V18|j+OknQN`=Eii)P)8@?jT$yWKDxy8ftoCQ5SxpaN{ngWJu
zeXosrMu>4m2_!QqwQ`!YLtG#~M_*OR?-BBbJDE;X`+mPxqHq`dvHo%wV<jy)b3q%^
zrsa?N3{jZ`U1@O;S9yUlXLprc!aMr_Mv_mH5P^>kW%YI=ySXQ0iP;%n>UH4)RiKKx
zs9=cmmv4eD!4IZnS4QSWHhd70W3&%w7imHpM{o;_T%BI*RVFgxS=dnUFUovAoU#0v
zItOyIvv0RRMd@d3c|<ZD3PxaS<Bx~h%Dn#i@NMDJ+<wu_15rJXFF78u9DYv(4|8VJ
z={enHE-jjY-4gxZ1P692*}c=uM&cNf$Y4*_1=n@EXHwuRMK#u&T%2ibo4^4pY{HzG
z$Uc9=ZaNOpR=U<wofzRu6=p40v{k&v)<Uu&D^8k+@bqG7)<{qBIc4FVS-Hpt$BUyK
z8Tl0c;SN`fXhB}ZqL8#q6W$e_%+-qPrSG%~Qh4@fytT^_eu|0Yr}FG&(mt>~*qv9J
z1MjDWZ9RHoZpc~(dQ98_eUkYYe5J8)l-e8?>9>Aa?WMMmT#*1ht3EESFCUJOT-&FO
z%MbRX45JY)Gnv#Q1?>SQO8QHzJhZhYv&Hg*@8)t`H|CJuWIgvjzIfTM^U7ay$luL;
zV$dZ2spUw2k+rr36))33T%U`5#E|76X?fh)%ZYa(Yw@qlZKvmseSZ!=rL=~!JMpo9
z#Xciu_EB0F)GtC7OEpSGFv(|%?O6sfckp%cbyt5b9!>Ut%k|{mlh|+=Yi_rkLB8z}
zqk2?XT~X`g+$YxMv=Z!ua?|us>C_s!GnFlaOgi1TACGu5HRNlc0w|x%fy`(oWO_^<
zPkj5n-09QOOkYBkneYvzPIm6z^@22Ax^l>1mda{wuF`1o)zv~OZeuC%lI0v*V`J{O
zI!1lMM>{eXzxmIscqGsJ9h7kD<e_va=-Qe$cqH3rjxXrMj^2?pSd3mnr8<lSyrzl7
zEev+w7wovc&~!~61-km>MT5dS)V<(+oFs7t*b~H+*L+yaK1|BS3!V&x-Wo-uC{hUt
zrtIqJI7MeKqz`Ty!?h~{Q3d2ER)jr!s=jGX<8x&tsAG~i!u0Xgwub9TY@F1`3d4Q7
z$Oo~dHC!5B%Y(^>j*XkPrP(!qd@_>$9u=Er3haRG^S;kDns%JVR0Y-89Y-W9tEY?d
zM>-W~q!vb>jA?(+2?#^xN;XCXGa|e@UYKVin>c0z9Fs;FL@l^QEt(s!snU^JjSWFH
z2X9|ImLnP@y;7%?E#*@P4yxj<;dK-}_Mbc4Kb*roc-_Z+G`liu#O=oix5IkR#R1|d
z=O`@y`o3iS?I>+khrUEnZn<k6qu7w?qV^(H`kbz^c1eV=B_vk-=hBeAlV|sJJ5j6q
zFET$4i1Ab=tQTC-dr;H|t}!+F(rCG?BS}nF57V(@qOuOv?`(_M_`Ig0zW|PksTjTO
z*;y#mYRaC!@L}1@`*07tQ+(iHHF8%|?1$0w@mS|bjVrtC@Mb~-d(y>hTZ13O{X_na
z_=UPNV~84R<|bY1hX}g?uhh_8e`~I|l_RHRX8U!Av0|ugS7$&E|I0;>qmjxv@}R<N
zWJkQSCieU>HI}{g#LSvC!g29rPu>@R(;D|rOC<5*KyF$K>kF^-PrbZOzHViBl`C%)
zY*I-ids*WbD>D}H8W0v|?&Gnzb|8!8u9|nrU5yq~#kpRS4vjQ+c1T4<+Jc(!3xrK=
zsz){kZk3jcneN$8=%{@4_lEPd5Ya`xbFaF&yO*Tswsk4Ls}XY7ZM&R#C@6v3P^fQ>
zQkxEz(0ak8YQIA4cuZR8%#0+6r{8zKSxGPgR&!<s;gqFFbs%wZOKQ(f3-?OH?TGBQ
zP=-|;)?TpU4Eir;e}5Y<(H2b}2;!g6(*NNSWdW@B^wRS5%79gq6Uhmum$gM%**MY5
zJQ5WKOu2dvHmv_R@BZ!M%kt<x1O(;=0mfZ04-^arse(Yv!0$gv0o(51R#X`)3td|s
zJuAT6`WL^Vu$3Ol25?w-gg&4E{+q&7%0l-X@Xs*`2P-{&zzhxsBk%xA@?S4{Za5qY
zr`M<djln=5H1#|BMQ{EW27(~aw)}G#6bwfbNu0xg3&0xw2L`w&{FxU7L82pdoy!YB
z@SqR5oWnpc2oyML{09a@2c$cX!GT6&U=XxG&ZmRH(Px+cz<7W#FMnddH3o*Zrk_s-
z963dY2t>dBBP+nBjV4<^hk+1aI7XWhP;?lN^XcF)bclpMFa+j(z#&kKK7zx55FY3A
z0)2#WZ$Lm4G%xzOIsiv>3^@bN<Y>Q=bLl{E1Qerv++Z~C%lULr9yFcd9~csS!2BEr
z$O=u)@dt)P`<9%;Ks+Gy(fdCzDEg$%c?=3gldS)N@xU<BA)%OhfP|wxa?a;PK+$Kc
z{=m2)808_jF=P)yLNM+Z#0^3FJDo3&8-_V%xM3K51jqo*4t72-j0b(d`40?<j!=0X
zL&DJnyyq}LcQN~j2MWRHbDsaA&%ltukO2=8&24<H4iFC#^WG3(=$~zd@W9cOvFGvv
z7@E8J90rC%kbl+z=0>0oN1RUw46Q$91%@Ip+6;jpG3F*<#zFj92LuYt(tp$efg&)*
z2QX$Z^c0wEG5QFQ0LJ(QL%Gp>@aOLX4nk9={eeNy%sl5X2po>lZxA>S^iN#_+JcUW
zbuKT^7PMQ=ISj~)F=ik<5Oe^g^XXvR7<~?bAu#855a2QNXI>BjQ_kq6_0M!XK(LrU
zFfbB{u~q;*gfZU0!1(=BuYkPhqlM@00}R6)(_mnk!H_)+fhMy*pBD~-W9SP6hB4m*
zx`r{=0Xl>x1UpwAkQGBNU>FQj$6?$Ua|s}O$bZodG~wL2I)E{QF^2)N!t7rV=9q>+
zV3>0GZ({)lfPd5r;o-(ur-5}JW2{0D2=o@`T$_Ql4DHeQ2gZ%2$UTpt2Q>zUxz0e5
zNR06YU>JP_?1%onF4$P10GBr_JhXqRypa=n+l8R#lC!V?d{BS+naZJgp!F>PkJMk}
z=1(30KB>IgJlrTAeHe;I8wAq<gaFhHg&|QW5CX{q(FY;<@&5NDzr75tZGc_nuYDX0
PXcr8Rkx@ii6z_ik^|D@p

literal 0
HcmV?d00001

diff --git a/book/tccs/shapestacks/SConstruct b/book/tccs/shapestacks/SConstruct
new file mode 100644
index 000000000..ef718ae21
--- /dev/null
+++ b/book/tccs/shapestacks/SConstruct
@@ -0,0 +1,8 @@
+from rsf.tex import *
+
+End(lclass='geophysics',options='reproduce',use='hyperref,listings,moreverb',hires='cmp2 nsnmo sstack wsstack', color='ALL')
+
+#lclass='segabs'
+#Paper('paper',lclass='article',use='pdfpages')
+
+#End(use='hyperref,listings,moreverb', color='ALL')
diff --git a/book/tccs/shapestacks/elf/SConstruct b/book/tccs/shapestacks/elf/SConstruct
new file mode 100644
index 000000000..b63059147
--- /dev/null
+++ b/book/tccs/shapestacks/elf/SConstruct
@@ -0,0 +1,422 @@
+from rsf.proj import *
+#import random, string, math
+#from rsf.recipes.beg import server as private
+from rsf.recipes.beg import server 
+
+#random.seed(2005)
+
+nr = 0
+
+def rnd(x):
+    global nr
+    r = str(random.randint(1,nr))
+    return r
+
+Fetch('elf0.H','elf',server)
+
+Flow('elf','elf0.H',
+     '''
+     dd form=native | cut n3=1 n2=1 n1=300 f3=663 f2=67 |
+     put unit1=s unit2=m unit3=m | bandpass fhi=60 | costaper nw3=50
+     ''')
+
+
+n1=800
+n2=128          # data dimensions
+o2=50
+d2=12.5         # lateral scale
+rect1=10
+rect2=10        # smoothing for dip estimation
+p0=0
+pmin=0          # initial and minimum dips
+clip=5          # clip percentile
+eps=0.1         # regularization
+nsp=200         # number of spikes
+ngath=800 	# gather number
+
+nd=1597
+dt=0.002
+
+# Data   
+Flow('gaths','elf','window n2=128')
+Flow('igaths','gaths','window j1=2')
+
+# One CMP gather
+Flow('gath','gaths','window n3=1 f3=%g' % ngath)
+Flow('igath','gath','window j1=2')
+
+# velocity scan for 4ms data
+
+Flow('vscand','gaths','vscan nv=120 dv=25 v0=1400 semblance=y',split=[3,1000])
+
+Flow('pick','vscand',
+     '''
+     mutter x0=2000 v0=1200 half=n |
+     scale axis=2 | pick rect1=70 rect2=70 | transp plane=23
+     ''')
+
+# NMO and stack
+
+Flow('nmo','gaths pick','nmo velocity=${SOURCES[1]}',split=[3,1000])
+Flow('stack1','nmo','stack')
+
+
+# Dense stack for one gather
+
+Flow('istack1','stack1','window f2=%g n2=1 | spline n1=%g d1=%g o1=0' % (ngath,nd,dt)) 
+
+# velocity scan for 8ms data
+Flow('vscan8','igaths','vscan nv=120 dv=25 v0=1400 semblance=y',split=[3,1000])
+
+Flow('picks8','vscan8',
+     '''
+     mutter x0=2000 v0=1200 half=n |
+     scale axis=2 | pick rect1=70 rect2=70 | transp plane=23
+     ''')
+
+#Result('pick0','pick',
+#       '''
+#       window max1=2.2 | scale dscale=0.001 | 
+#       put unit=km/s d2=0.0133333 label2=Midpoint unit2=km |
+#       grey color=j title="Stacking Velocity" scalebar=y barreverse=y
+#       allpos=y bias=1.8 clip=1.7
+#       ''')
+
+Flow('nmo8','igaths picks8','nmo velocity=${SOURCES[1]}',split=[3,1000])
+
+Result('nmo8',
+       '''
+       byte gainpanel=all |
+       grey3 title="Conventional NMO" frame1=400 frame2=50 frame3=500 point2=0.25
+       ''')
+
+Flow('istack8','nmo8','stack')
+Flow('istackn','istack8','window f2=%g n2=1 | spline n1=%g d1=%g o1=0' % (ngath,nd,dt)) 
+
+#####################################################################################
+############################## USE SHAPING ##########################################
+
+
+############################## Find dip of elf #################################
+# Define dip
+Flow('ipick','pick','spline n1=%g o1=0 d1=%g' % (nd,dt))
+
+# Predict dip from velocity
+#Flow('vdip','pick','remap1 n1=800 d1=0.004 o1=0 | v2d n=128 d=12.5 o=50 mute=y half=y v0=1880')
+#Result('vdip','grey unit1=s unit2=m  color=j title=Slope allpos=y scalebar=y')
+
+
+### Dip Method #1
+
+# t as a function of t0 and x
+maxvel = 3487 #3387
+minvel = 2500 #1859
+
+Flow('t','ipick',
+     '''
+     spray axis=2 n=128 d=12.5 o=50 label=Offset unit=m |
+     math output="sqrt(x1*x1+4*x2*x2*(1/(input*input)-%g))"
+     ''' % (1/(minvel*minvel)))
+
+Flow('td','t','window n4=1 f4=%g' % ngath)
+Plot('td','window j1=50 | window f1=1 | transp | graph wanttitle=n wantaxis=n yreverse=y min2=0 max2=3 pad=n')
+
+# dip as a function of t0 and x
+
+dt1=0.004
+dx=12.5
+
+# NMO with maximum velocity
+
+Flow('vmax','ipick','math output=%g' % minvel)
+
+Flow('nmo0',['gaths','vmax'],'nmo velocity=${SOURCES[1]} half=y | spline n1=%g d1=%g o1=0' % (nd,dt))
+
+#Flow('inmod0',['igaths','vmax'],'nmo velocity=${SOURCES[1]} half=y')
+#Flow('nmod0','gath vmax','nmo velocity=${SOURCES[1]} half=y | spline n1=%g d1=%g o1=0' % (nd,dt))
+
+Flow('nmod0','nmo0','window n3=1 f3=%g' % ngath)
+
+Plot('nmod0','grey title="NMO with max. velocity" ')
+Result('nmod0','nmod0 td','Overlay')
+
+Flow('p0','ipick t',
+     '''
+     spray axis=2 n=128 d=12.5 o=50 label=Offset unit=m |
+     math t=${SOURCES[1]} output="%g*4*x2/(t+0.001)*(1/(input*input)-%g)"
+     ''' % (dx/dt,1/(minvel*minvel)))
+
+# dip as a function of t
+
+Flow('vdip2','p0 t','iwarp warp=${SOURCES[1]} eps=1')
+Flow('taper','vdip2','math output=1 | mutter v0=1500 half=y | smooth rect1=10 rect2=10 rect3=20')
+Flow('vdip1','vdip2 taper','mul ${SOURCES[1]}')
+
+# One gather
+Flow('vdipd2','vdip2','window f4=%g n4=1' % ngath)
+Flow('taperd','vdipd2','math output=1 | mutter v0=1500 half=y | smooth rect1=10 rect2=10 rect3=20')
+Flow('vdipd1','vdipd2 taperd','mul ${SOURCES[1]}')
+
+Result('vdipd1','grey unit1=s unit2=m  color=j title=Slope scalebar=y')
+
+
+data4 = 'nmo0'  
+idip = data4+'idip2'
+
+# Test one dip
+#Flow(idip,[data4,'vdipd1'],'dip n4=0 rect1=10 rect2=5 rect3=20 idip=${SOURCES[1]} order=2')
+#Plot(idip,'grey unit1=s unit2=m  color=j title=Slope scalebar=y')
+#Result(idip,['nmod0',idip],'SideBySideAniso') 
+
+# All dips
+Flow(idip,[data4,'vdip1'],'dip n4=0 rect1=10 rect2=5 rect3=20 idip=${SOURCES[1]} order=2')
+
+Flow('ndip',idip,'window f3=%g n3=1' % ngath)
+Plot('ndip','grey unit1=s unit2=m  color=j title=Slope scalebar=y')
+Result('ndip','ndip nmod0','SideBySideAniso')
+
+# Seislet NMO
+def seislet(gather,dip,snmo):
+    nmos = []
+    for i2 in range(n2):
+    	traced = 'trace%d' % i2
+    	if i2 == 0:
+    		Flow(traced,gather,'cut f2=1')
+    	elif i2 == n2-1:
+    		Flow(traced,gather,'cut n2=%d' % i2)
+    	else:
+    		Flow(traced,gather,'cut n2=%d | cut f2=%d' % (i2,i2+1))
+    	
+    	nmo = 'nmodd%d' % i2
+    	nmos.append(nmo)
+    	Flow(nmo,[traced,dip],'''
+    	seislet dip=${SOURCES[1]} eps=0.1 adj=y inv=y unit=y type=haar order=2 
+    	| window n2=1 ''')
+    	Flow(snmo,nmos,''' cat axis=2 ${SOURCES[1:%d]} ''' % len(nmos))
+    	
+# Seislet Stack
+def seislet_stack(gather,dip,stack):
+    Flow(stack,[gather,dip],
+         '''
+         seislet dip=${SOURCES[1]} eps=0.1 adj=y inv=y type=haar order=2 |
+         window n2=1        
+         ''')
+
+#################################################################################
+############################### old dip #########################################
+
+#Flow('pat2','gaths','patch w=800,128,150')
+#Flow('dips2','pat2','dip n4=0 rect1=10 rect2=10 rect3=10 p0=0 pmin=0',split=[6,11])
+#Flow('dip4ms','dips2','patch inv=y weight=y')
+#Flow('idip4ms','dip4ms','spline o1=0 n1=1600 d1=0.002')
+#Flow('ndip4ms','idip4ms','window n3=1 f3=%g' % ngath)
+
+#################################################################################
+########################### Conventional stack ##################################
+Fetch('elf-stk.rsf','masha')
+
+Flow('stack','elf-stk.rsf',
+     '''
+     dd form=native | put d2=0.0133333 
+     ''')
+Plot('stack',
+	'''
+	grey label2=Midpoint unit2=km label1=Time 
+	unit1=s title="Dense"
+	''')
+Plot('istack8',
+	 '''     
+     window j1=2 |  put d2=0.0133333 | grey label2=Midpoint unit2=km label1=Time unit1=s
+     title="Conventional Stack"
+     ''') 
+Result('istack8',
+	 '''     
+     window j1=2 |  put d2=0.0133333 | grey label2=Midpoint unit2=km label1=Time unit1=s
+     title="Conventional Stack" screenratio=0.6
+     ''') 
+#Flow('nstack','stack','window n2=1 f2=%g | spline o1=0 n1=%g d1=%g' % (ngath, nd,dt))
+
+
+#################################################################################
+###################### Apply SShaping to one CMP gather #########################
+
+# Backward operator
+def backward(data,model,dip):
+    Flow(model,[data,dip,'vmax'],
+         '''
+	 spline pattern=${SOURCES[1]} |
+	 nmo velocity=${SOURCES[2]} half=y |
+         seislet dip=${SOURCES[1]} eps=0.1 adj=y inv=y type=linear order=1 |
+	 window n2=1
+	 ''')  
+ #	 pwstack dip=${SOURCES[1]} mode=1 order=2 eps=0.1 
+niter=10
+
+# Shaping NMO stack
+def shaping(cmp2,mod,plots,dip):
+    '''Put everything in a function'''
+
+    global niter
+    mod0 = mod+'0'
+
+    backward(cmp2,mod0,dip) 
+
+    m0 = mod0
+    old = m0
+    plots = []
+    for i in range(1,niter+1):
+    	new = '%s%d' % (mod,i)
+    	specn = 'spect%s%d' % (mod,i)
+    	Flow(new, [dip, old, m0,'vmax'],''' 
+    	pwpaint seed=${SOURCES[1]} eps=0.01 order=1   |	
+    	inmo velocity=${SOURCES[3]} half=y  	  |     	     
+    	window j1=4			 	  |	
+    	spline pattern=${SOURCES[1]}	 	  |
+    	nmo velocity=${SOURCES[3]} half=y   	  |
+    	seislet dip=${SOURCES[0]} eps=0.1 adj=y inv=y type=linear order=1 |
+    	window n2=1 |
+    	add ${SOURCES[1]} ${SOURCES[2]} scale=-1,1,1 |
+    	bandpass fhi=100 flo=2 ''')
+    	old = new
+    	
+    	if plots != None:
+    		plots.append(specn)
+    		
+    Flow(mod,new,'cp')
+
+#shaping('igath','shstacked',None,'ndip')
+
+Flow('gmres','igath ndip ipick','shpwstack flo=3 fhi=80 niter=3 half=y velocity=%g jump=4 dip=${SOURCES[1]} mode=1 order=2 eps=0.01 nmo=y' % minvel) 
+#Flow('gmres1','igath ndip','shpwstack flo=2 fhi=85 niter=3 half=y velocity=%g jump=4 dip=${SOURCES[1]} mode=1 order=2 eps=0.01 nmo=y' % minvel) 
+Flow('ugmres','gmres ipick','spray axis=2 n=1 d=12.5 o=50 | nmo velocity=${SOURCES[1]} | window n2=1') 
+
+# Spectral comparisons
+Plot('istack2','istack1','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" title="Dense vs. PWC Stack"')
+Plot('stackn2','istackn','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" title="Conventional vs. PWC Stack"')
+Plot('ugmres','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" plotcol=5 dash=1 title="" wantaxis2=n wantaxis1=n')
+
+Result('compspec','stackn2 ugmres','Overlay')
+
+Result('compspec2','istack2 ugmres','Overlay')
+
+#Result('compspec2','istack1 ugmres',
+#		'''
+#		cat axis=2 ${SOURCES[1]} | spectra | 
+#		window max1=125 |  graph label2=Amplitude max2=11000 
+#		title="PWC Stack vs Dense" max1=100
+#		''')
+#Result('compspec','istackn ugmres',
+#		'''
+#		cat axis=2 ${SOURCES[1]} | spectra | 
+#		window max1=125 | graph label2=Amplitude max2=11000
+#		title="PWC Stack vs Conventional" max1=100
+#		''')
+
+
+###################### Apply to all CMPs #############################
+
+#shaping('igaths','shstackV1',None,idip) 
+Flow('ngmres2',['igaths', idip],'shpwstack flo=3 fhi=110 niter=3 half=y velocity=%g jump=4 dip=${SOURCES[1]} mode=2 order=2 eps=0.01 nmo=y' % minvel)
+
+Flow('ungmres2','ngmres2 ipick','spray axis=2 n=1 d=12.5 o=50 | nmo velocity=${SOURCES[1]} | window n2=1') 
+
+# Apply lateral smoothing
+Flow('ungmres3','ungmres2','put d2=0.0133333 | transp | bandpass fhi=25 | transp')
+ 
+Plot('ungmres3', 
+	''' 
+        grey title="PWC Stack" 
+        unit1=s label2=Midpoint unit2=km 
+       ''')
+
+Result('ungmres3', 
+	'''
+        grey title="PWC Stack" 
+        unit1=s label2=Midpoint unit2=km 
+       ''')
+Plot('ngmres2', 
+	'''
+        grey title="PWC Stack 2" 
+        unit1=s label2=Midpoint unit2=km 
+       ''')
+
+Result('stacks','istack8 ungmres3','SideBySideAniso')
+
+## Window #1 
+Flow('wstack', 'istack8', 'put d2=0.0133333 | window n1=150 f1=50 n2=300 f2=600') 
+Plot('wstack', 'grey title="Conventional" unit1=s label2=Midpoint unit2=km')
+Result('wstack', 'grey title="Conventional" unit1=s label2=Midpoint unit2=km screenratio=2')
+Result('cwstack', 'wstack','grey title="Conventional" unit1=s label2=Midpoint unit2=km')
+
+Flow('wnstack', 'ngmres2', 'window n1=600 f1=197 n2=300 f2=600')  
+Plot('wnstack', 'grey title="PWC Stack" unit1=s label2=Midpoint unit2=km')
+Result('wnstack', 'grey title="PWC Stack" unit1=s label2=Midpoint unit2=km screenratio=2')
+Result('wpwstack','wnstack', 'grey title="PWC Stack" unit1=s label2=Midpoint unit2=km ')
+
+Result('wsstack','wstack wnstack','SideBySideAniso')
+
+## Window #2
+Flow('wstack1', 'istack8', 'put d2=0.0133333 | window n1=150 f1=25 n2=300 f2=300') 
+Plot('wstack1', 'grey title="Conventional" unit1=s label2=Midpoint unit2=km')
+Result('wstack1', 'grey title="Conventional" unit1=s label2=Midpoint unit2=km screenratio=2')
+Result('cwstack1', 'wstack1','grey title="Conventional" unit1=s label2=Midpoint unit2=km')
+
+Flow('wnstack1', 'ungmres3', ' window n1=600 f1=97 n2=300 f2=300')  
+Plot('wnstack1', 'grey title="PWC Stack" unit1=s label2=Midpoint unit2=km')
+Result('wnstack1', 'grey title="PWC Stack" unit1=s label2=Midpoint unit2=km screenratio=2')
+Result('wpwstack1', 'wnstack1','grey title="PWC Stack" unit1=s label2=Midpoint unit2=km')
+
+Result('wsstack1','wstack1 wnstack1','SideBySideAniso')
+
+
+#################################################################################
+############################## Shaping NMO Stack ################################
+
+## one CMP gather 
+
+Flow('shgmres','igath ipick','shstack niter=5 flo=2 fhi=100 jump=4 velocity=${SOURCES[1]}')
+
+Flow('shgmres2','igath ipick','shstack niter=5 flo=2 fhi=124 jump=4 velocity=${SOURCES[1]}')
+
+Result('mod','istack1 shgmres ugmres istackn',
+		'''
+		cat axis=2 ${SOURCES[1:4]} | window f1=300 n1=500 | scale axis=1 | 
+		dots labels=Dense:"SNMO Stack":"PWC Stack":Conventional
+		gaineach=n 
+		''')
+
+# Spectral comparisons
+Plot('istack1','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" title="Dense vs. Shaping NMO Stack"')
+Plot('istackn','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" title="Conventional vs. Shaping NMO Stack"')
+Plot('shgmres','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" plotcol=5 dash=1 title="" wantaxis2=n wantaxis1=n')
+Plot('shgmres3','shgmres','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" title="Shaping NMO Stack vs. PWC Stack"')
+Plot('shgmres2','spectra | window max1=100 | scale axis=1 | graph plotfat=3 label2=Amplitude unit2="" title="Shaping Operator Sensitivity"')
+
+Result('specn0','istackn shgmres','Overlay')
+
+Result('specn1','istack1 shgmres','Overlay')
+
+Result('spec8','shgmres2 shgmres','Overlay')
+
+Result('shpwspec','shgmres3 ugmres','Overlay')
+
+
+## All CMPs 
+
+Flow('shngmres','igaths ipick','shstack flo=2 fhi=100 niter=5 jump=4 velocity=${SOURCES[1]}')
+Result('shngmres', 
+	'''
+	put d2=0.0133333 |
+        grey title="Shaping NMO Stack" 
+        unit1=s label2=Midpoint unit2=km 
+       ''')
+
+Flow('shwnstack', 'shngmres', 'put d2=0.0133333 | window n1=600 f1=197 n2=300 f2=600')  
+Result('shwnstack', 'grey title="SNMO Stack" unit1=s label2=Midpoint unit2=km screenratio=2')
+Result('wshstack', 'shwnstack', 'grey title="SNMO Stack" unit1=s label2=Midpoint unit2=km ')
+
+Flow('shwnstack1', 'shngmres', 'put d2=0.0133333 | window n1=600 f1=97 n2=300 f2=300')  
+Result('shwnstack1', 'grey title="SNMO Stack" unit1=s label2=Midpoint unit2=km screenratio=2')
+Result('wshstack1', 'shwnstack1','grey title="SNMO Stack" unit1=s label2=Midpoint unit2=km')
+
+End()
diff --git a/book/tccs/shapestacks/paper.bib b/book/tccs/shapestacks/paper.bib
new file mode 100644
index 000000000..76ce6b38a
--- /dev/null
+++ b/book/tccs/shapestacks/paper.bib
@@ -0,0 +1,309 @@
+@Article{wisecup,
+    author={R[] D[] Wisecup},
+    title={Unambiguous signal recovery above the {N}yquist using random-sample-interval imaging},
+    year=1998,
+    journal={Geophysics},
+    volume=62,
+    issue=2,
+    pages={763-771}}
+@Article{stark,
+    author={T[] J[] Stark},
+    title={Signal recovery beyond conventional {N}yquist: {T}he sample rates used for seismic acquisition do not need to limit the maximum recoverable frequencies},
+    year=2013,
+    journal={The Leading Edge},
+    volume=32,
+    issue=11,
+    pages={1334-1339}}
+@Article{fomel,
+    author={S[] Fomel},
+    title={Shaping regularization in geophysical estimation problems},
+    year=2007,
+    journal={Geophysics},
+    volume=72,
+    issue=2,
+    pages={R29-R36}}
+@Article{fomel2,
+    author={S[] Fomel},
+    title={Nonlinear shaping regularization in geophysical inverse problems},
+    year=2008,
+    journal={78th Annual International Meeting, SEG, Expanded Abstracts},
+    pages={2046-2051}} 
+@Article{claerbout,
+    author={J[] Claerbout},
+    title={Earth soundings analysis: {P}rocessing versus inversion},
+    year=1992,
+    journal={Blackwell Science}}
+@Article{claerbout3,
+    author={J[] Claerbout},
+    title={Imaging the earth's interior},
+    year=1985,
+    journal={Blackwell Science}}
+@Article{sun,
+    author={Y[] Sun},
+    title={Inverse {NMO} stack in depth-variable velocity},
+    year=1997,
+    journal={Stanford Exploration Project (SEP) Report 94},
+    pages={213-222}}
+@Article{claerbout2,
+    author={J[] Claerbout},
+    title={A prospect for super resolution},
+    year=1996,
+    journal={Stanford Exploration Project (SEP) Report 93},
+    pages={133-136}}
+@Article{trickett,
+    author={S[] Trickett},
+    title={Stretch-free stacking},
+    year=2003,
+    journal={73rd Annual International Meeting, SEG, Expanded Abstracts},
+    pages={2008-2011}}
+@Article{kazemi,
+    author={N[] Kazemi and H[] R[] Siahkoohi},
+    title={Local stretch zeroing {NMO} correction},
+    year=2011,
+    journal={Geophysics},
+    volume=188,
+    issue=1,
+    pages={123-130}}
+@Article{fomel1,
+    author={S[] Fomel},
+    title={Transforming prestack seismic data by {G}ardner continuation},
+    year=2014,
+    journal={SEG Technical Program Expanded Abstracts},
+    pages={4643-4649}}
+@Article{fomel3,
+    author={S[] Fomel and Y[] Liu},
+    title={Seislet transform and frame},
+    year=2010,
+    journal={Geophysics},
+    volume=75,
+    issue=3,
+    pages={V25-V38}}  
+@Article{fomel4,
+    author={S[] Fomel},
+    title={Towards the seislet transform},
+    year=2006,
+    journal={76th Annual International Meeting: SEG},
+    pages={2847-2850}} 
+@Article{fomel5,
+    author={S[] Fomel},
+    title={Applications of plane-wave destruction filters},
+    year=2002,
+    journal={Geophysics},
+    volume=67,
+    pages={1946-1960}} 
+@Article{fomel6,
+    author={S[] Fomel},
+    title={Predictive painting of 3-{D} seismic volumes},
+    year=2010,
+    journal={Geophysics},
+    volume=75,
+    issue=4,
+    pages={A25-A30}}
+@Article{fomel7,
+    author={S[] Fomel and A[] Guitton},
+    title={Regularizing seismic inverse problems by model reparameterization using plane-wave construction},
+    year=2006,
+    journal={Geophysics},
+    volume=71,
+    issue=5,
+    pages={A43-A47}}  
+@Article{liu,
+    author={Y[] Liu and S[] Fomel and C[] Liu},
+    title={Signal and noise separation in prestack seismic data using velocity-dependent seislet transform},
+    year=2015,
+    journal={Geophysics},
+    volume=80,
+    issue=6,
+    pages={A25-A30}} 
+@Article{snieder,
+    author={R[] Snieder and J[] Trampert},
+    title={Inverse problems in geophysics},
+    year=2015,
+    journal={SEG Technical Program Expanded Abstracts},
+    pages={WD117-WD128}}
+@Article{silva,
+    author={M[] Silva and M[] Porsani and B[] Ursin},
+    title={Recursive stack to zero offset along local slopes},
+    year=2015,
+    journal={SEG Technical Program Expanded Abstracts},
+    pages={4323-4327}}
+@Article{mayne1,
+    author={W[] H[] Mayne},
+    title={Common reflection point horizontal data stacking techniques},
+    year=1962,
+    journal=Geophysics,
+    volume=27,
+    issue=6,
+    pages={927-938}}
+@Article{yilmaz,
+    author={O[] Yilmaz},
+    title={Seismic data analysis: Processing, inversion, and interpretation of seismic data},
+    year=2001,
+    journal={SEG}}
+@Article{shatilo,
+    author={A[] Shatilo and F[] Aminzadeh},
+    title={Constant normal-moveout ({CNMO}) correction: a technique and test results},
+    year=2000,
+    journal={Geophysical Prospecting},
+    volume=48,
+    pages={473-488}}
+@Article{rashed,
+    author={M[] Rashed},
+    title={Fifty years of stacking},
+    year=2014,
+    journal={Acta Geophysica},
+    volume=62,
+    issue=3,
+    pages={505-528}}
+@Article{ronen,
+    author={S[] Ronen and C[] Liner},
+    title={Least-sqaures {DMO} and migration},
+    year=2000,
+    journal={Geophysics},
+    volume=65,
+    issue=5,
+    pages={1364-1371}}
+@Article{ronen2,
+    author={S[] Ronen},
+    title={Wave-equation trace interpolation},
+    year=1987,
+    journal={Geophysics},
+    volume=52,
+    issue=7,
+    pages={973-984}}
+@Article{regimbal,
+    author={K[] Regimbal and S[] Fomel},
+    title={Improving resolution of {NMO} stack using shaping regularization},
+    year=2015,
+    journal={SEG},
+    pages={3859-3863}}
+@Article{kroode,
+    author={F[] Kroode and S[] Bergler and C[] Corsten and J[] Willem de Maag and F[] Strijbos and H[] Tijhof},
+    title={Broadband seismic data - The importance of low frequencies},
+    year=2013,
+    journal={Geophysics},
+    volume=78,
+    issue=2,
+    pages={WA3-WA14}}
+@Article{saad,
+    author={Y[] Saad and M[] Schultz},
+    title={{GMRES}: {A} generalized minimal residual algorithm for solving nonsymmetric linear systems},
+    year=1986,
+    journal={SIAM},
+    volume=7,
+    issue=3,
+    pages={856-869}}
+@Article{swan,
+    author={H[] W[] Swan},
+    title={Amplitude versus offset analysis in a finely layered media},
+    year=1988,
+    journal={SEG},
+    pages={1195-1198}}
+@Article{ma,
+    author={Y[] Ma and W[] Zhu and Y[] Luo and P[] Kelamis},
+    title={Super-resolution stacking based on Compressive Sensing},
+    year=2015,
+    journal={SEG},
+    pages={3502-3506}}
+@Article{miller,
+    author={R[] D[] Miller},
+    title={Normal moveout stretch mute on shallow-reflection data},
+    year=1992,
+    journal={Geophysics},
+    volume=57,
+    issue=11,
+    pages={1502-1507}}
+@Article{byun,
+    author={B[] S[] Byun and E[] S[] Nelan},
+    title={Method and system for correcting seismic traces for normal moveout stretch effects},
+    year=1997,
+    journal={U.S. Patent},
+    volume=5,
+    pages={684-754}}
+@Article{hicks,
+    author={G[] J[] Hicks},
+    title={Removing NMO stretch using the {R}adon and {F}ourier- {R}adon transforms},
+    year=2001,
+    journal={63rd Annual Conference and Exhibition, EAGE, Extended Abstracts},
+    pages={{A}-18}}
+@Article{hilterman,
+    author={F[] Hilterman and C[] Van Schuyver},
+    title={Seismic wide-angle processing to avoid NMO stretch},
+    year=2003,
+    journal={73rd Annual International Meeting, SEG, Expanded Abstracts},
+    pages={215-218}}
+@Article{rupert,
+    author={G[] B[] Rupert and J[] H[] Chun},
+    title={The block move sum normal moveout correction},
+    year=1975,
+    journal={Geophysics},
+    volume={40},
+    pages={17-24}}
+@Article{perroud,
+    author={H[] Perroud and M[] Tygel},
+    title={Nonstretch {NMO}},
+    year=2004,
+    journal={Geophysics},
+    volume={69},
+    pages={599-607}}
+@Article{masoomzadeh,
+    author={H[] Masoomzadeh and P[] J[] Barton and S[] C[] Singh},
+    title={Nonstretch moveout correction of long-offset multichannel seismic data for subbasalt imaging: Example from the North Atlantic},
+    year=2010,
+    journal={Geophysics},
+    volume={75},
+    issue={4},
+    pages={{R}83-{R}91}}
+@Article{zhang,
+    author={B[] Zhang and K[] Zhang and S[] Guo and K[] Marfurt},
+    title={Nonstretching NMO correction of prestack time-migrated gathers using a matching-pursuit algorithm},
+    year=2013,
+    journal={Geophysics},
+    volume={78},
+    issue={1},
+    pages={{U}9-{U}18}}
+@Article{dunkin,
+    author={J[] W[] Dunkin and F[] K[] Levin},
+    title={Effects of normal moveout on a seismic pulse},
+    year=1973,
+    journal={Geophysics},
+    volume={38},
+    pages={635-642}}
+@Article{barnes,
+    author={A[] E[] Barnes},
+    title={Another look at {NMO} stretch},
+    year=1992,
+    journal={Elsevier {S}cience {P}ublishing {C}o., {I}nc.}}
+@Article{zhdanov,
+    author={M[] S[] Zhdanov},
+    title={Geophysical inverse theory and regularization problems},
+    year=2002,
+    journal={Geophysics},
+    volume={57},
+    issue={5},
+    pages={749-751}}
+@Article{smith,
+    author={G[] C[] Smith and P[] M[] Gidlow},
+    title={Weighted stacking for rock property estimation and detection of gas},
+    year=1987,
+    journal={Geophysical {P}rospecting},
+    volume={35},
+    pages={993-1014}}
+@Article{bazelaire,
+    author={E[] de Bazelaire and J[] R[] Viallix},
+    title={Normal moveout in focus},
+    year=1994,
+    journal={Geophysical {P}rospecting},
+    volume={42},
+    pages={477-499}}
+@Article{keys,
+    author={R[] G[] Keys and D[] J[] Foster},
+    title={Comparison of seismic inversion methods on a single real data set},
+    year=1998,
+    journal={Tulsa: Society of Exploration Geophysicists}}
+@InProceedings{gardner,
+  author =       {G[] H[] F[] Gardner and S[] Y[] Wang and N[] D[] Pan and Z[] Zhang},
+  title =        {Dip moveout and prestack imaging},
+  booktitle = {18th Annual Offshore Technology Conference},
+  pages =        {OTC 5158},
+  year =         1986}
diff --git a/book/tccs/shapestacks/paper.tex b/book/tccs/shapestacks/paper.tex
new file mode 100644
index 000000000..36b28ad8b
--- /dev/null
+++ b/book/tccs/shapestacks/paper.tex
@@ -0,0 +1,529 @@
+\title{Increasing resolution of NMO stack using shaping regularization}
+\author{Kelly Regimbal$^{*1}$ and Sergey Fomel$^{1}$
+\\
+$^1$The University of Texas at Austin}
+
+\maketitle
+
+\lefthead{Regimbal and Fomel}
+\righthead{Shaping NMO stack}
+\footer{TCCS}
+
+\begin{abstract}
+Several problems arise with the assumptions and principles of
+conventional NMO and stack, resulting in lower amplitude and lower resolution stacks. We present two methods 
+that eliminate the effects of ``NMO stretch'' and restore a wider frequency band by replacing conventional 
+NMO and stack with a regularized iterative inversion to zero offset. The resulting stack is a model that 
+best fits the data using additional constraints imposed by the method of shaping 
+regularization. We use shaping regularization to achieve a 
+stack that has a denser time sampling and contains higher frequencies than the conventional 
+stack. In the first approach, we define the backward operator of shaping regularization
+using the principles of conventional NMO correction and stack. In the second approach, we introduce a recursive 
+stacking scheme using plane-wave construction in the backward operator of shaping regularization. The advantage of 
+using recursive stacking along local 
+slopes in the application to NMO and stack is that it avoids stretching effects caused by NMO correction and is insensitive to non-hyperbolic moveout in the data. Numerical tests confirm each algorithm's ability to attain a higher frequency stack with a denser temporal sampling interval compared to those of the conventional stack and to minimize stretching effects caused by NMO correction. We apply both methods to a marine dataset from the North Sea and achieve noticeable resolution improvements in the stacked sections compared with that of conventional NMO and stack. By treating NMO and stack
+as an iterative inversion using shaping regularization, resolution is enhanced by utilizing signal from different offsets and
+minimizing stretching effects to reconstruct a high resolution stack.
+\end{abstract}
+
+\section{Introduction}
+
+In conventional seismic data processing, common midpoint (CMP) stacking is one of the most fundamental processes \cite[]{yilmaz}.
+CMP stacking combines normal moveout (NMO)-corrected traces across a CMP gather to produce a single trace with a 
+higher signal-to-noise ratio \cite[]{rashed}. Several problems arise with the assumptions and principles that 
+set the foundation for conventional NMO and stack. Conventional NMO correction transforms traces at 
+non-zero offset into traces that are effectively at zero offset. NMO correction is, therefore, a time coordinate transformation,
+where time-variant nonlinear distortions of the seismic traces are introduced in the corrected traces \cite[]{barnes}.
+These undesirable distortions of signals on a seismic trace are known as NMO stretch and result in lower frequency 
+content of the corrected reflection event at far offsets \cite[]{dunkin}. 
+This violates the assumption of a uniform distribution of phase and frequency of 
+seismic reflections across the corrected gather. Common procedures to eliminate this stretching effect
+involve muting samples with severe distortions. This causes a decrease in fold and can destroy useful 
+far-offset information essential for amplitude versus offset (AVO) analysis \cite[]{swan}. Inaccuracy in stretch
+muting with residual stretching effects can also produce a lower-amplitude and lower-resolution stack \cite[]{miller}.  
+
+Conventional stacking is a simple process that combines the NMO-corrected traces by summing and normalizing.
+Stacking is based on the assumption that useful signal is coherent, whereas noise is random. In real seismic
+data, coherent noise and noise bursts are common, causing inaccuracy in the conventional stack \cite[]{rashed}.
+Traditional stacking also assumes that the NMO-corrected gather has perfectly aligned seismic reflections \cite[]{yilmaz}. 
+However, NMO correction is an approximation that assumes the travel-time as a function of offset follows a 
+hyperbolic trajectory in a CMP gather, which might fail in common geologic settings that involve lateral
+velocity variations or anisotropy \cite[]{bazelaire}.
+
+Several algorithms were developed to improve CMP stacking and enhance resolution 
+of stacked sections by reducing stretching effects. 
+\cite{claerbout} described inverse NMO stack, which recasts 
+NMO correction and stacking as an inversion process in the constant velocity case. This approach 
+combines conventional NMO and stack into one step by solving a set of simultaneous equations using
+iterative least-squares optimization. 
+\cite{sun} extended Claerbout's idea to the case of depth-variable 
+velocity. The inverse NMO stack operator applied depends on the hyperbolic moveout relation and can 
+be employed to remove non-hyperbolic events and random noise. 
+\cite{trickett} uses a variation of Claerbout's inverse NMO stack in his stretch-free stacking method 
+to avoid NMO stretch. Trickett's results tend to be higher frequency but noisier 
+than a conventional stack. Multiple other algorithms have been proposed that aim to reduce 
+NMO-stretching effects \cite[]{byun,hicks,hilterman,rupert,perroud,masoomzadeh,zhang,kazemi}.
+\cite{wisecup} introduced random sample interval imaging (RSI$^2$), which
+maps the CMP gather into the ``after NMO space'' using the exact moveout times and no 
+interpolation. The NMO-corrected values are collected in the stack, rather than 
+summed, where the input sample values are mapped to their correct time values in the stack. 
+\cite{shatilo} proposed a constant NMO-correction strategy, which applies a constant NMO shift 
+within a finite-time interval that is equal to the wavelet length of a trace. This approach eliminates wavelet 
+stretch and preserves higher frequencies than the conventional method, resulting in a higher resolution stack. However, 
+samples that exist in overlapping time windows are used twice during the correction, resulting in 
+an amplitude distortion.
+\cite{stark} discussed the idea of signal recovery beyond the conventional Nyquist 
+frequency using an approach similar to the RSI$^2$ algorithm. The method proposed is an output-driven process, 
+where the stack is defined as a merge trace and has a potentially higher sampling rate than the input 
+traces. Using this approach, the final stacked sections are not necessarily limited to the data-collected Nyquist frequencies.
+More recently, \cite{ma} proposed a stacking technique based 
+on a sparse inversion algorithm that computes the stack directly from a CMP gather by solving an optimization 
+problem using principles of compressive sensing. This method eliminates the stretch effect of conventional 
+CMP stacking and improves resolution in the stacked section. 
+\cite{silva} introduced a recursive stacking approach using local slopes to compute a stack without stretching effects. 
+This velocity-independent stacking technique starts at the farthest offset and accumulates the estimated stack until the near offset is reached. 
+
+In this paper, we propose two alternative approaches to NMO and stack using inversion by shaping regularization 
+\cite[]{fomel}. Shaping regularization implies a mapping of the input model to a space of acceptable models. 
+The shaping operator is integrated in an iterative inversion algorithm and provides an explicit
+control on the estimation result. In this application, the model is a stacked trace, and the data is a 
+CMP gather. We experiment with two different formulations of shaping regularization in the application of NMO and stack,
+ which we refer to as shaping NMO stack and plane-wave construction (PWC) stack, which both aim to optimize CMP stacking and 
+ enhance resolution of the final stacked seismic section. We start by reviewing shaping regularization in
+ the context of NMO and stack and define the operators for shaping NMO stack and PWC stack. We test these approaches 
+on synthetic examples to demonstrate each algorithm's ability to minimize 
+stretching effects and improve resolution. We then apply these methods to a field dataset from 
+the North Sea and, as a result, achieve noticeable resolution improvements in the stacked section in comparison 
+with conventional NMO and stack.
+
+
+\section{NMO and stack using shaping regularization}
+ 
+Geophysical estimation problems are often ill-posed due to insufficient data, 
+and small changes in the data may result in drastic changes in the model. Regularization is 
+used to solve ill-posed problems by providing additional constraints on the estimated model \cite[]{zhdanov}.
+Shaping regularization is defined as a mapping of the input model $\mathbf{m}$ to the space of acceptable functions. 
+The mapping is controlled by the shaping operator $\mathbf{S_m}$, which is integrated in an iterative
+inversion. Iterative inversion repeatedly uses the difference between the estimated model from the current step 
+and real data to update the model until convergence occurs \cite[]{ronen}. In the linear case, the solution of the estimation problem using 
+shaping regularization is defined by
+\begin{equation}
+\label{eq:inv}
+\mathbf{\hat{m}=[I+S_m(BF-I)]^{-1}S_mBd},
+\end{equation}
+where $\mathbf{F}$ is the forward operator, $\mathbf{B}$ is the backward operator, and $\mathbf{d}$
+is the input data. We implement the Generalized Minimum Residual (GMRES) algorithm \cite[]{saad} to solve the linear system in equation~\ref{eq:inv}. 
+GMRES is a general iterative method for inverting sparse, non-symmetric matrices.
+
+The main idea of shaping regularization in application to NMO stacking is to use signal from different 
+offsets to recover extra bandwidth in the stacked trace. In conventional NMO correction, estimates 
+of the zero-offset trace from the data recorded at different offsets may have different spectral bands 
+due to NMO stretch \cite[]{claerbout2}. At wide offsets, low frequencies exist, which extend the
+spectral bandwidth. We are interested in the recovery of both high and low frequencies. 
+We start by defining vectors and linear operators used in shaping NMO stack as follows:
+\begin{itemize}
+\item{$\mathbf{m}$ (model) is the desired seismic trace at zero-offset.}
+
+\item{$\mathbf{d}$ (data) is a CMP gather.}
+
+\item{$\mathbf{F}$ (forward operator) applies inverse moveout by spraying along hyperbolas using a known velocity model and subsamples in time.}
+
+\item{$\mathbf{B}$ (backward operator) interpolates the data to a denser temporal grid and applies NMO correction and stack.}
+
+\item{$\mathbf{S_m}$ (shaping operator) is a bandpass filter.}
+\end{itemize}
+The estimation result is controlled by the shaping operator. Therefore, the estimated stack is sensitive to how the 
+bandpass filter $\mathbf{S_m}$ is defined. The bounds of the bandpass operator depend
+on the range of frequencies that exist in the data and, hence, are defined by the dataset. 
+In the numerical experiments, we provide examples to illustrate how the 
+bandpass operator can be defined.   
+
+
+\section{Plane-wave construction stack using shaping regularization}
+
+Our second approach extends the method of shaping NMO stack further by introducing a recursive stacking scheme
+using plane-wave construction (PWC) \cite[]{fomel7} in the backward operator of shaping regularization to achieve a higher resolution stack. 
+The advantages of using recursive stacking along local slopes in the application to NMO 
+and stack is that it (1) avoids stretching effects caused by NMO correction and (2) is insensitive to 
+non-hyperbolic moveout in the data. We define the linear operators used in the shaping regularization 
+scheme in equation~\ref{eq:inv} as:
+\begin{itemize}
+\item{$\mathbf{F}$ (forward operator) applies predictive painting \cite[]{fomel6} to spread information using a known dip field 
+	and then subsamples in time.}
+\item{$\mathbf{B}$ (backward operator) interpolates the data in time to a denser grid and stacks in a recursive fashion using local slopes.}
+\item{$\mathbf{S_m}$ (shaping operator) is a bandpass filter that controls frequency content.}
+\end{itemize}
+
+% Recursive stacking
+We implement PWC stacking in the backward operator of shaping regularization which follows
+local-slope calculations from a CMP gather. The key idea of this stacking procedure is
+to start at the farthest offset trace of the gather and make a local slope prediction of the preceding trace using 
+PWC (Figure~\ref{fig:schematic}a). The partially corrected trace is then stacked with the uncorrected neighbor, 
+which is the input for the next local prediction (Figure~\ref{fig:schematic}b). The process is repeated in the 
+offset direction until the near-offset trace is reconstructed (Figure~\ref{fig:schematic}c). 
+The final step to this stacking scheme is to apply a near-offset NMO correction to reconstruct the zero-offset trace.
+This recursive stacking approach results in higher resolution stacks compared 
+to conventional NMO and stack. The procedure is equivalent to computing the zero scale of the seislet transform 
+\cite[]{fomel3}. Advantages of PWC stacking include eliminating 
+the effects of NMO stretch as well as the problem of non-hyperbolic moveout. 
+An approximate inverse of PWC stacking is defined by predictive painting \cite[]{fomel6}. 
+This algorithm is comprised of two main steps, namely estimating 
+local slopes of seismic events using plane-wave destruction (PWD) \cite[]{fomel5} and spreading information from a seed 
+trace inside a volume. To implement the forward operator $\mathbf{F}$, we use the updated model to spread information 
+across the CMP gather using the estimated dip field.
+
+\inputdir{.}
+\plot{schematic}{width=0.9\textwidth}{Schematic of PWC stacking algorithm. (a) Stack far offset trace T$_2$ with neighboring trace T$_1$, 
+	(b) stack updated trace T'$_1$ with neighboring trace T$_0$, (c) accumulated near-offset stack.} 
+
+In PWC stacking, each seismic trace is predicted from its neighbors that are shifted 
+along the event slopes. Slopes are estimated by PWD, which minimizes the prediction error to estimate optimal slopes.
+PWD can be sensitive to conflicting slopes at far offsets of a CMP gather where the dips are large.
+To account for this issue, we first apply a constant velocity NMO correction to the CMP gather, which results in 
+smoothly varying slopes without crossing events. 
+We then estimate the moveout $t(x)$ of the corrected seismic events at offset $x$ as follows:
+\begin{equation} 
+\label{eq:NMO2}
+t(x) = \sqrt{t_{0}^{2}+\frac{x^{2}}{v_{0}^{2}} + x^{2} \left(\frac{1}{v^{2}} - \frac{1}{v_{0}^{2}}\right)},
+\end{equation}
+where $t_0$ is the zero offset travel-time, $v$ is the NMO velocity estimated by a conventional method 
+and $v_0$ is a constant velocity. Adding the correction factor due to the constant velocity NMO correction 
+from equation~\ref{eq:NMO2}, the dip field becomes:
+\begin{equation} 
+\label{eq:dip2}
+p={\frac{x}{t}}\left(\frac{1}{v^{2}}-\frac{1}{v_{0}^{2}}\right).
+\end{equation}
+We use this estimated dip as the initial model for PWD. 
+This dip estimation scheme follows the velocity-dependent formulation of the seislet transform \cite[]{liu}
+and provides a better estimation of the dip field for CMP gathers with large dipping events and conflicting slopes at 
+far offsets. We incorporate this dip estimation method to compute the PWC stack in an iterative 
+fashion while using shaping regularization (equation~\ref{eq:inv}).
+
+\section{Examples}
+We evaluate the performance of shaping NMO stack and PWC stack by applying it to synthetic examples and a field dataset 
+from the North Sea. Numerical experiments test the ability of both algorithms to recover high frequencies and
+reduce the stretching effects common in conventional NMO correction. We additionally test the potential of the proposed
+methods to recover low frequencies and compare the results to conventional NMO and stack.
+
+\subsection{High frequency recovery}
+In our first experiment, we generated a synthetic trace with a sampling interval of 1 ms and 
+used it as a reference trace. This trace was then inverse NMO-corrected and subsampled to 4 ms 
+to produce a synthetic CMP gather, shown in Figure~\ref{fig:cmp2}, which became the input data for both 
+shaping methods and conventional NMO and stack. The result of applying a constant velocity NMO 
+correction is displayed in Figure~\ref{fig:nmo0} and is used to compute
+the dip field for the PWC stacking approach.
+
+The convergence of the GMRES algorithm in this example required only 3 iterations for PWC stack and 5 iterations for 
+shaping NMO stack to achieve a relative misfit tolerance of $10^{-5}$.
+The estimated shaping NMO and PWC stacks are displayed as a function of iteration in Figures~\ref{fig:shmod1} and~\ref{fig:mod1}, respectively.
+Iteration 0 of shaping NMO stack is the initial model and is comparable to a stack achieved by conventional NMO and stack 
+with lower amplitude and lower frequency content; iteration 5 is the estimation result. The PWC stack 
+exhibits a faster convergence, where iteration 3 shows the estimation result. The conventional stack shown in Figure~\ref{fig:mod4}
+results in lower amplitude and lower frequency content, while both PWC and shaping NMO stacks achieve results similar to the 
+reference trace.
+
+We next evaluate the spectral content recovered in the resulting stacks. The frequency content of the shaping NMO stack is compared
+with that of the PWC stack in Figure~\ref{fig:shpwspec1}. 
+In this synthetic example, the two shaping approaches reconstruct nearly the same amplitude scale and spectra. In Figure~\ref{fig:spec5,spec6},
+the spectral content of the PWC stack is compared with the conventional stack and the reference 
+trace. The conventional stack fails to recover frequencies ranging 
+from 110 Hz to 175 Hz, whereas the shaping NMO and PWC stacks recover useful frequencies up to 175 Hz. Thus, using inversion, we 
+accurately preserve the true amplitude scale and spectrum of the reference trace.
+Moreover, the high frequency information recovered is beyond
+the Nyquist frequency of the input data (125 Hz). \cite{ronen2} justifies the possibility of such recovery with the
+idea that a signal can be recovered with the combination of different aliased sequences. Since each offset
+corresponds to a different propagation path, it brings different information. Therefore,
+the combination of all offsets allows us to reconstruct a higher frequency signal that is not constrained
+by the Nyquist frequency of the input data. This is the fundamental principle of treating NMO and stack as an 
+inverse process using shaping regularization to recover higher resolution stacks. In this simple synthetic example, 
+we accurately preserve information in the recovered zero-offset trace with a 1-ms sampling 
+interval by using input data with only a 4-ms sampling interval. 
+
+
+To add complexity to the synthetic CMP gather, we incorporate random noise and artifical AVO effects, where amplitude
+is linearly decreasing with offest (Figure~\ref{fig:cmpavo}). The resulting shaping stacks and conventional stack are 
+compared with the reference trace in Figure~\ref{fig:mod5}. Both PWC stack and shaping NMO stack recover higher amplitude 
+in comparison with the conventional stack. We evaluate the difference in frequency content and amplitude recovered using 
+PWC stack and shaping NMO stack in Figure~\ref{fig:shpwspecavo}. Both approaches recover a similar bandwidth, but PWC
+stack recovers somewhat higher amplitude. A spectral comparison of the PWC stack versus conventional 
+NMO and stack is shown in Figure~\ref{fig:specavo2} and the reference trace in Figure~\ref{fig:specavo1}.
+The PWC stack still recovers higher frequency content in comparison
+to conventional NMO and stack. When comparing the spectral band of the PWC stack to that of the reference trace 
+in Figure~\ref{fig:specavo1}, the difference is minimal. Therefore, by adding the complexities of
+artificial AVO effects and random noise to the synthetic data, the ability of PWC stack to reconstruct the 
+zero-offset trace is not affected, whereas the frequency content recovered by the conventional stack is 
+limited to the sampling rate of the input data as well as other problems that arise with the assumptions of 
+conventional CMP stacking. The shaping NMO stack recovers approximately the same bandwidth as the reference 
+trace but has a mismatch in amplitude due to the decrease in energy across the input gather (Figure~\ref{fig:shspecavo1}).
+Conventional AVO methods are done through weighted 
+stacking \cite[]{smith} to display information about rock properties. Consequently, the principles of shaping NMO stack 
+and PWC stack can be extended to the application of AVO analysis to compute high resolution weighted stacks.  
+
+
+%Shaping NMO results
+\inputdir{synseis}
+\multiplot{2}{cmp2,nmo0}{width=0.6\textwidth}{(a) Synthetic CMP gather with a 4-ms sampling interval used as the input data for the two 
+	shaping stack methods and (b) constant velocity NMO-corrected gather to separate crossing events at far offsets used for PWC stacking.}
+
+\plot{shmod1}{width=0.8\textwidth}{Zoomed in portion of the estimated shaping NMO stack as a function of iteration.
+	Top: iteration 0 is similar to conventional NMO and stack, and bottom: iteration 5 is the estimation result, where frequency content 
+	and amplitude appear to be noticeably higher.}
+
+%PW stack
+\plot{mod1}{width=0.8\textwidth}{Zoomed in portion of the estimated PWC stack as function of iteration using shaping regularization.
+		Top: iteration 0 is the initial model and bottom: iteration 3 is the estimation result where convergence occurs.}
+
+%Both
+\plot{mod4}{width=0.8\textwidth}{Comparison of the resulting shaping stacks and conventional stack to the zero-offset reference trace.}
+
+\plot{shpwspec1}{width=0.7\textwidth}{Spectral comparison of the shaping NMO stack (blue) with the PWC stack (dashed red). In this synthetic experiment, 
+		both shaping stack approaches achieve similar results.} 
+
+%PW stack spectral comparisons
+\multiplot{2}{spec5,spec6}{width=0.6\textwidth}{Spectral comparison of the PWC stack (dashed red) with 
+	(a) conventional NMO and stack (blue) and (b) the reference trace (blue) with a 1-ms sampling interval. }
+
+
+%AVO example
+\plot{cmpavo}{width=0.8\textwidth}{Synthetic CMP gather with random noise and artificial AVO effects, where amplitude linearly decreases with offset.}
+
+\plot{mod5}{width=0.8\textwidth}{Comparison of the resulting shaping stacks and conventional stack with random noise and AVO effects added to the input gather to the zero-offset reference trace.}
+
+\plot{shpwspecavo}{width=0.7\textwidth}{Spectral comparison of resulting shaping stacks with random noise and artificial AVO effects applied 
+	to input data. The shaping NMO stack (blue) has a decrease in amplitude compared to the PWC stack (dashed red).}
+
+\multiplot{2}{specavo2,specavo1}{width=0.6\textwidth}{Spectral comparison of resulting stacks with random noise and artificial AVO effects applied 
+	to input data. Spectrum of the PWC stack (dashed red) with (a) conventional NMO and stack (blue) and (b) the reference 
+	trace (blue) with a 1-ms sampling interval.}
+
+\plot{shspecavo1}{width=0.7\textwidth}{Spectral comparison of resulting shaping NMO stack (dashed red) with random noise and AVO effects applied 
+	to input data with the reference trace (blue) with a 1-ms sampling interval.} 
+
+\plot{allspeclow}{width=0.7\textwidth}{Spectral comparison of the low frequencies recovered. Shaping NMO stack (dashed green) and PWC stack (dashed red)
+	recover lower frequencies in comparison with the conventional stack (dot-dash black) and are consistent with the reference trace (blue). Shaping 
+	NMO stack recovers lower frequencies than PWC stack.}
+
+\multiplot{2}{nmo3,nmo2}{width=0.4\textwidth}{(a) Conventional NMO correction and (b) NMO correction without stretch muting.}
+\multiplot{2}{shnsnmo,nsnmo}{width=0.4\textwidth}{Effective NMO correction using (a) shaping NMO stack and (b) PWC stack.}
+
+
+
+%\inputdir{synthetic}
+%\multiplot{2}{shspec5,shspec6}{width=0.7\textwidth}{Spectral comparison of the shaping NMO stack (dashed red) with 
+%		(a) conventional NMO and stack (blue) and (b) the reference trace (blue) with a 1-ms sampling interval.} 
+
+%\plot{shcmpnoise}{width=0.8\textwidth}{Synthetic CMP gather with random noise, which is used as the input data for shaping NMO stack
+%	and conventional NMO stack.}
+
+%\multiplot{2}{shspecnoise2,shspecnoise1}{width=0.7\textwidth}{Resulting spectral comparisons of stacks with random noise applied to input 
+%	data. Spectrum of the shaping NMO stack (dashed red) vs. (a) spectrum of conventional stack (blue) and (b) spectrum of reference trace (blue) 
+%	with 1-ms sampling interval.} 
+
+%\multiplot{2}{shspecavo2,shspecavo1}{width=0.7\textwidth}{Spectral comparison of resulting stacks with random noise and AVO effects applied 
+%	to input data. Spectrum of the shaping NMO stack (dashed red) vs. (a) conventional NMO and stack (blue) and (b) the reference 
+%	trace (blue) with a 1-ms sampling interval.} 
+
+
+% PW stack results
+%\inputdir{seislet}
+%\plot{nmo01}{width=0.8\textwidth}{Constant velocity NMO correction using minimum velocity. Estimated moveout (equation~\ref{eq:NMO2}) is overlain 
+%	to demonstrate how the crossing events now occur outside of the data. Using this approximation as the initial model for PWD, we are able to compute
+%	a more accurate dip field.}
+
+%\multiplot{2}{cmp2,nmo0}{width=0.7\textwidth}{(a) Synthetic CMP gather with a 4-ms sampling interval and (b) constant velocity NMO-corrected gather to separate crossing events at far offsets.}
+
+%\plot{cmpnoise}{width=0.8\textwidth}{Synthetic CMP gather with random noise, which is used as the input data for PWC stack and conventional NMO and stack.}
+
+%\multiplot{2}{specnoise2,specnoise1}{width=0.7\textwidth}{Resulting spectral comparisons of stacks with random noise applied to input 
+%	data. Spectrum of the PWC stack (dashed red) vs. (a) spectrum of conventional stack (blue) and (b) spectrum of reference trace (blue) 
+%	with 1-ms sampling interval.}
+
+
+%\plot{cmplow}{width=0.8\textwidth}{Synthetic CMP gather with a low-cut filter applied, cutting out all frequencies below 25 Hz.}
+
+
+
+\subsection{Low frequency recovery}
+Low frequencies play an important role in seismic inversion for velocity and impedance models \cite[]{kroode}.
+We next evaluate the ability of both algorithms to recover low frequency information in the 
+estimated stack. We apply a low cut filter to the synthetic CMP gather
+in order to remove all of the useful low frequency information below 25 Hz, which provides us with the input for both shaping
+methods and conventional NMO and stack. A spectral comparison of both approaches to conventional NMO and stack
+and the reference trace is shown in Figure~\ref{fig:allspeclow}. We zoom in to the low end of the frequency spectrum
+to see how well each method does in comparison to the reference trace. Conventional NMO stack fails
+to recover frequencies below 25 Hz, while shaping NMO stack and PWC stack recover 
+low frequency information that is consistent with the reference trace. In this experiment, shaping NMO stack 
+recovers more low frequencies compared to the PWC stack. 
+ 
+
+\subsection{``Non-stretch" NMO}
+To demonstrate how the shaping stacks reduce the effects of NMO stretch, 
+we seperately apply both methods to each trace of the synthetic CMP gather, setting other traces to zero. After repeating this process 
+for all traces in the gather, the output shaping stacks are concatenated to extract the effective NMO-corrected gathers.
+Conventional NMO correction with stretch muting fails to preserve far offset information (Figure~\ref{fig:nmo3}), 
+whereas NMO correction without stretch muting leads to prominent stretch at far offsets (Figure~\ref{fig:nmo2}). 
+Figure~\ref{fig:shnsnmo,nsnmo} compares the effective NMO-corrected gathers using PWC stacking and 
+shaping NMO stacking. The results indicate that the stretching effects that are prominent at far offsets and early times in the NMO-corrected gather 
+in Figure~\ref{fig:nmo2} are effectively reduced by implementing both shaping stack methods.
+
+
+\subsection{2-D North Sea data}
+
+\inputdir{elf}
+
+\plot{spec8}{width=0.7\textwidth}{Example of sensitivity to frequency bounds using shaping regularization. 
+		Spectral comparison of shaping NMO stack using a bandpass filter ranging from 2 Hz to 90 Hz (dashed red)
+		versus a bandpass filter ranging from 2 Hz to 124 Hz (blue). By using an upper bound that is too
+		large, spurious high frequencies are introduced.}
+
+\plot{mod}{width=0.8\textwidth}{Stacked result for one CMP gather. From top to bottom: 8-ms conventional stack, PWC stack,
+	shaping NMO stack, and dense 4-ms conventional stack}
+
+\plot{shpwspec}{width=0.7\textwidth}{Spectral comparison of the shaping NMO stack (blue) to the PWC stack (dashed red).
+	The PWC stack recovers slighly higher frequencies than the shaping NMO stack.}
+
+
+Our field dataset example is from the North Sea and was used previously by \cite{fomel4} and \cite{fomel3}. 
+This 2-D dataset was recorded over a complex salt dome region and has 1,000 CMP locations and 800 samples 
+per trace with a sampling interval of 4 ms. 
+We apply conventional NMO and stack to the 4-ms data to use as a reference for testing the accuracy of the stacking 
+approaches and subsample the data to 8 ms to provide the input data for shaping NMO stack, PWC stack and conventional NMO stack. 
+The shaping operator in this case is a bandpass filter ranging from 2 Hz to 100 Hz for shaping NMO stack and 
+3 Hz to 110 Hz for PWC stack. The parameters were selected based on the frequency content of the input data. 
+The estimation result using shaping regularization is sensitive to the defined bandpass filter. 
+Upper frequency bounds that are too high may lead to spurious high frequencies in the stack. 
+Figure~\ref{fig:spec8} demonstrates the effect of incorrect bandpass filter boundaries. 
+In this figure, the spectrum of the shaping NMO stack using a bandpass filter ranging from 2 Hz to 90 Hz is compared to 
+shaping NMO stack using a bandpass filter ranging from 2 Hz to 124 Hz. By using an upper bound that is too
+large in the latter case, spurious high frequencies are introduced. Therefore, the bounds of the bandpass filter
+are data-dependent and must be selected based on the frequency content of the data.
+The convergence of shaping NMO stack required 5 iterations, while PWC stack required only 3 iterations. 
+In Figure~\ref{fig:mod}, the result of applying the shaping methods to one CMP gather is compared to the resulting stacks from conventional
+NMO and stack and the densely-sampled reference stack.
+The frequency content and amplitude of the PWC stack appear to be the highest of the different stacking methods. We next evaluate the spectral content recovered for each approach by analyzing the amplitude spectra. We first compare the frequency content of the PWC stack to that of the shaping NMO stack, shown in Figure~\ref{fig:shpwspec}. 
+The PWC stack recovers slightly higher frequencies compared to the shaping NMO stack. 
+The spectrum of the shaping NMO stack is next compared to that of conventional NMO and stack in Figure~\ref{fig:specn0} 
+and the spectrum of the densely-sampled reference stack in Figure~\ref{fig:specn1}. Compared to the conventional stack, the 
+shaping NMO stack recovers extra bandwidth with higher frequencies present, which range to 80 Hz. 
+The accuracy of the recovered high frequencies is supported by the spectral comparison to the reference stack.
+We next analyze the ability of PWC stack to accurately recover higher frequencies, shown in Figure~\ref{fig:compspec,compspec2}.
+The PWC stack recovers higher frequency content compared to the conventional stack and also lower frequency content as 
+demonstrated in Figure~\ref{fig:compspec}. Compared to the dense stack (Figure~\ref{fig:compspec2}), PWC stack
+recovers slightly higher frequencies but is overall consistent with the spectral content.
+The resulting stacked sections using the 8-ms data as the input are displayed in Figures~\ref{fig:istack8},~\ref{fig:shngmres}, and~\ref{fig:ungmres3}. 
+For shallow reflectors, higher frequencies exist using shaping NMO stack and PWC stack in comparison with 
+the conventional NMO stack, which suggests that the effects of NMO stretch are reduced. This is seen clearly 
+in Figure~\ref{fig:wstack,shwnstack,wnstack}, which is a zoomed in section of a shallow region of the stacked sections.
+The PWC stacking approach recovers the most information out of the three approaches based on the thin layers that
+are resolved. Throughout the entire section, PWC stack and shaping NMO stack recover significantly higher frequencies compared to 
+the conventional stack. Thin layers that are unclear in the conventional stacked section become clearly resolved 
+and continuous in the PWC and shaping NMO stacked sections. Overall, events become more continuous and coherent throughout 
+the section using the shaping stack approaches, and resolution is noticeably improved. 
+
+\multiplot{2}{specn0,specn1}{width=0.6\textwidth}{Spectrum of the shaping NMO stack (dashed red) using the subsampled 8-ms
+		data versus (a) the conventional stack (blue) using the subsampled 8-ms 
+		data and (b) the conventional stack (blue) using the 4-ms data.}
+
+\multiplot{2}{compspec,compspec2}{width=0.6\textwidth}{Spectrum of PWC stack (dashed red) using the subsampled 8-ms 
+		data versus (a) the conventional stack (blue) using the subsampled 8-ms 
+		data and (b) the conventional stack (blue) using the 4-ms data.}
+
+\plot{istack8}{width=0.9\textwidth}{North Sea data example. NMO and stack using the conventional method with 8-ms data as input.}
+\plot{shngmres}{width=0.9\textwidth}{Stacked section using shaping NMO stack with 8-ms data as input.} 
+\plot{ungmres3}{width=0.9\textwidth}{Stacked section using PWC stack with 8-ms data as input.} 
+
+\multiplot{3}{wstack,shwnstack,wnstack}{width=0.25\textwidth}{Zoomed in section of the North Sea data using (a) conventional NMO and stack (b) shaping NMO stack and (c) PWC stack.} 
+
+
+%\subsection{Viking Graben data}
+%Our next dataset is the 2-D Mobil AVO Viking Graben Line 12, also from the North Sea \cite[]{keys}.
+%The dataset used was preprocessed, and most of the multiple reflections were attenuated. However, some residual multiple energy 
+%still remains in the dataset. 
+%This dataset contains 2,142 CMP locations and 1,001 samples per trace with a sampling interval of 4 ms. We subsample 
+%the data to 8 ms, which was the input for PWC stack, shaping NMO stack and conventional NMO stack. We apply conventional NMO and stack to the
+%4-ms data, which allowed us to have a densely-sampled reference stack in order to evaluate the accuracy of the shaping stack methods in recovering a broader frequency band. 
+%The shaping operator used is a bandpass filter ranging from 3 Hz to 90 Hz for shaping NMO stack and 2 Hz to 90 Hz for PWC stack. 
+%The convergence of shaping NMO stack required 6 iterations, while PWC stack required only 4 iterations. The result of applying
+%each method to one CMP location is shown in Figure~\ref{fig:mod3}. The conventional 8-ms stack appears to contain the lowest resolution, 
+%while the shaping stacks contain higher frequencies that are similar to the dense 4-ms stack. We next analyzed the frequency content
+%recovered by comparing the spectral bands of each stack. We first compare the spectrum of the PWC stack and shaping NMO stack in
+%Figure~\ref{fig:shspectest}. Overall, the frequency content recovered using both methods is similar, however, the shaping NMO stack recovers 
+%slightly lower frequencies than the PWC stack. The spectrum of the shaping NMO stack is compared to that of the conventional stack in Figure~\ref{fig:shspec1} 
+%and the spectrum of the densely-sampled reference stack in Figure~\ref{fig:shspec2}. The shaping NMO stack recovers significantly higher frequencies
+%compared to the conventional stack that are consistent with the dense stack. The shaping NMO stack also recovers lower frequencies than
+%the conventional stack. We next examine the frequency content of the PWC stack in comparison to the conventional stack in Figure~\ref{fig:pwspec2}
+%and the dense stack in Figure~\ref{fig:pwspec1}. The PWC stack is able to restore frequencies ranging to 80 Hz, while the conventional stack is
+%limited to approximately 60 Hz. Compared to the dense stack, the frequencies recovered using the PWC stack approach are accurate. 
+%The resulting stacked sections using the 8-ms data as the input to conventional NMO stack, shaping NMO stack, and PWC stack are displayed in Figures~\ref{fig:nmostack},~\ref{fig:shgmres1}, 
+%and~\ref{fig:pwgmres}, respectively. There are clear resolution improvements throughout the section using the shaping stack approaches in comparison with
+%the conventional stack. Shallow sections become more coherent and thin layers are resolved. We zoom in to a shallow portion of the stacked section, which
+%is displayed in Figure~\ref{fig:cwstack,shwstack,pwwstack}, to clearly detect the resolution improvements. In this windowed section, there appear to be dipping
+%beds that are truncated against a flat-lying upper layer. By using the shaping stack approaches, these dipping layers are coherent and continuous, whereas the
+%conventional stack appears low resolution and difficult to interpret. The flat layers above also have significant resolution improvements using the shaping
+%stack approaches, where thinner layers are detected in comparison with the conventional stack. Overall, by implementing the stacking schemes using
+%shaping regularization, we are able to reconstruct a higher resolution stack that is easier to interpret and not limited by the 8-ms sampling interval of the input gathers.  
+
+
+%\inputdir{viking}
+
+%\plot{mod3}{width=0.8\textwidth}{Viking Graben stacked results for one CMP gather. From top to bottom: dense 4-ms conventional stack, shaping NMO stack,
+%	PWC stack, and 8-ms conventional stack.}
+
+%\plot{shspectest}{width=0.7\textwidth}{Spectral comparison of the shaping NMO stack (blue) compared to the PWC stack (dashed red).
+%	Overall, the frequency content recovered using both methods is similar, however, the shaping NMO stack recovers slightly lower frequencies than the PWC stack.}
+
+%\multiplot{2}{shspec1,shspec2}{width=0.6\textwidth}{Viking Graben spectrum of the shaping NMO stack (dashed red) using the subsampled 8-ms
+%		data and 4-ms output stack versus (a) the conventional stack (blue) using the subsampled 8-ms 
+%		data and (b) the conventional stack (blue) using the 4-ms data.}
+
+%\multiplot{2}{pwspec2,pwspec1}{width=0.6\textwidth}{Viking Graben spectrum of PWC stack (dashed red) using the subsampled 8-ms 
+%		data and 4-ms output stack versus (a) the conventional stack (blue) using the subsampled 8-ms 
+%		data and (b) the conventional stack (blue) using the 4-ms data.}
+
+%\plot{nmostack}{width=0.9\textwidth}{Viking Graben data example. NMO and stack using the conventional method with the 8-ms data as input.}
+%\plot{shgmres1}{width=0.9\textwidth}{Viking Graben data example. NMO and stack using shaping NMO stack with the 8-ms data as input.} 
+%\plot{pwgmres}{width=0.9\textwidth}{Viking Graben data example. NMO and stack using PWC stack with the 8-ms data as input.} 
+
+%\multiplot{3}{cwstack,shwstack,pwwstack}{width=0.29\textwidth}{Zoomed in section of the Viking Graben data using (a) conventional NMO and stack (b) shaping NMO stack and (c) PWC stack.} 
+
+\section{Conclusions}
+% Compare to conventinal NMO and stack
+Conventional NMO stack can result in lower resolution stacked sections due to distortions caused by NMO correction and stretch
+muting. Treating the process of NMO and stack using iterative regularized inversion allows us to compute an optimal stack with
+higher frequency content and denser temporal sampling. As indicated by numerical experiments with synthetic data, 
+the recovered high frequency content is not limited
+by the sample rate of the input data and is consistent with the that of the reference trace for both shaping NMO and PWC stacking
+approaches. Furthermore, both stacking algorithms eliminate the effects of NMO stretch, which contributes to the 
+bandwidth of the reconstructed stack. Low frequency content also plays an important role in seismic data processing and imaging. As
+demonstrated by the synthetic and real data examples, PWC and shaping NMO stack both have the ability to recover both higher and
+lower frequencies compared to the conventional stack. The final stacked sections of the North Sea datasets have improved bandwidth and higher 
+resolution, which may aid in imaging and interpretation of small-scale features such as thin layers and channels. 
+
+% Advantages and disadvantages of each method
+The main difference between the two approaches described in this paper is that PWC stack follows local slopes of the data, while
+shaping NMO stack relies on velocity to compute the stack based on the conventional principles of NMO correction and stack.
+One benefit of shaping NMO stack, as demonstrated in the synthetic data example, is its ability to recover lower frequency content. 
+Advantages of PWC stack include its relative insensitivity to non-hyperbolic moveout and AVO effects. Furthermore,
+PWC stack recovered higher and lower frequencies in the real data example and converged in fewer iterations in comparison with shaping NMO stack.
+In the examples presented, both algorithms effectively reduce the effects of NMO stretch and improve the bandwidth of the
+final stack.
+
+% Broader impact
+By treating NMO and stack
+as an iterative inversion using shaping regularization, resolution is gained by utilizing signals from different offsets and
+minimizing stretching effects to reconstruct a high resolution stack that is not necessarily limited by the Nyquist
+frequency of the input data. In seismic processing, NMO correction and stack are simple processes applied to data. However,
+the same principles applied here can be extended to more complex processes, such as prestack migration, to achieve higher
+resolution results. With the increase in computational power, the extension of the inversion process that is proposed 
+in this paper to migration deserves further investigation and has the possibility to significantly enhance
+resolution in the resulting images.
+
+\section{Acknowledgments}
+We would like to thank the Texas Consortium for Computational Seismology (TCCS) members for helpful 
+discussions and TCCS sponsors for their financial support.
+
+\onecolumn
+\bibliographystyle{seg}
+\bibliography{paper}
+
+
diff --git a/book/tccs/shapestacks/synseis/SConstruct b/book/tccs/shapestacks/synseis/SConstruct
new file mode 100644
index 000000000..5b2b3c13a
--- /dev/null
+++ b/book/tccs/shapestacks/synseis/SConstruct
@@ -0,0 +1,476 @@
+from rsf.proj import *
+
+#################################################
+############# Synthetic CMP #####################
+n2=128
+d2=0.025 
+o2=0.05
+
+# 1. make a trace
+Result('test', None, 'spike n1=2001 d1=0.001  | noise rep=y seed=2014 | math output="input^3" | cut max1=0.25 | graph title=Spike')
+
+Flow('trace',None,
+	'''
+	spike n1=2001 d1=0.001 | noise rep=y seed=2014 | 
+	math output="input^3" | cut max1=0.25 | 
+	ricker1 frequency=75
+	''')
+Result('trace','graph title=trace')
+Result('spectrum','trace','spectra | graph title=Spectrum')
+
+# 2. make NMO velocity
+
+Flow('vel','trace','math output="1.5+x1" ')
+
+# 3. Apply inverse NMO 
+
+Flow('cmp','trace vel',
+	'''
+	spray axis=2 n=128 d=0.025 o=0.05 | 
+	inmo velocity=${SOURCES[1]} half=n 
+	''')
+Result('cmp','grey title="CMP Gather" label2=Offset unit2=km ')
+Flow('lowcmp','cmp','bandpass fhi=50') 
+Plot('lowcmp','cmp','bandpass fhi=50 | grey title="CMP Gather" label2=Offset unit2=km') 
+Result('lowcmp','cmp','bandpass fhi=50 | grey title="CMP Gather" label2=Offset unit2=km') 
+Result('mlowcmp','cmp','bandpass fhi=50 | grey wanttitle=n wantaxis=n') 
+
+Flow('cmpstack', 'cmp vel', 'nmo velocity=${SOURCES[1]} half=n | mutter v0=4.5 half=n | stack')
+Flow('difftest','cmpstack trace','add scale=1,-1 ${SOURCES[1]}')
+Result('modtest','difftest cmpstack trace',
+		'''
+		cat axis=2 ${SOURCES[1:3]} | 
+		dots labels=difference:interpolated:original 
+		gaineach=n
+		''')
+# Subsampled gather
+
+Flow('cmp2','cmp','window j1=4')
+
+Plot('cmp2','grey title="CMP Gather" label2=Offset unit2=km ')
+Result('cmp2','grey title="CMP Gather" label2=Offset unit2=km  ')
+Result('lowcmp2','cmp2','bandpass fhi=50 | grey title="CMP Gather" label2=Offset unit2=km') 
+Result('mlowcmp2','cmp','bandpass fhi=200 | window j1=4 | bandpass fhi=50 | grey wanttitle=n wantaxis=n')
+
+# Apply regular NMO
+Flow('nmo','cmp vel',
+	'''
+	nmo velocity=${SOURCES[1]} str=0.3 half=n | 
+	mutter v0=4.5 half=n 
+	''')
+Result('nmo','grey title="Conventional NMO Corrected Gather" label2=Offset unit2=km')
+Result('lownmo','nmo','bandpass fhi=50 | grey title="Conventional NMO Corrected Gather" label2=Offset unit2=km') 
+
+Flow('vel2','vel','window j1=4')
+Flow('nmo2','cmp2 vel2',
+	'''
+	nmo velocity=${SOURCES[1]} str=0.2 half=n |
+	mutter v0=4.5 half=n 	
+	''')
+Result('nmo2','grey title="With Stretch" label2=Offset unit2=km screenratio=2')
+Flow('nmo3','cmp2 vel2',
+	'''
+	nmo velocity=${SOURCES[1]} str=0.7 half=n |
+	mutter v0=4.5 half=n 	
+	''')
+Result('nmo3','grey title="NMO" label2=Offset unit2=km screenratio=2')
+Plot('nmo2','grey title="Conventional NMO" label2=Offset unit2=km')
+Result('lownmo2','nmo2','bandpass fhi=50 | grey title="Conventional NMO Corrected Gather" label2=Offset unit2=km') 
+
+# Regular stack
+
+Flow('stack2','nmo2','stack')
+Result('stack2','graph title="Subsampled Stack" ')
+Result('lowstack2','stack2','bandpass fhi=50 | graph title="Conventional Stack"')
+
+# Interpolate to dense sampling
+
+Flow('istack2','stack2 trace','bandpass fhi=124 | spline pattern=${SOURCES[1]}')
+Flow('diff2','istack2 trace','add scale=1,-1 ${SOURCES[1]}')
+Result('mod2','diff2 istack2 trace',
+		'''
+		cat axis=2 ${SOURCES[1:3]} | 
+		dots labels=difference:interpolated:original 
+		gaineach=n
+		''')
+
+Result('ispec2','trace istack2',
+		'''
+		cat axis=2 ${SOURCES[1]} | 
+		spectra 		 | 
+		window max1=75 	 | 
+		graph title=Spectrum
+		''')
+
+#########################################################
+################### Shaping Method ######################
+
+# Define dip
+
+data1 = 'cmp'
+
+# t as a function of t0 and x
+
+maxvel = 3.5
+minvel = 1.9
+
+Flow('t','vel',
+     '''
+     remap1 n1=2001 d1=0.001 o1=0 | 
+     spray axis=2 n=128 d=0.025 o=0.05 label=Offset unit=km |
+     math output="sqrt(x1*x1+x2*x2*(1/(input*input)-%g))"
+     ''' % (1/(minvel*minvel)))
+
+Plot('t','window j1=100 | window f1=3 | transp | graph wanttitle=n wantaxis=n yreverse=y min2=0 max2=2 pad=n')
+
+Flow('t1','vel',
+     '''
+     spray axis=2 n=128 d=0.025 o=0.05 label=Offset unit=km |
+     math output="x1*x1+(x2*x2)/(1/(input*input))"
+     ''')
+
+
+# dip as a function of t0 and x
+
+dt=0.001
+dx=0.025
+
+# NMO with maximum velocity
+
+Flow('vmax','vel','math output=%g' % minvel)
+Flow('nmo02','cmp2 vmax','nmo velocity=${SOURCES[1]} half=n')
+
+Flow('nmo0',[data1,'vmax'],'nmo velocity=${SOURCES[1]} half=n')
+Flow('inmo0','nmo0','bandpass fhi=125 | window j1=2')
+Plot('nmo0','grey title="NMO with Constant Velocity" label2=Offset unit2=km')
+Result('nmo0','grey title="NMO Corrected" label2=Offset unit2=km')
+Flow('lownmo0',['lowcmp','vmax'],'nmo velocity=${SOURCES[1]} half=n')
+Result('nmo01','nmo0 t','Overlay')
+
+Flow('p0','vel t',
+     '''
+     remap1 n1=2001 d1=0.001 o1=0 | 
+     spray axis=2 n=128 d=0.025 o=0.05 label=Offset unit=km |
+     math t=${SOURCES[1]} output="%g*x2/(t+0.001)*(1/(input*input)-%g)"
+     ''' % ((dx/dt),1/(minvel*minvel)))
+
+Flow('p1','vel t1',
+     ''' 
+     spray axis=2 n=128 d=0.025 o=0.05 label=Offset unit=km |
+     math t=${SOURCES[1]} output="%g*x2/(t+0.001)*(1/(input*input))"
+     ''' % (dx/dt))
+
+# dip as a function of t
+Result('p0','grey unit1=s unit2=km  color=j title=Slope scalebar=y')
+Flow('vdip2','p0 t','iwarp warp=${SOURCES[1]} eps=1')
+
+Flow('taper','vdip2','math output=1 | mutter v0=2.6 half=n | smooth rect1=10 rect2=10')
+
+Flow('vdip1','vdip2 taper','mul ${SOURCES[1]}')
+Flow('otaper','vdip2','math output=-7 | mutter v0=2.4 half=n inner=y | mutter v0=3.5 half=n | smooth rect1=15 rect2=15')
+Flow('ovdip1','vdip1 otaper','add ${SOURCES[1]}')
+Result('taper','grey title=Taper')
+Plot('vdip1','grey unit1=s unit2=km  color=j title=Slope scalebar=y')
+Result('ovdip1','grey unit1=s unit2=km  color=j title=Slope scalebar=y')
+
+data2 = 'nmo0'
+idip = data2+'idip2'
+Flow(idip,[data2,'ovdip1'],'dip nj1=1 both=n eps=1 rect1=15 rect2=15 idip=${SOURCES[1]} order=4')
+
+Plot(idip,'grey unit1=s unit2=km  color=j title=Slope scalebar=y label2=Offset')
+Result(idip,'grey unit1=s unit2=km  color=j title=Slope scalebar=y label2=Offset')
+
+Result('nmodip','cmp2 nmo0','SideBySideAniso')
+
+# Seislet NMO
+def seislet(gather,dip,snmo):
+    nmos = []
+    for i2 in range(n2):
+    	traced = 'trace%d' % i2
+    	if i2 == 0:
+    		Flow(traced,gather,'cut f2=1')
+    	elif i2 == n2-1:
+    		Flow(traced,gather,'cut n2=%d' % i2)
+    	else:
+    		Flow(traced,gather,'cut n2=%d | cut f2=%d' % (i2,i2+1))
+    		
+    	nmo = 'nmod%d' % i2
+    	nmos.append(nmo)
+    	Flow(nmo,[traced,dip],'''
+    	seislet dip=${SOURCES[1]} eps=0.1 adj=y inv=y type=haar order=4 |
+    	window n2=1 ''')
+    	Flow(snmo,nmos,'''cat axis=2 ${SOURCES[1:%d]} ''' % len(nmos))
+    	
+# Seislet Stack
+def seislet_stack(gather,dip,stack):
+    Flow(stack,[gather,dip],
+         '''
+         seislet dip=${SOURCES[1]} eps=0.1 adj=y inv=y type=haar order=4 |
+         window n2=1        
+         ''')
+
+#seislet('nmo0',idip,'dnmo2')
+#Plot('dnmo2','grey title="Seislet NMO"')
+#Result('dnmo2','grey title="Seislet NMO"')
+#Result('test2',['dnmo2','t'],'Overlay')
+
+############################## Shaping PW stack #################################
+
+# Backward operator
+def backward(data,model,mode):
+    Flow(model,[data,idip,'vmax'],
+         '''
+         spline pattern=${SOURCES[1]} | 
+	 nmo velocity=${SOURCES[2]} half=n |
+	 pwstack dip=${SOURCES[1]} mode=%d order=4 eps=0.01
+	 ''' % mode)
+    
+niter=5
+	 
+# Shaping PWC stack
+def shaping(cmp2,mod,plots):
+    '''Put everything in a function'''
+
+    global niter
+    mod0 = mod+'0'
+
+    backward(cmp2,mod0,1) 
+
+    m0 = mod0
+    old = m0
+    plots = []
+    for i in range(1,niter+1):
+    	new = '%s%d' % (mod,i)
+    	specn = 'spect%s%d' % (mod,i)
+    	Flow(new, [idip, old, m0, 'vmax'],'''
+    	pwpaint seed=${SOURCES[1]} eps=0.01 order=4 |
+    	inmo velocity=${SOURCES[3]} half=n  |
+    	window j1=4			 |	
+    	spline pattern=${SOURCES[0]} 	 |
+    	nmo velocity=${SOURCES[3]} half=n   |
+    	pwstack dip=${SOURCES[0]} order=4 mode=1 eps=0.01 | 
+    	add ${SOURCES[1]} ${SOURCES[2]} scale=-1,1,1 | 
+    	bandpass flo=1 fhi=300 ''')
+    	old = new
+    	
+    	Flow(specn, ['trace', new], '''
+    	cat axis=2 ${SOURCES[1]} | 
+    	spectra | 
+    	window max1=250 ''')
+    	Plot(specn, 'graph min2=0 max2=9 title=Spectrum%d'%i)
+    	if plots != None:
+    		plots.append(specn)
+    	
+    Flow(mod,new,'cp')
+    
+#shaping('cmp2','shstacked',None)
+
+
+############################## GMRES #################################
+backward('cmp2','gmres0',2)
+Flow('gmres1',['cmp2',idip],'shpwstack flo=1 fhi=300 niter=1 half=n jump=4 velocity=%g dip=${SOURCES[1]} mode=2 order=4 eps=0.01 nmo=y' % minvel)
+Flow('gmres2',['cmp2',idip],'shpwstack flo=1 fhi=300 niter=2 half=n jump=4 velocity=%g dip=${SOURCES[1]} mode=2 order=4 eps=0.01 nmo=y' % minvel)
+Flow('gmres3',['cmp2',idip],'shpwstack flo=1 fhi=300 niter=3 half=n jump=4 velocity=%g dip=${SOURCES[1]} mode=2 order=4 eps=0.01 nmo=y' % minvel)
+
+Result('niter','trace istack2 gmres3 gmres0',
+		'''
+		cat axis=2 ${SOURCES[1:4]} |
+		window f1=300 n1=500 | 
+		dots labels=reference:conventional:niter3:niter0 
+		gaineach=n
+		''')
+
+Result('mod1','gmres3 gmres2 gmres1 gmres0',
+		'''
+		cat axis=2 ${SOURCES[1:4]} | window f1=500 n1=500 |
+		dots labels=niter3:niter2:niter1:niter0
+		gaineach=n
+		''')
+
+# Spectral comparisons
+Plot('istack2','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude title="Conventional vs. PWC Stack"')
+Plot('trace','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude title="Reference vs. PWC Stack"')
+Plot('gmres3','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude plotcol=5 dash=1 title="" wantaxis2=n wantaxis1=n')
+
+Result('spec5','istack2 gmres3','Overlay')
+
+Result('spec6','trace gmres3','Overlay')
+
+
+
+Result('compspec','trace gmres3',
+		'''
+		cat axis=2 ${SOURCES[1]} | spectra | 
+		window max1=250 | graph label2=Amplitude 
+		title="Reference Trace vs. Seislet Stack"
+		''')
+
+Result('mod3','istack2 gmres3 trace',
+		'''
+		cat axis=2 ${SOURCES[1:3]} | 
+		dots labels=Conventional:Shaping:Reference 
+		gaineach=n
+		''')
+
+
+# With low cut filter
+Flow('cmplow','cmp','bandpass flo=25 | window j1=4 ')
+Flow('stacklow','cmplow vel2 vel','nmo velocity=${SOURCES[1]} half=n | stack | spline pattern=${SOURCES[2]} | bandpass flo=25')
+Flow('gmreslow',['cmplow',idip],'shpwstack flo=0 fhi=300 niter=4 half=n jump=4 velocity=%g dip=${SOURCES[1]} mode=2 order=4 eps=0.01 nmo=y' % minvel)
+#shaping('cmplow','shstacklow', None)
+Plot('gmreslow','spectra | window max1=35 | graph max2=3 label2=Amplitude plotcol=4 dash=4 title="" wantaxis2=n wantaxis1=n')
+Plot('traced','trace','spectra | window max1=35 | graph max2=3 label2=Amplitude title="" wantaxis2=n wantaxis1=n')
+Plot('stacklow','spectra | window max1=35 | graph max2=3 label2=Amplitude plotcol=7 dash=5 title="Low Frequency Spectrum"')
+Result('speclow','gmreslow traced stacklow','Overlay')
+
+
+Flow('widip',idip,'window j2=5 | window n1=200 f1=1000 n2=2 f2=10')
+Flow('pwdstack','cmp','window j2=5 | window n1=200 f1=1000 n2=2 f2=10')
+Flow('accstack','pwdstack widip','shpwstack flo=1 fhi=300 niter=2 half=n jump=1 dip=${SOURCES[1]} mode=2 order=4 eps=0.01 nmo=n')
+Plot('pwdstack','wiggle title=Data transp=y yreverse=y pclip=100')
+Plot('accstack','wiggle title="Accumulated stack" transp=y yreverse=y pclip=100 min2=0.1 max2=2.5')
+Result('schem','pwdstack accstack','SideBySideAniso')
+
+### Seislet Shaping NMO
+pwstacks = []
+for trace in range(n2):
+    pwonetrace = 'pwcmp-trace%d' % trace
+    if 0==trace:
+        Flow(pwonetrace,'cmp2','cut f2=%d' % (trace+1))
+    elif (n2-1)==trace:
+        Flow(pwonetrace,'cmp2','cut n2=%d' % trace)
+    else:
+        Flow(pwonetrace,'cmp2','cut n2=%d | cut f2=%d' % (trace,trace+1))
+
+    pwstack = 'pwcmp-%d-stack' % trace
+    pwstacks.append(pwstack)
+
+    # go through the shaping procedure: onetrace -> stack 
+    #Flow(stack,[onetrace,idip],'shpwstack flo=1 fhi=300 niter=4 half=n jump=4 velocity=%g dip=${SOURCES[1]} mode=1 order=4 eps=0.01 nmo=y' % minvel)
+    shaping(pwonetrace,pwstack,None)
+    
+Flow('nsnmo',pwstacks,'cat axis=2 ${SOURCES[1:%d]}' % n2)
+Result('nsnmo',
+       '''
+       transp | bandpass fhi=2 | transp | 
+       mutter v0=3.5 half=n | mutter slope0=0.4 half=n | grey title="PWC NMO" label2=Offset unit2=km screenratio=2
+       ''')
+
+############################## Shaping NMO stack #################################
+
+def shbackward(data,model):
+    Flow(model,[data,'vel'],
+         '''
+         spline pattern=${SOURCES[1]} |
+         nmo velocity=${SOURCES[1]} half=n str=0.0 |
+         stack
+         ''')
+shbackward('cmp2','shmod0')
+
+Flow('shgmres','cmp2 vel','shstack flo=1 fhi=350 niter=5 half=n jump=4 velocity=${SOURCES[1]}')
+Flow('shgmreslow','cmplow vel','shstack flo=1 fhi=350 niter=5 half=n jump=4 velocity=${SOURCES[1]}')
+
+Flow('shgmres1','cmp2 vel','shstack flo=1 fhi=300 niter=1 half=n jump=4 velocity=${SOURCES[1]}')
+Flow('shgmres2','cmp2 vel','shstack flo=1 fhi=300 niter=2 half=n jump=4 velocity=${SOURCES[1]}')
+
+Result('shmod1','shgmres shgmres2 shgmres1 shmod0',
+		'''
+		cat axis=2 ${SOURCES[1:4]} | window f1=500 n1=500 |
+		dots labels=niter5:niter2:niter1:niter0
+		gaineach=n
+		''')
+
+Plot('shgmres','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude title="SNMO Stack vs. PWC Stack"')
+
+Result('shpwspec1','shgmres gmres3','Overlay')
+
+
+Result('mod4','istack2 shgmres gmres3 trace',
+		'''
+		cat axis=2 ${SOURCES[1:4]} | 
+		dots labels=Conventional:SNMO:PWC:Reference 
+		gaineach=n
+		''')
+
+Plot('trace2','trace','spectra | window max1=35 | graph max2=3 plotfat=3 label2=Amplitude title=""')
+
+Plot('stacklow2','stacklow','spectra | window max1=35 | graph max2=3 plotfat=3 label2=Amplitude plotcol=7 dash=5 title="Low Frequency Recovery"')
+
+Plot('gmreslow2','gmreslow','spectra | window max1=35 | graph max2=3 plotfat=3 label2=Amplitude plotcol=5 dash=1 title="" wantaxis2=n wantaxis1=n')
+
+Plot('shgmreslow2','shgmreslow','spectra | window max1=35 | graph max2=3 plotfat=3 label2=Amplitude plotcol=3 dash=1 title="" wantaxis2=n wantaxis1=n')
+
+Result('allspeclow','trace2 gmreslow2 shgmreslow2 stacklow2','Overlay')
+
+
+# With AVO
+Flow('cmpavo1','trace vel',
+	'''
+	spray axis=2 n=%d d=%g o=%g  | 
+	inmo velocity=${SOURCES[1]} half=n |
+	noise seed=2015 var=0.1   	   |
+	math output="input*(1-x2*%g)" 	   
+	''' % (n2,d2,o2,2.0/5.0))
+Flow('cmpavo','trace vel',
+	'''
+	spray axis=2 n=%d d=%g o=%g  | 
+	inmo velocity=${SOURCES[1]} half=n |
+	noise seed=2015 var=0.1   	   |
+	math output="input*(1-x2*%g)" 	   |
+	window j1=4
+	''' % (n2,d2,o2,2.0/5.0))
+Result('cmpavo','grey title="CMP Gather with AVO" label2=Offset unit2=km ')
+Flow('nmo0avo','cmpavo1 vmax','nmo velocity=${SOURCES[1]} half=n')
+Flow('dipavo','nmo0avo ovdip1','dip nj1=1 both=n eps=1 rect1=15 rect2=15 idip=${SOURCES[1]} order=4')
+Flow('stackavo','cmpavo vel2 vel','nmo velocity=${SOURCES[1]} half=n  | mutter v0=4.5 half=n | stack | spline pattern=${SOURCES[2]}')
+Flow('pwavo','cmpavo dipavo','shpwstack flo=11 fhi=300 niter=4 half=n jump=4 velocity=%g dip=${SOURCES[1]} mode=2 order=4 eps=0.01 nmo=y | cut n1=255' % minvel) 
+
+Plot('stackavo','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude title="Conventional vs. PWC Stack"')
+Plot('pwavo','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude plotcol=5 dash=1 title="" wantaxis2=n wantaxis1=n')
+
+Result('specavo2','stackavo pwavo','Overlay')
+
+Result('specavo1','trace pwavo','Overlay')
+
+Flow('shgmresavo','cmpavo vel','shstack niter=6 flo=5 fhi=300 half=n jump=4 velocity=${SOURCES[1]} | cut f1=1 n1=200')
+Plot('shgmresavo','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude plotfat=3 plotcol=5 dash=1 title="" wantaxis2=n wantaxis1=n')
+Plot('shgmresavo1','shgmresavo','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude title="SNMO Stack vs. PWC Stack"')
+
+Result('shpwspecavo','shgmresavo1 pwavo','Overlay')
+
+Plot('trace1','trace','spectra | window max1=250 | graph max2=8 plotfat=3 label2=Amplitude title="Reference vs. Shaping NMO Stack"')
+Result('shspecavo1','trace1 shgmresavo','Overlay')
+
+Result('mod5','stackavo shgmresavo pwavo trace',
+		'''
+		cat axis=2 ${SOURCES[1:4]} | 
+		dots labels=Conventional:SNMO:PWC:Reference 
+		gaineach=n
+		''')
+# Try to extract effective NMO
+stacks = []
+for trace in range(128):
+    onetrace = 'cmp-trace%d' % trace
+    if 0==trace:
+        Flow(onetrace,'cmp2','add mode=m scale=128 |cut f2=%d' % (trace+1))  
+    elif 127==trace:
+        Flow(onetrace,'cmp2','add mode=m scale=128 |cut n2=%d' % trace)
+    else:
+        Flow(onetrace,'cmp2','add mode=m scale=128 |cut n2=%d | cut f2=%d' % (trace,trace+1))
+    stack = 'cmp-%d-stack' % trace
+    stacks.append(stack)
+
+    # go through the shaping procedure: onetrace -> stack 
+    Flow(stack,[onetrace, 'vel'],'shstack flo=2 fhi=200 niter=10 half=n jump=4 velocity=${SOURCES[1]}')
+
+Flow('shnsnmo',stacks,'cat axis=2 ${SOURCES[1:128]}')
+Result('shnsnmo', 'nsnmo',
+       '''
+       transp | bandpass fhi=2 | transp | 
+       mutter v0=4.5 half=n | mutter slope0=0.4 half=n | 
+       grey title="Shaping NMO" label2='Offset' unit2='km' screenratio=2
+       ''')
+
+Flow('nsstack','nsnmo','stack')
+
+End()