; Program to read the various intensity spectrum data from Cooray (2016)
; Data are text files and .dat files

; Data is integrated over wavelength to get W/m^2/sr, and then integrated over solid angle to get flux density
; Original units are in nW/m2/sr

; Andrew Kren
; July 15, 2016

begin

path = "/scratch4/BMC/shout/Andrew.Kren/EBL/"

; list files in directory to determine how many dropsondes to read
fils = systemfunc("csh -c 'cd " + path + " ; ls *.txt'")
fils2 = systemfunc("csh -c 'cd " + path + " ; ls *.dat'")

print(fils)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; end of user's settings
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


  ; read individual files, each one has its own separate format
  
  
  ; read text files first
  cobe_data_from_mats = readAsciiTable(path+fils(0),4,"double",(/2,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_cobe_data_from_mats = cobe_data_from_mats(:,0) ; wavelength in um
  bright_cobe_data_from_mats = cobe_data_from_mats(:,2)*1e-9 ; brightness W/m2/sr
  
  cobe_data_gorjian = readAsciiTable(path+fils(1),4,"double",(/2,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_cobe_data_gorjian = cobe_data_gorjian(:,0) ; wavelength um
  bright_cobe_data_gorjian = cobe_data_gorjian(:,2)*1e-9 ; brightness W/m2/sr
  
  cobe_data_wright_and_reese = readAsciiTable(path+fils(3),4,"double",(/1,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_cobe_data_wright_and_reese = cobe_data_wright_and_reese(:,0) ; wavelength in um
  bright_cobe_data_wright_and_reese = cobe_data_wright_and_reese(:,2)*1e-9 ; brightness W/m2/sr

  integrated_galaxies_madau = readAsciiTable(path+fils(4),4,"double",(/2,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_integrated_galaxies_madau = integrated_galaxies_madau(:,0) ; wavelength in um
  bright_integrated_galaxies_madau = integrated_galaxies_madau(:,1)*1e-9 ; brightness W/m2/sr
  
  integrated_galaxies_totani = readAsciiTable(path+fils(5),3,"double",(/2,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_integrated_galaxies_totani = integrated_galaxies_totani(:,0) ; wavelength in um
  bright_integrated_galaxies_totani = integrated_galaxies_totani(:,1)*1e-9 ; brightness W/m2/sr
  
  irts_data = readAsciiTable(path+fils(6),3,"double",(/2,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_irts_data = irts_data(:,0) ; wavelength in um
  wave_irts_data = wave_irts_data(::-1)
  bright_irts_data = irts_data(:,1)*1e-9 ; brightness W/m2/sr
  bright_irts_data = bright_irts_data(::-1)
  
  optical_ebl_bernstein = readAsciiTable(path+fils(7),3,"double",(/2,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_optical_ebl_bernstein = optical_ebl_bernstein(:,0) ; wavelength in um
  bright_optical_ebl_bernstein = optical_ebl_bernstein(:,1)*1e-9 ; brightness W/m2/sr
  
  
  ; read .dat files next
  cmb_ebl = readAsciiTable(path+fils2(0),2,"double",0) 
  wave_cmb_ebl = cmb_ebl(:,0)*1e6 ; wavelength in um
  bright_cmb_ebl = cmb_ebl(:,1)*1e-9 ; brightness W/m2/sr
  
  comptel_ebl = readAsciiTable(path+fils2(1),4,"double",0) 
  wave_comptel_ebl = comptel_ebl(:,0)*1e6 ; wavelength in um
  wave_comptel_ebl = wave_comptel_ebl(::-1) ; reverse
  bright_comptel_ebl = comptel_ebl(:,1)*1e-9 ; brightness W/m2/sr
  bright_comptel_ebl = bright_comptel_ebl(::-1)
  upper_comptel_ebl = comptel_ebl(:,2)*1e-9 ;brightness upper limit
  upper_comptel_ebl = upper_comptel_ebl(::-1)
  lower_comptel_ebl = comptel_ebl(:,3)*1e-9 ;brightness lower limit
  lower_comptel_ebl = lower_comptel_ebl(::-1)
  
  egret_ebl = readAsciiTable(path+fils2(2),2,"double",0) 
  wave_egret_ebl = egret_ebl(:,0)*1e6 ;wavelength in um
  wave_egret_ebl = wave_egret_ebl(::-1) ; reverse
  bright_egret_ebl = egret_ebl(:,1)*1e-9 ; brightness W/m2/sr
  bright_egret_ebl = bright_egret_ebl(::-1)
  
  fermilat_ebl = readAsciiTable(path+fils2(3),4,"double",0) 
  wave_fermilat_ebl = fermilat_ebl(:,0)*1e6 ; wavelength in um
  wave_fermilat_ebl = wave_fermilat_ebl(::-1) ; reverse
  bright_fermilat_ebl = fermilat_ebl(:,1)*1e-9 ; brightness W/m2/sr
  bright_fermilat_ebl = bright_fermilat_ebl(::-1)
  upper_fermilat_ebl = fermilat_ebl(:,2)*1e-9 ;brightness upper limit
  upper_fermilat_ebl = upper_fermilat_ebl(::-1)
  lower_fermilat_ebl = fermilat_ebl(:,3)*1e-9 ;brightness lower limit
  lower_fermilat_ebl = lower_fermilat_ebl(::-1)
  
  heao_ebl = readAsciiTable(path+fils2(4),4,"double",0) 
  wave_heao_ebl = heao_ebl(:,0)*1e6 ; wavelength in um
  wave_heao_ebl = wave_heao_ebl(::-1) ; reverse
  bright_heao_ebl = heao_ebl(:,1)*1e-9 ; brightness W/m2/sr
  bright_heao_ebl = bright_heao_ebl(::-1)
  upper_heao_ebl = heao_ebl(:,2)*1e-9 ;brightness upper limit
  upper_heao_ebl = upper_heao_ebl(::-1)
  lower_heao_ebl = heao_ebl(:,3)*1e-9 ;brightness lower limit
  lower_heao_ebl = lower_heao_ebl(::-1)
  
  heaomid_ebl = readAsciiTable(path+fils2(5),2,"double",0) 
  wave_heaomid_ebl = heaomid_ebl(:,0)*1e6 ;wavelength in um
  bright_heaomid_ebl = heaomid_ebl(:,1)*1e-9 ; brightness W/m2/sr
  
  integral_ebl = readAsciiTable(path+fils2(6),2,"double",0) 
  wave_integral_ebl = integral_ebl(:,0)*1e6 ;wavelength in um
  bright_integral_ebl = integral_ebl(:,1)*1e-9 ; brightness W/m2/sr
  
  optical_ebl = readAsciiTable(path+fils2(7),3,"double",(/2,0/)) ; skip two headers, no upper-lower limit on this dataset file
  wave_optical_ebl = optical_ebl(:,0) ; wavelength in um
  bright_low_optical_ebl = optical_ebl(:,1)*1e-9 ; brightness W/m2/sr
  bright_high_optical_ebl = optical_ebl(:,2)*1e-9 ; brightness W/m2/sr
  
  rxte_ebl = readAsciiTable(path+fils2(8),2,"double",0) 
  wave_rxte_ebl = rxte_ebl(:,0)*1e6 ;wavelength in um
  wave_rxte_ebl = wave_rxte_ebl(::-1)
  bright_rxte_ebl = rxte_ebl(:,1)*1e-9 ; brightness W/m2/sr
  bright_rxte_ebl = bright_rxte_ebl(::-1)
  
  smm_ebl = readAsciiTable(path+fils2(9),4,"double",0) 
  wave_smm_ebl = smm_ebl(:,0)*1e6 ; wavelength in um
  wave_smm_ebl = wave_smm_ebl(::-1) ; reverse
  bright_smm_ebl = smm_ebl(:,1)*1e-9 ; brightness W/m2/sr
  bright_smm_ebl = bright_smm_ebl(::-1)
  upper_smm_ebl = smm_ebl(:,2)*1e-9 ;brightness upper limit
  upper_smm_ebl = upper_smm_ebl(::-1)
  lower_smm_ebl = smm_ebl(:,3)*1e-9 ;brightness lower limit
  lower_smm_ebl = lower_smm_ebl(::-1)
  
  ; integrate over wavelength and solid angle (the integration over solid angle involves multiplying by 4pi)
  pi = 3.141592654
  
  test1 = wave_cobe_data_from_mats(0:1)
  test2 = bright_cobe_data_from_mats(0:1)
  int_cobe_data_from_mats = simpne(test1,test2)*4.0*pi ; W/m^2

  
  int_low_optical_ebl = simpne(wave_optical_ebl,bright_low_optical_ebl)*4.0*pi
  int_high_optical_ebl = simpne(wave_optical_ebl,bright_high_optical_ebl)*4.0*pi

  print(int_low_optical_ebl)
  print(int_high_optical_ebl)

end