fitacf_guess.f
88 lines
| 2.1 KiB
| text/x-fortran
|
FortranFixedLexer
| r1612 | 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 | ||||
