1. program medianfilter ! median filter for a lightcurve implicit none integer, parameter :: n=200 real, dimension(n) :: t, x, f1, f2 real :: m integer i, stat, npts ! read input data (time, value) i=0 do read(*,*, iostat=stat) t(i), x(i) if (stat /= 0) exit i=i+1 end do npts = i-1 ! nonrecursive filter do i=2, npts-1 f1(i) = median(x(i-1), x(i), x(i+1)) end do ! recursive filter f2=0 do i=2, npts-1 f2(i) = median(f2(i-1), x(i), x(i+1)) end do do i=1,npts write(*,'(4F8.3)'), t(i), x(i),f1(i),f2(i) end do contains real function median (a, b, c) real, intent(in) :: a,b,c median = (a+b+c)-max(a,b,c)-min(a,b,c) end function end program 4.878 6.515 0.000 0.000 4.795 6.506 6.506 6.492 4.949 6.492 6.503 6.492 5.063 6.503 6.497 6.497 5.176 6.497 6.503 6.497 5.445 6.504 6.504 6.504 5.570 6.519 6.519 6.519 5.725 6.531 6.531 6.531 5.861 6.555 6.555 6.555 5.996 6.569 6.569 6.569 6.028 6.575 6.575 6.575 6.049 6.578 6.578 6.578 6.122 6.600 6.600 6.600 6.289 6.617 6.617 6.617 6.413 6.620 6.619 6.620 6.547 6.623 6.620 6.620 6.609 6.615 6.615 6.615 6.711 6.605 6.605 6.605 6.968 6.583 6.583 6.583 6.946 6.575 6.574 6.575 7.255 6.565 6.565 6.565 7.823 6.559 6.559 6.559 7.936 6.553 6.558 6.558 8.060 6.558 6.558 6.558 8.050 6.562 6.560 6.560 8.215 6.560 6.560 6.560 8.472 6.539 6.539 6.539 8.564 6.522 6.521 6.521 8.615 6.515 6.519 6.519 8.842 6.519 6.515 6.519 8.904 6.515 6.515 6.515 9.089 6.502 6.502 6.502 9.305 6.497 6.502 6.502 9.812 6.517 6.517 6.517 9.864 6.525 6.525 6.525 9.999 6.538 6.538 6.538 10.063 6.556 6.556 6.556 10.137 6.586 6.586 6.586 10.345 6.604 6.604 6.604 10.408 6.623 6.604 6.604 10.726 6.598 6.598 6.598 10.870 6.589 6.589 6.589 11.014 6.587 6.587 6.587 11.117 6.578 6.578 6.578 11.199 6.567 6.567 6.567 11.446 6.557 6.557 6.557 11.611 6.550 6.556 6.556 11.694 6.556 6.553 6.556 11.786 6.553 6.553 6.553 11.869 6.550 6.553 6.553 12.024 6.556 6.554 6.554 12.127 6.554 0.000 0.000 The autocorrelation program cannot be easily used for non-equispaced data. The lightcurve is relatively smooth, and we'll use interpolation to create a new data set with equal spacing. The same dataset can also be used in the fft program. program psyint implicit none integer, parameter :: maxpoint=1000 integer :: n real, dimension(maxpoint) :: x, y, xequ, yequ integer i, stat, k real xx, yy, r, t ! read the data i=0 do read(*,*, iostat=stat) xx, yy if (stat /= 0) exit i=i+1 x(i) = xx ; y(i) = yy end do n=i ! interpolate equally spaced data ! an ugly solution (constants set manually ! according to the data file) i=0 t=4.6 do xequ(i)=t do k=1,n if (t > x(k) .and. t<= x(k+1)) exit end do r = y(k)+(t-x(k))*(y(k+1)-y(k))/(x(k+1)-x(k)) yequ(i)=r i=i+1 t=t+0.1 if (t > x(n-1)) exit end do do k=1,i write(*,'(F6.2,F9.3)') xequ(k), yequ(k) end do end program The new dataset is: 4.70 6.53 4.80 6.52 4.90 6.50 5.00 6.50 5.10 6.50 5.20 6.50 5.30 6.50 5.40 6.50 5.50 6.51 5.60 6.52 5.70 6.53 5.80 6.54 5.90 6.56 6.00 6.57 6.10 6.59 6.20 6.61 6.30 6.62 6.40 6.62 6.50 6.62 6.60 6.62 6.70 6.61 6.80 6.60 6.90 6.59 7.00 6.57 7.10 6.57 7.20 6.57 7.30 6.56 7.40 6.56 7.50 6.56 7.60 6.56 7.70 6.56 7.80 6.56 7.90 6.56 8.00 6.56 8.10 6.56 8.20 6.56 8.30 6.55 8.40 6.54 8.50 6.53 8.60 6.52 8.70 6.52 8.80 6.52 8.90 6.52 9.00 6.51 9.10 6.50 9.20 6.50 9.30 6.50 9.40 6.50 9.50 6.50 9.60 6.51 9.70 6.51 9.80 6.52 9.90 6.53 10.00 6.54 10.10 6.57 10.20 6.59 10.30 6.60 10.40 6.62 10.50 6.62 10.60 6.61 10.70 6.60 10.80 6.59 10.90 6.59 11.00 6.59 11.10 6.58 11.20 6.57 11.30 6.56 11.40 6.56 11.50 6.55 11.60 6.55 11.70 6.56 11.80 6.55 11.90 6.55 12.00 6.55 2. program fftlc ! fourier transform of a a lightcurve use fft implicit none integer, parameter :: m=7, n=2**m complex, dimension(n) :: a real, dimension(n) :: t, f, absval integer :: npoints, i, stat real :: xx, yy a = 0.0 t = 0.0 ; f=0.0 i=0 do read(*,*, iostat=stat) xx, yy if (stat /= 0) exit i=i+1 t(i) = xx ; f(i) = yy end do npoints=i write (*,*) 'original data' do i=1,npoints write(*,*) t(i), f(i) end do a = cmplx(f,0.0) call fourier(a, m, n, 1) absval = abs(a) write (*,*) 'Fourier transform' do i=1,npoints write(*,'(I3,F8.4)') i,absval(i) end do end program Fourier transform 1 484.9500 2 259.0823 3 63.2320 4 67.0025 5 54.7323 6 18.3917 7 43.3309 8 6.2032 9 30.8395 10 18.1149 11 17.2354 12 22.1367 13 4.4343 14 20.7122 15 5.7103 16 15.5283 17 12.1860 18 8.3230 19 14.6768 20 0.7592 21 13.5795 22 5.6662 23 9.8125 24 9.7901 25 4.4301 26 11.2699 27 1.1036 28 10.0283 29 5.7438 30 6.7396 31 8.5834 32 2.2974 33 9.2419 34 2.2192 35 7.7917 36 5.8140 37 4.7005 38 7.8255 39 0.8108 40 7.8823 41 3.0592 42 6.2252 43 5.9494 44 3.2623 45 7.2752 46 0.3462 47 6.9190 48 3.7059 49 5.0477 50 6.0093 51 2.0440 52 6.8757 53 1.3174 54 6.1385 55 4.2674 56 4.0135 57 6.1319 58 0.9743 59 6.5693 60 2.2091 61 5.4485 62 4.8533 63 3.0711 64 6.3528 65 0.0100 66 6.3528 67 3.0711 68 4.8533 69 5.4485 70 2.2091 71 6.5693 72 0.9743 73 6.1319 74 4.0135 Mainly noise; the light curve is too short to find the period 3-4. program autoc use plot implicit none integer, parameter :: n=1000 real, dimension(n) :: x, f, P integer i, npoints call getdata(npoints, x, f) ! do i=1,npoints ! write(*,*) x(i), f(i) ! end do do i=1,npoints/2 x(i)=0.1*i P(i)=ac(f, npoints, i) write(*,'(2F10.5)') x(i),P(i) end do call init("xx.ps", 3.0, 10.0, cm, cm) call frame(0.0, 0.0, 10.0, 5.0) call plotpoints(x, f, npoints/2) call plotcurve(x, 1000*P, npoints/2, 0.5, 0) call finish() contains real function ac(f, n, dt) ! autocorrelation for delay dt ! f(1:n)=observed values integer, intent(in) :: n real, dimension(n), intent(in) :: f integer, intent(in) :: dt real :: s, mean mean = sum(f)/real(n) s = sum((f(1:n-dt)-mean)*(f(1+dt:n)-mean)) ac = s/(n-dt-1) end function end program 0.10000 0.00138 0.20000 0.00129 0.30000 0.00114 0.40000 0.00094 0.50000 0.00072 0.60000 0.00047 0.70000 0.00023 0.80000 -0.00002 0.90000 -0.00024 1.00000 -0.00043 1.10000 -0.00059 1.20000 -0.00070 1.30000 -0.00077 1.40000 -0.00081 1.50000 -0.00082 1.60000 -0.00080 1.70000 -0.00078 1.80000 -0.00074 1.90000 -0.00070 2.00000 -0.00066 2.10000 -0.00062 2.20000 -0.00060 2.30000 -0.00060 2.40000 -0.00060 2.50000 -0.00060 2.60000 -0.00059 2.70000 -0.00058 2.80000 -0.00055 2.90000 -0.00050 3.00000 -0.00042 3.10000 -0.00032 3.20000 -0.00018 3.30000 -0.00000 3.40000 0.00020 3.50000 0.00044 3.60000 0.00068 3.70000 0.00092 The lightcurve is too short for finding the true period with the autocorrelation method. At the end the values are increasing, and the next local maximum would indicate the period.