pro cal_sun34, image, sky, quiet_sun, Range

;, frequency = frequency

;if n_elements(frequency) lt 0 then 
frequency = 34

imagemax=max(image, min=imagemin)

if n_elements(Range) le 0 then Range=2.e4

data=(smooth(float(image),3)-imagemin)/(imagemax-imagemin)*Range

hist=histogram(temporary(data), omin=omin,omax=omax)

N=n_elements(hist)

Nsm=fix(N/200. > 5 < 50)

h2=smooth(median(float(temporary(hist)),3), Nsm)

amax=max(h2)
area=total(h2)/amax

WM=wgtmax(h2)

Napp=[WM-3*area > 0, WM+3*area < (N-1)]

if n_elements(frequency) gt 0 then if frequency eq 34 then $
	Napp = [0, n_elements(h2)-1]

h2extr=h2(Napp(0):Napp(1))
Nextr=Napp(1)-Napp(0)+1
Npoints=128.
Npoints=512.*2

N0=16
;N0 = Npoints/8
for_approx=h2extr(findgen(Npoints)*Nextr/Npoints)

spec=fft(for_approx-mean(for_approx), -1)

spec(N0:Npoints-N0-1)=0

filtered=float(fft(spec, 1))

peaks=find_peaks(filtered)

peaks=(reverse(peaks(sort(filtered(peaks)))))([0,1])

peaks=peaks(sort(peaks))

amin=min(filtered(peaks(0):peaks(1)), imin)

border=peaks(0)+imin

peaks0=Napp(0)+peaks*Nextr/Npoints
;border0=Napp(0)+border*Nextr/Npoints

;thres=h2(border0)

;if thres lt amax*0.2 then thres = amax*0.3

;h3=h2*(h2 ge thres)

;sky=Napp(0)+wgtmax(h3(Napp(0):border0))+omin
;quiet_sun=border0+wgtmax(h3(border0:Napp(1)))+omin

sky = peaks0(0) + omin
quiet_sun = peaks0(1) +omin

sky=imagemin+sky/Range*(imagemax-imagemin)
quiet_sun=imagemin+quiet_sun/Range*(imagemax-imagemin)

	CASE frequency OF

17: T0 = 1e4
34: T0 = 1e4

	ELSE: T0 = 1.6e4
	ENDCASE
;stop
image = image - sky
quiet_sun = quiet_sun - sky
;;stop
image=image/quiet_sun*T0

;image=(image-sky)/quiet_sun*T0

	end