PROGRAM SLIDEFIT IMPLICIT NONE DOUBLE PRECISION PAR,YMODEL,OMEGA,D2,ST1,EPS,MINSTEP,MIN1A,C,PI, & YM,SD,SUMX,SUMXX,X,BP,TI1,TI2,MEAN,RSD INTEGER NPAR,I,J,K,N,NDATA,NFIT,NINT,WINDOW,RES,INP,DAT,MAXINT, & NTI1,NTI2 LOGICAL CONV CHARACTER DATAFILE*15 PARAMETER (NPAR=6,PI=3.14159265358979323846,RES=11,MAXINT=100) DIMENSION PAR(NPAR),C(3),BP(MAXINT),TI1(MAXINT),TI2(MAXINT) COMMON /SUNE/T(200000),Y(200000),DATA(2) DOUBLE PRECISION T,Y INTEGER DATA EXTERNAL CHISQ NDATA=200000 NINT=MAXINT EPS=1.0E-6 MINSTEP=1.0E-16 CALL SLAES(NDATA,T,Y,NINT,WINDOW,BP,TI1,TI2,WINDOW) OPEN(UNIT=RES,FILE='slidefit.out',STATUS='UNKNOWN') WRITE(RES,'(A)') & '# Bif. par. A^2 s.d.' * maybe the intervals given should be checked before starting... * OK, now do a sliding window fit for each of the NINT intervals specified DO 200,J=1,NINT DO 10,K=1,NDATA IF (T(K).GE.TI1(J)) THEN NTI1=K GOTO 20 END IF 10 CONTINUE 20 DO 30,K=NTI1,NDATA IF (T(K).GE.TI2(J)) THEN NTI2=K GOTO 40 END IF 30 CONTINUE 40 NFIT=NTI2-NTI1-WINDOW WRITE(*,*) & 'Now doing',NFIT,' fits from t =',TI1(J),' to t =',TI2(J) SUMX=0.0 SUMXX=0.0 * This loop moves a sliding window and calculates the period, the square * of the amplitude and the standard deviations of these. DO 100,I=NTI1,NTI1+NFIT DATA(1)=I DATA(2)=DATA(1)+WINDOW OMEGA=0.18 ST1=0.001 CALL BISEC(CHISQ,OMEGA,MIN1A,ST1,EPS,MINSTEP,D2) CALL GETC(OMEGA,C) IF (C(1).EQ.0.0) THEN PAR(1)=C(2) PAR(3)=ASIN(C(2)/PAR(1)) ELSE PAR(3)=ATAN2(C(2),C(1)) PAR(1)=C(2)/SIN(PAR(3)) END IF PAR(2)=OMEGA PAR(6)=C(3) * X=2.0*PI/OMEGA X=PAR(1)**2.0 SUMX=SUMX+X SUMXX=SUMXX+X**2.0 100 CONTINUE N=NFIT+1 MEAN=SUMX/N RSD=SD/MEAN SD=SQRT((SUMXX-(SUMX**2.0)/N) / (N-1)) WRITE(RES,'(3(XG15.8))') BP(J),MEAN,SD WRITE(*,'(3(XG15.8))') BP(J),MEAN,SD 200 CONTINUE CLOSE(RES) END ************************************************************************* DOUBLE PRECISION FUNCTION YMODEL(T,P) IMPLICIT NONE DOUBLE PRECISION T,P DIMENSION P(6) YMODEL=P(1)*SIN(P(2)*T+P(3))+P(6) * YMODEL=P(1)*SIN(P(2)*T+P(3))+P(4)*SIN(2.0*P(2)*T+P(5)) + P(6) RETURN END ************************************************************************* SUBROUTINE GNU(T1,T2) IMPLICIT NONE DOUBLE PRECISION T1,T2 INTEGER IUNIT PARAMETER (IUNIT=10) OPEN(UNIT=IUNIT,FILE='gnufile',STATUS='UNKNOWN',FORM='FORMATTED') WRITE(IUNIT,'(A)') 'set term X11' WRITE(IUNIT,'(A,F8.1,A,F8.1,A)') 'set xrange[',T1,':',T2,']' WRITE(IUNIT,'(2A)') "plot 'fitplt' title 'data' w l,", & "'fitplt' u 1:3 title 'fit' w l" WRITE(IUNIT,'(A)') 'pause -1 "Press to continue."' CLOSE(IUNIT) CALL SYSTEM("gnuplot gnufile") OPEN(UNIT=IUNIT,FILE='gnufile',STATUS='OLD',FORM='FORMATTED') CLOSE(IUNIT,STATUS='DELETE') RETURN END