;program for bremsstrahlung 
;for known B(z)
PRO TB_CS3D_B_vect_all_RL_24, path, x,y, height, temp, nne, nhi, np, b, bz, NN, ConFun=ConFun, CFOnly=CFOnly

;input prameters are from corona down to photosphere

COMMON TEM;, TB, T_XO, P
COMMON GAUNT, GAUNT_ATM
COMMON COFU;, CF
common cfxy, CF_R, CF_L


;ALMA frequency bands + longer and shorter lambda  24 lambda

FREQ=[35, 37, 39, 41, 67, 71, 79, 83, 90.5, 92.5, 102.5, 104.5,$
      224, 226, 240, 242 ,336.5, 338.5, 348.5, 360.5, 632, 634, 636, 638]*1e9

XL=reform(transpose(3e10/freq)*1e8)  ;angstrem


SW=SIZE(XL)

sz=size(height)
height00=height
s=sz(1)

;INTERPOLATION OF ORIGINAL HEIGHT INTO FINER SCALE (NN=10000 WAS COMMONLY USED)
INTERP='Yes'
;INTERP='no'
nn=5000

TEMP_IN=DBLARR(NN)
NNE_IN=DBLARR(NN)
NP_IN=DBLARR(NN)
NHI_IN=DBLARR(NN)
b_in=DBLARR(NN)
bz_in=DBLARR(NN)
HEIGHT_IN=DBLARR(NN)
h2=DBLARR(NN)

TAU=DBLARR(NN,SW(1))
STEP=DBLARR(NN)
STEP1=DBLARR(NN)

NNE_UP=DBLARR(NN)

;stop

STEP(0)=0
;;;;;;;; HEIGHT SCALE IS IN CM!!!!!!!!!!!!!!!!!!
IF INTERP EQ 'Yes' THEN $
   BEGIN

    STEP_IN=ABS(HEIGHT(0)-HEIGHT(S-1))/(NN-1)
    FOR I=0,NN-1 DO HEIGHT_IN(I)=HEIGHT(0)-I*STEP_IN
    TEMP_IN = INTERPOL(TEMP, HEIGHT, HEIGHT_IN)
    NNE_IN = INTERPOL(NNE, HEIGHT, HEIGHT_IN)
    NP_IN = INTERPOL(NP, HEIGHT, HEIGHT_IN)
    NHI_IN = INTERPOL(NHI, HEIGHT, HEIGHT_IN)
    b_in = INTERPOL(b, HEIGHT, HEIGHT_IN)
    bz_in = INTERPOL(bz, HEIGHT, HEIGHT_IN)
    STEP(*)=STEP_IN(*);*1E5
   END $
 ELSE $
   BEGIN
    ;step array btw original layers of the atmosphere
    NN=S
    TEMP_IN=DBLARR(NN)
    NNE_IN=DBLARR(NN)
    NP_IN=DBLARR(NN)
    NHI_IN=DBLARR(NN)
    b_in=DBLARR(NN)
    bz_in=DBLARR(NN)
    HEIGHT_IN=DBLARR(NN)
    h2=DBLARR(NN)

    TAU=DBLARR(NN,SW(1))
    STEP=DBLARR(NN)
    STEP1=DBLARR(NN)

    NNE_UP=DBLARR(NN)
    FOR I=1,NN-1 DO STEP(I)=(HEIGHT(I-1)-HEIGHT(I));*1E5
    ;PRINT,'STEP:',STEP
    TEMP_IN(*)=TEMP(*)
    NNE_IN(*)=NNE(*)
    NP_IN(*)=NP(*)
    NHI_IN(*)=NHI(*)
    b_in(*)=b(*)
    bz_in(*)=bz(*)
    HEIGHT_IN(*)=HEIGHT(*)
    step_in=step

   END

;stop
TTB=FLTARR(4,SW(1))

INTERPOLATE_gaunt,XL,TEMP_IN
GAUNT_ATM_R=reverse(GAUNT_ATM,2)

;opacity as function of height and wavelength
TAU_X=DBLARR(NN,SW(1))
TAU_X(*,*)=0
TAU_O=DBLARR(NN,SW(1))
TAU_O(*,*)=0
TAU_X_R=DBLARR(NN,SW(1))
TAU_O_R=DBLARR(NN,SW(1))
TAU_R=DBLARR(NN,SW(1))
TAU_R(*,*)=0
TAU_L=DBLARR(NN,SW(1))
TAU_L(*,*)=0
TAU_R_R=DBLARR(NN,SW(1))
TAU_L_R=DBLARR(NN,SW(1))


;contribution functions as function of height and wavelength

CF_R=DBLARR(NN,SW(1))
CF_R(*,*)=0
CF_L=DBLARR(NN,SW(1))
CF_L(*,*)=0
cf_R_in=DBLARR(n_elements(height00),SW(1))
cf_L_in=DBLARR(n_elements(height00),SW(1))


RADIUS=6.96E10 ;cM
photo=where(abs(height_in) eq min(abs(height_in)))  ;determine photosphere layer

;in CGS units
C=3E10 ;LIGHT VELOCITY [CM/SEC]
E= 4.8032068E-10  ;ELECTRON CHARGE [ESU]
M=9.1093897E-28 ;GRAMM
PI=3.14159265
k = 1.380658E-16;


coefB=E/2./PI/M/C ;2.8e6
COEF=E^2/M*E^2/sqrt(M)*E^2/k/sqrt(k)/c*8./3./sqrt(2.)/sqrt(!pi)
coef_ei=4*sqrt(2.*!pi)/3.*e/k*e/sqrt(k)*e^2/sqrt(m)


  TAU_X(*,*)=0
  TAU_O(*,*)=0
  TAU_X_R(*,*)=0
  TAU_O_R(*,*)=0
  TAU_R(*,*)=0
  TAU_L(*,*)=0
  TAU_R_R(*,*)=0
  TAU_L_R(*,*)=0


  TB_X=0
  TB_O=0
  TB_R=0
  TB_L=0


  For K=0,SW(1)-1 Do $  ;wavelength cycle
   Begin
   CFreq= 3E10/(XL(K)/1E8)
   CFreq=CFreq^2

   NNE_IN_R=reverse(NNE_IN)
   TEMP_IN_R=reverse(TEMP_IN)
   NP_IN_R=reverse(NP_IN)
   NHI_IN_R=reverse(NHI_IN)
   STEP_IN_R=reverse(STEP_IN)
   B_IN_R=reverse(B_IN)

   NHe_IN_R=0.0846*(NHI_IN_R+NP_IN_R)

   Bz_IN_R=reverse(bz_IN)
   HEIGHT_IN_R=reverse(HEIGHT_IN)
   STEP_R=reverse(STEP)

    LFreq = 8.97E3*Sqrt(NNE_IN_R)
    LFreq=LFreq^2

    EI_OPACITY=fltarr(2,nn)
    EA_OPACITY=fltarr(2,nn)
    EI_OPACITY(*,*)=0
    EA_OPACITY(*,*)=0
    NNE_UP_R=NNE_IN_R
    res=where(LFreq GE CFreq, count)
    if count ne 0 then begin  NNE_UP_R(res)=0 &
                              TAU_X_R(res,K)=0 ; {wave with omega < lengmur omega can not propagate}
                              TAU_O_R(res,K)=0
                              TAU_R_R(res,K)=0 ; {wave with omega < lengmur omega can not propagate}
                              TAU_L_R(res,K)=0
                           END

    ;*****************************************************************************

    NUB_R=COEFB*B_IN_R    ;FREQUENCY IN HZ  2.8^1E6*B
    CS_IN_R=BZ_IN_R/B_IN_R
    sub_Bneg=where(Bz_in_R lt 0.)
    ;stop
       ;
    ;*****************************************************************************
    ;TAU CALCULATION
    ;IF NNE_UP(j) NE 0 and j le photo THEN BEGIN
   ; IF NNE_UP(j) NE 0 THEN BEGIN
    res=where(NNE_UP_R NE 0, count, complement=res2, ncomplement=count2)
     ; EI_OPACITY = OPACITY_EI(XL(K),NP_IN(I),NNe_UP(I),TEMP_IN(I),STEP(I))
    ; EI_OPACITY = OPACITY_EI_3(XL(K),NP_IN(I),NNe_IN(I),TEMP_IN(I),STEP(I))
     ;EI_OPACITY = OPACITY_EI_G(XL(K),NP_IN(j),NNe_IN(j),TEMP_IN(j),STEP(j),GAUNT_ATM(K,j))
    ;;;;;;;;;;;; EI_OPACITY = OPACITY_EI_G_B(XL(K),NP_IN(j),NNe_IN(j),TEMP_IN(j),STEP(j),GAUNT_ATM(K,j), NUB, CS)
    ; print, XL(K),NNe_IN(j)/1E9,TEMP_IN(j)/1E3,STEP(j)/1E9,NUB, ALPHA
     ;stop
    ;EI_OPACITY = OPACITY_EI_ZLOTNIK(XL(K),NNe_IN(j)/1E9,TEMP_IN(j)/1E3,STEP(j)/1E9,NUB, ALPHA)
    ;EI_OPACITY = OPACITY_EI_WHITE(XL(K),NNe_IN(j),TEMP_IN(j),STEP(j),NUB, ALPHA)

    if count ne 0 then $
    begin
    ;EI_OPACITY[*,res] = OPACITY_EI_G_B_vect(XL(K),NP_IN_R[res],NNe_IN_R[res],TEMP_IN_R[res],STEP_R[res],GAUNT_ATM_R(K,res), NUB_R[res], CS_IN_R[res])
    FREQ_EI=EFFFREQ_EI(Nne_IN_R[res],TEMP_IN_R[res], GAUNT_ATM_R(K,res), coef_ei)
    ;FREQ_EI=EFFFREQ_EI_LOG(Nne_IN_R[res],TEMP_IN_R[res], 3E18/XL(K), coef_ei)
    EI_OPACITY[*,res] = OPACITY_TOT_B_ZH_VECT(XL(K),NNe_IN_R[res],TEMP_IN_R[res],STEP_R[res],NUB_R[res], CS_IN_R[res], FREQ_EI, 0)
    ;EI_OPACITY[*,res] = OPACITY_TOT_B_ZH_VECT(XL(K),NNe_IN_R[res],TEMP_IN_R[res],STEP_R[res],NUB_R[res], CS_IN_R[res], FREQ_EI, 0,/isotropy)

    ;EA_OPACITY[*,res] = OPACITY_EA_Q_FL_VECT(XL(K),NNe_IN_R[res],NHI_IN_R[res],TEMP_IN_R[res],STEP_R[res],NUB_R[res], CS_IN_R[res])
    ;HM + HeM
    EA_OPACITY[*,res] = OPACITY_EA_Q_FL_VECT_TOT(XL(K),NNe_IN_R[res],NHI_IN_R[res],NHe_IN_R[res],TEMP_IN_R[res],STEP_R[res],NUB_R[res], CS_IN_R[res],/nohe)

    ;EA_OPACITY[*,res] = OPACITY_EA_Q_B_vect(XL(K),NNe_UP_R[res],NHI_IN_R[res],TEMP_IN_R[res],STEP_R[res], NUB_R[res], CS_R[res]);
;
    ;OPTICAL DEPTH INCLUDES BOTH OPACITY SOURCES - ELECTRON-ION AND ELECTRON-NEUTRAL ATOM
    ;X and O modes
    ;EI_OPACITY(*,*)=0
    Tau_X_R(*,K) = EI_OPACITY(0,*) + EA_OPACITY(0,*)
    Tau_O_R(*,K) = EI_OPACITY(1,*) + EA_OPACITY(1,*)

    ;R and L circular polarized emission
    Tau_R_R(*,K) = EI_OPACITY(0,*) + EA_OPACITY(0,*)  ;x-mode
    ;R is O and L is X where B<0
    Tau_R_R(sub_Bneg,K) = EI_OPACITY(1,sub_Bneg) + EA_OPACITY(1,sub_Bneg)
    Tau_L_R(*,K) = EI_OPACITY(1,*) + EA_OPACITY(1,*)  ;o-mode
    Tau_L_R(sub_Bneg,K) = EI_OPACITY(0,sub_Bneg) + EA_OPACITY(0,sub_Bneg)

   END ;internal cycle

  TB_X=0
  TB_O=0
  TB_R=0
  TB_L=0

;stop

 HEIGHT_IN=reverse(HEIGHT_IN_R)
 TEMP_IN=reverse(TEMP_IN_R)
 NNE_IN=reverse(NNE_IN_R)
 TAU_X(*,K)=reverse(TAU_X_R(*,K),1)
 TAU_O(*,K)=reverse(TAU_O_R(*,K),1)
 TAU_R(*,K)=reverse(TAU_R_R(*,K),1)
 TAU_L(*,K)=reverse(TAU_L_R(*,K),1)


  ;BRIGHTNESS CALCULATION FOR 2 WAVE MODES
  TbCalc,XL(K),NN,TAU_X(*,K),HEIGHT_IN,TEMP_IN,NNE_IN,C,path,1,'X';,/file ;$
  TB_X=TB


  TbCalc,XL(K),NN,TAU_O(*,K),HEIGHT_IN,TEMP_IN,NNE_IN,C,path,1,'O';,/file ;$
  TB_O=TB


  ;BRIGHTNESS CALCULATION FOR 2 polarizations - RIGHT and LEFT
  TbCalc,XL(K),NN,TAU_R(*,K),HEIGHT_IN,TEMP_IN,NNE_IN,C,path,1,'R';,/file ;$
  TB_R=TB
  CF_R(*,K)=CF   ;CF_X[height,lambda]

  TbCalc,XL(K),NN,TAU_L(*,K),HEIGHT_IN,TEMP_IN,NNE_IN,C,path,1,'L';,/file ;$
  TB_L=TB
  CF_L(*,K)=CF


  TTB(*, K)=[TB_X,TB_O,TB_R,TB_L]
  Freq = 30*1E8/XL(K);{ nu=c/E9/lambda  frequency is in GHz}

 END ;end of lambda cycle


 T_tot=total(TTB(0:1,*),1)/2.  ;TOTAL BRIGHTNESS
 P=(TTB(2,*)-TTB(3,*))/(TTB(2,*)+TTB(3,*))   ;POLARIZATION DEGREE
 T_XO=TTB(0:1,*)  ;(2,24)
 T_RL=TTB(2:3,*)

nxl=n_elements(xl)
if not KEYWORD_SET(CFOnly) THEN BEGIN
 tb=fltarr(2,nxl)
 tb(0,*)=xl/1e7
 tb(1,*)=t_tot
end


;SAVE CONTRIBUTION FUNCTION INTO SAV FILE
IF KEYWORD_SET(ConFun) or keyword_set(CFOnly) THEN BEGIN
 for i=0,nxl-1 do begin $
  cf_R_in(*,i)=interpol(cf_R(*,i), height_in, height00)
  cf_L_in(*,i)=interpol(cf_L(*,i), height_in, height00)
 end
cf_R=transpose(cf_R_in) ;lambda,height
cf_L=transpose(cf_L_in)

 END

;stop
END