FUNCTION _sin2, x, p, nx, np ; p0 - background ; p1 - amplitude ; p2 - period ; p3 - phase p0 = transpose(rebin(reform(p[0, *]), np, nx)) p1 = transpose(rebin(reform(p[1, *]), np, nx)) p2 = transpose(rebin(reform(p[2, *]), np, nx)) p3 = transpose(rebin(reform(p[3, *]), np, nx)) f = p0+p1*sin(2*!pi*x/p2+p3) return, f END FUNCTION _sin2_chi2, p, np, funa = funa nx = n_elements(funa.x) xx = rebin(funa.x, nx, np) yy = rebin(funa.y, nx, np) ee = rebin(funa.e, nx, np) f = _sin2(xx, p, nx, np) chi = total((yy-f)^2/ee^2, 1) return, chi END col0 = 180 th = 5 prange = [2e-2, 1] ; DATA nx = 15 p = [37.2, 7.15, 0.333, 0.5*!pi] ndim = n_elements(p) x = randomu(seed, nx) y = _sin2(x, p, nx, ndim) e = 1.*sqrt(abs(y))*randomn(seed, nx) y = y+e e = abs(e) ; FUNA funa = {x:x, y:y, e:e} ; MODEL nm = 1000 xm = findgen(nm)/nm ym = _sin2(xm, p, nm, 1) ; PERIODOGRAM lnp = lnp_test(x, y, wk1 = freq, wk2 = amp, ofac = 16, hifac = 32) ; SOLBER sol = solber('_sin2_chi2', ndim, funa = funa, npop = 300, lim = [[1e-6, 100], [1e-6, 10], [0, 1], [0, 2*!pi]], plot_flag = 2, ngen_max = 500, term_fit = 11, status = status, ngen_tot = ngen_tot, gfit_best = gfit_best) snapshot = TVRD(True = 1) write_png, 'solber_sin2_stat.png', snapshot ; PLOT window, 1, retain = 2, xs = 300, ys = 400 device, decomposed = 0 !p.multi = [0, 1, 3] !y.margin = [3, 2] !p.charsize = 1.6 ; usym a = findgen(17) * (!PI*2/16.) usersym, cos(a), sin(a), /FILL ; plot data plot, [0], [0], background = 255, color = 0, xr = [0, 1], yr = minmax([y-e, y+e]), title = 'Model+Data' oplot, xm, ym, color = col0, thick = th oplot, x, y, psym = 8, color = 0 errplot, x, y-e, y+e, color = 0 ; plot periodogram yr = [0, 1.1*max(amp)] plot, [0], [0], xr = prange, xs = 1, color = 0, yr = yr, ys = 1, title = 'Periodogram' oplot, p[2]*[1, 1], yr, color = col0, thick = th oplot, 1./freq, amp, color = 0 ; plot solber solution oplot, sol[2]*[1, 1], yr, color =0 ; compare solutions plot, [0], [0], background = 255, color = 0, xr = [0, 1], yr = minmax([y-e, y+e]), title = 'Model+Data vs Solber' oplot, xm, ym, color = col0, thick = 5 oplot, x, y, psym = 8, color = 0 errplot, x, y-e, y+e, color = 0 ysol = _sin2(xm, sol, nm, 1) oplot, xm, ysol, color = 0 !p.multi = 0 snapshot = TVRD(True = 1) write_png, 'solber_sin2.png', snapshot finish: END