;To run program:
; .r spectrum_plotting.pro

;-----------------------------------------------------------------------
;This program reads in a specific fits file defined in the code
;It reads in the wavelength and flux for each pixel and for each
;order.  There are 5 orders, starting with 3, thus wave3 is the
;wavelength of the pixels in the 3rd order and flux3 is the flux value
;of those pixels in the 3rd order.  

;The notation [*,0,0] refers to the matrix of data of the FITS
;files. It is a 3-D array of data consisting of sets of triplet arrays
;of data for each aperture and each order.  Each triplet is composed
;of an array for wavelength, flux, and error.

;If no orders have been skipped or deslected in the extraction
;process, the contents of aperture Y in order X can be found as
;follows:
;lambda = array[*,0,( X - (smallest order number))*naps + (Y-1)]
;flux   = array[*,1,( X - (smallest order number))*naps + (Y-1)]
;error  = array[*,2,( X - (smallest order number))*naps + (Y-1)]

;For example, for an SXD file with two apertures, the wavelengtharray
;for aperture 1 in order 3 is located in array [*,0,0], the flux is in
;array [*,1,0] and the error is in array [*,2,0]. For aperture 2, the
;wavelength is in array [*,0,1], the flux is in array [*,1,1] and the
;error is in array [*,2,1]. For order 4, the flux for aperture 1 is in
;array [*,1,2], while the flux for aperture 2 is in array [*,1,3]. For
;order 5, the fluxes for the two apertures are in arrays [*,1,4] and [*,1,5]

;I am interested in the wavelength of 1.27um which is in the 5th
;order (wave5, flux5) and is the only one plotted.
;---------------------------------------------------------------------------

path='/tharsis3/candaceg/Venus/Dec2010/Q4NM05/UT101221/tspec/Reduced2/'
file='Venus.0013.REDUCED.16aprt.sky0019.fits'
name=' '
name=strcompress(path+file,/remove_all)
spec=readfits(name)

wave3=spec[*,0,0]
flux3=spec[*,1,0]

wave4=spec[*,0,1]
flux4=spec[*,1,1]

wave5=spec[*,0,2]
flux5=spec[*,1,2]

wave6=spec[*,0,3]
flux6=spec[*,1,3]

wave7=spec[*,0,4]
flux7=spec[*,1,4]

;-----------------------------------------------------------------------
; Used 18" of the 45" shown.  But APO says the slit is 43".  If that's
; true and xspextool display is wrong, then the actual length is
; 18*(43/45)=17.2"
; 1.1*17.2=18.92 sqare arc sec
; 18.92/(3600^2) = 1.33e-6 square degrees
; converting to radians:
; 1.33e-6 / (180/pi)^2 = 4.04e-10
;-----------------------------------------------------------------------


set_plot,'X'
plot,wave5,flux5/4.1e-10,xstyle=1,xrange=[1.26,1.29],yrange=[4,8],ytit='Flux (W m!E-2 !Nsr!E-1 !Nµm!E-1!N)',xtit='Wavelength (µm)',tit='O!I2 !NEmission from Venus Night Side Dec 27, 2010 (Image 0011)',ystyle=1 

set_plot,'PS'
device, filename='Venus.0013.REDUCED.16aprt.sky0019.ps'
plot,wave5,flux5/4.1e-10,xstyle=1,xrange=[1.26, 1.29],yrange=[4,8],ytit='Flux (W m!E-2 !Nsr!E-1 !Nµm!E-1!N)',xtit='Wavelength (µm)',tit='O!I2 !NEmission from Night Sky Dec 21, 2010 (Image 0013 16" aperture',ystyle=1 
device, /close


;------------------------------------------------------------------------
;plotting wavelength versus flux. The wavelength range is from 1.26 -
;1.29um.  The flux units read out from xspextool should be W m-2 um-1
;but we want sr-1 also.  To get the division by steridian, you have to 
;find how many steridians the slit is on the sky.  The slit used was 
;1.1x43" which is 47.3 square arcsec on the sky converting to degrees,
;this is 3.65e-6 square degrees.  A radian is 180/pi degrees, or about 57.296°.
;A steradian is (180/pi)2 square degrees or about 3282.8 square degrees.
;Therefore 3.65e-6/[(180/pi)2] is 1.11e-9 steridians on the sky.
;Divided the flux values by 4e-6 and you should have values of order 
;0.1 W m-2 sr-1 µm-1.
;------------------------------------------------------------------------


;---------------------------------------------------------------------
;Irrelevent code, but possible good commands to learn in the future.

;ystyle=1,$
;yrange=[0,9e-12]
;plot,wave3,flux3,xstyle=1,xrange=[4.5,5.4],ystyle=1,$
;yrange=[0,1e-6]
;oplot,1.e4/wave3_sb,flux3_sb,linestyle=1
;oplot,1.e4/wave3_sz,flux3_sz,linestyle=2
;oplot,1.e4/wave3_ssb,flux3_ssb,linestyle=3
;oplot,1.e4/wave3_pol,flux3_pol,linestyle=4

;legend,linestyle=[0,1,2,3,4],['EZ','SB','SZ','SSB','POL'],$
;   pos=[0.7,0.9],/norm,box=no

;readcol,'iso/swaa38602541.asc',num,lam,flux,err,a,b
;oplot,1.e4/lam,flux+300

;readcol,'/home/phoebe/nchanove/photom/atmostrans1-6.dat',wave,tran
;oplot,1.e4/wave,(tran*200)+700.

;wait,5

;plot,wave4,flux4,ystyle=1,yrange=[0,1e-12],$
;   xstyle=1,xrange=[3.9,4.2]
;oplot,wave4_sb,flux4_sb,linestyle=1
;oplot,wave4_sz,flux4_sz,linestyle=2
;oplot,wave4_ssb,flux4_ssb,linestyle=3
;oplot,wave4_pol,flux4_pol,linestyle=4
;
;wait,5

;plot,wave5,flux5,ystyle=1,yrange=[0,1e-12]
;;   xstyle=1,xrange=[3.9,4.2]
;oplot,wave5_sb,flux5_sb,linestyle=1
;oplot,wave5_sz,flux5_sz,linestyle=2
;oplot,wave5_ssb,flux5_ssb,linestyle=3
;oplot,wave5_pol,flux5_pol,linestyle=4

;wait,5

;plot,wave6,fluz6
;;   xstyle=1,xrange=[3.9,4.2]
;oplot,wave6_sb,flux6_sb,linestyle=1
;oplot,wave6_sz,flux6_sz,linestyle=2
;oplot,wave6_ssb,flux6_ssb,linestyle=3
;oplot,wave6_pol,flux6_pol,linestyle=4

;wait,5

;plot,wave7,flux7,ystyle=1,yrange=[0,2e-10],$
;   xstyle=1,xrange=[2.6,2.75]
;oplot,wave7_sb,flux7_sb,linestyle=1
;oplot,wave7_sz,flux7_sz,linestyle=2
;oplot,wave7_ssb,flux7_ssb,linestyle=3
;oplot,wave7_pol,flux7_pol,linestyle=4

;wait,2

;openw,1,'junk.dat'
;printf,1,wave3
;printf,1,' '
;printf,1,wave3_pol
;close,1
;---------------------------------------------------------------------------

end