readcol, '/Users/colkin/Work04/Pluto_3um/ATRAN/atran_z29.3.dat', tmp, wav29, tran29
plot, wav29, tran29
readcol, '/Users/colkin/Work04/Pluto_3um/ATRAN/atran_z52.4.dat', tmp, wav52, tran52
oplot, wav52, tran52

nontelluric = where((tran52 gt .97) and (tran29 gt .97))
help, nontelluric

clearWav = wav52[nontelluric]

dirS = '/Users/colkin/Work04/Pluto_3um/Spectra3/'

s22 = readfits(dirS+'12aus0022.fits')
s23 = readfits(dirS+'12aus0023.fits')
s24 = readfits(dirS+'12aus0024.fits')
s25 = readfits(dirS+'12aus0025.fits')
s26 = readfits(dirS+'12aus0026.fits')
s27 = readfits(dirS+'12aus0027.fits')
s28 = readfits(dirS+'12aus0028.fits')
s29 = readfits(dirS+'12aus0029.fits')
s30 = readfits(dirS+'12aus0030.fits')
s31 = readfits(dirS+'12aus0031.fits')
s32 = readfits(dirS+'12aus0032.fits')
s33 = readfits(dirS+'12aus0033.fits')
s66 = readfits(dirS+'12aus0066.fits')
s67 = readfits(dirS+'12aus0067.fits')
s68 = readfits(dirS+'12aus0068.fits')
s69 = readfits(dirS+'12aus0069.fits')
s82 = readfits(dirS+'12aus0082.fits')
s83 = readfits(dirS+'12aus0083.fits')
s84 = readfits(dirS+'12aus0084.fits')
s85 = readfits(dirS+'12aus0085.fits')
 
 w = readfits("/Users/colkin/Work04/Pluto_3um/Rectified6/WavelengthScale.fits")

n22 = s22/total( interpol(s22, w, clearWav) )
n23 = s23/total( interpol(s23, w, clearWav) )
n24 = s24/total( interpol(s24, w, clearWav) )
n25 = s25/total( interpol(s25, w, clearWav) )

n26 = s26/total( interpol(s26, w, clearWav) )
n27 = s27/total( interpol(s27, w, clearWav) )
n28 = s28/total( interpol(s28, w, clearWav) )
n29 = s29/total( interpol(s29, w, clearWav) )

n30 = s30/total( interpol(s30, w, clearWav) )
n31 = s31/total( interpol(s31, w, clearWav) )
n32 = s32/total( interpol(s32, w, clearWav) )
n33 = s33/total( interpol(s33, w, clearWav) )

n66 = s66/total( interpol(s66, w, clearWav) )
n67 = s67/total( interpol(s67, w, clearWav) )
n68 = s68/total( interpol(s68, w, clearWav) )
n69 = s69/total( interpol(s69, w, clearWav) )

n82 = s82/total( interpol(s82, w, clearWav) )
n83 = s83/total( interpol(s83, w, clearWav) )
n84 = s84/total( interpol(s84, w, clearWav) )
n85 = s85/total( interpol(s85, w, clearWav) )

all = [n23, n24, n25, n26, n27, n28, n29, n30, n31, n32, $
n33, n66, n67, n68, n69, n82, n83, n84, n85]

airmass = [1.14401, 1.1445, 1.14688, 1.15003, 1.15107, $
1.15173, 1.15262, 1.15379, 1.15505, 1.15579, 1.1569, 1.4114, 1.41607, $
1.41953, 1.42404, 1.61382, 1.62244, 1.62786, 1.63487]

sz = size(n85,/DIM)

all2 = reform( all, [sz,19])
 print, max( all2[*,0]-n23) 
 
 ctrs = [90, 90, 225, 225, 90, 90,  226, 225, 90, 90, 225, $
 	223, 87, 87, 223, 227, 91, 91, 227]

index1 = where(ctrs lt 100)
index2 = where(ctrs gt 100)
 	
for i=0, 20 do begin plot, airmass[index2], all2[i+50,index2], psym=1, color=5 &$
oplot, airmass[index1], all2[i+50,index1], psym=2, color=4 & wait, 1.0

for i=0, 20 do begin plot, airmass[index2], all2[i+100,index2], psym=1, color=5 &$
oplot, airmass[index1], all2[i+100,index1], psym=2, color=4 & wait, 1.0

for i=0, 20 do begin plot, airmass[index2], all2[i+500,index2], psym=1, color=5 &$
oplot, airmass[index1], all2[i+500,index1], psym=2, color=4 & wait, 1.0

for i=0, 20 do begin plot, airmass[index2], all2[i+700,index2], psym=1, color=5 &$
oplot, airmass[index1], all2[i+700,index1], psym=2, color=4 & wait, 1.0

spec225 = (n25 + n26 + n29 + n30 + n33 + n66 + n69 + n82 + n85)/9.
spec90 = (n23 + n24 + n27 + n28 + n31 + n32 + n67 + n68 + n83 + n84)/10. 

plot, w, spec225, xrange=[3.3, 3.4], background=255, color=0
oplot, w, spec90, color=3

length = sz

offset = findgen(length)
slope = findgen(length)

for i=0, length-1 do begin $&
	res = GOODPOLY(airmass, all2[i,*], 1, 3.0)  &$
	offset[i] = res[0]  &$
	slope[i] = res[1]  &$
end

posS = where(slope gt 0)


plot, w, n32, background=255, color=0, xtitle='Wavelength, microns', ytitle='Normalized Flux', title='+ marks positive slope extinction, slit position 90'
oplot, w[poss], n32[poss], psym=1, color=0
plot, w, n33, background=255, color=0, xtitle='Wavelength, microns', ytitle='Normalized Flux', title='+ marks positive slope extinction, slit position 225
oplot, w[poss], n33[poss], psym=1, color=0


	
a = (findgen(8)+10.0)/10.

model = fltarr(length, 8) 
for i=0, length-1 do model[i,*] = offset[i] + slope[i]*a

;movie of extinction
for i=0, 20 do begin plot, a, model[i+100,*] & $
oplot, airmass, all2[i+100, *], psym=1 & wait, 1.0


dirD = '/Users/colkin/Work04/Pluto_3um/PlutoSpectra3/'
p34 = readfits(dirD+'12aus0034.fits')
p35 = readfits(dirD+'12aus0035.fits')
p36 = readfits(dirD+'12aus0036.fits')
p37 = readfits(dirD+'12aus0037.fits')
p38 = readfits(dirD+'12aus0038.fits')
p39 = readfits(dirD+'12aus0039.fits')
p40 = readfits(dirD+'12aus0040.fits')
p41 = readfits(dirD+'12aus0041.fits')
p42 = readfits(dirD+'12aus0042.fits')
p43 = readfits(dirD+'12aus0043.fits')
p44 = readfits(dirD+'12aus0044.fits')
p45 = readfits(dirD+'12aus0045.fits')
p46 = readfits(dirD+'12aus0046.fits')
p47 = readfits(dirD+'12aus0047.fits')
p48 = readfits(dirD+'12aus0048.fits')
p49 = readfits(dirD+'12aus0049.fits')
p50 = readfits(dirD+'12aus0050.fits')
p51 = readfits(dirD+'12aus0051.fits')
p52 = readfits(dirD+'12aus0052.fits')
p53 = readfits(dirD+'12aus0053.fits')
p54 = readfits(dirD+'12aus0054.fits')
p55 = readfits(dirD+'12aus0055.fits')
p56 = readfits(dirD+'12aus0056.fits')
p57 = readfits(dirD+'12aus0057.fits')
p58 = readfits(dirD+'12aus0058.fits')
p59 = readfits(dirD+'12aus0059.fits')
p60 = readfits(dirD+'12aus0060.fits')
p61 = readfits(dirD+'12aus0061.fits')
p62 = readfits(dirD+'12aus0062.fits')
p63 = readfits(dirD+'12aus0063.fits')
p64 = readfits(dirD+'12aus0064.fits')
p65 = readfits(dirD+'12aus0065.fits')
p70 = readfits(dirD+'12aus0070.fits')
p71 = readfits(dirD+'12aus0071.fits')
p72 = readfits(dirD+'12aus0072.fits')
p73 = readfits(dirD+'12aus0073.fits')
p74 = readfits(dirD+'12aus0074.fits')
p75 = readfits(dirD+'12aus0075.fits')



readcol, '/Users/colkin/Work04/Pluto_3um/Reduction/airmassAll', fn, obj, am, format="A, A, F"
print, obj

readcol, '/Users/colkin/Work04/Pluto_3um/Reduction/fileNameMaxRow2Pluto.txt', t1, t2, pctr, format="A, A, F"

pindex = where(obj eq "PLUTO")
help, pindex
pairmass = am[pindex]
pfn = fn[pindex]

pfn = fn[pindex[0:37]]
pairmass = am[pindex[0:37]]

dirExt = '/Users/colkin/Work04/Pluto_3um/Extinction3/'

ext35 = fltarr(length)
	for i=0, length-1 do ext35[i] = offset[i] + slope[i]*pairmass[0] 
	
	
for j=0, n_elements(pairmass)-1 do begin &$
	ext = fltarr(length) &$
	for i=0, length-1 do ext[i] = offset[i] + slope[i]*pairmass[j] &$
	writefits, dirExt+t1[j], ext &$
end

ext35 = readfits(dirExt+'12aus0035.fits')
ext36 = readfits(dirExt+'12aus0036.fits')
ext37 = readfits(dirExt+'12aus0037.fits')
ext38 = readfits(dirExt+'12aus0038.fits')
ext39 = readfits(dirExt+'12aus0039.fits')
ext40 = readfits(dirExt+'12aus0040.fits')
ext41 = readfits(dirExt+'12aus0041.fits')
ext42 = readfits(dirExt+'12aus0042.fits')
ext43 = readfits(dirExt+'12aus0043.fits')
ext44 = readfits(dirExt+'12aus0044.fits')
ext45 = readfits(dirExt+'12aus0045.fits')
ext46 = readfits(dirExt+'12aus0046.fits')
ext47 = readfits(dirExt+'12aus0047.fits')
ext48 = readfits(dirExt+'12aus0048.fits')
ext49 = readfits(dirExt+'12aus0049.fits')
ext50 = readfits(dirExt+'12aus0050.fits')
ext51 = readfits(dirExt+'12aus0051.fits')
ext52 = readfits(dirExt+'12aus0052.fits')
ext53 = readfits(dirExt+'12aus0053.fits')
ext54 = readfits(dirExt+'12aus0054.fits')
ext55 = readfits(dirExt+'12aus0055.fits')
ext56 = readfits(dirExt+'12aus0056.fits')
ext57 = readfits(dirExt+'12aus0057.fits')
ext58 = readfits(dirExt+'12aus0058.fits')
ext59 = readfits(dirExt+'12aus0059.fits')
ext60 = readfits(dirExt+'12aus0060.fits')
ext61 = readfits(dirExt+'12aus0061.fits')
ext62 = readfits(dirExt+'12aus0062.fits')
ext63 = readfits(dirExt+'12aus0063.fits')
ext64 = readfits(dirExt+'12aus0064.fits')
ext65 = readfits(dirExt+'12aus0065.fits')
ext70 = readfits(dirExt+'12aus0070.fits')
ext71 = readfits(dirExt+'12aus0071.fits')
ext72 = readfits(dirExt+'12aus0072.fits')
ext73 = readfits(dirExt+'12aus0073.fits')
ext74 = readfits(dirExt+'12aus0074.fits')
ext75 = readfits(dirExt+'12aus0075.fits')

px35 = p35/ext35
px36 = p36/ext36
px37 = p37/ext37
px38 = p38/ext38
px39 = p39/ext39
px40 = p40/ext40
px41 = p41/ext41
px42 = p42/ext42
px43 = p43/ext43
px44 = p44/ext44
px45 = p45/ext45
px46 = p46/ext46
px47 = p47/ext47
px48 = p48/ext48
px49 = p49/ext49
px50 = p50/ext50
px51 = p51/ext51
px52 = p52/ext52
px53 = p53/ext53
px54 = p54/ext54
px55 = p55/ext55
px56 = p56/ext56
px57 = p57/ext57
px58 = p58/ext58
px59 = p59/ext59
px60 = p60/ext60
px61 = p61/ext61
px62 = p62/ext62
px63 = p63/ext63
px64 = p64/ext64
px65 = p65/ext65
px70 = p70/ext70
px71 = p71/ext71
px72 = p72/ext72
px73 = p73/ext73
px74 = p74/ext74
px75 = p75/ext75

pxall = (px37 + px38+ px41+ px42+ px45+ px46+ px49+ px50+ px53+ px54+ px57+ px58+ px61+ px62+ px65+ px70+ px73+ px74 +$
px35+ px36+ px39+ px40+ px43+ px44+ px47+ px48+ px51+ px52+ px55+ px56+ px59+ px60+ px63+ px64+ px71+ px72+ px75)/(18+19)

plot, w, pxall
plot, w, smooth(pxall,20), yrange=[-1000,60000.]
oplot, w, smooth(pxall, 20), color=4

extArr1 =reform([ext37, ext38, ext41, ext42, ext45, ext46, ext49, ext50, ext53, ext54, ext57, ext58, $
ext61, ext62, ext65, ext70, ext73, ext74], [length, 18])
extArr2 =reform([ext35, ext36, ext39, ext40, ext43, ext44, ext47, ext48, ext51, ext52, ext55, ext56, $
ext59, ext60, ext63, ext64, ext71, ext72, ext75], [length, 19])

for i=0, 17 do begin plot, extArr1[*,i] & wait, 1.0
for i=0, 18 do begin plot, extArr2[*,i] & wait, 1.0

writefits, '/Users/colkin/Work04/Pluto_3um/PlutoSpectra3/pAve.fits', pxall