;+
; NAME:
;       rebin
;
; PURPOSE:
;       To put input solar spectra to a selected binning scheme in
;       production mode, or to get both rebinned spectra and 
;       major species' cross section and branching ratio in 
;       development mode.
;
; CATEGORY:
;       core procedure called by the main program.
;
; CALLING SEQUENCE:
;       rebin,mode,bin_file,in_spectra,out_spectra,sigma_o,sigma_o2,sigma_n,sigma_n2
;
; INPUTS:
;       mode:         set to 'dev' to run under development mode; 
;                     set to 'prod' to run under production mode.
;                     Development mode means to calculate both
;                     tabulated spectra and tabulated cross section
;                     and branching ratio. Production mode means to
;                     calculate tabulated spectra only.
;       bin_file:     the selected binning scheme file.
;       in_spectra:   a two dimensional array that holds the input
;                     spectra. It has at least 3 columns:
;                     wave short boundary in A, long boundary in A, 
;                     solar flux in photon cm^-2 s^-1, more fluxes if any.
;                     For example, for Hinteregger solar minimum
;                     spectrum, in_spectra has 5 columns:
;                     wave short, wave long, solar flux, flux from
;                     chromsphere, flux from corona. For Woods 
;                     reference spectrum, it also has 5 columns:
;                     wave short, wave long, solar flux, 
;                     27-day variation, 11-year variation.
;
; OUTPUTS:
;       out_spectra:  a two dimensional array that holds the output
;                     spectra. It has same data structure as in_spectra.
;       sigma_o:      a two dimensional array that holds cross section
;                     and branching ratios for atomic Oxygen. It has 
;                     8 columns.Its first 6 columns are branching ratios,
;                     last two columns are ionization and absorption 
;                     cross section. Cross section is in Megabarns.
;       sigma_o2:     has same data structure as sigma_o but for 
;                     molecular Oxygen.
;       sigma_n:     has same data structure as sigma_o but for 
;                     atomic Nitrogen.
;       sigma_n2:     has same data structure as sigma_o but for
;                     molecular Nitrogen.
;
; COMMON BLOCKS:
;       None.
;
; PROCEDURE:
;       branching ratios are calculated using Conway data (Availability: 
;       N2 and O2: 18.62 A to 1050.01 A, for O: 18.62A to 910.01 A). In the
;       wave length range without Conway data, data are assigned according
;       to the data in the GLOW model.
;
;       Major species' cross section data are calculated using Fennelly data
;       (Availabilty: 23.70 A - 1027 A). In the wave length range less than 50 A,
;       Henke data are used. In the wave length range greater than 1027, data
;       are assigned according to the data in the GLOW model.
;
; ROUTINES CALLED:
;       getbins, read_gen
;
; MODIFICATION HISTORY:
;       12/12/02, Liying Qian, Initial Version.
;
;+

pro rebin, mode,bin_file,in_spectra,out_spectra,sigma_o, sigma_o2, sigma_n, $
                sigma_n2

; get binning scheme

getbins, bin_file, wave1,wave2

n_bins=n_elements(wave1)

; The following is commented out so that the output wave bins are the same 
; as GLOW and T*GCM

;wave=(wave1+wave2)/2.0
;index=sort(wave)
;wave1=wave1(index)
;wave2=wave2(index)
;wave=wave(index)

iflines=(wave1 eq wave2)  ; 1 where bins are lines, 0 for continuum bins

; unpack input spectra

refwvln=0.5*(in_spectra[0,*]+in_spectra[1,*])
n_rows=n_elements(refwvln)
result=size(in_spectra)
n_columns=result[1]
n_cind=n_columns-2
refflux=fltarr(n_cind,n_rows)
refflux[0:n_cind-1,*]=in_spectra[2:n_columns-1,*]

; check lines in the input spectra

lines=wave1*iflines
ifreflines=fltarr(n_rows)
for i=0,n_rows-1 do begin
  for j=0,n_bins-1 do begin
     if ((lines(j) eq refwvln(i)) and (iflines(j) eq 1)) then $
     ifreflines(i)=1
  endfor
endfor

; read cross section data 

at_n2= 986.30   ; N2 absorption threshhold
at_o2= 2000.00  ; O2 absorption threshhold
at_o = 913.00   ; O absorption threshhold
it_n2= 798.00   ; N2 ionization threshold
it_o2= 1025.72  ; O2 ionization threshold,
it_o = at_o     ; O ionization threshold

loc='/hao/ada1/lqian/x/'
fenn=read_gen(loc+'phfenn.tab' ,1,10,1945)
fenno=read_gen(loc+'phfenno.tab',1,8,1945)
henkeo=read_gen(loc+'henkeo.tab',1,3,343)
henken=read_gen(loc+'henken.tab',1,3,342)

; Compile Fennely, Henke, and GLOW cross section data

ind_h1=where(henkeo(0,*) le 50)
ind_h2=where(henken(0,*) le 50)
ind_f1=where(fenno(0,*) gt 50)
ind_f2=where(fenn(0,*) gt 50)
lo=[transpose(henkeo(0,ind_h1)),transpose(fenno(0,ind_f1))]
lo2=[transpose(henkeo(0,ind_h1)),transpose(fenn(0,ind_f2))]
ln1=[transpose(henken(0,ind_h2)),transpose(fenn(0,ind_f2))]
ln2=[transpose(henken(0,ind_h2)),transpose(fenn(0,ind_f2))]

ao_ext=[transpose(henkeo(1,ind_h1)),transpose(fenno(7,ind_f1))]
io_ext=ao_ext
an_ext=[transpose(henken(1,ind_h2)),transpose(fenn(9,ind_f2))]
in_ext=an_ext
ao2_ext=[transpose(henkeo(2,ind_h1)),transpose(fenn(5,ind_f2))]
io2_ext=[transpose(henkeo(2,ind_h1)),transpose(fenn(8,ind_f2))]
an2_ext=[transpose(henken(2,ind_h2)),transpose(fenn(1,ind_f2))]
in2_ext=[transpose(henken(2,ind_h2)),transpose(fenn(4,ind_f2))]

n_fenn=n_elements(fenn(0,*))
print,'fenn(0,n_fenn-1)=',fenn(0,n_fenn-1)
ind=where(refwvln gt fenn(0,n_fenn-1),count)
if (count ne 0) then begin
  lo2=[lo2,refwvln(ind)]
  ao2_ext=[ao2_ext,refwvln(ind)]
  io2_ext=[io2_ext,refwvln(ind)]
endif

wave_arr=[fenn(0,n_fenn-1),1050.00,1100.00,1150.00,1200.00,1250.00,1300.00,1350.00,$
	  1400.00,1450.00,1500.00,1550.00,1600.00,1650.00,1700.00,1750]
sigma_arr=[1.15,1.00,0.40,1.40,13.00,0.40,$
	       2.20,12.00,15.00,13.00,10.00,6.00,3.40,1.50,0.50]
n_fenn=n_elements(wave_arr)
for i=0,n_fenn-2 do begin
  ind=where((lo2 ge wave_arr[i]) and (lo2 lt wave_arr[i+1]),count)
  if count ne 0 then begin
    ao2_ext(ind)=sigma_arr[i]
    io2_ext(ind)=0
  endif
endfor
ind=where(lo2 ge wave_arr(n_fenn-1),count)
if (count ne 0) then begin
  ao2_ext(ind)=0
  io2_ext(ind)=0
endif

; Fennallly file has Ly beta cross section = 0 for ionization of O2,
; so assign a correct value

indx=where(lo2 eq it_o2)
io2_ext(indx(0))=1.0
	
fennsigin2=interpol(in2_ext, ln2,refwvln)
fennsigin =interpol(in_ext,ln1,refwvln)
fennsigio2=interpol( io2_ext,lo2,refwvln)
fennsigio =interpol(io_ext,lo,refwvln)
fennsigan2=interpol( an2_ext,ln2,refwvln)
fennsigan =interpol(an_ext,ln1,refwvln)
fennsigao2=interpol( ao2_ext,lo2,refwvln)
fennsigao =interpol(ao_ext,lo,refwvln)


; read Conway cross section data for development mode

if mode eq 'dev' then begin
  conwayn2_file=read_gen(loc+'photon2.tab',2,15,808)
  conwayo2_file=read_gen(loc+'photoo2.tab',2,14,808)
  conwayo_file=read_gen(loc+'photoo.tab' ,2,8,645)

  ; map cross section data to solar spectra data  

  conwayn2=fltarr(15,n_rows)
  conwayo2=fltarr(14,n_rows)
  conwayo =fltarr(8,n_rows)

  conwayn2(0,*)=refwvln
  conwayo2(0,*)=refwvln
  conwayo (0,*)=refwvln

  for i=1,14 do $
    conwayn2(i,*)=interpol(conwayn2_file(i,*),conwayn2_file(0,*),refwvln)

  for i=1,13 do $
    conwayo2(i,*)=interpol(conwayo2_file(i,*),conwayo2_file(0,*),refwvln)

  for i=1,7 do $
    conwayo(i,*)=interpol(conwayo_file(i,*),conwayo_file(0,*),refwvln)

  ; zero out the cross section data where there is no original data

  n_con=n_elements(conwayn2_file(0,*))
  index=where((conwayn2(0,*) gt conwayn2_file(0,n_con-1)) $
  	and (conwayn2(0,*) lt conwayn2_file(0,0)))
  if (index[0] ne -1) then $
    for i=1,14 do conwayn2(i,index)=0.

  n_con=n_elements(conwayo2_file(0,*))
  index=where((conwayo2(0,*) gt conwayo2_file(0,n_con-1)) $
      	and (conwayo2(0,*) lt conwayo2_file(0,0)))
  if (index[0] ne -1) then $
    for i=1,13 do conwayo2(i,index)=0.

  n_con=n_elements(conwayo_file(0,*))
  index=where((conwayo(0,*) gt conwayo_file(0,n_con-1)) $
  	and (conwayo(0,*) lt conwayo_file(0,0)))
  if (index[0] ne -1) then $
    for i=1,7 do conwayo(i,index)=0.
endif

; Get solar flux into model bins.
; Calculate flux weighted cross sections,branching ratio
; as well if in development mode

modflux=fltarr(n_cind,n_bins)

if mode eq 'dev' then begin
  modsigin2=fltarr(n_bins)
  modsigin=fltarr(n_bins)
  modsigio2=fltarr(n_bins)
  modsigio =fltarr(n_bins)
  modsigan2=fltarr(n_bins)
  modsigan=fltarr(n_bins)
  modsigao2=fltarr(n_bins)
  modsigao =fltarr(n_bins)
  probn2=fltarr(6,n_bins)
  probo2=fltarr(4,n_bins)
  probo =fltarr(5,n_bins)
endif

flag_low=(n_bins lt 45)
count_1=0
count_2=0
count_3=0

for ibin=0,n_bins-1 do begin

  ; the following if statement should handle lines or continua, it does
  ; assume though that the line is one of the Hinteregger lines

  if iflines(ibin) eq 1 then begin

    ind=where((refwvln ge wave1(ibin)) and (refwvln le wave2(ibin)))

  endif else begin

    if flag_low eq 1 then begin

       case wave1[ibin] of

	 650.0: begin
	   case count_1 of
	     0: begin
	       ind=where( (refwvln ge wave1(ibin)) and $
		     (refwvln lt wave2(ibin)) and $
		     (fennsigan2 lt 31) )
	       count_1=count_1+1
             end

             1: ind=where( (refwvln ge wave1(ibin)) and $
                     (refwvln lt wave2(ibin)) and $
                     (fennsigan2 ge 31) )
           endcase
	 end

	 798.0: begin
	   case count_2 of
	     0: begin
	       ind=where( (refwvln ge wave1(ibin)) and $
		     (refwvln lt wave2(ibin)) and $
		     (fennsigan2 lt 4) ) 
               count_2=count_2+1
             end

             1: begin
	       ind=where( (refwvln ge wave1(ibin)) and $
                     (refwvln lt wave2(ibin)) and $
                     (fennsigan2 ge 4)  and (fennsigan2 lt 31) )
               count_2=count_2+1
             end

             2: ind=where( (refwvln ge wave1(ibin)) and $
                     (refwvln lt wave2(ibin)) and $
                     (fennsigan2 ge 31) )
           endcase
	 end

	 913.0: begin
	   case count_3 of
	     0: begin
	       ind=where( (refwvln ge wave1(ibin)) and $
		     (refwvln lt wave2(ibin)) and $
		     (fennsigan2 lt 4) ) 
               count_3=count_3+1
             end

             1: begin
	       ind=where( (refwvln ge wave1(ibin)) and $
                     (refwvln lt wave2(ibin)) and $
                     (fennsigan2 ge 4)  and (fennsigan2 lt 31) )
               count_3=count_3+1
             end
	     2: begin
	       ind=where( (refwvln ge wave1(ibin)) and $
		     (refwvln lt wave2(ibin)) and $
		     (fennsigan2 ge 31) ) 
             end

           endcase
	 end

	 else: begin
           ind=where((refwvln ge wave1(ibin)) and (refwvln lt wave2(ibin)))
	 end

       endcase

    endif else ind=where((refwvln ge wave1(ibin)) and (refwvln lt wave2(ibin)))

    if (ind(0) ne -1) then begin
      indcont=(where(ifreflines(ind) eq 0))
      oldind=ind
      if indcont(0) ne -1 then begin 
	indcont=ind(indcont)
	ind=indcont
      endif
    endif

  endelse

  if (ind(0) ne -1) then begin

;    if wave1(ibin) eq 1050 then begin
;       print,[refwvln,refflux]
;    endif
    for cind=0,n_cind-1 do begin
      modflux(cind,ibin)=total(refflux(cind,ind))
    endfor
	
    if mode eq 'dev' then begin
      totflux=modflux[0,ibin]
      if totflux eq 0.0 then totflux=1.0
		
      modsigin2(ibin)=total(fennsigin2(ind)*refflux(0,ind))/totflux
      modsigin(ibin)=total(fennsigin(ind)*refflux(0,ind))/totflux
      modsigio2(ibin)=total(fennsigio2(ind)*refflux(0,ind))/totflux
      modsigio (ibin)=total(fennsigio (ind)*refflux(0,ind))/totflux
      modsigan2(ibin)=total(fennsigan2(ind)*refflux(0,ind))/totflux
      modsigan(ibin)=total(fennsigan(ind)*refflux(0,ind))/totflux
      modsigao2(ibin)=total(fennsigao2(ind)*refflux(0,ind))/totflux
      modsigao (ibin)=total(fennsigao (ind)*refflux(0,ind))/totflux

      ; calculate branching ratio
      ; O+ states: 4S,2Do,2Po,4P4,2Pe
      ; O2+ states: X,a+A,b,dissociative ionization  
      ; N2+ states: X,A,B,C,F,dissociative ionization (Conway thinks C and F 
      ; states are dissociative states, but Kirby disagree. Here we treat 
      ; C and F states as what Kirby thinks to be consistent with GLOW model
      ; For N2+, we also need to subtract N2++ column from the FRAG column

      num_r=n_elements(ind)
      sum_sigi=fltarr(num_r)

      ; N2+, Conway array 5,6,7,8,9,3


      conwayn2(3,ind)=conwayn2(3,ind)-conwayn2(8,ind)- $
  		    conwayn2(9,ind)-conwayn2(13,ind)

      ind_arr=[5,6,7,8,9,3]
      n_ind=n_elements(ind_arr)
      term=fltarr(n_ind)

      for ii=0,n_ind-1 do term[ii]=total(conwayn2(ind_arr[ii],ind)*refflux(0,ind))/totflux
      sum_term=total(term)

      if (sum_term ne 0) then begin
        for ii=0,n_ind-1 do probn2(ii,ibin)=term[ii]/sum_term
      endif

      ; O2, Conway array 5,6,7,3

      ind_arr=[5,6,7,3]
      n_ind=n_elements(ind_arr)
      term=fltarr(n_ind)

      for ii=0,n_ind-1 do term[ii]=total(conwayo2(ind_arr[ii],ind)*refflux(0,ind))/totflux
      sum_term=total(term)

      if (sum_term ne 0) then begin
        for ii=0,n_ind-1 do probo2(ii,ibin)=term[ii]/sum_term
      endif

      ; O, Conway array 2,3,4,5,6

      ind_arr=[2,3,4,5,6]
      n_ind=n_elements(ind_arr)
      term=fltarr(n_ind)

      for ii=0,n_ind-1 do term[ii]=total(conwayo(ind_arr[ii],ind)*refflux(0,ind))/totflux
      sum_term=total(term)

      if (sum_term ne 0) then begin
        for ii=0,n_ind-1 do probo(ii,ibin)=term[ii]/sum_term
      endif

    endif
  endif

endfor

; manully set N2 absorption cross section for 1nm bin schme at wave length 970 A

if mode eq 'dev' then begin
  if (strpos (bin_file,'nm_bins') ne -1) then begin
    ind=where(wave1 eq 970)
    modsigan2(ind(0))=2.24
  endif

  ; Assign a correct O2 absorption cross section data at Ly_alpha
  ; or for the bin that include Ly-alpha (1215.67)

  indx=where(refwvln eq 1215.67)
  if (indx(0) ne -1) then fennsigao2(indx(0))=0.01

  indx=where(wave1 eq 1215.67)
  if (indx(0) ne -1) then modsigao2(indx(0))=0.01

  indx=where((wave1 lt 1215.67) and (wave2 gt 1215.67))
  ind_temp=where(wave1 eq 1215.67)
  if ((indx(0) ne -1) and (ind_temp[0] eq -1)) then begin
    modsigao2(indx(0))=0.01
  endif

  ; Assign a correct cross section O2 for the bin that 
  ; include Ly-beta (1025.72)

  indx=where((wave1 lt 1025.72) and (wave2 gt 1025.72))
  ind_temp=where(wave1 eq 1025.72)
  if ((indx(0) ne -1) and (ind_temp[0] eq -1)) then begin
    modsigao2(indx(0))=1.632
    modsigio2(indx(0))=1.0
  endif

  ; Zero out long wavelength bins

  indx=where(wave1 gt it_n2) & modsigin2(indx)=0.0 & probn2(*,indx)=0.0
  indx=where(wave1 gt it_o2) & modsigio2(indx)=0.0 & probo2(*,indx)=0.0
  indx=where(wave1 gt it_o ) & modsigio (indx)=0.0 & probo (*,indx)=0.0


  ; Conway does not have data for the shortest wavelength bins, therefore 
  ; assign data for branching ratio according to GLOW model

  ind=where(wave1 lt conwayo_file(0,0)) 
  for i=0,n_elements(ind)-1 do probo[0:4,ind[i]]=[0.30,0.32,0.21,0.09,0.08]
  ind=where(wave1 lt conwayo2_file(0,0)) 
  for i=0,n_elements(ind)-1 do probo2[0:3,ind[i]]=[0.00,0.00,0.00,1.00]
  ind=where(wave1 lt conwayn2_file(0,0)) 
  for i=0,n_elements(ind)-1 do probn2[0:5,ind[i]]=[0.02,0.02,0.00,0.00,0.00,0.96]
endif

; write solar flux/cross section data with the original solar spectra resolution
; for diagnosis in the development mode

if mode eq 'dev' then begin
  file_name='orig'+'_'+strcompress(n_rows,/remove_all)+'.dat'
  openw,1,file_name
  printf,1,n_rows
  printf,1,'Output of rebin.pro, wave length in A, cross section in megabarns'
  printf,1,'  wave ', '   flux  ','    O-I ','       O-A ','       O2-I ',$
   '      O2-A ','      N2-I ','    N2-A '
  format='$(f7.2,2x,e10.4,6(2x,f9.5))'
  for i=0,n_rows-1 do begin
    printf,1,format,refwvln[i],refflux[0,i],fennsigio[i],fennsigao[i],fennsigio2[i],$
  	     fennsigao2[i],fennsigin2[i],fennsigan2[i]
  endfor
  close,1
endif

; pack output 

out_spectra=fltarr(n_columns,n_bins)
out_spectra[0,*]=wave1
out_spectra[1,*]=wave2
out_spectra[2:n_columns-1,*]=modflux[0:n_cind-1,*]

if mode eq 'dev' then begin
  bso2=fltarr(n_bins)
  if n_bins eq 59 then begin
     bso2[0:36]=0.
     bso2[37:46]=[0.01,0.03,0.07,0.07,0.09,0.10,0.09,0.10,0.03,0.01]
     bso2[47:58]=0.
  endif else if n_bins eq 37 then begin
     bso2[0:15]=0.
     bso2[16:24]=[0.04,0.04,0.04,0.07,0.08,0.09,0.10,0.03,0.01]
     bso2[25:36]=0.
  endif else if n_bins eq 132 then begin
     bso2[0:83]=0.
     bso2[84:115]=[0.002,0.004,0.006,0.008,0.01,0.01,0.02,0.02,0.03,0.03, $
		   0.04,0.05,0.05,0.06,0.07,0.07,0.07,0.07,0.07,0.08,0.09,0.09, $
		   0.09,0.10,0.010,0.10,0.09,0.07,0.06,0.04,0.03,0.01]
     bso2[116:131]=0.
  endif

  sigma_o=fltarr(5,n_bins)
  sigma_n=fltarr(1,n_bins)
  sigma_o2=fltarr(7,n_bins)
  sigma_n2=fltarr(6,n_bins)

  sigma_o[0,*]=modsigao
  sigma_n[0,*]=modsigan
  sigma_o2[0,*]=modsigao2
  sigma_n2[0,*]=modsigan2

  sigma_o[2,*]=probo[0,*]+probo[3,*]
  sigma_o[3,*]=probo[1,*]+probo[4,*]*0.72
  sigma_o[4,*]=probo[2,*]+probo[4,*]*0.18

  sigma_o2[2,*]=probo2[0,*]
  sigma_o2[3,*]=probo2[1,*]+probo2[2,*]
  sigma_o2[4,*]=probo2[3,*]

  sigma_n2[2,*]=probn2[0,*]
  sigma_n2[3,*]=probn2[1,*]+probn2[2,*]+probn2[3,*]+probn2[4,*]
  sigma_n2[4,*]=probn2[5,*]

  for k=0,n_bins-1 do begin
    if sigma_o[0,k] eq 0 then sigma_o[1,k]=0 else $
			sigma_o[1,k]=modsigio[k]/modsigao[k]
    if sigma_o2[0,k] eq 0 then sigma_o2[1,k]=0 else $
			sigma_o2[1,k]=modsigio2[k]/modsigao2[k]
    if sigma_n2[0,k] eq 0 then sigma_n2[1,k]=0 else $
			sigma_n2[1,k]=modsigin2[k]/modsigan2[k]

    sigma_o[2:4,k]=sigma_o[2:4,k]*sigma_o[1,k]
    sigma_o2[2:4,k]=sigma_o2[2:4,k]*sigma_o2[1,k]
    sigma_n2[2:4,k]=sigma_n2[2:4,k]*sigma_n2[1,k]
  endfor

  ; changes made on 02/20/04: for a version of cross section and branching for WACCM, which
  ; covers wavelength from 0.5A to 1210A, we consider that dissociation of O2 leads 
  ; to O(3P) only. If the wavelength extends to 1750A, then use the code for O(1D) and O(1S)

  ; O2 dissociation yields (O1D)

;  index=where( (wave1 gt 1349) and (wave2 lt 1751),count)
;  if (count ne 0) then begin
;     sigma_o2[5,*]=0
;     sigma_o2[5,index]=1
;     sigma_o2[5,index]=sigma_o2[5,index]*(1-sigma_o2[1,index])
;  endif

  ; O2 dissociation yields O(1S)

;  sigma_o2[6,*]=bso2
;  sigma_o2[6,*]=sigma_o2[6,*]*(1-sigma_o2[1,*])

  ; O2 dissociation yields O(3P): O2+hv-->O(3P)+O(3P), a good approximation for EUV

  for k=0,n_bins-1 do begin
    if sigma_o2[0,k] eq 0. then sigma_o2[5,k]=0. else $
        sigma_o2[5,k]=1-sigma_o2[1,k]
  endfor

 ; dissociation of N2 (yields N(4S) and N(2D))

  for k=0,n_bins-1 do begin
    if sigma_n2[0,k] eq 0. then sigma_n2[5,k]=0. else $
        sigma_n2[5,k]=1-sigma_n2[1,k]
  endfor

endif

; energy conservation check

ind1=where(in_spectra[0,*] lt wave2[n_bins-1])
ind2=where(refwvln lt wave2[n_bins-1])

if n_columns eq 3 then begin
   print,'in rebin, input:',total(in_spectra[2,ind1])
   print,'in rebin, output:',total(out_spectra[2,*])
endif else if n_columns gt 3 then begin
   print,'in rebin, input:',total(in_spectra[2,ind1]),total(in_spectra[3,ind1])
   print,'in rebin, output:',total(out_spectra[2,*]),total(out_spectra[3,*])
endif

return
end