#!/fs/pkgs/share/perl/bin/perl -w 

$arch_las = "/fs/proyectos/epep/geopetro/Trabajo/CI_135.las";
$arch_sal = "/fs/proyectos/epep/geopetro/Trabajo/borrame";

open (LAS, "<$arch_las") or die "No puede abrir $arch_las: $!";

for ($i=1; $i<=7; ++$i) # y en caso de ser mas o menos columnas, hay que estar sustituyendo aca ????????????
  {
  $j[$i] = 0; $j1[$i] = 0;
  }
while (<LAS>)  # Lectura de todas las lineas  
    {
    $linea = $_;
    if ($linea=~ /^\s*[1-9]/)
       {
       @tmp = split (/\s+/, $linea);  # Recorta la linea en el # de columnas existentes
    
       $rt1[$j1[2]]   = $tmp[2]; ++$j1[2]; ###########################################
       $vcl1[$j1[3]]  = $tmp[3]; ++$j1[3]; #
       $pign1[$j1[4]] = $tmp[4]; ++$j1[4]; #   
       $suwi1[$j1[5]] = $tmp[5]; ++$j1[5]; # Guarda los valores originales en vectores
       $kint1[$j1[6]] = $tmp[6]; ++$j1[6]; #
       $r351[$j1[7]]  = $tmp[7]; ++$j1[7]; ############################################
       ++$j1;
       #printf ("%10.4f    %10.4f    %10.4f   %10.4f   %10.4f   %10.4f    %10.4f \n", $j1, $rt1[$j1[2]], $vcl1[$j1[3]], $pign1[$j1[4]], $suwi1[$j1[5]], $kint1[$j1[6]], $r351[$j1[7]]);
       
################################ OMISION DE REGISTRO FALTANTES #####################################
        
        if (int($tmp[2]) != -999) #estoy obligando a rt ser la columna 2 ???????????????????????
           {
           $rt[$j[2]] = $tmp[2];                   #  prof: Profundidad 
           $prof_rt[$j[2]]  = $tmp[1];             #  rt:   True Resistivity
           ++$j[2];
           }
        if (int($tmp[3]) != -999) 
           {
           $vcl[$j[3]]      = $tmp[3];             #  vcl:  Volume of Clay Type Material Relative to Total Volume
           $prof_vcl[$j[3]] = $tmp[1];
           ++$j[3];
           }
        if (int($tmp[4]) != -999) 
           {
           $pign[$j[4]]      = $tmp[4];            #  pign: Effective Porosity, VISO, VPAR, and Bound Water Removed 
           $prof_pign[$j[4]] = $tmp[1];
           ++$j[4];
           }
        if (int($tmp[5]) != -999) 
           {
           $suwi[$j[5]]      = $tmp[5];            #  suwi: Non-Clay Water Intergranular Water Saturation for Undisturbed Zone
           $prof_suwi[$j[5]] = $tmp[1];
           ++$j[5];
           }
        if (int($tmp[6]) != -999)        # Es necesario colocar el [6] en $kint[$j[6]] ????????????????????
           {
           $kint[$j[6]]      = $tmp[6];            #  kint: Intrinsic Permeability
           $prof_kint[$j[6]] = $tmp[1];   
           ++$j[6];
           }
        if (int($tmp[7]) != -999) 
           {
           $r35[$j[7]]      = $tmp[7];             #  r35:  RADIO DE GARGANTES DE POROS
           $prof_r35[$j[7]] = $tmp[1];
           ++$j[7];
           }   # OJO: No habra problemas cuando haya mas o menos de 7 columnas??????????????????????

############################### CIERRE OMISION DE REGISTRO FALTANTES ##############################
      # printf ("%15.4f    %15.4f    %15.4f   %15.4f   %15.4f   %15.4f    %15.4f \n", $j1, $rt1[$j1], $vcl1[$j1], $pign1[$j1], $suwi1[$j1], $kint1[$j1], $r35[$j1]);
       }  
    }   

$long_rt = $#rt; print("Tam_rt   = $long_rt\n"); print("Tam_rt1   = $#rt1\n"); 
$long_vcl = $#vcl; print("Tam_vcl  = $long_vcl\n"); print("Tam_vcl1  = $#vcl1\n");
$long_pign = $#pign; print("Tam_pign = $long_pign\n"); print("Tam_pign1 = $#pign1\n");
$long_suwi = $#suwi; print("Tam_suwi = $long_suwi\n"); print("Tam_suwi1 = $#suwi1\n");
$long_kint = $#kint; print("Tam_kint = $long_kint\n"); print("Tam_kint1 = $#kint1\n");
$long_r35 = $#r35; print("Tam_r35  = $long_r35\n"); print("Tam_r351  = $#r351\n"); #$long_r35 = $#r35;

close(LAS);
open (SAL, ">$arch_sal") or die ("No puede abrir $arch_las: $!");
print ("Introduzca el ancho de banda (b): ");
$b = <STDIN>;
chop ($b);

############################ LLAMADO A RUTINA QUE HACE LA SUVIZACION E IMPRESION DE RESULTADOS ###########################

if ($b > 0)
  {
  printf (SAL "       Prof Obj         Intensidad RT\n");
  guardar(\@prof_rt, \@rt, \@rt1, 2, $long_rt);                                
  printf (SAL "       Prof Obj         Intensidad VCL\n");
  guardar(\@prof_vcl,\@vcl, \@vcl1, 3, $long_vcl);                              
  printf (SAL "       Prof Obj         Intensidad PIGN\n");
  guardar(\@prof_pign, \@pign, \@pign1, 4, $long_pign);                            
  printf (SAL "       Prof Obj         Intensidad SUWI\n");
  guardar(\@prof_suwi, \@suwi, \@suwi1, 5, $long_suwi);                            
  printf (SAL "       Prof Obj         Intensidad KINT\n");
  guardar(\@prof_kint, \@kint, \@kint1, 6, $long_kint);                            
  printf (SAL "       Prof Obj         Intensidad R35\n");
  guardar(\@prof_r35, \@r35, \@r351, 7, $long_r35);
  }
else
  {
  print("El ancho de banda debe ser mayor a 0\n"); 
  }

############################ SUBRUTINA QUE HACE LA SUAVIZACION PARA UNA PROFUNDIDAD OBJETIVO K #######################

sub suavizar # el ($$@@) es opcional
   {
   #$prof_obj    = $_[0]; 
   #$b           = $_[1];
   #@$prof_val   = $_[2];
   #@$valor      = $_[3];
   #$long        = $_[4];
   #($prof_obj, $b, @$prof_val, @$valor, $long) = (shift, shift, shift, shift); # Es equivalente a lo de arriba
   ($prof_obj, $b, @$prof_val, @$valor, $long) = @_; #se puede sustituir por la linea arriba
   local ($u, $j2, $denominador, $numerador);
   $denominador = 0;
   #$long_val = .scalar(@$prof_val); # como se guarda el valor de .scalar(@$prof_var).
   printf($long);  # No esta leyendo la longitud del vector @$prof_val
     for ($j2 = 0; $j2 < ($long + 1); ++$j2) 
        {
        $u = ($prof_obj - $prof_val[$j2])/$b;
        $denominador += kernel($u);
        print("j = $j2  ");
        print("Denominador = $denominador  ");
        print("Kernel =", kernel($u)); print("\n");
        }
     $numerador = 0;
     for ($j2 = 0; $j2 <= $long; ++$j2)
        {
        $u = ($prof_obj - $prof_val[$j2])/$b;
        $numerador += $valor[$j2]* kernel($u); 
        print("j1 = $j2");
        print("numerador = $numerador\n");
        }
     return ($numerador / $denominador);
   }

############################# SUBRUTINA QUE EJECUTA LA FUNCION KERNEL ############################################

sub kernel  # El ($u) es opcional
   {
   $u = $_[0]; 
   $a = 0.399;
   return ($a * exp (-0.5 * $u * $u ));
   }

############################# SUBRUTINA QUE GUARDA LAS SUAVIZACIONES ######################

sub guardar 
   {
   (@$prof_val, @$valor, @$valor1, $p, $long) = ($_[0], $_[1], $_[2], $_[3]); 
   local $j3  = 0;
   for ($prof_obj = 622; $prof_obj <= 2590; $prof_obj += .5)
      { 
      $intens[$j3] = suavizar ($prof_obj, $b, \@prof_val, \@valor, $long); 
      printf (SAL "%15.4f    %15.4f    %15.4f\n", $prof_obj, $valor1[$j3[$p]], $intens[$j3]);  
      ++$j3;
      }
   }   
############################## CLAUSURA DEL ARCHIVO DE SALIDA ############################################

close (SAL);

