##// END OF EJS Templates
schain
RSS
add fitacf_guess.f
rflores -
r1612:3e57003ca20e
parent child
Show More
Side by Side Unified
Context file:
r1612:3e57003ca20e Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,88
1 subroutine guess(acf,tau,npts,zero,amin,te,tr)
2 c
3 c find zero crossing (zero), depth of minimum (amin), height of maximum
4 c
5
6 real acf(npts),tau(npts)
7 c write(*,*) 'acf: ',acf
8 c write(*,*) 'tau: ',tau
9 c read(*,*) xx
10
11 zero=0.0
12 amin=1.0
13 tmin=0.0
14 jmin=0
15
16 do i=npts,2,-1
17 if(acf(i)*acf(i-1).lt.0.0) then
18 zero=(tau(i-1)*acf(i)-tau(i)*acf(i-1))/(acf(i)-acf(i-1))
19 end if
20 if(acf(i).lt.amin) then
21 amin=acf(i)
22 jmin=i
23 end if
24 end do
25
26 if(jmin.gt.0) then
27 call parab1(tau(jmin-1),acf(jmin-1),a,b,c)
28 tmin=-b/(2.0*a)
29 amin=c+tmin*(b+tmin*a)
30 end if
31
32 tr=cdtr1(-amin)
33 te=czte1(zero*1000.0,tr)
34 c write(*,*) 'zero: ',zero
35 c write(*,*) 'amin: ',amin
36 c write(*,*) 'te: ',te
37 c write(*,*) 'tr: ',tr
38 c read(*,*) xx
39
40 return
41 end
42
43 subroutine parab1(x,y,a,b,c)
44 C-----
45 dimension x(3),y(3)
46 delta=x(1)-x(2)
47 a=(y(1)-2.*y(2)+y(3))/(2.*delta*delta)
48 b=(y(1)-y(2))/delta - a*(x(1)+x(2))
49 c=y(1)-a*x(1)*x(1)-b*x(1)
50 return
51 end
52
53 real function cdtr1(depth)
54 C-----convert depth to te/ti ratio
55 dimension tr(4)
56 C according to the curve published in farley et al 1967
57 c modified for 2004 conditions on axis
58 c data tr/7.31081,3.53286,5.92271,.174/
59 data tr/9.5,4.0,8.5,.3/
60 data nt/4/
61 cdtr1=tr(1)
62 do i=2,nt
63 cdtr1=cdtr1*depth + tr(i)
64 end do
65 return
66 end
67
68 real function czte1(zlag,tr)
69 C-----convert zero crossing point to te
70 C according to the curve published in farley et al 1967
71 c modified for 2004 conditions on axis
72 dimension dt(4)
73 c data dt/0.00945025,-0.0774338,.203626,.812397/,nd/4/
74 data dt/0.00945025,-0.0774338,.2,0.9/,nd/4/
75 data t0/1000./
76 tr1=min(abs(tr),5.)
77 if(zlag .eq. 0)then
78 czte1=1000000.
79 else
80 dt0=dt(1)
81 do i=2,nd
82 dt0=dt0*tr1 + dt(i)
83 end do
84 czte1=t0*(dt0/zlag)**2
85 end if
86 return
87 end
88
Show file before | Show file after | Hide comments
@@ -114,7 +114,8 setup(name='schainpy',
114 extra_f77_compile_args=["-fallow-argument-mismatch"]),
114 extra_f77_compile_args=["-fallow-argument-mismatch"]),
115 Extension("schainpy.model.proc.fitacf_guess",
115 Extension("schainpy.model.proc.fitacf_guess",
116 sources=[
116 sources=[
117 "schainf/Ffiles/fitacf_guess.pyf"],
117 "schainf/Ffiles/fitacf_guess.pyf",
118 "schainf/Ffiles/fitacf_guess.f",],
118 extra_f77_compile_args=["-fallow-argument-mismatch"]),
119 extra_f77_compile_args=["-fallow-argument-mismatch"]),
119 Extension("schainpy.model.proc.fitacf_acf2",
120 Extension("schainpy.model.proc.fitacf_acf2",
120 sources = [
121 sources = [
@@ -227,6 +228,7 setup(name='schainpy',
227 "schainf/Ffiles/lmdif1.f",
228 "schainf/Ffiles/lmdif1.f",
228 "schainf/Ffiles/reader.c",
229 "schainf/Ffiles/reader.c",
229 "schainf/Ffiles/cbesi.f",
230 "schainf/Ffiles/cbesi.f",
231 "schainf/Ffiles/lagp.f",
230 "schainf/Ffiles/i1mach.f",
232 "schainf/Ffiles/i1mach.f",
231 "schainf/Ffiles/zeta.f",
233 "schainf/Ffiles/zeta.f",
232 "schainf/Ffiles/qc25f.f",
234 "schainf/Ffiles/qc25f.f",
General Comments 0
You need to be logged in to leave comments. Login now