pro specmod_counts

; 7 MAY 2004 - DHC
; 
; MODIFIED FROM SPECMOD.PRO (6MAY04 VERSION) 
; PLOTS COUNTS/BIN RATHER THAN FLUX UNITS
; NO OPTION FOR MODEL PLOTTING 
; NO OPTION FOR RESIDUALS WINDOW
; CUSTOMIZED FOR ZETA ORI - RELEVANCE IS TWO SEPARATE OBS ID
; 	FILES, BUT -1 +1 ALREADY COADDED
; 

; PRIMARILY, THIS PROGRAM IS USED TO MAKE (LINE) LABELLED PLOTS


; Parameters to change between runs:

ps ='specmod_counts'	;Name of the ps plot

; Cosmetics:

name = 'zeta Ori MEG +/-1, coadded ObsID 1524 + 610'


; zeta Ori 

exptime1=59631.   ;obs610
exptime2=13783.   ;obs1524


; CALCULATE POISSON ERRORS FROM TOTAL COUNTS PER BIN?

ipoisson = 1 	; if 1, then calculates poisson errors from counts


;NOTE: CURRENTLY (7/8/03) THE ONLY WAY TO GET THE 
;RESIDUALS PLOT IS TO CHOOSE PLOTMOD=1

ploterr=1	;Plot uncertainty? 0 if no, 1 if yes

titleloc=0	;Include title? 0 if no, 1 if the title 
		;should go at the top, 2 if at the 
		;bottom (because the line is in the way)

; RANGE OF X-AXIS

first=24.2
last=25.3
xrdata=[first,last]



;SETTING INPUT FILES

;TWO INPUT DATA FILES (MEG+1 AND MEG-1, TYPICALLY)
filename1a = '/home/cohen/chandra/idl/zori/*.dat'
file1a = dialog_pickfile(/read, filter = filename1a)
filename1b = '/home/cohen/chandra/idl/zori/*.dat'
file1b = dialog_pickfile(/read, filter = filename1b)


;READING IN DATA FROM FILES1a and 1b
width=6
filler=fltarr(width)
empty=' '

;To determine the length of File1a:
length = 0
openr,unit,file1a,/get_lun,/compress
while not eof(unit) do begin
	length = length+1
	if (length LE 1) then readf,unit,empty $
		else readf,unit,filler
endwhile
free_lun,unit
length = length-1
print,'Length of file1a is ',length

;Read in the data:
array = fltarr(width,length)
openr,unit,file1a,/get_lun
for tag=1,1 do readf,unit,empty
readf,unit,array
free_lun,unit

energya = reform(array[0,*])
countsa = reform(array[1,*])
uncerta = reform(array[4,*])

;IMPLEMENTING POISSON ERROR OPTION

;8/4/03 DHC: GETTING BIN SIZE OF DATASET 1
;SPOT CHECK (ASSUMES IMPLICITLY THE BINNING IS CONSTANT)
;AND THE BINNING BETWEEN THE TWO DATASETS IS CHECKED BELOW

binwidth = energya[1]-energya[2]
print, 'binwidth =', binwidth

if (ipoisson EQ 1) then begin
countsa = binwidth*exptime1*countsa
;uncerta = sqrt(countsa)
endif

;To determine the length of File1b:
length = 0
openr,unit,file1b,/get_lun,/compress
while not eof(unit) do begin
	length = length+1
	if (length LE 1) then readf,unit,empty $
		else readf,unit,filler
endwhile
free_lun,unit
length = length-1
print,'Length of file1b is ',length

;Read in the data:
array = fltarr(width,length)
openr,unit,file1b,/get_lun
for tag=1,1 do readf,unit,empty
readf,unit,array
free_lun,unit

energy = reform(array[0,*])
counts = reform(array[1,*]) 
uncert = sqrt(reform(array[4,*])^2 + uncerta^2)


if (ipoisson EQ 1) then begin
counts = binwidth*exptime2*counts
;uncert = sqrt(counts)
endif

counts = counts + countsa
phcounts = counts
uncert_phcounts = sqrt(counts)

; VARIABLE PHCOUNTS DEFINED TO KEEP POSSIBILITY OF PLOTTING
; DATA AS COUNTS/BIN

if (ipoisson EQ 1) then begin
uncert = sqrt(counts)
counts = counts/binwidth/(exptime1+exptime2)
uncert = uncert/binwidth/(exptime1+exptime2)

if (exptime1 EQ exptime2) then begin
print, 'exptime1 eq exptime2'
counts = counts*2.
uncert = uncert*2.
endif

endif

;NOTE THAT THE UNCERTAINTIES ARE ADDED IN QUADRATURE; THE COUNT RATE
;DATA ARE SIMPLY ADDED, AND THE WAVELENGTHS (STORED IN A VARIABLE CURRENTLY
;CALLED 'ENERGY') ARE TAKEN FROM THE SECOND FILE
;THOUGH WE DO A SPOT CHECK HERE TO MAKE SURE THAT THE FIRST AND LAST 
;ELEMENTS IN THESE ARRAYS ARE THE SAME

;THIS IS ALL OK AS LONG AS WE'RE ADDING TWO ORDERS FROM THE SAME 
;OBSERVATION; BUT IF OBSERVATIONS OF DIFFERENT LENGTH ARE BEING ADDED
;OR WE WANT POISSON ERRORS CALCULATED DIRECTLY FROM THE COUNTS, THEN 
;THE OBSERVATION EXPOSURE TIMES MUST BE USED - ipoisson=1 
;NOTE: IF YOU WANT TO EXPLICITLY SPECIFY EXP TIMES (I.E. ADD 2 SEP OBS)
;THEN YOU MUST USE THE POISSON OPTION

if (energy[0] NE energya[0]) then begin
	print, 'WARNING: two data wavelength arrays are not equal'
endif

if (energy[length-1] NE energya[length-1]) then begin
	print, 'WARNING: two data wavelength arrays are not equal'
endif



;/////////////////////////////////////////


;to determine the x range:

;	first=energy[length-1]
;        last=energy[0]


;NOTE: THIS IS THE ENTIRE RANGE; AND THE ABOVE LINE OF CODE
;REVERSES THE ORDER-IT ASSUMES THAT THE WAVELENGTHS ARE 
;TABULATED IN DECREASING ORDER IN THE INPUT FILES
;
; 7MAY04 - DHC - FIRST AND LAST ARE NOW MANUALLY SET AT THE TOP 
; OF THE PROGRAM

;xrdata=[first,last]

yrdata = [min(phcounts),max(phcounts)]

;THE USER CAN EASILY HARDWIRE CHANGES TO THE X- OR Y-RANGES
; AND GENERALLY WILL - AT THE TOP OF THE PROGRAM (3AUG03)

;THESE ARE AXIS LABELS
	xtit = 'Wavelength ('+STRING("305B)+')'
;	ytit = 'Count Rate (counts s!u-1!n '+string("305B)+'!u-1!n)'

; 7MAY04 - DHC - ASSUMING THAT Y-VALS ARE COUNTS/BIN

	ytit = 'Counts'


; FIRST SEND PLOT TO IDL DISPLAY WINDOW

window

plot,energy,phcounts,ps=10,xr=xrdata, $
	ytit=ytit,chars=1.5,thick=2,xs=1,xtit=xtit


xcoor = (last-first)/2 + first

if (titleloc EQ 1) then xyouts,xcoor, $
	(yrdata[1]*11./10.),name,alignment=0.5,chars=1.5

if (titleloc EQ 2) then xyouts,xcoor, $
	-0.1,name,alignment=0.5,chars=1.5



if (ploterr EQ 1) then oploterr,energy,phcounts,uncert_phcounts,1


; 7MAY04 - DHC - WRITING LINE LABELS


name1 = 'O VII'
name2 = 'O VIII'
name3 = 'Fe XVIII'
name4 = 'Fe XVII'
name5 = 'Ne X'
name6 = 'Ne IX'
name7 = 'Mg XI'
name8 = 'Mg XII'
name9 = 'Si XIII'
name10 = 'Fe XIX'
name11 = 'N VII'

; FOR PANEL 2 OF THE ATLAS (13:25 A)

;xyouts,13.5,38,name6,alignment=0.5,chars=1.25
;xyouts,13.45,36.,'|',alignment=0.5,chars=1
;xyouts,13.55,36.,'|',alignment=0.5,chars=1

;xyouts,15.13,39.5,name4,alignment=0.5,chars=1.25
;xyouts,15.01,37.5,'|',alignment=0.5,chars=1
;xyouts,15.261,37.5,'|',alignment=0.5,chars=1

;xyouts,16.01,28.5,name2,alignment=0.5,chars=1.25
;xyouts,16.01,26.5,'|',alignment=0.5,chars=1

;xyouts,16.45,24.,name4,alignment=0.5,chars=1.25
;xyouts,16.787,22.,'|',alignment=0.5,chars=1
;xyouts,16.07,22.,'|',alignment=0.5,chars=1

;xyouts,17.0735,41.,name4,alignment=0.5,chars=1.25
;xyouts,17.05,39.,'|',alignment=0.5,chars=1
;xyouts,17.095,39.,'|',alignment=0.5,chars=1

;xyouts,18.4,12.,name1,alignment=0.5,chars=1.25
;xyouts,18.627,10.,'|',alignment=0.5,chars=1

;xyouts,19.5,41.5,name2,alignment=0,chars=1.25
;xyouts,19.15,41.6,'-',alignment=0,chars=1

;xyouts,21.7,24.,name1,alignment=0.5,chars=1.25
;xyouts,21.6,22.,'|',alignment=0.5,chars=1
;xyouts,21.8,22.,'|',alignment=0.5,chars=1

;xyouts,24.78,10.,name11,alignment=0.5,chars=1.25
;xyouts,24.78,8.,'|',alignment=0.5,chars=1

; FOR PANEL 1 OF THE ATLAS (6:13 A)

xyouts,6.68,20.,name9,alignment=0.5,chars=1.25
xyouts,6.62,18.,'|',alignment=0.5,chars=1
xyouts,6.74,18.,'|',alignment=0.5,chars=1

xyouts,8.42,15.5,name8,alignment=0.5,chars=1.25
xyouts,8.42,13.5,'|',alignment=0.5,chars=1

xyouts,9.26,32.5,name7,alignment=0.5,chars=1.25
xyouts,9.17,30.5,'|',alignment=0.5,chars=1
xyouts,9.34,30.5,'|',alignment=0.5,chars=1

xyouts,11.544,15.,name6,alignment=0.5,chars=1.25
xyouts,11.544,13.25,'|',alignment=0.5,chars=1

xyouts,12.134,47.,name5,alignment=0.5,chars=1.25
xyouts,12.134,45.,'|',alignment=0.5,chars=1



;NEXT, SEND SAME PLOT TO PS FILE
;NOTE SLIGHTLY DIFFERENT CHOICES ABOUT LINE THICKNESS

psopen,ps
!p.font=0
device,/times

plot,energy,phcounts,ps=10,xr=xrdata, $
	ytit=ytit,chars=1.5,thick=4,xs=1,xtit=xtit
if (ploterr EQ 1) then oploterr,energy,phcounts,uncert_phcounts,1

xcoor = (last-first)/2 + first

if (titleloc EQ 1) then xyouts,xcoor, $
	(yrdata[1]*11./10.),name,alignment=0.5,chars=1.5

if (titleloc EQ 2) then xyouts,xcoor, $
	-0.1,name,alignment=0.5,chars=1.5

; 7MAY04 - DHC - WRITING LINE LABELS


name1 = 'O VII'
name2 = 'O VIII'
name3 = 'Fe XVIII'
name4 = 'Fe XVII'
name5 = 'Ne X'
name6 = 'Ne IX'
name7 = 'Mg XI'
name8 = 'Mg XII'
name9 = 'Si XIII'
name10 = 'Fe XIX'
name11 = 'N VII'

; FOR PANEL 2 OF THE ATLAS (13:25 A)

;xyouts,13.5,38,name6,alignment=0.5,chars=1.25
;xyouts,13.45,36.,'|',alignment=0.5,chars=1
;xyouts,13.55,36.,'|',alignment=0.5,chars=1

;xyouts,15.13,39.5,name4,alignment=0.5,chars=1.25
;xyouts,15.01,37.5,'|',alignment=0.5,chars=1
;xyouts,15.261,37.5,'|',alignment=0.5,chars=1

;xyouts,16.01,28.5,name2,alignment=0.5,chars=1.25
;xyouts,16.01,26.5,'|',alignment=0.5,chars=1

;xyouts,16.45,24.,name4,alignment=0.5,chars=1.25
;xyouts,16.787,22.,'|',alignment=0.5,chars=1
;xyouts,16.07,22.,'|',alignment=0.5,chars=1

;xyouts,17.0735,41.,name4,alignment=0.5,chars=1.25
;xyouts,17.05,39.,'|',alignment=0.5,chars=1
;xyouts,17.095,39.,'|',alignment=0.5,chars=1

;xyouts,18.4,12.,name1,alignment=0.5,chars=1.25
;xyouts,18.627,10.,'|',alignment=0.5,chars=1

;xyouts,19.5,41.5,name2,alignment=0,chars=1.25
;xyouts,19.15,41.6,'-',alignment=0,chars=1

;xyouts,21.7,24.,name1,alignment=0.5,chars=1.25
;xyouts,21.6,22.,'|',alignment=0.5,chars=1
;xyouts,21.8,22.,'|',alignment=0.5,chars=1

;xyouts,24.78,10.,name11,alignment=0.5,chars=1.25
;xyouts,24.78,8.,'|',alignment=0.5,chars=1

; FOR PANEL 1 OF THE ATLAS (6:13 A)

xyouts,6.68,20.,name9,alignment=0.5,chars=1.25
xyouts,6.62,18.,'|',alignment=0.5,chars=1
xyouts,6.74,18.,'|',alignment=0.5,chars=1

xyouts,8.42,15.5,name8,alignment=0.5,chars=1.25
xyouts,8.42,13.5,'|',alignment=0.5,chars=1

xyouts,9.26,32.5,name7,alignment=0.5,chars=1.25
xyouts,9.17,30.5,'|',alignment=0.5,chars=1
xyouts,9.34,30.5,'|',alignment=0.5,chars=1

xyouts,11.544,15.,name6,alignment=0.5,chars=1.25
xyouts,11.544,13.25,'|',alignment=0.5,chars=1

xyouts,12.134,47.,name5,alignment=0.5,chars=1.25
xyouts,12.134,45.,'|',alignment=0.5,chars=1




psclose

!P.MULTI = 0	;restores the default: 1 plot per page

end