SUBROUTINE BISEC(F,X,MIN1A,ST1,EPS,MINSTEP,D2) IMPLICIT NONE DOUBLE PRECISION X,D2,ST1,EPS,MINSTEP,Q0,Q1,Q2,X0,X1,Q, & ST0,MIN1A,F LOGICAL TEST ST0=ST1 IF (ABS(ST0).LT.2.0*MINSTEP) THEN ST0=2.0*MINSTEP ENDIF X0=X Q0=F(X0) X1=X0+ST0 Q1=F(X1) IF (Q1.GT.Q0) THEN ST0=-1.0*ST0 Q2=Q1 X1=X0+ST0 Q1=F(X1) ENDIF IF(Q1.LT.Q0)THEN X0=X1 Q2=Q0 TEST=.FALSE. 10 Q0=Q1 X1=X0+ST0 Q1=F(X1) IF(Q1.LT.Q0)THEN X0=X1 ST0=2*ST0 ELSE TEST=.TRUE. ENDIF IF(.NOT.TEST)THEN GOTO 10 ENDIF ENDIF ST0=ST0*0.5 20 IF ((ABS(ST0).GE.MINSTEP).AND.(3.0*(Q1 + Q2 - 2.0*Q0).GT. & (ABS(Q1) + ABS(Q2) + ABS(Q0))*EPS)) THEN IF (Q1.GT.Q2) THEN ST0=-1.0*ST0 Q=Q1 Q1=Q2 Q2=Q ENDIF X1=X0+ST0 Q=F(X1) IF (Q.GT.Q0) THEN ST0=-1.0*ST0 Q1=Q2 Q2=Q X1=X0+ST0 Q=F(X1) ENDIF IF (Q.LT.Q0) THEN X0=X1 Q2=Q0 Q0=Q ELSE Q1=Q ENDIF ST0=ST0*0.5 GOTO 20 ENDIF MIN1A=Q0 X=X0 ST1=ST0 D2=(Q2+Q1-2.0*Q0)/(2.0*ST0)**2.0 END *********************************************************************** SUBROUTINE GETC(OMEGA,C) IMPLICIT NONE INTEGER I,J,N,M,NMAX DOUBLE PRECISION F,YFIT,TFIT,OMEGA,C PARAMETER(M=3,NMAX=10000) DIMENSION F(NMAX,M),YFIT(NMAX),TFIT(NMAX),C(M) COMMON /SUNE/T(200000),Y(200000),DATA(2) DOUBLE PRECISION T,Y INTEGER DATA * Prepare vectors and matrix N=DATA(2)-DATA(1)+1 DO 10,I=1,N YFIT(I)=Y(DATA(1)+I-1) TFIT(I)=T(DATA(1)+I-1) F(I,1)=SIN(OMEGA*TFIT(I)) F(I,2)=COS(OMEGA*TFIT(I)) F(I,3)=1.0 10 CONTINUE * Do a least square fit to obtain the vector c CALL LSF(F,C,YFIT,N,M,NMAX,M) RETURN END *********************************************************************** FUNCTION CHISQ(OMEGA) IMPLICIT NONE INTEGER I,J,N,M,NMAX,MMAX DOUBLE PRECISION CHISQ,F,YFIT,TFIT,OMEGA,C,SUM PARAMETER(M=3,NMAX=10000) DIMENSION F(NMAX,M),YFIT(NMAX),TFIT(NMAX),C(M) COMMON /SUNE/T(200000),Y(200000),DATA(2) DOUBLE PRECISION T,Y INTEGER DATA * Check whether OMEGA has a sensible value. (5s < T < 1000s) * If not, then put CHISQ = 1.0E37 IF ((OMEGA.LT.0.006283).OR.(OMEGA.GT.1.257)) THEN CHISQ=1.0E37 RETURN END IF * Prepare vectors and matrix N=DATA(2)-DATA(1)+1 DO 10,I=1,N YFIT(I)=Y(DATA(1)+I-1) TFIT(I)=T(DATA(1)+I-1) F(I,1)=SIN(OMEGA*TFIT(I)) F(I,2)=COS(OMEGA*TFIT(I)) F(I,3)=1.0 10 CONTINUE * Do a least square fit to obtain the vector c CALL LSF(F,C,YFIT,N,M,NMAX,M) * Calculate CHISQ CHISQ=0.0 SUM=0.0 DO 30,I=1,N DO 20,J=1,M SUM=SUM+C(J)*F(I,J) 20 CONTINUE CHISQ=CHISQ+(YFIT(I)-SUM)**2.0 SUM=0.0 30 CONTINUE RETURN END ********************************************************************** * * Subroutine for least square fit of the linear system Ax=b * ********************************************************************** SUBROUTINE LSF(A,X,B,N,M,NMAX,MMAX) IMPLICIT NONE INTEGER N,M,NMAX,MMAX,I,J,K DOUBLE PRECISION A,X,B,ATA,ATB,P,SUM DIMENSION A(NMAX,MMAX),X(MMAX),B(NMAX),ATA(MMAX,MMAX),ATB(MMAX), & P(MMAX) * Ax=b, x is found by solving AT*Ax=AT*b * where AT is A transposed * Calculate ATA, that is A(transposed)*A SUM=0.0 DO 30,I=1,M DO 20,J=I,M DO 10,K=1,N SUM=SUM+A(K,I)*A(K,J) 10 CONTINUE ATA(I,J)=SUM ATA(J,I)=SUM SUM=0.0 20 CONTINUE 30 CONTINUE * Calculate ATB, that is A(transposed)*B DO 50,I=1,M DO 40,K=1,N SUM=SUM+A(K,I)*B(K) 40 CONTINUE ATB(I)=SUM SUM=0.0 50 CONTINUE CALL CHOLDC(ATA,M,MMAX,P) CALL CHOLSL(ATA,M,MMAX,P,ATB,X) RETURN END ************************************************************************ * * Subroutine for Cholesky decomposition. * Matrix should be symmetric and positive definite. * ************************************************************************ SUBROUTINE choldc(a,n,np,p) INTEGER n,np DOUBLE PRECISION a(np,np),p(np) INTEGER i,j,k DOUBLE PRECISION sum do 13 i=1,n do 12 j=i,n sum=a(i,j) do 11 k=i-1,1,-1 sum=sum-a(i,k)*a(j,k) 11 continue if(i.eq.j)then if(sum.le.0.)pause 'choldc failed' p(i)=sqrt(sum) else a(j,i)=sum/p(i) endif 12 continue 13 continue return END *********************************************************************** * * Subroutine for Cholesky backsubstitution * *********************************************************************** SUBROUTINE cholsl(a,n,np,p,b,x) INTEGER n,np DOUBLE PRECISION a(np,np),b(n),p(n),x(n) INTEGER i,k REAL sum do 12 i=1,n sum=b(i) do 11 k=i-1,1,-1 sum=sum-a(i,k)*x(k) 11 continue x(i)=sum/p(i) 12 continue do 14 i=n,1,-1 sum=x(i) do 13 k=i+1,n sum=sum-a(k,i)*x(k) 13 continue x(i)=sum/p(i) 14 continue return END