#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')


# In[2]:


from pkprocess import *


# In[3]:


from scipy.ndimage.filters import laplace


# In[4]:


sd=read_su('data/marmousi.su')


# # Kirchhoff prestack depth migration using the true velocity model

# In[5]:


vel=np.fromfile("data/marm16km.drt",dtype=np.float32)


# In[6]:


print(vel)


# In[7]:


nx=576
nz=188
h=0.016
vel.shape=(nx,nz)


# In[8]:


plot_vel(vel,h)


# In[9]:


times=np.empty((nx,nx,nz),dtype=float32)


# In[10]:


for isrc in range(nx):
    if isrc %50 ==0:
        print(isrc,nx)
    times[isrc]=traveltime(vel,h,isrc,0)
#times.tofile("ttime_marm_org.bin")
#times=np.fromfile("ttime_marm_org.bin",dtype=float32)
#times.shape=(nx,nx,nz)


# In[11]:


image=kirchhoff(sd,h,times,0.1)


# In[12]:


dimg=zdiff2(image)
dimg.min(),dimg.max()


# In[13]:


plot_mig(perc_clip(dimg,97),h)


# # Kirchhoff prestack depth migration using a smooth velocity model

# In[15]:


svel=moving_average2d(vel,10,10)
plot_vel(svel,h)


# In[16]:


stimes=np.empty((nx,nx,nz),dtype=float32)
for isrc in range(nx):
    if isrc %50 ==0:
        print(isrc,nx)
    stimes[isrc]=traveltime(svel,h,isrc,0)
#stimes.tofile("ttime_marm_sm10.bin")
#stimes=np.fromfile("ttime_marm_sm10.bin",dtype=np.float32)
#stimes.shape=(nx,nx,nz)


# In[17]:


simage=kirchhoff(sd,h,stimes,0.1)


# In[18]:


dsimg=zdiff2(simage)
dsimg.min(),dsimg.max()


# In[19]:


plot_mig(perc_clip(dsimg,97),h)


# In[ ]: