datestr = datestr6_pwd()
echo=1
NCP = nsxng_load_ncp()
redo = 0
c = nsxng_constants() ; load NSX contants (were globals)
ncp = nsxng_load_ncp()          ; load NSX edge of slit polynomials
nsxng_load_WSC, c.nsxdir, WSC

; read the shortlog file...
nsets = rd_2_shortlog(datestr, sets, /nocal)

;for iset = 4, nsets-1 do begin
for iset = 2,2 do begin
   set = sets[iset]
   count = rd_2_position(datestr, set.setno, seqno, object, dpatname, beam)
   if strtrim(dpatname[0],2) eq 'ABBA' then begin
      beamarr = ['ABEAM','BBEAM']
   endif else begin
      beamarr = ['CENTER','SKY']
   endelse
   for ibeam = 0, 1 do begin
      beammatch = beamarr[ibeam]
      indx = where(beam eq beammatch, nim)
      if nim gt 0 then begin
         SPXIM = replicate(nsxng_SPXtype(),9, nim)
         ; The main loop over images

         for iim = 0, nim-1 do begin
            seq = seqno[indx[iim]]
            fn = fn1s( datestr, seqno[indx[iim]], ok)
            print, fn
            if not ok then stop
            nsxng_read_image, IA, echo, file=fn[0], logfu=-1
            IA.slitoffset = slitoffset(IA.root)
            
            ; THIS IS WHERE THE MAGIC HAPPENS
            nsxng_arcsec_profile_mash, IA, SPX, AVP, NCP, WSC
            SPXIM[*,iim] = SPX

         endfor                 ; loop over valid images in the beam
         ; After loop, save and plot the profiles
         ; are all the x's the same?
         for nso = 3, 7 do begin
            spx0 = spxim[nso, 0]
            for iim = 1, nim-1 do begin
               spxn = spxim[nso, iim]
               diff = max(abs((spxn.pro_apx)[0:spxn.pro_apn-1] - (spx0.pro_apx)[0:spxn.pro_apn-1]))
;               print, nso, iim, diff
               if diff gt 1e-6 then stop
            endfor
         endfor
         
         fnout = 'sav/'+'3_profile_'+string(set.setno,for='(I03)')+'_'+beammatch+'.sav'
         
         colors = [x11rgb('black'), x11rgb('DarkGray'), x11rgb('sienna'), x11rgb('red'), $
                   x11rgb('orange'), x11rgb('DarkGreen'), x11rgb('DarkBlue'), x11rgb('purple'), $
                   x11rgb('black'), x11rgb('DarkGray'), x11rgb('sienna'), x11rgb('red') ]
         line = [replicate(0,8), replicate(2,4) ]
         window, xs=840, ys=1000
         !p.multi=[0,1,5]
         for nso = 3, 7 do begin
            nmax = max( (spxim[nso,*].pro_apn) )
            xmax = max ( (spxim[nso,*].pro_apx)[nmax-1,*,*] )
            ymin = min(spxim[nso,*].pro_apymed)
            ymax = max(spxim[nso,*].pro_apymed)
            for iim = 0, nim-1 do begin
               spxn = spxim[nso, iim]
               if iim eq 0 then begin
                  plot, [0, xmax], [ymin, ymax], /xs, /ys, /nodata
                  xyouts, 0, (ymin+ymax)/2, string(nso)
               endif
               oplot, (spxn.pro_apx)[0:spxn.pro_apn-1], (spxn.pro_apymed)[0:spxn.pro_apn-1], $
                      color=colors[iim], line=line[iim]
               xyouts, xmax-5*ch_size(/x), ymax-(iim+1)*ch_size(/y), seqno[indx[iim]], color=colors[iim], /align
            endfor
         endfor
         fnout = 'fig/'+'3_profile_'+string(set.setno,for='(I03)')+'_'+beammatch+'.png'
         tv2png, fnout, /co, /tr
      endif                     ; any valid images
   endfor                       ; loop over beams
   print, 'end of set ', set.setno
endfor                  ; loop over the sets

end