PRO lfit_2, y, s2, x, nch, ax, fx chmn = 1.0e20 nch1 = nch-1 FOR n = 1, nch-3 DO BEGIN ;at least 2 points below or above ebr n1 = n+1 y1 = y(0:n) ;below break energy x1 = x(0:n) s21 = s2(0:n) lfit, y1, s21, x1, ax1 ;fit a power law y2 = y(n1:nch1) ;above break energy x2 = x(n1:nch1) s22 = s2(n1:nch1) lfit, y2, s22, x2, ax2 ;fit a power law dgamma = -ax2(1)+ax1(1) ebr = (ax2(0)-ax1(0))/dgamma ;the break energy ;will be restricted ;to x(n)<ebr<x(n1) IF(ebr LT x(n)) THEN ebr = x(n) ELSE $ ;reset ebr IF(ebr GT x(n1)) THEN ebr = x(n1) ax2(0) = ax1(0)+dgamma*ebr ;reset k2 ys1 = ax1(0)+ax1(1)*x1 ;values of y from fit ys2 = ax2(0)+ax2(1)*x2 ys = [ys1, ys2] ;get chi^2 ch2 = total((ys-y)^2/s2) IF (ch2 LT chmn) THEN BEGIN chmn = ch2 ax = [ax1, ax2] ;[ax1(0),ax1(1),ax2(0),ax2(1)] ;[k1,gm1,k2,gm2] fx = ys ENDIF ENDFOR RETURN END