C C The Kwee-Woerden's method. C C References: C C K K Kwee, H van Woerden : A Method for Computing accurately C the epoch of minimum of an eclipsing variable, BAN, Vol 12, No. 464, C p 327 (1956) C C * F.Hroch, hroch@physics.muni.cz, October 1996, October 1997, 1999 Feb., C SUBROUTINE KW(METODA,MINMAX,POCET,T,M,T0,ER,POCSUM,KOD) C C Interface to the Fast KW method. C C INPUT PARAMETERS: C C METODA. method for estimation of equidistant C data from raw data. Only two value (size sensitive!): C 'L'..linear interpolation C 'S'..spline interpolation (used procedures Spline and Seval C from Forsythe, Malcom, Moler:Computer methods for C mathematical computations, Prentice-Hall Inc., 1977, C Packages FMM on netlib library. C C MINMAX. find the maximum or the minimum? C 'X'..maximum C 'N'..minimum (size sensitive!) C C POCET..number of points. C C T..array of input x-coordinates (time) (sorted by decreasing order!) C C M..array of input y-coordinates (magnitude) C C OUTPUT PARAMETERS: C C T0..Estimated time of the minimum C C ERR..standart deviation of T0 C C KOD.. exit code: C 10..METODA isn't 'L' or 'S' C 13..Minimum is nearly limit of interval or no extreme found C 14..Number of interations during finding extreme was exceeded. C Time of exterme is probably wrong. C 31..Error during evaluating of linear interpolation C (probably any points at same time or points not sort) C 41..MINMAX isn't 'X' or 'N' C 99..wrong declared dimensions of MAXDIM. C INTEGER MAXDIM PARAMETER( MAXDIM = 1000 ) CHARACTER METODA,MINMAX INTEGER POCET,KOD,J,MIN0,POCSUM REAL DT,T0,ER,T(POCET),M(POCET),E(MAXDIM),F(MAXDIM),G(MAXDIM) REAL M1(MAXDIM),T1(MAXDIM) REAL SEVAL,INTERPOL INTEGER ODHAD C C Test for correct dimensions of declared arrays. IF( POCET.GT.MAXDIM )THEN KOD = 99 RETURN ENDIF KOD=0 C Times's step. DT=(T(POCET)-T(1))/(POCET-1.0) C Interpolation. IF(METODA.EQ.'S')THEN CALL SPLINE(POCET,T,M,E,F,G) DO J=1,POCET T1(J)=(J-1)*DT+T(1) M1(J)=SEVAL(POCET,T1(J),T,M,E,F,G) ENDDO ELSEIF(METODA.EQ.'L')THEN DO J=1,POCET T1(J)=(J-1)*DT+T(1) M1(J)=INTERPOL(T1(J),J,POCET,T,M,KOD) IF(KOD.NE.0) THEN KOD=31 RETURN ENDIF ENDDO ELSE KOD=10 RETURN ENDIF C First estimation of time of extrem. IF((MINMAX.EQ.'X').OR.(MINMAX.EQ.'N'))THEN MIN0 = ODHAD(MINMAX,POCET,M) ELSE KOD=41 RETURN ENDIF CALL KWMIN(MIN0,POCET,T1,M1,T0,ER,DT,POCSUM,KOD) END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE KWMIN(IMIN0,POCET,T,M,T0,ER,DT,POCSUM,KOD) C C The core of this package. The KW method. C C INTEGER I,IMIN,IMIN0,POCET,KOD,N,L,POCSUM REAL T(POCET),M(POCET),T0,ER DOUBLE PRECISION SUM(-1:1),A,B,CAS0 C IMIN = IMIN0 DO I = 1,POCET C ************************************************************************** C P = \min( K-1 , N-K ) - 1 C*************************************************************************** POCSUM = MIN(IMIN-1,POCET-IMIN) - 1 IF(POCSUM.LT.2)THEN KOD=13 RETURN ENDIF DO N=-1,1 SUM(N)=0.0 ENDDO C*************************************************************************** C S_{-} = \sum_{i=1}^P ( M_{K-1+i} - M_{K-1-i} )^2 C S_{0} = \sum_{i=1}^P ( M_{K+i} - M_{K-i} )^2 C S_{+} = \sum_{i=1}^P ( M_{K+1+i} - M_{K+1-i} )^2 C*************************************************************************** DO N=-1,1 DO L=1,POCSUM SUM(N) = SUM(N) + (M(IMIN+N+L) - M(IMIN+N-L))**2 ENDDO ENDDO IF( (SUM(-1).GE.SUM(0)).AND.(SUM(1).GE.SUM(0)) ) GOTO 20 C Direction of iterations. IF( SUM(-1).GT.SUM(1) ) IMIN = IMIN+1 IF( SUM(-1).LT.SUM(1) ) IMIN = IMIN-1 ENDDO 20 CONTINUE IF( I.GE.POCET ) KOD = 14 C CAS0=T(IMIN) C*************************************************************************** C A = \frac{ S_{+} + S_{-} - 2 S_{0} }{ 2 D^2 } C B = \frac{ S_{0} - (S_{-} + S_{+})/2 }{D} - A ( 2 T_0 - D ) C*************************************************************************** A=(SUM(-1) + SUM(1) - 2.0*SUM(0))/(2.0*DT**2) B=(SUM(0) - 0.5*(SUM(-1)+SUM(1)))/DT - A*(2.0*CAS0 - DT) C*************************************************************************** C T_{min} = - \frac{B}{ 2 A } C \sigma^2 = \frac{ S_{0}/A - ( T_{min} - T_{0} )^2 }{ 2 P } C************************************************************************** T0=-B/2.0/A ER=(SUM(0)/A - (T0 - CAS0)**2)/2.0/POCSUM IF(ER.LT.0.0) ER=0.0 ER=SQRT(ER) write(*,*) sum,a,b,dt C END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL FUNCTION INTERPOL(U,KK,POCET,T,M,USPECH) C Code for linear interpolation in arrays C T (time coordination) and M (intensity, magnitude). C C U ... evaluated time (value will be returned INTERPOL's) C KK... initial estimate index's position U between T(i) C REAL U INTEGER K,KK,L,POCET,USPECH REAL T(POCET),M(POCET) C USPECH = 0 K = KK C Searching in tabulated points. 10 continue IF( (K.LT.1) .OR. (K.GE.POCET) ) GOTO 30 IF( T(K+1).GE.U .AND. U.GT.T(K) ) GOTO 20 IF( T(K+1)-T(K).LE.0.0 )THEN C Error. Times aren't sort or you are (machine) nearly. USPECH = 10 RETURN ENDIF IF( U.LE.T(K) )THEN K = K - 1 ELSEIF( U.GT.T(K+1) ) THEN K = K + 1 ENDIF GOTO 10 20 CONTINUE C Value U between boundary points T(1) and T(POCET). Linear interpolation. INTERPOL = M(K)+(U-T(K))*((M(K+1)-M(K))/(T(K+1)-T(K))) RETURN C Value U outside the interval. Linear extrapolations. 30 IF( K.LT.1 )THEN INTERPOL = M(1)+(U-T(1))*((M(2)-M(1))/(T(2)-T(1))) ELSEIF( K.GE.POCET )THEN L = POCET-1 INTERPOL = M(L)+(U-T(L))*((M(L+1)-M(L))/(T(L+1)-T(L))) ENDIF END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC INTEGER FUNCTION ODHAD(MINMAX,POCET,M) C First estimation of time of extreme. INTEGER I,J,POCET CHARACTER MINMAX REAL X,M(POCET) X = M(1) J = 1 IF( MINMAX.EQ.'X' )THEN DO I=2,POCET IF( M(I).GT.X )THEN X = M(I) J = I ENDIF ENDDO ELSEIF( MINMAX.EQ.'N' )THEN DO I=2,POCET IF( M(I).LT.X )THEN X = M(I) J = I ENDIF ENDDO ENDIF ODHAD = J END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC REAL function seval(n, u, x, y, b, c, d) integer n REAL u, x(n), y(n), b(n), c(n), d(n) c c this subroutine evaluates the cubic spline function c c seval = y(i) + b(i)*(u-x(i)) + c(i)*(u-x(i))**2 + d(i)*(u-x(i))**3 c c where x(i) .lt. u .lt. x(i+1), using horner's rule c c if u .lt. x(1) then i = 1 is used. c if u .ge. x(n) then i = n is used. c c input.. c c n = the number of data points c u = the abscissa at which the spline is to be evaluated c x,y = the arrays of data abscissas and ordinates c b,c,d = arrays of spline coefficients computed by spline c c if u is not in the same interval as the previous call, then a c binary search is performed to determine the proper interval. c integer i, j, k REAL dx data i/1/ if ( i .ge. n ) i = 1 if ( u .lt. x(i) ) go to 10 if ( u .le. x(i+1) ) go to 30 c c binary search c 10 i = 1 j = n+1 20 k = (i+j)/2 if ( u .lt. x(k) ) j = k if ( u .ge. x(k) ) i = k if ( j .gt. i+1 ) go to 20 c c evaluate spline c 30 dx = u - x(i) seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i))) return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine spline (n, x, y, b, c, d) integer n REAL x(n), y(n), b(n), c(n), d(n) c c the coefficients b(i), c(i), and d(i), i=1,2,...,n are computed c for a cubic interpolating spline c c s(x) = y(i) + b(i)*(x-x(i)) + c(i)*(x-x(i))**2 + d(i)*(x-x(i))**3 c c for x(i) .le. x .le. x(i+1) c c input.. c c n = the number of data points or knots (n.ge.2) c x = the abscissas of the knots in strictly increasing order c y = the ordinates of the knots c c output.. c c b, c, d = arrays of spline coefficients as defined above. c c using p to denote differentiation, c c y(i) = s(x(i)) c b(i) = sp(x(i)) c c(i) = spp(x(i))/2 c d(i) = sppp(x(i))/6 (derivative from the right) c c the accompanying function subprogram seval can be used c to evaluate the spline. c c integer nm1, ib, i REAL t c nm1 = n-1 if ( n .lt. 2 ) return if ( n .lt. 3 ) go to 50 c c set up tridiagonal system c c b = diagonal, d = offdiagonal, c = right hand side. c d(1) = x(2) - x(1) c(2) = (y(2) - y(1))/d(1) do 10 i = 2, nm1 d(i) = x(i+1) - x(i) b(i) = 2.*(d(i-1) + d(i)) c(i+1) = (y(i+1) - y(i))/d(i) c(i) = c(i+1) - c(i) 10 continue c c end conditions. third derivatives at x(1) and x(n) c obtained from divided differences c b(1) = -d(1) b(n) = -d(n-1) c(1) = 0. c(n) = 0. if ( n .eq. 3 ) go to 15 c(1) = c(3)/(x(4)-x(2)) - c(2)/(x(3)-x(1)) c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3)) c(1) = c(1)*d(1)**2/(x(4)-x(1)) c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3)) c c forward elimination c 15 do 20 i = 2, n t = d(i-1)/b(i-1) b(i) = b(i) - t*d(i-1) c(i) = c(i) - t*c(i-1) 20 continue c c back substitution c c(n) = c(n)/b(n) do 30 ib = 1, nm1 i = n-ib c(i) = (c(i) - d(i)*c(i+1))/b(i) 30 continue c c c(i) is now the sigma(i) of the text c c compute polynomial coefficients c b(n) = (y(n) - y(nm1))/d(nm1) + d(nm1)*(c(nm1) + 2.*c(n)) do 40 i = 1, nm1 b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.*c(i)) d(i) = (c(i+1) - c(i))/d(i) c(i) = 3.*c(i) 40 continue c(n) = 3.*c(n) d(n) = d(n-1) return c 50 b(1) = (y(2)-y(1))/(x(2)-x(1)) c(1) = 0. d(1) = 0. b(2) = b(1) c(2) = 0. d(2) = 0. return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC