;+
; NAME:
;	get_spectra
;
; PURPOSE:
;       Given an solar spectra file, read the solar spectra,
;       convert wave length unit to Angstrom and solar flux 
;       unit to photon cm^-2 s^-1 if necessary. For model
;       spectra that carry solar activity scaling indices,
;       the indices are converted to corresponding fluxes
;       to carry along implicitly as fluxes.
;
; CATEGORY:
;       utility function called by main program.
;
; CALLING SEQUENCE:
;       out_spectra=get_spectra(solar_file,group)
;
; INPUTS:
;       solar_file: (='dir_name/file_name')
;       group:      (=1, =2, =3)
;
; OUTPUTS:
;       out_spectra: two dimensional array that holds 
;                     output solar spectra. The array 
;                     has at least two columns: wave 
;                     length in A, solar flux in
;                     photon cm^-2 s^-1, any number of
;                     extra flux columns. For example,
;                     Hinteregger solar minimum spectrum
;                     has 4 columns: wave length in A,
;                     solar flux in photon cm^-2 s^-1,
;                     flux from chromsphere, flux from
;                     corona, Woods reference 
;                     spectrum has four columns:
;                     wave length in A, solar flux in
;                     photon cm^-2 s^-1, 27-day variation
;                     11-year variation.
; COMMON BLOCKS:
;	None.
;
; PROCEDURE:
;       some usual/useful conversion of units:
;       1 nm = 10 A
;       1 mW m^-2 = 1 erg cm^-2 s^-1
;       flux(photon cm^-2 s^-1) = (6.242e11/12398)*wave(A) *flux(erg cm^-2 s^-1) 
;
;       a list of cases and their explainations:
;       number=1: default case: read input data given the following parameters:
;         n_comm: number of comment lines
;         n_columns: number of columns of data
;         n_lines: number of lines of data
;         unit_w: factor of unit transfer of wave length from others to Angstrom.
;         is_e: =1 if solar flux unit is in energy, =0 if not
;         unit_f: factor of unit transfer of solar flux from others to photon cm^-2 s^-1
;       number=2: read Hinteregger reference spectrum
;       number=3: read Woods reference spectrum
;       number=4: read SEE spectra in either ASCII or NetCDF format
;
; ROUTINES CALLED:
;       read_gen, read_hs, read_ws, read_netcdf, get_boundary
;
; MODIFICATION HISTORY:
;       12/12/02, Liying Qian, Initial Version
;
;+
;
function get_spectra, solar_file, group

number=1

if (group eq 1) then number=2 $
else if (strpos(solar_file,'ref_min_27day_11yr') ne -1) then number=3 $
else if (strpos(solar_file,'see.dat') ne -1) then number=4 $
else if (strpos(solar_file,'see_flare_2002111_v7b.dat') ne -1) then number=4 $
else if (strpos(solar_file,'L3_merged_2002343_006.ncdf') ne -1) then number=4 $
else if (strpos(solar_file,'see__L3_merged_2003062_006.ncdf') ne -1) then number=4 $
else if (strpos(solar_file,'see__L3_merged_2003269_006.ncdf') ne -1) then number=4 $
else if (strpos(solar_file,'see__L3_merged_2003308_007.ncdf') ne -1) then number=4 $
else if (strpos(solar_file,'see__L3A_merged_007_2003308.ncdf') ne -1) then number=4 $
else if (strpos(solar_file,'see__L3_merged_2003320_007.ncdf') ne -1) then number=4 $
else if (strpos(solar_file,'see__L3A_merged_007_2003320.ncdf') ne -1) then number=4 $
else if (strpos(solar_file,'EGS_1nm_Rocket_8Feb2002_v5.dat') ne -1) then number=5 $
else if (strpos(solar_file,'rocket_1nm_8Feb2002_v5.dat') ne -1) then number=6     $
else if (strpos(solar_file,'see__L3A_merged_007_2004047.ncdf') ne -1) then number=7 $
else if (strpos(solar_file,'see__L3_merged_2004040_007.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2005007_007.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2005184_008.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2005286_008.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3A_merged_2005286_008.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3A_merged_2006187_008.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2006187_008.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2007019_008.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2007120_009.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2007238_009.ncdf.repaired') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2008180_009.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2009003_009.ncdf') ne -1) then number=7 $ 
else if (strpos(solar_file,'see__L3_merged_2009316_010.ncdf') ne -1) then number=7 $ 
;else if (strpos(solar_file,'FISM') ne -1) then number=7 $ 
else if (strpos(solar_file,'woods.dat') ne -1) then number=8                       $
else if (strpos(solar_file,'2003301.dat') ne -1) then number=8   $
else if (strpos(solar_file,'2003308.dat') ne -1) then number=8   $
else if (strpos(solar_file,'2003301_spectra_1nm.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2006339flare_kent_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2006339flare_fism_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2003301flare_fism_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2004198flare_fism_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2005251flare_fism_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2001_fism_daily_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2001104-106_fism_flare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2001104-106_fism_noflare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2005_fism_daily_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2005_0906-0910_fism_flare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2005_0906-0910_fism_noflare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2011_fism_daily_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2011_0212-0215_eve_flare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2011_0212-0215_eve_noflare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2012_fism_daily_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2012_0122-0124_fism_flare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2012_0122-0124_fism_noflare_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2011001_365_fism_0_195nm_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'2003storm_fism_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X17_cent_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X17_noflare.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X17_limb_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_year2_control_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X10_preflare_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X10_long_gp_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X10_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X10_quick_ip_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_X10_afterflare_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'fism_20100503_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'eve_20100503_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2003301-305_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2003300-305_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2003300-305_noflare_orig.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'2005250-254_orig.dat') ne -1) then number=8 $
else if (strpos(solar_file,'fism2008_min.dat') ne -1) then number=8 $
else if (strpos(solar_file,'eve2008_min.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'Woods2008_min.dat') ne -1) then number=8 $
else if (strpos(solar_file,'Woods1996_min.dat') ne -1) then number=8 $
else if (strpos(solar_file,'Woods2008_min_ver2.dat') ne -1) then number=8 $
else if (strpos(solar_file,'Woods1996_min_ver2.dat') ne -1) then number=8 $ 
else if (strpos(solar_file,'HaoBinsEUV0Fmt.txt') ne -1) then number=9 

case number of

   1:begin

; for ssflux.dat
;     n_comm=0
;     n_columns=3
;     n_lines=59
;     w_column=2
;     unit_w=1
;     unit_f=1
;     is_e=0

; for euvac.dat
     n_comm=0
     n_columns=4
     n_lines=61
     w_column=2
     unit_w=1
     unit_f=1
     is_e=0

; for euvac_spec.dat
;     n_comm=0
;     n_columns=3
;     n_lines=22
;     w_column=2
;     unit_w=1
;     unit_f=1
;     is_e=0

     out_spectra=fltarr(n_columns,n_lines)
     data=read_gen(solar_file,n_comm,n_columns,n_lines)
     out_spectra[0:w_column-1,*]=data[0:w_column-1,*]*unit_w
     if is_e eq 0 then begin
       out_spectra[w_column:n_columns-1,*]=data[w_column:n_columns-1,*]*unit_f
     endif else begin
       out_spectra[1,*]=data[f_column]*unit_f *out_spectra[0,*]
     endelse 
   end

   2:begin
     if strpos(solar_file, 'sc21refw') ne -1 then begin
	solar_activity='low'
	n_lines=1661
     endif else begin
	solar_activity='high'
	n_lines=1664
     endelse

     out_spectra=read_hs(solar_file,solar_activity,n_lines)

   end

   3:begin
     
     data=read_ws(solar_file,/silent)
     out_spectra=fltarr(4,n_elements(data[0,*]))
     out_spectra[0,*]=data[0,*]*10   ; convert from nm to Angstrom
     out_spectra[1,*]=data[1,*]
     out_spectra[2,*]=data[1,*]*(data[2,*]-1) ; convert 27-day variation 
					      ; index to flux
     out_spectra[3,*]=data[1,*]*(data[3,*]-1) ; convert 11-year solar cycle 
					      ; variation index to flux
   end

   4:begin
     
     if strpos(solar_file,'see.dat') ne -1 then begin
       n_comm=2
       n_columns=3
       n_lines=373
       out_spectra=read_gen(solar_file,n_comm,n_columns,n_lines)
     endif else if strpos(solar_file,'see_flare_2002111_v7b.dat') ne -1 then begin
       n_comm=21
       n_columns=5
       n_lines=346
       out_spectra=read_gen(solar_file,n_comm,n_columns,n_lines)
     endif else if strpos(solar_file,'.ncdf') ne -1 then begin
       read_netcdf,solar_file,data,attributes
       days=data.date
       n_days=n_elements(days)

       if strpos(solar_file,'L3A') ne -1 then begin
	  w=data[0].sp.wave
	  n_lines=n_elements(w)
	  f=fltarr(n_lines,n_days)
	  for i=0,n_days-1 do begin
	     f[*,i]=data[i].sp.FLUX_FILLED
          endfor
; the following part of code only applies to see__L3A_merged_007_2003320.ncdf since the 
; part of the data the following code is trying to fix was missing in the
; above original data file
	  print,'Printing the problem day:'
	  print,days[9139]
	  print,'  '
	  print,f[0:20,9139]
	  print,'  '
	  print,f[0:20,9140]
	  f[7:16,9139]=10.*f[7:16,9140]
	  print,'  '
	  print,f[0:20,9139]
	  print,'  '
	  print,f[0:20,9140]
       endif else begin
          w=data.sp_wave
          f=data.sp_flux_filled
          n_lines=n_elements(w)
       endelse

       out_spectra=dblarr(n_days+1,n_lines)
       out_spectra[0,*]=double(w)
       for i=1,n_days do begin
         out_spectra[i,*]=double(f[*,i-1])
       endfor
     endif  
      
     result=size(out_spectra)
     n_columns=result[1]

     ; see flare data has unit of W m^-2 nm^-1 for solar flux
     ; but the actual wave interval is less than 1 nm (about 0.5 nm)
     ; thus calculate actural solar flux for each interval

     temp_spectra=get_boundary(out_spectra)
     for i=2,n_columns do begin
       temp_spectra[i,*]=temp_spectra[i,*]*(temp_spectra[1,*]-temp_spectra[0,*])
     endfor 
     out_spectra[1:n_columns-1,*]=temp_spectra[2:n_columns,*]

     ; convert wave unit from nm to Angstrom, solar flux unit
     ; from W M^-2 to photon cm^-2 s^-1 (from w m^-2 to erg cm^-2 s^-1 (1e3),
     ; then from erg cm^-2 s^-1 to photon cm^-2 s^-1)

     out_spectra[0,*]=10*out_spectra[0,*]
     for i=1,n_columns-1 do begin
       out_spectra[i,*]=(6.242e11/12398)*out_spectra[0,*]*out_spectra[i,*] *1e3
     endfor
   end


   5:begin

     n_comm=3
     n_columns=3
     n_lines=164
     temp_spectra=read_gen(solar_file,n_comm,n_columns,n_lines)

     out_spectra=fltarr(n_columns-1,n_elements(temp_spectra[0,*]))
     out_spectra[0:1,*]=temp_spectra[0:1,*]

     ; Convert wave length unit from nm to Angstrom

     out_spectra[0,*]=10*out_spectra[0,*]

   end

   6:begin

     n_comm=3
     n_columns=2
     n_lines=191
     out_spectra=read_gen(solar_file,n_comm,n_columns,n_lines)

     ; convert wave unit from nm to Angstrom, solar flux unit
     ; from mW M^-2 to photon cm^-2 s^-1 

     out_spectra[0,*]=10*out_spectra[0,*]
     out_spectra[1,*]=(6.242e11/12398)*out_spectra[0,*]*out_spectra[1,*]
   end

   7:begin
     
     read_netcdf,solar_file,data,attributes
     days=data.date
     n_days=n_elements(days)

;     if strpos(solar_file,'L3A') ne -1 then begin
;        w=data[0].sp.wave
;        n_lines=n_elements(w)
;        f=fltarr(n_lines,n_days)
;        for i=0,n_days-1 do begin
;            f[*,i]=data[i].sp.FLUX
;        endfor
;     endif else begin
        w=data.sp_wave
        f=data.sp_flux
        n_lines=n_elements(w)
;     endelse

     out_spectra=dblarr(n_days+1,n_lines)
     out_spectra[0,*]=double(w)
     for i=1,n_days do begin
        out_spectra[i,*]=double(f[*,i-1])
     endfor
      
     result=size(out_spectra)
     n_columns=result[1]

     ; see flare data has unit of W m^-2 nm^-1 for solar flux
     ; but the actual wave interval is less than 1 nm (about 0.5 nm)
     ; thus calculate actural solar flux for each interval
     ;
     ; new SEE data wave interval is 1nm, but the following code 
     ; won't cause problem-----lqian, 07/07/2005

     temp_spectra=get_boundary(out_spectra)
     for i=2,n_columns do begin
       temp_spectra[i,*]=temp_spectra[i,*]*(temp_spectra[1,*]-temp_spectra[0,*])
     endfor 
     out_spectra[1:n_columns-1,*]=temp_spectra[2:n_columns,*]

     ; convert wave unit from nm to Angstrom, solar flux unit
     ; from W M^-2 to photon cm^-2 s^-1 (from w m^-2 to erg cm^-2 s^-1 (1e3),
     ; then from erg cm^-2 s^-1 to photon cm^-2 s^-1)

     out_spectra[0,*]=10*out_spectra[0,*]
     for i=1,n_columns-1 do begin
       out_spectra[i,*]=(6.242e11/12398)*out_spectra[0,*]*out_spectra[i,*] *1e3
     endfor

   end

   8:begin
     n_comm=0
     ; need to change n_columns based on how many days in this flare study
;     n_columns=5761
;for 2003storm_fism_orig     
;     n_columns=23041
; for 2009 LWS project (fism_X17_cent_orig.dat, 9 days on 1 minute time resolution)
;     n_columns=12961
;     n_lines=195
; for 05/03/2010 EVE rocket, 22 day daily data(fism_20100503_orig.dat)
;     n_columns=23
;     n_lines=195
; for 2010 LWS project (fism_X10_preflare_orig.dat, 20 days on 1 minute time resolution)
;     n_columns=28801
;     n_lines=195
; for 2010 LWS project (fism_X10_preflare_orig.dat, 22 days on 1 minute time resolution)
;     n_columns=31681
;     n_lines=195
; for 2010 LWS project (fism_control_orig.dat, 25 hours on 1 minute time resolution)
;     n_columns=1501
;     n_lines=195
; for 2010 LWS project (fism_afterflare_orig.dat, 4 days on 1 minute time resolution)
;     n_columns=5761
;     n_lines=195
; for 2009 LWS project (2003301-305_orig.dat)
; and (2005250-254_orig.dat)
;     n_columns=7201
;     n_lines=195
; for 2009 LWS project (2003300-305_orig.dat)
;     n_columns=8641
;     n_lines=195
; for 2012 LWS eastern michgan proposal (2001_fism_daily_orig.dat)
;     n_columns=366
;     n_lines=195
; for 2012 LWS eastern michgan proposal (2001104-106_fism_flare_orig.dat)
; for 2012 LWS eastern michgan proposal (2001104-106_fism_noflare_orig.dat)
     n_columns=525601L
     n_lines=195
; for 2012 LWS team meeting  (2005_fism_daily_orig.dat and 2011_fism_daily_orig.dat)
;     n_columns=366
;     n_lines=195
; for 2012 LWS team meeting (2005_0906-0910_fism_flare_orig.dat,7days)
;     n_columns=10081
;     n_lines=195
; for 2012 LWS team meeting (2011_0212-0215_eve_flare_orig.dat,4days)
;     n_columns=5761
;     n_lines=195
; for 2012 flare for Phil (2012_0122-0124_fism_flare_orig.dat,3days)
;     n_columns=4321
;     n_lines=195
; for 2012 flare for Phil (2012_fism_daily_orig.dat,70days)
;     n_columns=71
;     n_lines=195
; for fism2008 min spectrum (50 days run)
;     n_columns=51
;     n_lines=195
; for eve2008 min spectrum (50 days run)
;     n_columns=52
;     n_lines=1221
; for Woods2008 min spectrum  and Woods1996 min spectrum (50 days run)
;     n_columns=51
;     n_lines=195
; for Woods2008 min spectrum  and Woods1996 min spectrum (366 days run)
;     n_columns=367
;     n_lines=195
     out_spectra=read_gen(solar_file,n_comm,n_columns,n_lines)
     print,' in get_spectra'
     print,out_spectra[0,0:10]
     ; convert wavelength from nm to A
     ; need to check whether the input data use wavelength in Angstrom or nm
     out_spectra[0,*]=out_spectra[0,*]*10.
   end
   ;
   ; read Juan Fontela's data. 
   9:begin
      nrows=2435
      ncolumns=34
      d=dblarr(ncolumns,nrows)
      a=' '
      nbins=ncolumns-2
      out_spectra=dblarr(nrows+2,nbins)
   ;Juan uses 37 low-resolution bins but does not have virtual bins,
   ; therefore, bin boundaries are hard-coded here.
      out_spectra[0,*]=[0.5,4,8,18,32,70,155,224,290,320,540,650,798,913,975,987,1027,1050,1100, $
                        1150,1215.67,1200,1250,1300,1350,1400,1450,1500,1550,1600,1650,1700]
      out_spectra[1,*]=[4,8,18,32,70,155,224,290,320,540,650,798,913,975,987,1027,1050,1100, $
                        1150,1200,1215.67,1250,1300,1350,1400,1450,1500,1550,1600,1650,1700,1750]
      close,1
      openr,1,solar_file
      readf,1,a
      readf,1,d
      close,1
      for j=0,nrows-1 do begin
         out_spectra[2+j,*]=d[2:ncolumns-1,j]
      endfor
   end

endcase

return, out_spectra
end