pro bigplot_v9b

; dhc - 8jun14 - based on bigplot_v8.pro in, for example,
;       xspec/2014/t1oric/

; v.9b reads in two files and adds them together

;
; reads in two four-column ascii files: 

; col 1 = wavelength
; col 2 = bin half-width 
; col 3 = data
; col 4 = uncertainty
;
; this version of the code reads in TWO files, assuming no header
;    -  see skiplines in transread call 
;

; AXIS RANGES 

xrdata=[4.5,22.8]
yrdata=[0.0,0.23]

; INPUT FILE NAMES - these two spectra will be coadded; they must have
;                    exactly the same binning

fname1 = 'zpup_meg-1_dataonly.dat'
fname2 = 'zpup_meg+1_dataonly.dat'

; observation-specific exposure time

exptime=67636. ; for zeta pup

; exposure time not currently in use 

; list of lab wavelengths of lines - user can edit

lines=[5.0387,5.0648,5.1015,6.1821,6.6479,6.685,6.7403,7.85,8.421,9.1688,9.2308,9.3143,9.7082,10.239,11.544,12.1339,15.014,16.780,17.051,17.096,18.9689,20.9099,21.602,21.804,22.098]

num_lines=n_elements(lines)

print, 'number of lines =',num_lines
print, 'wavelengths =',lines

; SWITCH FOR DRAWING THREE VERTICAL LINES

ivert=0

; terminal velocity in km/s, used only when drawing two extreme
; velocity vertical lines around each marked line 

vinf=2250.
c=3.e5

;
; READ IN THE FIRST FILE AND POPULATE ARRAYS

wave1=0.0d
width1=0.0d
data1=0.0d
error1=0.0d

wave2=0.0d
width2=0.0d
data2=0.0d
error2=0.0d

transread,unit,wave1,width1,data1,error1,count=count1,filename=fname1
transread,unit,wave2,width2,data2,error2,count=coun2,filename=fname2

; this block of code and the variable naming convention allow for
; intuitive extension in a version that coadds two or more spectra

wave_check = where((wave1-wave2),count)
print, 'NUMBER OF ELEMENTS FOR WHICH WAVELENGTH ARRAYS DISAGREE = ',count

wave=wave1
width = width1

data = data1+data2
error = error1+error2

; data and error (and eventual model) are in counts/s/A (it is
; assumed)

; note that native bin width is 2.5 mA for HEG and 5 mA for MEG
; the program can just as easily accept (unevenly) grouped data

; MAKE THE PLOT

; code to automatically choose x and y axis ranged 
; turned off/commented out currently

;xrdata=[min(wave),max(wave)]
;yrdata=[0.,(max(data)+max(error))]

xtit = 'Wavelength ('+STRING("305B)+')'
ytit = 'Count Rate (cnt s!u-1!n '+string("305B)+'!u-1!n)'

; loop through screen plotting and ps-file writing

for i=1,2 do begin

if (i eq 1) then begin
window

; setting up colors for .ps using coyote graphics

blue = cgColor("blue")
red = cgColor("red")
green = cgColor("green")

endif else begin
psopen,ps,/color

; setting up colors for .ps

blue = cgColor("blue")
red = cgColor("red")
green = cgColor("green")

!p.font=0
device,/times
device,/inches,xsize=35.0,ysize=8.0

endelse

plot,wave,data,ps=10,thick=4,xr=xrdata,yr=yrdata,xstyle=1,ystyle=1,ytit=ytit,chars=3.8,xtit=xtit,xtickinterval=1

; plot error bars 
;oploterr,wave,data,error,10

; plot model
;oplot,wave,model,ps=10,linestyle=1,thick=4,color=red

; PLOTTING VERTICAL LINES AT LINE CENTER AND AT THE TWO EDGES

if ivert EQ 1 then begin

for j=0,num_lines-1 do begin

lam_o=lines(j)

x_center=[lam_o,lam_o]
y_center=[0.,1.]

lam_left=lam_o*(1.-vinf/c)
lam_right=lam_o*(1.+vinf/c)

x_left=[lam_left,lam_left]
y_left=y_center

x_right=[lam_right,lam_right]
y_right=y_center

   oplot,x_center,y_center,linestyle=2,thick=0.75

; next two lines should be uncommented if you want flanking v_inf
; lines too

;   oplot,x_left,y_left,linestyle=1,thick=0.5
;   oplot,x_right,y_right,linestyle=1,thick=0.5

endfor

endif

endfor

if (i EQ 1) then begin
print, i
endif else begin
psclose
endelse

end