      Program PRSONDA
C                                                                      *
C     Programa para el proceso de los registros brutos de sonda sonica
C facilitados en formato LAS.
C     El fichero de entrada contiene los datos correspondientes a los
C a los tres receptores de la sonda. Primero se dan los del receptor
C cercano, luego los del receptor lejano y finalmente los del receptor
C intermedio.
C     Cada registro contiene 513 valores. El primer valor es la
C profundidad. Los siguientes 512 valores corresponden a la amplitud
C registrada cada 4 microsegundos.
C     El programa recupera los tres registros correspondientes a cada
C profundidad y realiza las siguientes operaciones:                    *
C         (1) Filtrado "high pass" para eliminar el DC offset de los
C             registros,
C         (2) Determinar los FAT-First Arrival Time de cada registro 
C             y las amplitudes correspondientes,
C         (3) Obtener los valores de Vp asociados a los FAT, por 
C             diferencia entre los FAT, 
C         (4) Obtener los parametros alpha de atenuacion con la 
C             distancia,
C         (5) Determinar las amplitudes de los picos 2 y 3 del registro
C             del receptor cercano y calcular el parametro de forma
C             (1-P2/P3)
C         (6) Calcular el cociente de Poisson por correlacion con la
C             velocidad Vp y el factor de forma
C         (7) Obtener la velocidad Vs a partir de Vp y el cociente de
C             Poisson: Vs=Vp SQRT((1-2nu)/(2-2nu))
C         (8) Escribir en el fichero de output los resultados a la
C             profundidad de los registros
C                                                                      *
C              Programado por F. Beltran. Belgar Ingenieros.  
C                          2019.02.14
C              Cambio en la matriz de factores de forma, segun
C  resultados de ESTEYCO. F. Beltran. Belgar Ingenieros.
C                          2019.02.28
C
C **********************************************************************
C                          Declaracion de variables - Apertura ficheros*
C **********************************************************************
C
      Implicit none
C                   Maximo numero de profundidades con registro sonico *
      Integer maxprof
      Parameter (maxprof=14000)
C                   Parametro del filtro high-pass (frecuencia corte)  *
      Real fc, pi
      Parameter (fc=1000.0, pi=3.1415927) 
      Real rc, alpha
C
      Real dist, dt, amax, pct
      Real prof
      Real recep1(513,maxprof),recep2(513,maxprof),recep3(513,maxprof)
      Real rec1hp(513,maxprof),rec2hp(513,maxprof),rec3hp(513,maxprof)
C                                                                      *
      Real Vp1, Vp2, Vp3, Vp, a1, a2, a3, nu1,nu2,nu3, Vs1,Vs2,Vs3
      Real P2r1, tP2r1, P3r1, tP3r1
      Real P2r2, tP2r2, P3r2, tP3r2
      Real P2r3, tP2r3, P3r3, tP3r3
      Real Poisson
C
      Integer  i, iprof, j, k, icount, icnu
      Integer  iflgP2, iflgP3, iflbad
      Character linea*80
C                             Distancia entre receptores               *
C                             Delta t de los registros                 *
C                             Amplitud maxima de los registros         *
C                 Fraccion de la amplitud maxima umbral del ruido      *
      Parameter (dist=0.20, dt=4e-6, amax=32768.0, pct=0.05)
C
      Open (25,file='input',status='unknown',form='formatted')
      Open (26,file='output',status='unknown',form='formatted')
C
C **********************************************************************
C                    Lectura de los registros de la sonda sonica       *
C **********************************************************************
C
C                                                  Encabezado
      Do i=1,11
         read (25,'(a80)') linea
         write (26,'(a80)') linea
      End do
C                                             Receptor cercano         *
C                                             Busca comienzo registros
  100 Continue
         Read (25,'(a80)') linea
         If (linea(2:9).ne.'LOG_DATA') goto 100
C                                                                      *
      iprof= 0
  200 Continue
         Read (25,'(a80)') linea
         If (linea(2:9).ne.'--------') then
            Backspace 25
            iprof=iprof+1
            if (iprof.gt.maxprof) STOP 
     *              '***** Excedido el numero maximo de profundidades'
            Read (25,*) (recep1(i,iprof),i=1,513)
            Write (*,*) ' Leyendo receptor cercano, profundidad ',
     *                                              recep1(1,iprof)
         else
            Goto 300
         end if
C
         Goto 200
C                                             Receptor lejano
C                                             Busca comienzo registros
  300 Continue
         Read (25,'(a80)') linea
         If (linea(2:9).ne.'LOG_DATA') goto 300
C                                                                      *
      iprof= 0
  400 Continue
         Read (25,'(a80)') linea
         If (linea(2:9).ne.'--------') then
            Backspace 25
            iprof=iprof+1
            if (iprof.gt.maxprof) STOP 
     *              '***** Excedido el numero maximo de profundidades'
            Read (25,*) (recep3(i,iprof),i=1,513)
            Write (*,*) ' Leyendo receptor lejano, profundidad ',
     *                                              recep3(1,iprof)
         else
            Goto 500
         end if
C
         Goto 400
C                                             Receptor intermedio
C                                             Busca comienzo registros
  500 Continue
         Read (25,'(a80)') linea
         If (linea(2:9).ne.'LOG_DATA') goto 500
C                                                                      *
      iprof= 0
  600 Continue
         iprof=iprof+1
         if (iprof.gt.maxprof) STOP 
     *              '***** Excedido el numero maximo de profundidades'
         Read (25,*,END=700) (recep2(i,iprof),i=1,513)
         Write (*,*) ' Leyendo receptor intermedio, profundidad ',
     *                                              recep2(1,iprof)
         Goto 600
C
  700 Continue
C
      iprof=iprof-1
C
C ******************************************************************** *
C                      Filtro high-pass para eliminar el DC off-set    *
C ******************************************************************** *
C
      rc=1.0/(2.0*pi*fc)
      alpha=rc/(rc+dt)
C                                                                      *
      do i=1,iprof
C                                      copia la profundidad
         rec1hp(1,i)=recep1(1,i)
         rec2hp(1,i)=recep2(1,i)
         rec3hp(1,i)=recep3(1,i)
C                                     pone a cero el primer punto      *
         rec1hp(2,i)=0.0
         rec2hp(2,i)=0.0
         rec3hp(2,i)=0.0
C                                     ejecuta el filtro
         do j=3,513
            rec1hp(j,i)=alpha*(rec1hp(j-1,i)+recep1(j,i)-recep1(j-1,i))
            rec2hp(j,i)=alpha*(rec2hp(j-1,i)+recep2(j,i)-recep2(j-1,i))
            rec3hp(j,i)=alpha*(rec3hp(j-1,i)+recep3(j,i)-recep3(j-1,i))
         end do
      end do
C
C **********************************************************************
C                     Proceso principal para cada profundidad          *
C **********************************************************************
C 
      write (26,1000)
 1000 Format (//'Prof  P2r1 tP2r1 P3r1 tP3r1   P2r2 tP2r2 P3r2 tP3r2 ',
     *          ' P2r3 tP2r3 P3r3 tP3r3  Vp1 a1  Vp2 a2  Vp3 a3 ',
     *          ' Vpmed  fformr1  fformr2  fformr3  nu1 nu2 nu3 ',
     *          ' Vs1  Vs2  Vs3   Vsmed'//)
C
C                      Inicializa contador de registros defectuosos
      icount=0
      icnu=0
C      
      do i=1, iprof
C      
         Write (*,*) ' Procesando registros profundidad ',
     *                                              rec1hp(1,i)
C         
C                   ****************************************************
C                              Determinar los picos 2 y 3 de cada
C                              registro y los tiempos en que ocurren   *
C                       (los registros nulos vienen dados por -99999.0)
C
      if (recep1(2,i).gt.-99999.0) then
C
        iflbad=0
C                      Receptor cercano
C                                             Inicializa
        iflgP2=0
        iflgP3=0
        P2r1=-amax
        P3r1=0.0
C                   Salta 20 primeros puntos (blanking window)       
        do j=22,513
C                                                 Pico 3 (negativo)
           if (iflgP3.eq.0) then
               if (rec1hp(j,i).lt.P3r1) then
                  P3r1=rec1hp(j,i)
                  if ((P3r1.lt.rec1hp(j+1,i)).and.
     *                          (abs(P3r1).gt.(pct*amax))) then
                    iflgP3=1
                    tP3r1=dt*(j-2)
                  end if
               end if
           end if
C
           if ((iflgP3.eq.1) .and. (iflgP2.eq.0)) then
C                       Vuelta atras para calcular Pico 2 (positivo)   *
             k=0
 5001        continue
             k=k+1
             if ((k+1).eq.j) then
                 iflbad=1
                 iflgP2=1
             else if (rec1hp(j-k,i).gt.P2r1) then
                  P2r1=rec1hp(j-k,i)
                  if (P2r1.gt.rec1hp(j-(k+1),i)) then
                     iflgP2=1
                     tP2r1=dt*(j-k-2)
                  end if
             end if
             if (iflgP2.eq.0) go to 5001
           end if
C                                                                      *
        end do
C                                Comprobacion 
C                      Receptor intermedio                      ****** *
C                                             Inicializa
        iflgP2=0
        iflgP3=0
        P2r2=-amax
        P3r2=0.0
C                   Salta 20 primeros puntos (blanking window)          
        do j=22,513
C                                                 Pico 3 (negativo)
           if (iflgP3.eq.0) then
               if (rec2hp(j,i).lt.P3r2) then
                  P3r2=rec2hp(j,i)
                  if ((P3r2.lt.rec2hp(j+1,i)).and.
     *                          (abs(P3r2).gt.(pct*amax))) then
                    iflgP3=1
                    tP3r2=dt*(j-2)
                  end if
               end if
           end if
C
           if ((iflgP3.eq.1) .and. (iflgP2.eq.0)) then
C                       Vuelta atras para calcular Pico 2 (positivo)   *
             k=0
 5002        continue
             k=k+1
             if ((k+1).eq.j) then
                 iflbad=1
                 iflgP2=1
             else if (rec2hp(j-k,i).gt.P2r2) then
                  P2r2=rec2hp(j-k,i)
                  if (P2r2.gt.rec2hp(j-(k+1),i)) then
                     iflgP2=1
                     tP2r2=dt*(j-k-2)
                  end if
             end if
             if (iflgP2.eq.0) go to 5002
           end if
C
        end do
C                                Comprobacion 
C                      Receptor lejano                          ****** *
C                                             Inicializa
        iflgP2=0
        iflgP3=0
        P2r3=-amax
        P3r3=0.0
C                   Salta 20 primeros puntos (blanking window)          
        do j=22,513
C                                                 Pico 3 (negativo)
           if (iflgP3.eq.0) then
               if (rec3hp(j,i).lt.P3r3) then
                  P3r3=rec3hp(j,i)
                  if ((P3r3.lt.rec3hp(j+1,i)).and.
     *                          (abs(P3r3).gt.(pct*amax))) then
                    iflgP3=1
                    tP3r3=dt*(j-2)
                  end if
               end if
           end if
C
           if ((iflgP3.eq.1) .and. (iflgP2.eq.0)) then
C                       Vuelta atras para calcular Pico 2 (positivo)   *
               k=0
 5003          continue
               k=k+1
               if ((k+1).ge.j) then
                  iflbad=1
                  iflgP2=1
               else if (rec3hp(j-k,i).gt.P2r3) then
                  P2r3=rec3hp(j-k,i)
                  if (P2r3.gt.rec3hp(j-(k+1),i)) then
                     iflgP2=1
                     tP2r3=dt*(j-k-2)
                  end if
               end if
               if (iflgP2.eq.0) go to 5003
           end if
C
        end do
C                                Comprobacion 
C
C                   ****************************************************
C                      Comprobaciones de calidad y calculos de Vp y Vs
C                                                                      *
        if (abs(tP2r2-tP2r1).le.(4.0*dt)) then
C                                            receptor 1 y 2 repetidos
           icount=icount+1
        else if ((P3r1/P3r2).lt.(0.95)) then
C               amplificacion con la distancia - pico quedo recortado
           icount=icount+1
        else if (iflbad .eq. 1) then
C               Pico 2 muy cerca del origen en alguno de los receptores
           icount=icount+1
        else if ((P2r1.eq.(-amax)).or.(P2r2.eq.(-amax))) then
C    *                .or.(P2r3.eq.(-amax))) then
C               Pico 2 no identificado en alguno de los receptores       
           icount=icount+1
        else
C                      Registros validos - Calculos
           Vp1=dist/(tP2r2-tP2r1)
C                                                                      *
           a1=LOG(P3r1/P3r2)/dist
C                                                                      *
           nu1=Poisson(Vp1,1.0-(P2r1/P3r1))
C                     Poisson = 0.0 si interpolacion fuera de rango
           Vs1=Vp1*sqrt((1-2.0*nu1)/(2.0-2.0*nu1)) 
C                                                Poisson dentro de rango 
           if ((nu1.gt.0.01)) then
     *                          
C
          write (26,1010) rec1hp(1,i),P2r1,tP2r1,P3r1,tP3r1,P2r2,tP2r2,
     *           P3r2,tP3r2,P2r3,tP2r3,P3r3,tP3r3,Vp1,a1,Vp2,a2,Vp3,a3,
     *           (Vp1+Vp2+Vp3)/3, 
     *           1.0-(P2r1/P3r1),1.0-(P2r2/P3r2),1.0-(P2r3/P3r3),
     *           nu1,nu2,nu3, Vs1,Vs2,Vs3, (Vs1+Vs2+Vs3)/3
 1010     format (30(E12.5,1X))
           else
              icnu=icnu+1
           end if
C                                                             **********
C                                     Final condicional calidad
        end if
C                                     Final condicional sobre registros
C                                     nulos (-99999)                   *
      end if
C                                     Final bucle sobre profundidades  
      end do
C
      write (*,*) ' Registros defectuosos ', icount
      write (*,*) ' Poisson fuera de rango de interpolacion ', icnu
      write (*,*)
C      
      stop '********* Terminacion normal'
C                                                                      *
      end

      Real function Poisson (Vp,ff)
C
C     Funcion que interpola el cociente de Poisson a partir del factor *
C de forma ff obtenido de los registros de sonda sonica y de la Vp
C velocidad de las ondas P obtenida tambien de las "primeras llegadas" *
C observadas en los registros de los receptores cercano, intermedio y
C lejano.
C            Programado por F.Beltran. Belgar Ingenieros.
C                   2019.02.18
C                   2019.02.28                                         *
C
C **********************************************************************
C
      Implicit none
C
      Real Vp, ff
C
      Integer nvp, nnu
      Parameter (nvp=14,nnu=13)
      Real Vps(nvp), nus(nnu), Fforms(nvp,nnu)
      Real finterp(nnu)
      Integer i, j
C                       Datos tabla de interpolacion                   *
      Data Vps /
     *          1500.0, 1600.0, 1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 
     *          2200.0, 2300.0, 2400.0, 2500.0, 2600.0, 3000.0, 9000.0/
      Data nus /
     *          0.100, 0.125, 0.150, 0.175, 0.200, 0.225, 0.250, 
     *          0.275, 0.300, 0.325, 0.350, 0.375, 0.400 /

      Data Fforms / 
     * 0.5774, 0.6571, 0.7368, 0.8165, 0.8741, 0.8877, 0.9416, 0.9369, 
     *         0.9346, 0.9315, 0.9287, 0.9263, 0.9201, 0.9551,
     * 0.6838, 0.7501, 0.8163, 0.8826, 0.9325, 0.9503, 0.9850, 0.9880, 
     *         0.9903, 0.9906, 0.9899, 0.9885, 0.9796, 0.9850,
     * 0.7903, 0.8431, 0.8959, 0.9487, 0.9910, 1.0129, 1.0284, 1.0391, 
     *         1.0460, 1.0497, 1.0511, 1.0507, 1.0391, 1.0149,
     * 0.9002, 0.9362, 0.9721, 1.0148, 1.0495, 1.0755, 1.0953, 1.1100, 
     *         1.1200, 1.1261, 1.1290, 1.1297, 1.1198, 1.0946,
     * 1.0299, 1.0494, 1.0690, 1.0960, 1.1232, 1.1473, 1.1675, 1.1831, 
     *         1.1940, 1.2005, 1.2037, 1.2046, 1.1959, 1.1735,
     * 1.1481, 1.1563, 1.1645, 1.1787, 1.1977, 1.2178, 1.2361, 1.2506, 
     *         1.2604, 1.2661, 1.2686, 1.2689, 1.2599, 1.2403,
     * 1.2588, 1.2547, 1.2506, 1.2550, 1.2672, 1.2829, 1.2983, 1.3105, 
     *         1.3183, 1.3223, 1.3234, 1.3227, 1.3119, 1.2938,
     * 1.3648, 1.3450, 1.3252, 1.3222, 1.3301, 1.3423, 1.3543, 1.3636, 
     *         1.3690, 1.3710, 1.3706, 1.3685, 1.3549, 1.3370,
     * 1.4041, 1.4117, 1.3812, 1.3780, 1.3861, 1.3967, 1.4061, 1.4127, 
     *         1.4158, 1.4159, 1.4138, 1.4105, 1.3930, 1.3733,
     * 1.4435, 1.4301, 1.4167, 1.4245, 1.4373, 1.4483, 1.4563, 1.4611, 
     *         1.4628, 1.4617, 1.4583, 1.4536, 1.4311, 1.4074,
     * 1.4571, 1.4601, 1.4632, 1.4768, 1.4914, 1.5027, 1.5105, 1.5153, 
     *         1.5171, 1.5157, 1.5115, 1.5054, 1.4761, 1.4447,
     * 1.5857, 1.5585, 1.5313, 1.5319, 1.5441, 1.5596, 1.5741, 1.5850, 
     *         1.5909, 1.5915, 1.5873, 1.5801, 1.5400, 1.4951,
     * 1.7142, 1.6568, 1.5994, 1.5437, 1.5812, 1.6241, 1.6633, 1.6934, 
     *         1.7120, 1.7197, 1.7185, 1.7110, 1.6531, 1.5803
     *      /
C
C ******************************************************************** *
C                       Comprobar que Vp esta en rango de datos        *
      If ((Vp.lt.Vps(1)).or.(Vp.gt.Vps(nvp))) then
         write (*,*) 'Aviso:  Vp fuera de rango en funcion Poisson'
         Poisson=0.0
      Else
C                       Interpolar factor forma linealmente en Vp      *
         do i=1,nvp-1
            if ((Vp.ge.Vps(i)).and.(Vp.le.Vps(i+1))) then
               do j=1,nnu
                  finterp(j)=Fforms(i,j)+(Fforms(i+1,j)-Fforms(i,j))*
     *                         (Vp-Vps(i))/(Vps(i+1)-Vps(i))
               end do
               go to 100
            end if        
         end do
 100     continue
C         
C                       Comprobar que FF esta en rango de datos        *
         if ((ff.lt.finterp(1)).or.(ff.gt.finterp(nnu))) then
            write (*,*) 'Aviso: ff fuera de rango en funcion Poisson'
            Poisson=0.0
         else
C                       Interpolar factor forma linealmente en nu      *
            do i=1,nnu-1
               if ((ff.ge.finterp(i)).and.(ff.le.finterp(i+1))) then
                  Poisson=nus(i)+(nus(i+1)-nus(i))*
     *                    (ff-finterp(i))/(finterp(i+1)-finterp(i))
                  go to 200
               end if
            end do
 200        continue
         end if
C
      End if
C                                                                      *
      Return
      End
