From 3e57003ca20ece95eb52c2a5590f04cb2c75f617 2023-06-02 19:29:43
From: Roberto Flores <rflores@igp.gob.pe>
Date: 2023-06-02 19:29:43
Subject: [PATCH] add fitacf_guess.f

---

diff --git a/schainf/Ffiles/fitacf_guess.f b/schainf/Ffiles/fitacf_guess.f
new file mode 100644
index 0000000..3ac16f5
--- /dev/null
+++ b/schainf/Ffiles/fitacf_guess.f
@@ -0,0 +1,88 @@
+      subroutine guess(acf,tau,npts,zero,amin,te,tr)
+c
+c     find zero crossing (zero), depth of minimum (amin), height of maximum
+c
+
+      real acf(npts),tau(npts)
+c      write(*,*) 'acf: ',acf
+c      write(*,*) 'tau: ',tau
+c      read(*,*) xx
+
+      zero=0.0
+      amin=1.0
+      tmin=0.0
+      jmin=0
+
+      do i=npts,2,-1
+         if(acf(i)*acf(i-1).lt.0.0) then
+            zero=(tau(i-1)*acf(i)-tau(i)*acf(i-1))/(acf(i)-acf(i-1))
+         end if
+         if(acf(i).lt.amin) then
+            amin=acf(i)
+            jmin=i
+         end if
+      end do
+
+      if(jmin.gt.0) then
+         call parab1(tau(jmin-1),acf(jmin-1),a,b,c)
+         tmin=-b/(2.0*a)
+         amin=c+tmin*(b+tmin*a)
+      end if
+      
+      tr=cdtr1(-amin)
+      te=czte1(zero*1000.0,tr)
+c      write(*,*) 'zero: ',zero
+c      write(*,*) 'amin: ',amin
+c      write(*,*) 'te: ',te
+c      write(*,*) 'tr: ',tr
+c      read(*,*) xx
+
+      return
+      end
+      
+      subroutine parab1(x,y,a,b,c)
+C-----
+      dimension x(3),y(3)
+      delta=x(1)-x(2)
+      a=(y(1)-2.*y(2)+y(3))/(2.*delta*delta)
+      b=(y(1)-y(2))/delta - a*(x(1)+x(2))
+      c=y(1)-a*x(1)*x(1)-b*x(1)
+      return
+      end
+
+      real function cdtr1(depth)
+C-----convert depth to te/ti ratio
+      dimension tr(4)
+C     according to the curve published in farley et al 1967
+c     modified for 2004 conditions on axis
+c      data tr/7.31081,3.53286,5.92271,.174/
+      data tr/9.5,4.0,8.5,.3/
+      data nt/4/
+      cdtr1=tr(1)
+      do i=2,nt
+        cdtr1=cdtr1*depth + tr(i)
+      end do
+      return
+      end
+
+      real function czte1(zlag,tr)
+C-----convert zero crossing point to te
+C     according to the curve published in farley et al 1967
+c     modified for 2004 conditions on axis
+      dimension dt(4)
+c      data dt/0.00945025,-0.0774338,.203626,.812397/,nd/4/
+      data dt/0.00945025,-0.0774338,.2,0.9/,nd/4/
+      data t0/1000./
+      tr1=min(abs(tr),5.)
+      if(zlag .eq. 0)then
+        czte1=1000000.
+      else
+        dt0=dt(1)
+        do  i=2,nd
+          dt0=dt0*tr1 + dt(i)
+        end do
+        czte1=t0*(dt0/zlag)**2
+      end if
+      return
+      end
+
diff --git a/setup.py b/setup.py
index e7de4c9..4dda05e 100644
--- a/setup.py
+++ b/setup.py
@@ -114,7 +114,8 @@ setup(name='schainpy',
             extra_f77_compile_args=["-fallow-argument-mismatch"]),
         Extension("schainpy.model.proc.fitacf_guess",
             sources=[
-                "schainf/Ffiles/fitacf_guess.pyf"],
+                "schainf/Ffiles/fitacf_guess.pyf",
+                "schainf/Ffiles/fitacf_guess.f",],
             extra_f77_compile_args=["-fallow-argument-mismatch"]),
         Extension("schainpy.model.proc.fitacf_acf2",
             sources = [
@@ -227,6 +228,7 @@ setup(name='schainpy',
                 "schainf/Ffiles/lmdif1.f",
                 "schainf/Ffiles/reader.c",
                 "schainf/Ffiles/cbesi.f",
+                "schainf/Ffiles/lagp.f",
                 "schainf/Ffiles/i1mach.f",
                 "schainf/Ffiles/zeta.f",
                 "schainf/Ffiles/qc25f.f",