;code to read in the level 2 OMI data and plot the ozone for
;each scan
@/home/labow/nimbus/level2/make_coltabtv.pro

;**********put omi reader here*******

;restore,'oct05.sav'
restore,'oct24.sav'

;**************now enter in OMI scan sizes in Kilometers***********
lonstretch=fltarr(5,61)
latstretch=lonstretch
lonstretch(*,1)= [64,64,-64,-64,64]
lonstretch(*,2)= [62.21,62.21,-62.21,-62.21,62.21]
lonstretch(*,3)= [60.41,60.41,-60.41,-60.41,60.41]
lonstretch(*,4)= [58.62,58.62,-58.62,-58.62,58.62]
lonstretch(*,5)= [56.83,56.83,-56.83,-56.83,56.83]
lonstretch(*,6)= [55.03,55.03,-55.03,-55.03,55.03]
lonstretch(*,7)= [53.24,53.24,-53.24,-53.24,53.24]
lonstretch(*,8)= [51.45,51.45,-51.45,-51.45,51.45]
lonstretch(*,9)= [49.65,49.65,-49.65,-49.65,49.65]
lonstretch(*,10)= [47.86,47.86,-47.86,-47.86,47.86]
lonstretch(*,11)= [46.07,46.07,-46.07,-46.07,46.07]
lonstretch(*,12)= [44.27,44.27,-44.27,-44.27,44.27]
lonstretch(*,13)= [42.48,42.48,-42.48,-42.48,42.48]
lonstretch(*,14)= [40.69,40.69,-40.69,-40.69,40.69]
lonstretch(*,15)= [38.89,38.89,-38.89,-38.89,38.89]
lonstretch(*,16)= [37.1,37.1,-37.1,-37.1,37.1]
lonstretch(*,17)= [35.31,35.31,-35.31,-35.31,35.31]
lonstretch(*,18)= [33.52,33.52,-33.52,-33.52,33.52]
lonstretch(*,19)= [31.72,31.72,-31.72,-31.72,31.72]
lonstretch(*,20)= [29.93,29.93,-29.93,-29.93,29.93]
lonstretch(*,21)= [28.14,28.14,-28.14,-28.14,28.14]
lonstretch(*,22)= [26.34,26.34,-26.34,-26.34,26.34]
lonstretch(*,23)= [24.55,24.55,-24.55,-24.55,24.55]
lonstretch(*,24)= [22.76,22.76,-22.76,-22.76,22.76]
lonstretch(*,25)= [20.96,20.96,-20.96,-20.96,20.96]
lonstretch(*,26)= [19.17,19.17,-19.17,-19.17,19.17]
lonstretch(*,27)= [17.38,17.38,-17.38,-17.38,17.38]
lonstretch(*,28)= [15.59,15.59,-15.59,-15.59,15.59]
lonstretch(*,29)= [13.79,13.79,-13.79,-13.79,13.79]
lonstretch(*,30)= [12,12,-12,-12,12]	;;***NADIR***
for n=31,60 do begin
  lonstretch(*,n)=lonstretch(*,61-n)		;make sp 30 eq 31
endfor

;*********************size doesn't change in latitude??*****
x=[-6.5,6.5,6.5,-6.5,-6.5]
for n=0,60 do begin
  latstretch(*,n)=x
endfor


device,pseudo=8
window, 0, XSIZE=600, YSIZE=480,colors=248
make_coltabtv

read,'enter latitude min and max (-90,90) ',lat1,lat2
read,'enter longitude min and max (-180,180) ',lon1,lon2

w=where((lat_all gt lat1) and (lat_all lt lat2) and (lon_all gt lon1) and (lon_all lt lon2) and (so2_all gt -100))

   latit=lat_all(w)
   longit=lon_all(w)
   spp=is_all(w)
   soi=so2_all(w)
;   ozz=bytscl(soi,-4,4)
   ozz=bytscl(soi,0,32)
   orb=SECDAY_ALL(w)

  ;now define each orbit
    s=shift(orb,1)
    d=s-orb
    ww=where(abs(d) gt 3000)            ;find orbital breaks
    ww=[ww,n_elements(orb)-1]           ;add endpoint
    x=fltarr(600,480,n_elements(ww)-1)

  if (n_elements(ww) eq 2) then begin		;if only one orbit
     orblat=latit
     orblon=longit
     orbsp=spp
     pixcol=ozz
     map_set,limit=[lat1,lon1,lat2,lon2],/merc,/cont,/grid,color=244
    
   for c=0,n_elements(orblat)-1 do begin 
    distance=110.*cos(orblat(c)/180.*!pi)
    ptlon=orblon(c)+lonstretch(*,orbsp(c))/distance
    ptlat=orblat(c)+latstretch(*,orbsp(c))/110.
    w=where (ptlon gt 180)
     if (w(0) ge 0) then begin
       ptlon(w)=180.
     endif
    w=where (ptlon lt -180)
     if (w(0) ge 0) then begin
       ptlon(w)=-180.
     endif 
    polyfill,ptlon,ptlat,color=pixcol(c)
    plots,ptlon,ptlat,/noclip
   endfor

   ;plots,-122.2,46.2,psym=2,thick=2                ;mt st helens
   plots,15.00,37.75,psym=2,thick=2               ;mt etna

   map_continents,/coasts,/hires,mlinethick=2
   x(*,*,0)=tvrd()

  goto,skip
  endif				;end 1 orbit loop

  for cc=0,n_elements(ww)-2 do begin            ;for each orbit if more than one
    orblat=latit(ww(cc):ww(cc+1))
    orblon=longit(ww(cc):ww(cc+1))
    orbsp=spp(ww(cc):ww(cc+1))
    pixcol=ozz(ww(cc):ww(cc+1))
    map_set,limit=[lat1,lon1,lat2,lon2],/merc,/cont,/grid,color=244

   for c=0,n_elements(orblat)-1 do begin 
    distance=110.*cos(orblat(c)/180.*!pi)
    ptlon=orblon(c)+lonstretch(*,orbsp(c))/distance
    ptlat=orblat(c)+latstretch(*,orbsp(c))/110.
    w=where (ptlon gt 180)
     if (w(0) ge 0) then begin
       ptlon(w)=180.
     endif
    w=where (ptlon lt -180)
     if (w(0) ge 0) then begin
       ptlon(w)=-180.
     endif 
    polyfill,ptlon,ptlat,color=pixcol(c)
    plots,ptlon,ptlat,/noclip
   endfor


   plots,-122.2,46.2,psym=2,thick=2                ;mt st helens
   ;plots,15.00,37.75,psym=2,thick=2               ;mt etna

   map_continents,/coasts,/hires,mlinethick=2
   x(*,*,cc)=tvrd()

  endfor

skip:
erase

for n=0,n_elements(ww)-2 do begin
  tv,x(*,*,n)
  xyouts,.01,.175,'orbit',size=2.5,charthick=3,/normal
  xyouts,-.4,.05,n,size=4,charthick=4,/normal
endfor

num=0
;read,'Choose Orbit Number to Plot ',num

date='October 24, 2004'
window, 0 , XSIZE=750, YSIZE=670,colors=248,retain=2
tv,x,30,70

 XYOUTS, .145, 0.08, 'SO2 Index', /normal, SIZE=1.2,  color=255, alignment=0.5
 XYOUTS, .14, 0.04, 'NASA/GSFC', /normal, SIZE=1.2,   color=255, alignment=0.5
 XYOUTS, .42, 0.80, 'AURA-OMI', /normal, SIZE=3.0,  charthick=2, color=220, alignment=0.5
 XYOUTS, .42, 0.91, date, /normal, SIZE=2.3, charthick=2, $
  color=90, alignment=0.5

a=intarr(30,28,16)
for n=0,15 do a(*,*,n)=109+(25*n)

back=intarr(34,480)+245
tv,back,693,68
                                                      ;make colorbar
for n=0,15 do tv,bytscl(a(*,*,n),100,500),695,n*30+70


;bar=['-4','-3','-2','-1',' 0',' 1',' 2',' 3',' 4']
bar=[' 0',' 4',' 8','12','16','20','24','28','32']
for n=0,8 do xyouts,655,n*59+65,bar(n),color=244,/device            ;make colorbar scale
  XYOUTS, .84, 0.50, 'SO2 Index', /normal, SIZE=1.2, color=244, alignment=0.5, orientation=270


omi=tvrd()
write_gif,'omi_manam_soi.gif',omi,rr,gg,bb
; set_plot,'PS' & device,filename='omi_manam_ai.ps',/landscape

; device,/close


end