% ------------------------------------------------------------------------- %
 %                                                                           %
 %          Elimination of Radon artifacts using the KL transform            %
 %                                                                           %
 %                M.D.Sacchi - for Geoph 426 - UofA                          %
 %                                                                           %
 % I have used a large trade-off parameter in the Parabolic Radon            %
 % transform. This is to guarantee stability. The problem is                 %
 % that a large trade-off will produce artifacts in the estimated            %
 % primaries (residual energy from the multiples is left                     %
 % in the panel of primaries). The artifacts are attenuated                  %
 % with the KL transform.                                                    %
 %                                                                           %
 % From SeismicLab: pradon_demultiple.m                                      %
 %                  readsegy.m                                               %
 %                  kl.m                                                     %
 %                  pimage.m                                                 %
 %                  clip.m                                                   %
 %                  seismic.m                                                %
 %                                                                           %
 % ------------------------------------------------------------------------- %



  clear;
  close all;


 % Read data, a nmo-ed corrected marine gather
 % from the Gulf of Mexico (./SeismicLab_data)

  [D,H] = readsegy('gom_cdp_nmo.su');

 % We will use a subset of the data

  D = D(800:1200,:);

 % Extract areas that were muted

  Mutes = ones(size(D)); 
  I = find(D==0); 
  Mutes(I) = 0;

 % Extract offet and sampling interval

  h = [H.offset];
  dt = H(1).dt/1000/1000;

 % Define parameters for Radon Demultiple

  qmin = -0.3;
  qmax = .8;
  nq = 120;
  flow = 2.;
  fhigh = 90;
  mu = 1;              % Trade-off parameter
  q_cut = 0.05;


  [P,M,tau,q] = pradon_demultiple(D,dt,h,qmin,qmax,nq,flow,fhigh,mu,q_cut);

 % Apply data mutes back to the estimated primaries

  P = P.*Mutes;

 % KL transform filter using a subspace of 15 eigenvectors

  P_kl = kl(P,15);
  P_kl = P_kl.*Mutes;
 
 % Plot data, Radon panel, Primaries and Primaries after
 % fitering with the kl transform

 Big =  [D,10*M,P,P_kl,P-P_kl];
 [nt,nx]=size(Big);

 pimage([1:1:nx],tau,clip(Big,50,50));

 colormap(seismic(2));