%------------------------------------------------------------------------------------
%           Migracao Kirchhoff 2,5D em angulo comum
%           (por empilhamento, em tempo e em amplitude verdadeira)
%------------------------------------------------------------------------------------
%------------------------------------------------------------------------------------
% EXECUCAO do programa:
%------------------------------------------------------------------------------------
% mig_ca1(sf1,x,h,t,v,xm,tm,ar,daf,dxm,ab,opm,opd);
%------------------------------------------------------------------------------------
%------------------------------------------------------------------------------------
% Parametros de ENTRADA:
%------------------------------------------------------------------------------------
% sf1    dado sismico em paineis de afastamento comum    Matriz   nx by nt by nh
%        - traco analitico com 1/2 derivada(tempo) -
%        (ver OBS6)
%        (formato de uma matriz 3D)
% dip    arquivo c/ mergulhos estimados (grau)           Matriz    nxm by ntm
% x      Coord. ponto medio (entrada) (m)                Vetor     1 by nx
% h      Meio afastamento (m)                            Vetor     1 by nh
% t      Coordenada de tempo (s)                         Vetor     1 by nt (PAR!!!!)
% v      Velocidade p/ migracao (m/s) - se v=cte      => Escalar   1 by 1  or
%                                     - se v=v(xm)    => Vetor     1 by nxm or
%                                     - se v=v(tm)    => Vetor     1 by ntm or        
%                                     - se v=v(xm,tm) => Matriz    nxm by ntm       
% xm     Coordenada de saida (m)                         Vetor     1 by nxm
% tm     Coordenada tempo migrado (s)                    Vetor     1 by ntm
% ar     Angulos de reflexao (grau)                      Vetor     1 by nar
% daf    Intevalo entre dips de mig (grau)               Escalar   1 by 1
% dxm    Intervalo entre coordenada de saida (m)         Esclar    1 by 1
% tw     estimativa do tamanho do pulso sismico          Escalar   1 by 1 
% opm => Opcao da migracao
%        opm=0 => migracao sem taper e abertura = comprimento da linha
%        opm=2 => migracao com abertura e taper no operador
%                 taper em +- 15% do angulo de mergulho (alfa) maximo (ver OBS7) 
% opd => Opcao do angulo de mergulho (alfa) de migracao (ver OBS8)
%        opd=1 => alfa max. e min. variavel(depende do ang.reflex de cada painel)
%                 c/ incremento cte. (nr.alfas varia p/ cada painel)
%        opd=3 => alfa max e min cte (depende do ang.reflex. max de saida
%                 segundo a relacao alfa+gamamax < 90)
%                 c/ incremento cte (nr.de alfas igual. p/ todos os paineis)
%
%------------------------------------------------------------------------------------
%------------------------------------------------------------------------------------
% Arquivos de SAIDA:
%------------------------------------------------------------------------------------
% mig_ca_arnn.mat => Arquivos binarios em disco no formato Matlab (.mat)
%                    para cada painel de ang. de reflexao comum (nn)
% logmig_ca.txt   => Arquivo com relatorio da migracao (tempo de execucao, etc)  
%
%------------------------------------------------------------------------------------
%------------------------------------------------------------------------------------
% Funcao externa
%  - cos_taper (calcula taper cossenoidal)
%------------------------------------------------------------------------------------
%------------------------------------------------------------------------------------
% agosto/2004
% Autora: Tânia M. Godinho Santiago (antares@petrobras.com.br)
%
% OBS:
%  1- Este programa foi realizado durante o mestrado em Ciências e Engenharia
%     do Petróleo (UNICAMP/FEM e IG - Campinas-SP-Brasil) cuja dissertação foi
%     defendida em 12/2004 com o título "Migração Kirchhoff 2,5D em tempo no
%     domínio de ângulo comum e em amplitude verdadeira".
%
%  2- Otimizacao do programa com a retirada do "loop" no tempo migrado 
%     deve-se à sugestão do prof. Dr. Lucio Tunes dos Santos - UNICAMP/IMECC.
%
%  3- A funcao externa "cos_taper.m" é de Xinxiang Li, CREWES project,
%     U of C. FEB. 1996.
%
%  4- Equações da parte cinemática da migração: "Fomel et al., 1999 -
%     Angle gather time migration - Stanford Exploration Project -
%     Report 100, p.141-151".
%
%  5- Equação do peso p/ migração em amplitude verdadeira: "Tygel, 2004 -
%     Migração em ângulo comum: tempo de trânsito e pesos de verdadeira amplitude -
%     Campinas: DMA/IMECC/Unicamp".
%
%  6- A geração do traco analítico e aplicação da 1/2 derivada em relação
%     ao tempo pode ser realizada usando a funcao "cria_hilb.m" (utiliza
%     a função "half.m" do prof. Dr. Lucio Tunes dos Santos).
%
%  7- O taper aplicado é cossenoidal iniciando (multiplica por 1) no ang. de mergulho
%     cujo valor é 15% abaixo (ou acima) do ang.mergulho máximo
%     e terminando (multiplica por 0) no ang.mergulho máximo
%     Ex: alfa max = +- 60 graus
%         início do taper (1) => alfa = +- 51 graus
%         fim do taper (0)    => alfa = +- 60 graus
%
%  8- opd => Opcao do angulo de mergulho (alfa) de migracao 
%      opd=1 => alfa max. e min. variavel(depende do ang.reflex-gama de cada painel)
%               c/ incremento cte. (nr.alfas varia p/ cada painel) =>
%                   alfa max. e min = +- (90 - gama - daf)
%      opd=3 => alfa max e min cte (depende do ang.reflex. max de saida
%               segundo a relacao alfa+gamamax < 90)
%               c/ incremento cte (nr.de alfas igual. p/ todos os paineis) =>
%                   alfa max. e min = +- (90 - gama max - daf)
%------------------------------------------------------------------------------------
%lembretes
%1- tirar ab da entrada e correlatos
%mig_cae1(ofallh_inc10,dip,xmid,h,t,3000,xm,tm,10,1,50,0.125,2,3);
%mig_cae1(inc10h_parc,dip,xmid,hn,tn,3000,xm,tm,10,1,50,0.125,2,3);

function mig_cae2(sf1,dip,x,h,t,v,xm,tm,ar,daf,dxm,tw,opm,opd);
diary logmig_ca.txt %grava neste arquivo o relatorio da migracao
if opm == 0  | opm == 2
 %
 %controle de tempo
 cput1=cputime;
 cli1=date; cli2=clock;
 disp(['....']);
 disp([num2str(cli1),'   ',num2str(cli2(4)),...
   ':',num2str(cli2(5)),' h']); 
 %
 ti = t .* 1000; %passa p/ ms
 arr = ar*pi/180; %passa para radianos
 dipr = dip*pi/180; %passa para radianos
 %dipr=dip; %jah estah em radianos
 
 nx = length(x); %nr. ptos.medios/secao de entrada
 nh = length(h); %nr. secoes de entrada
 nt = length(t); %nr. amostras/traco de entrada
 nxm = length(xm); %nr. coord.migradas
 nar = length(ar); %nr. secoes de saida (paineis ang.reflex.cte)
 ntm = length(tm); %nr. amostras migradas/traco
 ns = nh * nx; %nr. tracos de entrada
 nm = nar * nxm; %nr. tracos saida

 dafr = daf*pi/180; %passa para radianos
 dx = x(2) - x(1); %intervalo entre ptos.medios
 dh = h(2) - h(1); %intervalo entre offsets
 dt = t(2) - t(1); %razao de amostragem
 dtm = tm(2) - tm(1); %razao de amostragem migrada
 if nar ~= 1
  dar = ar(2) - ar(1); %intervalo entre ang.reflexao
 else
  dar = 0
 end
  %
 %nab = ab/dx; %nt traços na meia abertura (s/ contar o tr.central)
 
 %prepara campo de velocidade (dominio de saida)
 [m,n] = size(v); %(nr.linha,nr.coluna de v)
 if (m == 1) & (n == 1); v = v * ones(nxm,ntm); end; %p/ v=cte
 if (m == 1) & (n == nxm); v = v' * ones(1,ntm); end; %p/ v=v(xm)
 if (m == 1) & (n == ntm); v = ones(nxm,1) * v; end; %p/ v= v(tm)
 %p/ v=v(xm,tm) a matriz de entrada ja deve ser completa
  
 %parte final do peso (constante)
 fac= dafr*4*pi/sqrt(2*pi);
 
 disp(['....']);
 if opm == 0
  disp(['Migracao estah sendo feita sem restricao:']);
  disp(['- abertura = ',num2str((nx*dx)-dx),'m']);
  disp(['- sem taper (suave) no operador']);
 elseif opm == 2
  disp(['Migracao estah sendo feita com restricoes:']);
  %disp(['- abertura = ',num2str(2*ab),'m']);
  disp(['- taper (suave) no operador']);
  
 end
 disp(['....'])
 
 %sobre opcao do angulo de mergulho
 if opd == 1
  disp(['....'])
  disp(['Angulo maximo de mergulho limitado pelo angulo de reflexao de cada painel'])
  disp(['e com numero de amostras variavel'])
  disp(['....'])
 elseif opd == 3
  %-----------------cria vetor de ang. mergulho
  arrmax = max(arr);
  alfa = arrmax-(pi/2)+dafr : dafr : (pi/2)-arrmax-dafr;
  naf = length(alfa); %nr. elementos de alfa
  iaf = 1 : 1 :naf; %indices de alfa
  %--------------TAPER do operador em fcao do ang.mergulho
  %cria vetor de taper cossenoidal do operador (+-15% do alfa max)
  tpi1 = alfa(1) - alfa(1)*0.15;
  tpi2 = alfa(naf) - alfa(naf)*0.15;
  ind = find(alfa <= tpi1);
  itpi1 = length(ind);
  ind = find(alfa >= tpi2);
  itpi2 = ind(1);   
  %fcao externa CREWES
  coef = cos_taper(itpi2,naf); %coeficientes do taper
  ncoef = length(coef);
  fat1 = ones(1,naf);
  fat1(1:itpi1) = fat1(1:itpi1) .* coef(ncoef:-1:1);
  fat1(itpi2:naf) = fat1(itpi2:naf) .* coef;
  
  %auxiliares
  salfa = sin(alfa);       calfa = cos(alfa);    calfa2 = calfa.^2;
  
  disp(['....']);
  disp(['Angulo maximo de mergulho = +-',num2str(max(alfa)*180/pi),'^o']);
  disp(['Nr. de amostras igual a = ',num2str(naf)]);
  disp(['....']);
 end
 %-----------------------------  
  %informacoes p/ aparecer na tela
 disp(['....']);
 disp(['ENTRADA:']);
 disp(['nr.offsets-h = ',num2str(nh)]);
 disp(['int.h = ',num2str(dh),' m']);
 disp(['nr.CMP-x/secao = ',num2str(nx)]);
 disp(['int.x = ',num2str(dx),' m']);
 disp(['tini = ',num2str(t(1)),' s  -   tfin = ',num2str(t(nt)),' s']);
 disp(['R.A. = ',num2str(dt),' s']);
 disp(['total de tracos = ',num2str(ns)]);
 disp(['total de amostras = ',num2str(nt*ns),]);
 disp(['....']);
 disp(['SAIDA:']);
 disp(['nr.ang.reflexao-gm = ',num2str(nar)]);
 disp(['int.gm = ',num2str(dar),' grau']);
 disp(['nr.coor.mig.-xm/secao = ',num2str(nxm)]);
 disp(['int.xm = ',num2str(dxm),' m']);
 disp(['amostras/traco = ',num2str(ntm)]);
 disp(['R.A. = ',num2str(dtm),' s']);
 disp(['total de tracos = ',num2str(nm)]);
 disp(['....']);
  
 %----------------  inicio dos loops  --------------------
  
 for iar = 1 : nar %loop no angulo reflex 
  m1 = zeros(nxm,ntm); %cria matriz zerada de saida (nxm by ntm)
  %--------------------escreve no arquivo de log
  disp([' ---------- secao de ang.reflexao = ',num2str(ar(iar))])
  cli3=date; cli4=clock;
  disp([num2str(cli3),'   ',num2str(cli4(4)),...
    ':',num2str(cli4(5)),' h']);
  
  %-----------------cria vetor de ang. mergulho, taper e auxliares
  %--- p/ opd == 1
  if opd == 1   
   dipmax=88*pi/180;
   afmax=min(dipmax,(pi/2)-arr(iar)-dafr);
   alfa = -afmax : dafr : afmax;
   naf = length(alfa); %nr. elementos de alfa
   iaf = 1 : 1 :naf; %indices de alfa
   %--------------------------------------------
      %--------------TAPER do operador em fcao do ang.mergulho
   %cria vetor de taper cossenoidal do operador (+-15% do alfa max)
   tpi1 = alfa(1) - alfa(1)*0.15;
   tpi2 = alfa(naf) - alfa(naf)*0.15;
   ind = find(alfa <= tpi1);
   itpi1 = length(ind);
   ind = find(alfa >= tpi2);
   itpi2 = ind(1);   
   %fcao externa CREWES
   coef = cos_taper(itpi2,naf); %coeficientes do taper
   ncoef = length(coef);
   fat1 = ones(1,naf);
   fat1(1:itpi1) = fat1(1:itpi1) .* coef(ncoef:-1:1);
   fat1(itpi2:naf) = fat1(itpi2:naf) .* coef;
   
   %auxiliares
   salfa = sin(alfa);   calfa = cos(alfa);    calfa2 = calfa.^2;
  end
    
  %-----------------------------------------------------------
  %auxiliares
  sgm = sin(arr(iar));  sgm2 = sgm^2;  cgm = cos(arr(iar)); dcs = calfa2 - sgm2;
  
  %aux do dip estimado (matrizes nxm X ntm)
  sdip = sin(dipr);       cdip = cos(dipr);    cdip2 = cdip.^2;
  dcsg = cdip2 - sgm2;
  
  %calcula  parte 1 do peso Tygel e multiplica pelo taper (fat1)
  w1 = (fat1) .* sqrt( calfa * cgm ) ./ dcs; %parte do w multipl.atenuacao
  
    
  for ixm = 1 : nxm %loop na coord. migrada
  %--------------------escreve no arquivo de log
  disp([' ---------- x de saida = ',num2str(xm(ixm)/dxm)])
   %auxiliares
   v2tm = v(ixm,:) .* tm/2;

   %calcula xe = coord. pto estacionario p/ todos os tm's
   xe = xm(ixm) + v2tm .* (sdip(ixm,:) .* cdip(ixm,:) ./ dcsg(ixm,:)); %vetor 1Xntm 
   
   %calcula, em torno de cada xe, ZFresnel_1 projetada p/todos os tm's
   fz1 = v(ixm,:) .* sqrt(tw.*tm.*cdip(ixm,:).*cgm) ./ (2*dcsg(ixm,:)); %vetor 1Xntm
   
   %calcula rab (dist entre xm e xe) e acha a max abertura
   rab = abs(xm(ixm) - xe); %vetor 1Xntm
   ab = max(fz1 + rab) %escalar
   
   %definindo abertura
   abi = xm(ixm) - ab; abf = xm(ixm) + ab;
   if abi <= x(1); abi = x(1); end;
   if abf >= x(nx); abf= x(nx); end;
   %indices de x limites da abertura   
   iabi = floor(((abi - x(1))/dx)) +1;
   iabf = ceil(((abf - x(1))/dx)) +1;
   
   %se a opcao de migracao for sem taper e "sem" abertura (opm=0)
   if opm == 0; abi = x(1); abf= x(nx); iabi = 1; iabf = nx; end
   
   %calcula  parte 2 do peso Tygel
   w2 = v(ixm,:).* sqrt(tm.^3);
   
   
   %calcula curvas de difracao para todos os tm's da coordenada xm(ix),
   %jogando nas matrizes (ntm X naf) de tc (tempo nao migrado), hc (meio offset)
   %e x (coord.nao migrada)
   tc = tm' * ( calfa * cgm ./ dcs);
   hc = v2tm' * (sgm * cgm ./ dcs);
   xc = xm(ixm) + v2tm' * (salfa .* calfa ./ dcs);
   
   % passa tc p/ ms (diminuir tempo de exec da interpolacao)
   tc = tc .* 1000;
   
   
   %kizumba p/ zerar xc's fora da abertura centrada no pto estacionario
   %1-calcula xi e xf (coord ini e fim de cada ab centrada) p/ cada tm
   xi = xe - fz1; %vetor 1Xntm
   xf = xe + fz1; %vetor 1Xntm
   
   %2-
   ci = ones(1,naf); cf = ci;
   xiaux = xi'* ci; xfaux = xf'*cf; %matrizes ntm X naf
   
   %3-acha os valores de xc que sao menores que xi ou maiores que xf
   %  e coloca 9999999 nestas posicoes => sairao fora das coord da linha
   lli = find(xc < xiaux);
   xc(lli) = 999999999.0;
   llf = find(xc > xfaux);
   xc(llf) = 999999999.0;
   
   %interpola as amplitudes
   %gera matriz (ntm X naf) com amplitudes interpoladas
   ampi = interp3(ti,x(iabi:iabf),h,sf1(iabi:iabf,:,:),tc,xc,hc,'cubic');
   
   %tira nan
   indn = find(isnan(ampi));
   ampi(indn) = 0;
   
   %aplica peso e soma as amplitudes interpoladas p/cada tm
   %gera vetor (1 X ntm)
   amps = real( w2 .* (ampi * w1')' );
   
   %soma as amplitudes de cada tm na matriz de saida e multiplica por    % fac=dafr*4*pi/sqrt(2*pi)
   m1(ixm,:) = (m1(ixm,:) + amps).*fac;   
  end % fim do loop na coord. migrada
  
  if opm == 0 | opm == 2
   %grava painel de ang.comum (ar(iar)) como arquivo .mat (mig.sem restricao ou semi)
   saida1 = ['mig_ca_ar',num2str(ar(iar))];
   save(saida1,'m1');
  end
 end %fim do loop no angulo reflex
 
 %----------------  fim dos loops  -----------------------
 
 disp(['....']);
 %controle de tempo
 disp(['Inicio:  ',num2str(cli1),'   ',num2str(cli2(4)),...
   ':',num2str(cli2(5)),' h']);
 clf1=date; clf2=clock;
 disp(['Fim:     ',num2str(clf1),'   ',num2str(clf2(4)),...
   ':',num2str(clf2(5)),' h']);
 disp(['....']);
 et=etime(clf2,cli2);
 disp(['elapsed_time:   ',num2str(et/60),' min'])
 cput2=cputime;
 disp(['cpu_time:       ',num2str((cput2-cput1)/60),' min']) 
 
elseif opm ~=0 & opm ~= 2
 disp(['....']);
 disp(['Parametro "opm" (opcao de migracao) estah incorreto']);
 disp(['....']);
end 
diary off %fecha arquivo com o relatorio da migracao
