; ----------------------------------------------------------------------
;  Extract object and background spectrum.
;  Also apply wavelength scale.  Note use of EPERDN value.
;
; IMGSKY used only to find the wavelength scale
pro nsx_extract_spectrum, IMG, IMGSKY, AVP, SPX, nsxdir, nsxout
  common nsx

  in_nsx_extract_spectrum = 1

  
  print, "Extracting spectrum."

; Load sky spectrum(.spsky[]), load (shifted) wavelength scale(.spwav[]), and dispersion(.spdsp[]). 
  nsx_compute_skyline_shift, IMGSKY, SPX, nsxdir, nsxout 

; Allocate, create background subtracted image.
  image = PTR_NEW( fltarr(IMG.nc, IMG.nr) )
;> Only used if ATMABS_CORRECTION
;mktwave = (double *)calloc((NUMMKT),sizeof(double)) 
;MKT = (MKTtype *)calloc((NUMMKT),sizeof(MKTtype)) 

; Create background subtracted image. 
  for ii=0L,IMG.nc-1 do begin
     for jj=0L,IMG.nr-1 do begin
        pixno = ii + (jj * IMG.nc) 
        (*image)[pixno] = (*IMG.corimg)[pixno] - (*IMG.bckimg)[pixno] 
     endfor
  endfor


; Construction of atmospheric absorption data */
;if (CONSTRUCT_ATMABS) {
;  printf("construct atmospheric absorption data..\n") 
;  for (airmass=1.0  airmass<2.01  airmass=airmass+0.1) {
;    nsx_construct_atmabs_data( airmass, IMG, SPX ) 
;  }
;  printf("early exit\n") 
;  exit(0) 
;}

; Extract object spectrum from corrected, background subtracted object image for each order. */
; sets SPX.spobj and SPX.sperr
  for nso=3,7 do begin
     if (nso eq 7) then ecol=IMG.nc/2 else ecol=IMG.nc
     jjoff = ((7-nso) * 200) 
     rblo= double(jjoff) 
     rbhi= double(jjoff + SPX[nso].numpro) 
     rb1 = double(jjoff) + (SPX[0].asp1[0] / ARCSEC_PER_PIXEL) 
     if (rb1 lt rblo) then rb1 = rblo 
     if (rb1 gt rbhi) then rb1 = rbhi 
     rb2 = double(jjoff) + (SPX[0].asp2[0] / ARCSEC_PER_PIXEL) 
     if (rb2 lt rblo) then rb2 = rblo 
     if (rb2 gt rbhi) then rb2 = rbhi 
     as  = (SPX[0].asp1[0] + SPX[0].asp2[0]) / 2. 
     for icol=0L,ecol-1 do begin
        SPX[nso].spobj[icol] = nsx_fractional_pixel_rb( IMG.nc, image, rb1, rb2, icol ) 
        sum_corimg           = nsx_fractional_pixel_rb( IMG.nc, IMG.corimg, rb1, rb2, icol ) 
        sum_varimg           = nsx_fractional_pixel_rb( IMG.nc, IMG.varimg, rb1, rb2, icol ) 
        if (sum_varimg gt 0.) then SPX[nso].sperr[icol] = sqrt(sum_varimg/EPERDN) else SPX[nso].sperr[icol] = 0.
        edge = nsx_find_real_image_row( 1, icol, nso ) 
        SPX[nso].sprow[icol] = edge + nsx_AVPinv(AVP,nso,icol,as,IMG) ; Row position in original image. 
; Scale to per second. */
        SPX[nso].spobj[icol] = SPX[nso].spobj[icol] / IMG.exptime 
        SPX[nso].sperr[icol] = SPX[nso].sperr[icol] / IMG.exptime 
     endfor
     SPX[nso].numsp = ecol 
  endfor
;/* -------------------------  -tab 09feb2018
; I noticed that the AVP polynomials have slight 'notches' in them when you plot offset
; vs. column.. but the notches are small (~<0.1pixels) so probably not significant..
; you can see the affect in the following outputs: 
;  sprintf(wrd,"jaf%d.dat",nso) 
;  outfu = fopen_write(wrd) 
;  as = (SPX[0].asp1[0] + SPX[0].asp2[0]) / 2. 
;  rowoff = as / ARCSEC_PER_PIXEL   /x this is approx. since AVP refers to the uncorrected image x/
;  for (icol=0  icol<ecol  ++icol) {
;    rj1 = nsx_AVPinv(AVP,nso,icol,as,IMG) 
;    rj2 = nsx_AVP(AVP,nso,icol,rowoff,IMG) 
;    fprintf(outfu,"%d %f %f\n",icol,rj1,rj2) 
;  }
;  fclose(outfu) 
;  -------------------------- */
;
;
;
;/* Decide on best atmospheric transmission data set. */
;if (ATMABS_CORRECTION) {
;  nsx_find_best_mktfile( IMG, SPX, MKT, mktwave, nsxdir ) 
;}
;
;/* Vega comparison */
;if (VEGA_FLUX_CALIBRATION) {
;  nsx_vega_flux_calibration( IMG, SPX, nsxdir ) 
;  if (nso != 43234) exit(0) 
;}


; Extract background from corrected object image for each order (in background window(s)). */
  for nso=3,7 do begin
     if (nso eq 7) then ecol=IMG.nc/2 else ecol=IMG.nc 
     jjoff = ((7-nso) * 200) 
     rblo = double(jjoff )
     rbhi = double(jjoff + SPX[nso].numpro) 
     for icol=0L,ecol-1 do begin
        sum=0.
        npx=0. 
        for kk=0L,SPX[0].nbk-1 do begin
           rb1 = double(jjoff) + (SPX[0].abk1[kk] / ARCSEC_PER_PIXEL) 
           if (rb1 lt rblo) then rb1 = rblo 
           if (rb1 gt rbhi) then rb1 = rbhi 
           rb2 = double(jjoff) + (SPX[0].abk2[kk] / ARCSEC_PER_PIXEL) 
           if (rb2 lt rblo) then rb2 = rblo 
           if (rb2 gt rbhi) then rb2 = rbhi 
           sum = sum + nsx_fractional_pixel_rb( IMG.nc, IMG.corimg, rb1, rb2, icol ) 
           npx = npx + ABS((rb2 - rb1)) 
        endfor
        if (npx gt 0.) then SPX[nso].spbck[icol] = sum/npx else SPX[nso].spbck[icol] = 0.
        SPX[nso].spbck[icol] = SPX[nso].spbck[icol] / IMG.exptime 
     endfor 
  endfor

;  stop

; Free. 
  free,image
;free(MKT) 
;free(mktwave) 
;

  print, "returning from nsx_extract_spectrum to nsx_main"  ;> DEBUG
  return

end