pro surfer,qlin,qline,time,zen,line,timlin,numbin,atap,rev,day,year,mins,maxs,ig,ir,lclnoon,nsec,zbin

answer = ' '
ang = '!6 !SA!R!U!9 % !N!6'
a: read, 'Do you want the ' + line + ' intensities plotted ? ', answer
if (answer ne 'n') and (answer ne 'N') and $
   (answer ne 'y') and (answer ne 'Y') then begin
   print, ' Invalid response, try again bozo ! '
   goto, a
end
if (answer ne 'y') and (answer ne 'Y') then goto, b

print, ' '
numbns = numbin
print, 'NUMBNS = ', numbns
print, 'Please enter an integer(n) to multiply NUMBNS by.'
read, 'NUMBNS * n must be less than or equal to 512.  ', n
print, ' '

set_plot, 'ps'
device, /color, bits = 8, /landscape
loadct, 5

for inc=0,1 do begin
!p.position = [.025,.05,.675,.7]
 
 psub = '!6' + line + ang + ' Intensity (R)'
 q = qlin(*,*,inc)
 ;q = smooth(q,3)

 ;minq = min(q)
 ;maxq = max(q)

 ;if (fix(minq) eq 0) then minq = abs(minq)

 surface, q, time, zen, ax=45, az=15, charsize=2, subtitle=psub

overplot, rebin(q, numbns*n, 36*7), q, time, zen, [.675,.75,.925,.99], 1, $
   line + ang + ' Intensity (R)', 'Universal Time (hours)', $
   'Zenith Angle (!Uo!N)', 1
 
scl = max(q) 
xdm = 30
ydm = 255
bar = fltarr(ydm, xdm)
ymx = scl
ybar = (findgen(xdm) * ymx) / (xdm - 1)
for i = 0, 255-1 do begin
   bar(i,*) = ybar
endfor

!p.title = ' '
!y.title = ' '
!x.title = ' '
!y.range = 0
!x.style = 0
!y.style = 0
bar = (bar*255)/scl

overplotbar, bar, bar, ybar, [.94,.75,.975,.99], 1, '(R)', 1

!x.title = '!6Universal Time (hours)'
!y.title = '!6Zenith Angle (!Uo!N)'

magnet, atap, numbns, inc, lclnoon
weekap, day

;azimuth, atap, numbns, zbin, inc

xtmp = intarr(1)
ytmp = intarr(1)

ig = 3.0

!p.position = [0,0,1,1]
plot, xtmp, ytmp, xthick=0, ythick=0, xtitle=' ', ytitle=' ', $
   ticklen=0, xtickname=[' ', ' ', ' ', ' ', ' ', ' ', ' '], $
   ytickname=[' ', ' ', ' ', ' ', ' ', ' '], xstyle = 4, ystyle = 4, $
   /nodata, /noerase

timel = '!6'+strtrim(string(mins/3600.),2)
timeh = '!6'+strtrim(string(maxs/3600.),2)
timbin = '!6'+strtrim(string(fix(nsec)),2)
zenbin = '!6'+strtrim(string(fix(zbin)),2)

;xyouts, 0, .725, "!6 In Geomagnetic Coordinates'
xyouts, .725, .600, '!6Tape: '+ atap
xyouts, .725, .550, '!6Orbit: '+ rev
xyouts, .725, .500, '!6Day: '+ day
xyouts, .725, .450, '!6Year: ' +year
xyouts, .725, .400, '!6Grating Step: '+ strtrim(string(fix(ig)),2)
xyouts, .725, .350, '!6Resolution Mode: '+ strtrim(string(fix(ir)),2)
xyouts, .725, .300, '!6Start Time: '+ timel +' h'
xyouts, .725, .250, '!6End Time: '+ timeh + ' h'
xyouts, .725, .200, '!6Time Bin: '+ timbin + ' s'
xyouts, .725, .150, '!6Zenith Angle Bin: '+ zenbin +'!Uo!N'

endfor 
set_plot, 'x'

b: exit = 1
return
end


pro twozapf,dummy
print, 'WARNING - calculating counts not rayleighs!'
;set_plot,'ps'
;device,/landscape
atap = 'tapeidno'
scprsp = 0.192  ; seconds per spectra
read, ' What is the name of the tape to be plotted ? ', atap
openr, 1,   atap + '.two'
rv = 'rev#'
dy = 'day#'
yr = 'year'
tml = 0.0
tmh = 1.0
numbns = fix(0)
nzbns = fix(0)
numlns = fix(0)
ig = fix(0)
ir = fix(0)
izbnw = 0
nsec = 0

readf, 1, atap
readf, 1, rv
readf, 1, dy
readf, 1, yr
readf, 1, tml,tmh,numbns,nzbns,numlns,ig,ir,izbnw,nsec

print, ' atap = ',atap,', rv = ',rv,', dy = ',dy,', yr = ',yr
;if (dy lt 10) then info=string(format='(a8," Day:",i1," Year:",i4)', $
;atap,dy,yr)
;if (dy ge 10) then info=string(format='(a8," Day:",i2," Year:",i4)', $
;atap,dy,yr)
print, ' tml = ',tml/3600.,', tmh = ',tmh/3600.,', numbns = ',numbns,'
print, ' nzbns = ',nzbns, ' izbnw = ',izbnw,', nsec = ',nsec
lclnoon = fix((tml+tmh)/(3600*2))
  wline = 0.0
  wl1 = fix(0)
  wh1 = fix(0)
  wl = fix(0)
  wh = fix(0)
  wl2 = fix(0)
  wh2 = fix(0)
  sb1 = 0.0
  sb = 0.0
  sb2 = 0.0
  dn1 = 0.0
  dn2 = 0.0

   readf, 1, wline,wl1,wh1,wl,wh,wl2,wh2,sb1,sb,sb2,dn1,dn2

print,'sb = ',sb,' sb1 = ',sb1,' sb2 = ',sb2
pttime = fltarr(numbns)
ptglon = fltarr(numbns,2)
ptglat = fltarr(numbns,2)
ptblat = fltarr(numbns,2)
ptszen = fltarr(numbns,2)
sdglon = fltarr(numbns,2)
sdglat = fltarr(numbns,2)
sdblat = fltarr(numbns,2)
sdszen = fltarr(numbns,2)
timlin = fltarr(numbns,nzbns,2)
qlin = fltarr(numbns,nzbns,2)
qline = fltarr(numbns,nzbns,2)
ptzen = fltarr(numbns,nzbns,2)
sdzen = fltarr(numbns,nzbns,2)
ptazim = fltarr(numbns,nzbns,2)
sdazim = fltarr(numbns,nzbns,2)
ptsaz = fltarr(numbns,nzbns,2)
sdsaz = fltarr(numbns,nzbns,2)
for inc=0,1 do begin
  for itbn = 0, numbns - 1 do begin
    pttimex = 0.0
    ptglonx = 0.0
    ptglatx = 0.0
    ptblatx = 0.0
    ptszenx = 0.0
    sdglonx = 0.0
    sdglatx = 0.0
    sdblatx = 0.0
    sdszenx = 0.0

    readf, 1, pttimex,ptglonx,ptglatx,ptblatx,ptszenx
    readf, 1, sdglonx,sdglatx,sdblatx,sdszenx

    pttime(itbn) = pttimex
    ptglon(itbn,inc) = ptglonx
    ptglat(itbn,inc) = ptglatx
    ptblat(itbn,inc) = ptblatx
    ptszen(itbn,inc) = ptszenx
    sdglon(itbn,inc) = sdglonx
    sdglat(itbn,inc) = sdglatx
    sdblat(itbn,inc) = sdblatx
    sdszen(itbn,inc) = sdszenx
    for izbn = 0, nzbns-1 do begin
      pcntx = 0.0
      pbk1x = 0.0
      plinx = 0.0
      pbk2x = 0.0
      ptzenx = 0.0
      sdzenx = 0.0
      ptazimx = 0.0
      sdazimx = 0.0
      ptsazx = 0.0
      sdsazx = 0.0

      readf, 1, pcntx,pbk1x,plinx,pbk2x
      readf, 1, ptzenx,sdzenx,ptazimx,sdazimx,ptsazx,sdsazx

      qlinx = plinx / sb - pbk1x * dn1 - pbk2x * dn2
      ptzen(itbn,izbn,inc)=ptzenx
      sdzen(itbn,izbn,inc)=sdzenx
      ptazim(itbn,izbn,inc)=ptazimx
      sdazim(itbn,izbn,inc)=sdazimx
      ptsaz(itbn,izbn,inc)=ptsazx
      sdsaz(itbn,izbn,inc)=sdsazx
      qlinxe=plinx/sb^2+pbk1x*dn1^2+pbk2x*dn2^2
      if (qlinxe gt 0.0) then qlinxe=sqrt(qlinxe) else qlinxe=0.0
      timlin(itbn,izbn,inc) = pcntx
      if (pcntx gt 0) then begin
        qlin(itbn,izbn,inc) = plinx - (pbk1x*dn1 + pbk2x*dn2)*sb
        qline(itbn,izbn,inc) = sqrt(plinx + (pbk1x*dn1 + pbk2x*dn2)*sb)
        qlin(itbn,izbn,inc) = qlinx / pcntx / scprsp
        qline(itbn,izbn,inc) = qlinxe / pcntx / scprsp
      endif else begin
        qlin(itbn,izbn,inc) = 0.0
        qline(itbn,izbn,inc) = 0.0
      endelse
    end
  end
end
 
close,1
free_lun, 1

nsec = float(nsec)
zbin = float(izbnw)
time = findgen(numbns) * nsec + nsec / 2.0 + tml
time = time / 3600.0
zen = findgen(nzbns) * zbin + zbin / 2
!x.title = '!6Universal Time (hours)'
!y.title = '!6Zenith Angle (!Uo!N)'
!p.title = ' '
wq=where(qlin lt 0.0)
if (wq(0) ne -1) then qlin(wq) = 0.0
if (wline ge 1000.0) then chline=string(format='(i4)',wline)
if (wline lt 1000.0) then chline=string(format='(i3)',wline)

surfer, qlin, qline, time, zen, chline, timlin,numbns,atap,rv,dy,yr,tml,tmh,ig,ir,lclnoon,nsec,zbin

end