pro comp_spec, ramaty=ramaty,therm1=therm1,$
        label=label,brokep = brokep,tek=tek,simple=simple,bkyellow=bkyellow,$
        nofill = nofill,composite=composite,powercomp=powercomp,pinot=pinot,$
        therm2=therm2
;+
; NAME:  
;	COMP_SPEC (Main Program and a Subprogram called from HESI_MENU)
; PURPOSE:
;       To plot various combinations of solar flare spectra, from 
;       Thermal Bremsstahlung to pi not meson decay spectra.
;       Used currently for the HESI Instrument Proposal
; Currently Available:  Composite or Gamma-Ray Solar Spectrum for
; R. Ramaty's Positron/Nulcear lines, 2 Thermal Bremsstrahlung spectra, 
; nonthermal Bremsstrahlung (power law or Schwartz's broken
; power law), pi not meson decay, and a composite spectrum
;
; The program has options to plot all the spectra on one page,
; to plot only the gamma-ray portion of the spectra (and to not
; include the thermal spectra), to plot only ramaty's lines, or
; to plot the aforementioned with solid colors
; INSTRUCTIONS:
; Call Comp_spec with the desired keyword parameters, the default
; is to plot only the Fluence vs. Energy contained in the file
; perm$data:ramaty_lines_no_Co.dat  
;
; CALLing PROCEDURES:
;	COMP_SPEC,[/ramaty,/label,...]
; 
; 	Called Subprograms:
;       INIT_LINECOLORS
;       COMP_SPEC_PRO      
;       F_COMPOSITE
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
; file - the file containing one column of energies in MeV and
; one column of Fluence values in photons cm^-2 Mev^-1
; over_file - same as file but contains Ramaty positron and nuclear
; lines, if set without specifying a filename
; user_disk0:[eric.hesi_read]hesi2.dat is used, the file's plot will be
; overlayed onto the file keyword's plot unless file is not set
;
; OPTIONAL OUTPUT KEYWORDS:
; POWER - if the keyword is set (using /power or power=1), then the
; composite of the thermal, nonthermal (powerlaw), and ramaty lines
; will be plotted
; OVERPOWER - if set produces the composite spectrum with its
; constituent spectra overlayed in different colors
; ALLPOWER - if set produces a plot of only the constituent spectra
; without the composite plot
; RAMATY - if set produces a plot of the powerlaw(nonthermal), ramaty lines, and
; the composite spectrum in the gamma-ray range
; BIC_RAM - if set produces the ramaty keyword's plot with solid colors for
; each spectrum
; SIMPLE - replaces the lables "positron and nuclear..." with 'gamma Rays' and
; "nonthermal bremmstrahlung" with X-rays
; BKYELLOW - for color plots eliminates use of yellow backround
; SCHWARTZ - if set produces the overpower keyword's plot using R. Schwartz's
; PLOT_COMPOSITE program values instead of the values calculated in this 
; program
; LABEL- if set produces labels for plots using bic_Ram, schwartz, ramaty, and
; overpower keywords
; BROKEP - if set substitutes R. Schwartz's PLOT_COMPOSITE program values for the
; nonthermal (powerlaw) spectrum
; TEK - if set set_plot,'tek' instead of set_plot,'X'
; NOFILL - if set, the schwartz Composite Flare Spectrum will not plot solid 
; color blocks under the component spectra lines
; COMPOSITE - if set only the schwartz composite spectrum will be plotted
; POWERCOMP - if sets plots gamma ray composite spectrum with power law spec.
;      only works with ramaty keyword
;
; COMMON BLOCKS:
;  none
; SIDE EFFECTS
;  depending on your device and previous plotting activity, color tables
; on your screen may be rearranged when linecolors is called by the program
; RESTRICTIONS
;  Currently the program is not in the path, must be run from
;  USER_DISK0:[eric.hesi_read] on SDAC
; MODIFICATION HISTORY
;  -Written June 1994 by Eric Carzon, Hughes/STX for Brian Dennis
;  - Modifed October 1994 by EC to use corrected Schwartz program
;  - EAC 1-3 95 made program compatible with HESI_MENU, a user interface
;    routine
;  - RAS made modifications to subprogram comp_spec_plot.pro, to 
;    convolve Bremsstralung spectra in 1 keV bandwidths     3/10/95
;  - EAC added a pi not meson decay spectra to comp_spec_plot 3/15/95
;  - EAC modified calling keywords to simplify program use 3/15/95 added
;    use of init_linecolors for color values
;______________________________________________________________
;-
hd_cp = 'n'
do_color = 1
nxsize = 21
nysize = 14
xtitle = 'Photon Energy (MeV)'
ytitle = 'Fluence (photons MeV!u-1!n cm!u-2!n)
m2title = 'Composite Solar Flare Spectrum'
m1title = 'Solar Flare Gamma-Ray Spectrum'
mtitle = 'Fluence vs. Energy'
over_plot=0
over_power=0
all_power = 0
ramaty_plot = 0
bicolor_ram = 0
schw_pl = 0
powerlabs = 0
ramatylabs = 0
broken_power = 0
use_schwartz = 0
simp = 0
yell = 0
nfill = 0 
comp = 0
powcow = 0
pi_not = 1
prfont = -1

Print,'REMEMBER to CLEAN the FILE in an EDITOR, SO ONLY TWO COLUMNS OF NUMBERS REMAIN!'
if (keyword_set(tek)) then set_plot,'tek' else set_plot,'X'

plot_it:
if (keyword_set(composite)) then comp = 1
if (keyword_set(bkyellow)) then yell = 1
if (keyword_set(simple)) then simp = 1
if (keyword_set(nofill)) then nfill = 1
if (not keyword_Set(filename)) then filename='[sdac.hesi.menu]hesi.dat'
if (not keyword_set(over_filename)) then over_plot = 0 else over_plot = 1
if (keyword_set(power)) then begin 
          over_plot = 2
          over_filename = '[eric.hesi_read]hesi2.dat'
endif
if (keyword_set(overpower)) then begin
          over_plot = 2
          over_filename = '[eric.hesi_Read]hesi2.dat'
          over_power = 1
          use_schwartz = 1     ; use schwartz equations labels, etc.
          mtitle = m2title
          if (keyword_set(allpower)) then all_power = 1
          if (keyword_set(label)) then powerlabs = 1
endif
if (keyword_set(ramaty)) then begin
          over_plot=2
          over_filename = '[eric.hesi_read]hesi2.dat'
          over_power = 1
          mtitle = m1title
          ramaty_plot = 1
          if (keyword_set(powercomp)) then powcow = 1
          if (keyword_set(label)) then ramatylabs = 1
endif else ramaty = 1         
if (keyword_set(bic_ram)) then  begin
          over_plot=2
          over_filename = '[eric.hesi_Read]hesi2.dat'
          bicolor_ram = 1
          mtitle = m1title
          if (keyword_set(label)) then ramatylabs = 1
endif
if (keyword_set(schwartz)) then begin
          if (keyword_set(label)) then powerlabs = 1
          mtitle = m2title
	  schw_pl = 1
          over_filename = '[eric.hesi_Read]hesi2.dat'
endif
if (keyword_set(brokep)) then broken_power = 1 

read_loop:
; read in file hesi.dat

if (over_plot ne 2) then begin

  hesi_plot = dblarr(2,5000)
  i = 0
  openr,lun,filename,/get_lun

  while (not eof(lun)) do begin
 	readf,lun, hesx, hesy
        hesi_plot(0,i) = hesx
        hesi_plot(1,i) = hesy
        if (i lt 4999) then i=i+1 else goto,error_exit
  endwhile

  hesi_plot = hesi_plot(*,0:i)
  close,lun
  free_lun,lun
;  help, hesi_plot
;  print,hesi_plot(0,0:10)  
;  print,hesi_plot(1,0:10)  

endif
; read in file hesi2.dat, contains ramaty positron + nuclear lines

if (over_plot gt 0 or schw_pl eq 1) then begin

    hesi_oplot = dblarr(2,5000)
    i = 0
    openr,lun,over_filename,/get_lun

    while (not eof(lun)) do begin
	readf,lun, hesox, hesoy
        hesi_oplot(0,i) = hesox
        hesi_oplot(1,i) = hesoy
        if (i lt 4999) then i=i+1 else goto,error_exit
    endwhile

    hesi_oplot = hesi_oplot(*,0:i-1)
    close,lun
    free_lun,lun
;    help, hesi_oplot
;    print,hesi_oplot(0,0:10)  
;    print,hesi_oplot(1,0:10)  

endif

power_plot:

; write arrays for nonthermal bremsstrahlung (powerlaw), thermal
; bremsstrahlung (thermal), and composite (addtiion of powerlaw/thermal
; and ramaty lines)

if (over_plot eq 2 or schw_pl eq 1) then begin

    nelem_hesi2 = n_elements(hesi_oplot(0,*))
    en_all = fltarr(nelem_hesi2)
    en_all(*) = hesi_oplot(0,*)
;    help,en_all
    en1_index = where(en_all lt .05, ct_en1)
    en1 = en_all(en1_index)
    en1a = en_all(en1_index(ct_en1-1)+1:*)
    Flux_1 = 5.94e11*En1^(-1.4)*(exp(-en1/.002))  ; Thermal Bremss.
    Flux_1a = fltarr((nelem_hesi2-ct_en1))        ; EMISS. = 10^48 cm^-3
    Flux_1a(*) = 0.0				; T =2*10^7 K t = 100s
    flux_1tot = [flux_1,flux_1a]
    Flux_4 = 400.*En_all^(-2.4)

    pl_flux = flux_1tot + Flux_4 + hesi_oplot(1,*)

    if (use_schwartz eq 1) then begin
      EM = en_all*1000.
      FTHERM = Flux_1/1000.
      FCOMP = pl_flux/1000.
      FNUCL = hesi_oplot(1,*)/1000.
      FBPOW = flux_4/1000.
    endif

endif

; create smaller arrays for plots which use polyfill to fill
; in solid color under each line

if (bicolor_Ram eq 1) then begin
    bic_index = where(en_all gt .1 and en_all le 8.496, bic_ct)
    en_bic = en_all(bic_index)
    fl_bic = fltarr(bic_ct)
    fl_bic(*) = hesi_oplot(1,bic_index)
    en_bic = [.1,en_bic]
    fl_bic = [164.,fl_bic]
endif

;plot_it:

; initialize colors
if do_color eq 2 then do_col=0 else do_col = do_color
init_linecolors,col,do_color=do_col

;if (do_color ne 2) then orange = 4 else orange = 135
;if (do_color ne 2) then red = 2 else red = 90
;if (do_color ne 2) then blue = 11 else blue = 120
;if (do_color ne 2) then green = 8 else green = 50
;if (do_color ne 2) then pink = 3 else pink = 99
;if (do_color ne 2) then purple = 13 else purple = 185
;if (do_color ne 2) then olive = 6 else olive = 78
;if (do_color ne 2) then yellow = 5 else yellow = 195
;if (do_color ne 2) then maroon = 1 else maroon = 25
;if (do_color ne 2) then cyan = 9 else cyan = 130
;if (do_color ne 2) then magenta = 12 else magenta = 105
;if (do_color ne 2) then black = 0 else black = 5
;if (do_color ne 2) then dark_bl = 11 else dark_bl = 15

; initialize xmins/maxs and ymins/maxs
if (over_plot eq 0 or ramaty_plot eq 1 or bicolor_ram eq 1) then begin
  ymin = 1.
  ymax = 10.^5.
  if (ramaty_plot eq 0 and bicolor_Ram ne 1) then xmax= 100. else xmax = 10.
  xmin = .1
endif else begin
  print,'using default mins for overplotting'
  ymin = .01 ;  min(hesi_plot(1,*))
  xmin = .001    ;min(hesi_plot(0,*))
  if(over_plot eq 2) then ymax = (max(pl_flux) < 1.e20) else ymax=10e5
  xmax = 100.    ; max(hesi_plot(0,*))
endelse

; Plot Richard Schwartz's calculated Thermal and Powerlaw Spectrums with
; Ramaty's positron/nuclear lines
if (schw_pl eq 1 or use_schwartz eq 1) then begin
    ; overwrites any previous values for title, calculations, etc.
    Comp_spec_plot,EM,FCOMP,FTHERM,FBPOW,FNUCL,$
    mtitle,xtitle,ytitle,col,$
    use_schwartz,$
    broken_power,en1,do_color,yell,nfill,comp,ramaty,pi_not,prfont,$
    therm1,therm2,power,nuke
;yellow,orange,purple,red,green,cyan,
    goto,labeling
endif

; plot the spectrum given by hesi.dat
if (over_plot ne 2 and schw_pl eq 0) then $
    plot_oo,(hesi_plot(0,*)>0.),(hesi_plot(1,*)>0.),$
    yrange=[ymin,ymax],xrange=[xmin,xmax],charsize=2,$
    charthick=4,thick=4,xthick=3,ythick=3

; overlay the Ramaty Positron/Nuclear lines from hesi2.dat
if (over_plot eq 1) then oplot,(hesi_oplot(0,*)>0.),(hesi_oplot(1,*)>0.),$
    color=col.red,thick=4

; For the Ramaty lines overlayed by a powerlaw, composite (of powerlaw, ramaty,
; and/or thermal spectrum), and/or thermal spectrum set up the plot window
; with appropriate titles and margins
if ((over_plot eq 2 and all_power ne 1) or bicolor_ram eq 1) then begin 
    plot_oo,[.00001],[.000001], font=prfont,$
    title=mtitle,xtitle=xtitle,ytitle=ytitle,yrange=[ymin,ymax],$
    xrange=[xmin,xmax],ytick_get=ytick,xtick_get=xtick,$
    charsize=2,charthick=4,thick=4,xthick=3,ythick=3 ;ystyle=1
    if (bicolor_ram ne 1) then oplot,en_all, pl_flux,thick=4   ;,color=col.cyan

endif
   ; if the all_power option is on, then set up to plot only the 
   ; Ramaty lines, thermal spectrum, and powerlaw (otherwise
   ; the composite line which adds the others together will be
   ; included)

if (over_power eq 1 and ramaty_plot ne 1) then begin
   if(all_power eq 1) then plot_oo,[.00001],[.00001],title=mtitle,$
       xtitle=xtitle,ytitle=ytitle,$
       yrange=[ymin,ymax],xrange=[xmin,xmax],charsize=1.5,charthick=3,$
       ystyle=1,thick=4,xthick=3,ythick=3
   oplot,en_all,hesi_oplot(1,*),color=col.red,thick=5
   oplot,[en1,.11],[flux_1,1.e-5],color=col.orange,thick=4
   if (bicolor_ram ne 1) then oplot,en_all, pl_flux,thick=4   ;,color=col.cyan
endif

; if ramaty option is on, then plot the Ramaty, power, and composite along
; with labels for the positron/nuclear lines
if (ramaty_plot eq 1) then begin
   if(powcow ne 1) then oplot,en_all,hesi_oplot(1,*),color=col.red,thick=5
   if(broken_power eq 1) then oplot,em,fbpow,color=col.green,thick=4 else $
      oplot,en_all,flux_4,color=col.green,thick=4
endif                           

; create the ramaty plot immediately above without labels and 
; blocked in solid colors underneath each plotted line
if (bicolor_ram eq 1) then begin
   xfill1 = [.1,en_bic,en_bic(bic_ct-1),.1]
   yfill1 = [1.,fl_bic,1.,1.]
   xfill2 = [.1,.1,10.,10.]
   yfill2 = [1.,1.e5,1.e5,1.]
   gr_index = where(en_all ge .1 and en_all le 10.,gr_ct)
   yf_plf = pl_flux(gr_index)
   xfill3 = [.1,.1,en_all(gr_index),10.,.1]
   yfill3 = [1.,1.e5,flux_4(gr_index),1.,1.]
   xfill4 = [.1,.1,en_all(gr_index),10.0,.1]
   yfill4 = [1.,1.e5,yf_plf,1.,1.]

   if (yell eq 1) then polyfill,/normal,[0,0,1,1,0],[0,1,1,0,0],color=col.yellow

   polyfill,xfill4,yfill4,color=col.purple,/data,clip=[.1,1,10.,1.e5]
   polyfill,xfill3,yfill3,color=col.green,/data,clip=[.1,1,10.,1.e5]
   polyfill,xfill1,yfill1,color=col.cyan,/data,clip=[.1,1,10.,1.e5]
   oplot,en_bic,fl_bic,thick=2,color=col.red
   oplot,en_all,flux_4,thick=2,color=col.green
   oplot,en_all,pl_flux,thick=2

    plot_oo,[.00001],[.000001],/noerase, font=prfont,$
    title=mtitle,xtitle=xtitle,ytitle=ytitle,yrange=[ymin,ymax],$
    xrange=[xmin,xmax],ytick_get=ytick,xtick_get=xtick,$
    charsize=2,charthick=4,ystyle=1,thick=4,xthick=3,ythick=3

endif

Labeling:
if (ramatylabs eq 1 or powerlabs eq 1) then begin
        if (bicolor_ram eq 1 or schw_pl eq 1) then color1 = col.red else $
            color1 = col.red
	if (ramatylabs eq 1) then begin
           if (simp eq 1) then begin 
            poslab = 'Gammay-Rays' 
            plabx = 1.0
           endif else begin
            poslab = 'Positron and Nuclear Gamma-Ray Lines'
            plabx = .3
           endelse
	   xyouts,plabx,4,poslab,color=col.red,charthick=4,charsize=1.2
           if (simp eq 1) then nonlab = 'X-Rays' else $
            nonlab = 'Nonthermal Bremsstrahlung!c(400E!u-2.4!n)'
	   xyouts,.2,5.e4,nonlab,color=col.green,charthick=4,charsize=1.2
	   if (simp ne 1 and powcow ne 1) then xyouts,.2,1.3e2,'Positronium',charthick=4,charsize=1.2,color=color1
           if (powcow eq 1) then xyouts,.17,1.7e3,'Positronium',charthick=4,charsize=1.2,color=color1
	   !p.font = 1
           if (powcow ne 1) then begin
            xyouts,.135,7.e2,'!17!u56!nCo',charthick=3,color=col.red,charsize=1.5,font=-1
  	    xyouts,.28,1.0e3,'!u59!nNi!X',charthick=3,color=col.red,charsize=1.5,font=-1
            xyouts,.405,1.25e3,'!4a!x!4a!X',charthick=4,color=col.red,charsize=1.5,font=-1
	   endif else xyouts,.405,1.e4,'!4a!x!4a!X',charthick=4,color=col.red,charsize=1.5,font=-1
	   xyouts,.509,2.e4,'!17e!u+!n',charthick=3,color=col.red,charsize=1.5,font=-1
	   xyouts,.7,2.e3,'!u56!nFe',charthick=3,color=col.red,charsize=1.5,font=-1
           xyouts,1.4,1.1e3,'!u20!nNe',charthick=3,color=col.red,charsize=1.5,font=-1
	   xyouts,2.2,3.e4,'n',charthick=3,color=col.red,charsize=1.5,font=-1
	   xyouts,4,3.e2,'!u12!nC',charthick=3,color=col.red,charsize=1.5,font=-1
	   xyouts,5.3,1.7e2,'!u16!nO!X',charthick=3,color=col.red,charsize=1.5,font=-1
	   !p.font = -1
	endif
	if (powerlabs eq 1) then begin
	   if (schw_pl ne 1 and use_schwartz ne 1) then $
               xyouts,.04,1.e11, color=col.orange,charthick=4,$
	      'Thermal Bremsstrahlung!c6.0x10!u11!nE!u-1.4!nexp(-E/2.0keV)!cT = 2x10!u7!nK !cEM = 10!u50!ncm!u-3!n !c(Duration = 100s)'
           if (schw_pl eq 1 or use_schwartz eq 1) then begin
               xyouts,30.,5.e9,charthick=4,charsize=1.2,$
                'Thermal Bremsstrahlung',color=col.purple
               xyouts,30.,5.e8, color = purple,charthick=4,charsize=1.2,$
	     'T = 2x10!u7!nK, EM = 10!u50!ncm!u-3!n'
              ; xyouts,30.,1.5e7,color=col.purple,charthick=4,charsize=1.2,$
	     ;'T = 3.4x10!u7!nK !cEM = 10!u50!ncm!u-3!n'
           endif
	   if (schw_pl eq 1 or use_schwartz eq 1) then $
            cnon = 'Nonthermal Bremsstrahlung' else $
            cnon = 'Nonthermal Bremsstrahlung!c(400!nE!u-2.4!n)'
           cpos = 'Positron and Nuclear!c Gamma-Ray Lines'
	   if (schw_pl eq 1 or use_schwartz eq 1) then xpos = 200. else xpos = .09
           if (schw_pl eq 1 or use_schwartz eq 1) then ypos = 1.e-3 else ypos = .125
           xyouts,xpos,ypos,''+cpos,color=color1,charthick=4,charsize=1.2 
	   if (schw_pl eq 1 or use_schwartz eq 1) then xnon = 200. else xnon = .15
           if (schw_pl eq 1 or use_schwartz eq 1) then ynon = 1.e3 else ynon = 1.e6
	   xyouts,xnon,ynon,''+cnon,color=col.green,charthick=4,charsize=1.2
	endif
endif

hard_copy:

if (hd_cp eq 'y' or hd_cp eq 'Y') then begin
    ; print thru net$print if do_color=2
    if (do_color eq 2) then begin
        tek_end
        tek_print
        hd_cp = 'n'
        if keyword_set(tek) then alpha_page
        goto,good_End
     endif
     psplot,color=do_color 
     if keyword_set(tek) then alpha_page
     goto,good_end
endif else begin
    read,'hard copy (y/n)?  ',hd_cp
	if (hd_cp eq 'y' or hd_cp eq 'Y') then begin
          if (over_plot gt 0 or schw_pl eq 1) then begin
             print,'Color Printing Key:  0 = B/W only, 1 = Color, 2 = Psuedogray'
             read,'color added (0/1/2)?  ',do_color
          endif
          prfont = 0
          linecolors

          ; print through net$print w/out using psplot if do_color=2
          if (do_color eq 2) then begin
              set_graphics,printer='ps'
              set_plot,'ps',/copy
              device,/landscape,xoff=3,yoff=25,xsize=22,ysize=16,/times
              tek_init
              goto,plot_it
          endif

          set_plot,'ps',/copy
          device,color=do_color,/landscape,bits=8,xoff=3,yoff=25,$
          xsize=nxsize,ysize=nysize,/times,/bold
	  goto,plot_it
	endif else begin
		if keyword_set(tek) then alpha_page
		goto,good_end
	endelse
endelse

error_exit:
print, 'Error in Comp_spec'

good_end:
over_plot=0
over_power=0
all_power = 0
ramaty_plot = 0
bicolor_ram = 0
schw_pl = 0
powerlabs = 0
ramatylabs = 0
return
end