path2='c:\!work\sswdb\CS3D\tomo\
;save, tbr, tbl, effhr, effhl, temp, height, freq, xl, filename=path2+'data_385_142_270.sav'
restore, path2+'data_385_142_270.sav'

s=n_elements(freq)
dh_all=fltarr(s)
slope_all=fltarr(s)
coef_all=fltarr(s)
sigma_all=fltarr(s)
dh_all(*)=0.
sigma_all(*)=0.
slope_all(*)=0.

tb_all=tbr
tt2=shift(tb_all,-1)
dt_all=tt2-tb_all
h_all=effhr
tt2=shift(h_all,-1)
deffh=tt2-h_all
teff_all=interpol(temp,height/1e5, h_all)
  
  
  
for l=0,5 do $
 begin
  ;Tb within 1 band
  tb_b3=tbr[4*l:4*l+3]
  tt2=shift(tb_b3,-1)
  dt=tt2-tb_b3
  dt=dt[0:2]
  ;frequency
  fr_b3=freq[4*l:4*l+3]/1e9
  tt2=shift(fr_b3,-1)
  df=tt2-fr_b3
  df=df[0:2]
  ;effective heights 
  h_b3=effhr[4*l:4*l+3]
  tt2=shift(h_b3,-1)
  deffh=tt2-h_b3
  deffh=deffh[0:2]
  teff_b3=interpol(temp,height/1e5, h_b3)
  
  n=3
  slope=fltarr(n)
  sigma=fltarr(n)
  coef=fltarr(n)
  for i=0,n-1 do $
   begin
    tb=tb_b3[i:i+1];+tbl[8+i:8+i+1])/2.
    x=alog(fr_b3[i:i+1])
    y=alog(tb)  ;check definitions of Tb in paper Wang et al. (2015)
    res=linfit(x,y,yfit=yfit,CHISQR=CHISQR, sigma=sigma0)
    slope(i)=res(1)
    coef(i)=res(0)
    sigma(i)=sigma0(1)
   end 

  hrange=effhr(4*l+3)-effhr(4*l) 
  h0=effhr(4*l)

  a=fltarr(n-1)
  ;slope=-slope
  for i=0, n-2 do a(i)=slope(i)*dt(i+1)/slope(i+1)/dt(i)

  dheight=get_dheight(dt, df, hrange, a)  
 ; nh_b3=fltarr(4)
 ; nh_b3[0]=h0
 ; for i=1,3 do nh_b3[i]=nh_b3[i-1]+dheight[i-1]
  
  dh_all[4*l:4*l+2]=dheight
  slope_all[4*l:4*l+2]=slope
  coef_all[4*l:4*l+2]=coef
  sigma_all[4*l:4*l+2]=sigma
end
;stop
;slopes and dheights between the bands
;;method1 - connect bands with a line and determine its slope
;for l=0,4 do $
;  begin
;   ind=[3+4*l,4+4*l]
;   tb=tbr[3+4*l:4+4*l];+tbl[8+i:8+i+1])/2.
;   x=alog(freq[3+4*l:4+4*l]/1e9)
;   y=alog(tb)  ;check definitions of Tb in paper Wang et al. (2015)
;   res=linfit(x,y,yfit=yfit,CHISQR=CHISQR, sigma=sigma0)
;   slope=res(1)
;   sigma=sigma0(1)
;   slope_all(3+4*l)=slope
;   sigma_all(3+4*l)=sigma
;   dh=slope_all[2+4*l]*dt_all[3+4*l]*dh_all[2+4*l]/slope/dt_all[2+4*l]
;   dh_all[3+4*l]=dh
;  end
  
 ;method2  - extrapolate fit in each band, find crosssection of the two fits and 2 corresponding dh 
for l=0,4 do $
  begin
   ind=[3+4*l,4+4*l]
   tb=tbr[3+4*l:4+4*l];+tbl[8+i:8+i+1])/2.
   x=alog(freq[3+4*l:4+4*l]/1e9)
   y=alog(tb)  ;check definitions of Tb in paper Wang et al. (2015)
   ;cross point
   x0=(coef_all(4+4*l)-coef_all(2+4*l))/(slope_all(2+4*l)-slope_all(4+4*l))
   y0=coef_all(4+4*l)+slope_all(4+4*l)*x0
   print, y0
   tb0=exp(y0)
   ;print, coef_all(2+4*l)+slope_all(2+4*l)*x0
   if ((x0 gt x(0)) and (x0 lt x(1))) then $
    begin
    dh1=dh_all(2+4*l)*(tb0-tb(0))/(dt_all(2+4*l))
    dh2=dh_all(4+4*l)*(tb(1)-tb0)/(dt_all(4+4*l))
    dh=dh1+dh2
    end else $
      begin
       res=linfit(x,y,yfit=yfit,CHISQR=CHISQR, sigma=sigma0)
       slope=res(1)
       sigma=sigma0(1)
       slope_all(3+4*l)=slope
       sigma_all(3+4*l)=sigma
       dh=slope_all[2+4*l]*dt_all[3+4*l]*dh_all[2+4*l]/slope/dt_all[2+4*l]
      end 
   dh_all[3+4*l]=dh
  end

nh_all=fltarr(s)
h0=effhr(0)
nh_all[0]=h0
for i=1,s-1 do nh_all[i]=nh_all[i-1]+dh_all[i-1]
  

stop  
  
  
  

window,0,xs=900,ys=400
TVLCT, rgb, /Get
TVLCT, Reverse(rgb,1)
!p.multi=[0,2,1]

plot, freq/1e9, tb_all/1e3, psym=-2, xtitle='freq, GHz', ytitle='Tb, 10^3K', /ystyle, yrange=[5, 13], /xstyle, xrange=[30,700]
plot, h_all, teff_all/1e3, psym=-2, xtitle='height, km', ytitle='Te, 10^3K',  yrange=[5, 13]
set_line_color
oplot, h_all, teff_all/1e3, psym=-2, color=5
oplot, nh_all, tb_all/1e3, psym=-2, color=3
xyouts, 1200, 5.4,/data, 'Te(heff)', color=5, charsize=1.5
xyouts, 1200, 4.4,/data, 'Tb(heff(0)+dh)', color=3, charsize=1.5

WRITE_PNG,path2+'temp_allbands_385_142_270_m2.png', TVRD(/TRUE)
 end