;program to calculate Tb in ALMA bands
pro calc_sp_cs3d_B_all_RL_24_all, path, filenum, x, y, ALL=ALL, ConFun=Confun

common tem;, tb
common sp, sp_tb, sp_tb_X, sp_tb_O, sp_tb_R, sp_tb_L, sp_P
common cfxy

nlambda=24
path='/data/slam/OSLO_simulations/processed/'
flstr=strcompress(string(filenum), /remove_all)
fl='cb24bih_'+flstr+'.sav'

path2='/data/slam/OSLO_simulations/processed/'+flstr+'/'
file_mkdir,path2

restore,'/data/slam/OSLO_simulations/cb24bih_mesh.idlsave'
height=reverse(z)*1e3
file='/data/slam/OSLO_simulations/en024048_hion/atmos/BIFROST_en024048_hion_bz_'+flstr+'.fits'
print,file
b_z=mrdfits(file, 0,hdr)
file='/data/slam/OSLO_simulations/en024048_hion/atmos/BIFROST_en024048_hion_bx_'+flstr+'.fits'
b_x=mrdfits(file, 0,hdr1)
file='/data/slam/OSLO_simulations/en024048_hion/atmos/BIFROST_en024048_hion_by_'+flstr+'.fits'
b_y=mrdfits(file, 0,hdr2)

sub=where((height ge 0) and (height le 5000),count)
h=height[sub]
;print, h(0),h(250)
b_x=b_x(*,*,sub)
b_y=b_y(*,*,sub)
b_z=b_z(*,*,sub)
;print, b_z(0,0,0),b_z(0,0,250)
b=sqrt(b_x^2+b_y^2+b_z^2)
b=reverse(b*1e4,3)
b_z=reverse(b_z*1e4,3)

fileinp=path+fl
print, fileinp
restore, fileinp
s=size(temp)




IF KEYWORD_SET(ALL) THEN BEGIN

 nelem=504*504.
 sp_tb=fltarr(504, 504, nlambda)
 sp_tb(*,*,*)=0.
 sp_tb_X=fltarr(504, 504, nlambda)
 sp_tb_X(*,*,*)=0.
 sp_tb_O=fltarr(504, 504, nlambda)
 sp_tb_O(*,*,*)=0.
 sp_tb_R=fltarr(504, 504, nlambda)
 sp_tb_R(*,*,*)=0.
 sp_tb_L=fltarr(504, 504, nlambda)
 sp_tb_L(*,*,*)=0.

 sp_P=fltarr(504, 504, nlambda)
 sp_P(*,*,*)=0.


confun_R=fltarr(504, nlambda, 251)
confun_R(*,*,*)=0.
confun_L=fltarr(504, nlambda, 251)
confun_L(*,*,*)=0.

; for i=0, s(1)-1 do $
; for j=0,s(2)-1 do $

  for i=0,s(1)-1 do $
   begin
    confun_r(*,*,*)=0.
    confun_l(*,*,*)=0.
    for j=0,s(2)-1 do $
     begin
      IF KEYWORD_SET(ConFun) THEN $
       TB_CS3D_B_vect_all_RL_24, path2, i, j, height, reform(temp(i,j,*)), reform(nne(i,j,*)), reform(nhi(i,j,*)), reform(np(i,j,*)), reform(b(i,j,*)), reform(b_z(i,j,*)), 10000,  /ConFun $
        Else TB_CS3D_B_vect_all_RL_24, path2, i, j, height, reform(temp(i,j,*)), reform(nne(i,j,*)), reform(nhi(i,j,*)), reform(np(i,j,*)), reform(b(i,j,*)), reform(b_z(i,j,*)),10000
     sp_tb(i,j,*)=tb(1,*)
  
     sp_tb_X(i,j,*)=T_XO(0,*)
     sp_tb_O(i,j,*)=T_XO(1,*)
     confun_R(j,*,*)=cf_R(*,*)
     confun_L(j,*,*)=cf_L(*,*)
     sp_tb_R(i,j,*)=T_RL(0,*)
     sp_tb_L(i,j,*)=T_RL(1,*)

     sp_P(i,j,*)=P(*)
     end
     IF KEYWORD_SET(ConFun) THEN save, confun_R, confun_L, filename=path2+'cf_'+flstr+'_'+strtrim(string(i),2)+'_all_RL.sav'
    end
   save, sp_tb, filename=path2+'sp_tb_'+flstr+'_all.sav'
   save, sp_tb_X, sp_tb_O, sp_P, filename=path2+'sp_tb_XO_'+flstr+'_all.sav'
   save, sp_tb_R, sp_tb_L, sp_P, filename=path2+'sp_tb_RL_'+flstr+'_all.sav'


END $
ELSE begin
       i=x
       j=y
       IF KEYWORD_SET(ConFun) THEN $
       begin
       TB_CS3D_B_vect_all_RL_24, path+fl, i, j, height, reform(temp(i,j,*)), reform(nne(i,j,*)), reform(nhi(i,j,*)), reform(np(i,j,*)), reform(b(i,j,*)), reform(b_z(i,j,*)),10000, /ConFun
       save, cf_x, cf_o, cf_r, cf_l, filename=path2+'cf_'+fl+'_'+strcompress(string(i),/remove_all)+'_'+strcompress(string(j),/remove_all)+'_all.sav'
       end Else TB_CS3D_B_vect_all_RL_24, path+fl, i, j, height, reform(temp(i,j,*)), reform(nne(i,j,*)), reform(nhi(i,j,*)), reform(np(i,j,*)), reform(b(i,j,*)), reform(b_z(i,j,*)),10000
       save, tb, t_xo, t_rl, p, filename=path2+'tb_'+fl+'_'+strcompress(string(i),/remove_all)+'_'+strcompress(string(j),/remove_all)+'_all.sav'
     end

end;program