pro computeDxDy, im_final, im, old_xref, old_yref, xref, yref, old_xcen, old_ycen, xcen, ycen, dx, dy
COMPILE_OPT IDL2, HIDDEN

  cntrd, im_final, old_xref, old_yref, xref, yref, 5, /silent
  cntrd, im, old_xcen, old_ycen, xcen, ycen, 5, /silent
  if xcen eq -1 then begin
    print, 'Unable to determine centroid at (' + strtrim(xref, 2) + ', ' + strtrim(yref, 2) + ')'
    print, 'Using old values...'
    xcen = old_xcen
    ycen = old_ycen
  endif

  dx = xcen-xref
  dy = ycen-yref
end


pro findStars, im, xstar, ystar, xpos, ypos, plot
COMPILE_OPT IDL2, HIDDEN

  cntrd, im, xstar, ystar, xc, yc, 5, /silent

  if plot eq 1 then begin
    test = where(xc ne -1)
    if test[0] ne -1 then tvcircle, 10, xc, yc
  endif

  test = where(xc eq -1)
  if test[0] ne -1 then begin
    xc[test] = xstar[test]
    yc[test] = ystar[test]
  endif

  if n_elements(xpos) eq 0 then begin
    xpos = xc
    ypos = yc
  endif else begin
    xpos = [[xpos], [xc]]
    ypos = [[ypos], [yc]]
  endelse

  xstar = xc
  ystar = yc

end

function reorderAdaChronogically, ada
  cycles = strmid(ada, 4, 2)
  channels = strmid(ada, 6, 2)

  for i=1, max(cycles) do begin
    if (i mod 2) eq 1 then continue
    
    idx = where(cycles eq i)
    if idx[0] eq -1 then continue
    
    ada[idx] = ada[reverse(idx)]
  endfor

  return, ada
end

function integrateAdaAndCorrectGuiding, directory, prb, cal_cen, nobreak, xstars, ystars, flux, p, ignorecycle=ignorecycle, correct_continuum=correct_continuum, stars_file=stars_file
COMPILE_OPT IDL2

  if n_elements(p) eq 0 then p = 1
  if (p eq 0) and ((n_elements(xstars) + n_elements(ystars)) lt 2) then begin
    print, 'Error! Unable to correct guiding with plot disabled and no reference star given!'
    return, -1
  endif

  if n_elements(nobreak) eq 0 then nobreak=0

  spawn, 'pwd', current_wd

  pushd, directory

  sz = readadt(adtch, adtcy, adteff)

  xsize = sz[1]
  ysize = sz[2]
  zsize = sz[0]

;  if n_elements(correctNum) eq 0 then correctNum = 5  ; Correct position every 5 ada

  ada = file_search('*.[aA][dD][aA]')

  ada = reorderAdaChronogically(ada)

  cycles = strmid(ada, 4, 2)
  channels = strmid(ada, 6, 2)

  if channels[0] eq 0 then begin
    ada = ada[1:n_elements(ada)-1]
    channels = channels[1:n_elements(channels)-1]
    cycles = cycles[1:n_elements(cycles)-1]
  endif

  ; On retire les cycles incomplets...
  nchannels = max(channels)
  ncycles = max(cycles)
  bad = 0
  for i=1, ncycles do begin
    elem = where(cycles eq i)
    if n_elements(elem) lt nchannels then begin
      bad = [bad, elem]
      print, 'Removing cycle ' + strtrim(i,2)
    endif
  endfor

  if n_elements(bad) ne 1 then begin
    bad = bad[1:n_elements(bad)-1]
    remove, bad, ada
    remove, bad, channels
    remove, bad, cycles
  endif

  if n_elements(ignorecycle) ne 0 then begin
    for i=0, n_elements(ignorecycle)-1 do begin
      idx2 = where(cycles eq ignorecycle[i])
      idx = where(cycles ne ignorecycle[i])
      if (idx[0] ne -1) and (idx2[0] ne -1) then begin
        print, 'Ignoring cycle ' + strtrim(ignorecycle[i],2)
        ada = ada[idx]
        cycles = cycles[idx]
        channels = channels[idx]
      endif
    endfor
  endif

  nphotons = fltarr(n_elements(ada))

  if p eq 1 then begin
    window, 0, xsize=xsize, ysize=ysize
    loadct, 0, /silent
  endif

;  goto, do_correct_continuum

  if n_elements(xstars) eq 0 then begin
    ; On commence par trouver les etoiles de reference
    im = readada(ada[where(cycles eq min(cycles))], xsize, ysize, silent=1)

    tvset, 0
    im[*, ysize-1] = 0 ; Pour eviter la ligne de photons dans le bas...
    im[*, 0] = 0       ; probleme avec FaNTOmM
    tvscl, rescale(im, 0.001, 0.999)
    print, 'Click on reference regions...'
    print, 'Double click to add no more!'

    while 1 do begin
      cursor, x, y, /device, /up
      cntrd, im, x, y, xc, yc, 5
    
      if xc eq -1 then begin
        print, 'Unable to find a reference star... Try again...'
        continue
      endif

      print, 'Star is at ', xc, yc

      tvcircle, 10, xc, yc
      if n_elements(xstars) eq 0 then begin
        xstars = xc
        ystars = yc
      endif else begin
        if (floor(xc) eq floor(xstars[n_elements(xstars)-1])) and (floor(yc) eq floor(ystars[n_elements(ystars)-1])) then break

        xstars = [xstars, xc]
        ystars = [ystars, yc]
      endelse
    endwhile
    print, 'Writing ' + stars_file
    toto = writead1(stars_file, [xstars, ystars])
  endif
 
  print, 'Ok, got ' + strtrim(n_elements(xstars),2) + ' stars!'

  im_total = fltarr(xsize, ysize)

  for i=1, max(cycles) do begin ;ceil(float(n_elements(ada))/factor) do begin
;    im = readada(ada[i*factor:((i+1)*factor-1)<(n_elements(ada)-1)], xsize, ysize)

    test = where(cycles eq i)
    if test[0] ne -1 then begin
      im = readada(ada[test], xsize, ysize, silent=1)
    endif else continue

    if p eq 1 then begin
      tvscl, rescale(im, 0.0001, 0.9999)
      xyouts, 10, 10, 'Cycle ' + strtrim(i,2), /device
    endif

    findStars, im, xstars, ystars, xpos, ypos, p
  endfor

  ; On trouve le barycentre de chaque cycle...
  baryx = total(xpos, 1)/n_elements(xpos[*,0])
  baryy = total(ypos, 1)/n_elements(ypos[*,0])

  if nobreak eq 0 then begin
    ; On trouve le moment de la cassure principale...
    maxx = max(deriv(baryx), idx_maxx)
    minx = min(deriv(baryx), idx_minx)
    maxy = max(deriv(baryy), idx_maxy)
    miny = min(deriv(baryy), idx_miny)

    max = max(abs([maxx, minx, maxy, miny]), idx)
    idx = ([idx_maxx, idx_minx, idx_maxy, idx_miny])[idx]

    ; Maintenant, on sait ou est la cassure principale...
    ; Il faut retourner dans ce cycle et integrer a plus
    ; petits pas pour trouver le moment exact de la cassure...
    test = where((cycles eq idx) or (cycles eq (idx-1)) or (cycles eq (idx+1)))
    if test[0] eq -1 then begin
      print, 'Big fuckin problem man!'
      return, -1
    endif

    xstar = xpos[*,idx-1]
    ystar = ypos[*,idx-1]

    ada_fine = ada[test]
    factor = 8
    for i=0, n_elements(ada_fine)/factor-1 do begin
      im = readada(ada_fine[i*factor:((i+1)*factor-1)<(n_elements(ada_fine)-1)], xsize, ysize, silent=1)

      if p eq 1 then begin
        tvscl, rescale(im, 0.0001, 0.9999)
        xyouts, 10, 10, strtrim(i,2), /device
      endif

      findStars, im, xstar, ystar, xposfine, yposfine, p
      print, ada_fine[i*factor:((i+1)*factor-1)<(n_elements(ada_fine)-1)]
      print, '  ', xposfine[n_elements(xposfine)-1], yposfine[n_elements(yposfine)-1]
    endfor

    baryfinex = total(xposfine, 1)/n_elements(xposfine[*,0])
    baryfiney = total(yposfine, 1)/n_elements(yposfine[*,0])

    maxx = max(deriv(baryfinex), idx_maxx)
    minx = min(deriv(baryfinex), idx_minx)
    maxy = max(deriv(baryfiney), idx_maxy)
    miny = min(deriv(baryfiney), idx_miny)

    max = max(abs([maxx, minx, maxy, miny]), idxfine)
    idxfine = ([idx_maxx, idx_minx, idx_maxy, idx_miny])[idxfine]

    print, 'Idxfine: ', idxfine

    fitxBegin = ladfit(findgen(2)+(idx-2), baryx[idx-2:idx-1])
    fitxEnd = ladfit(findgen(2)+idx, baryx[idx:idx+1])

    fityBegin = ladfit(findgen(2)+(idx-2), baryy[idx-2:idx-1])
    fityEnd = ladfit(findgen(2)+idx, baryy[idx:idx+1])

    xcenBegin = findgen(n_elements(baryx)*max(channels))/max(channels)
    xcenEnd = findgen(n_elements(baryx)*max(channels))/max(channels)
    ycenBegin = findgen(n_elements(baryy)*max(channels))/max(channels)
    ycenEnd = findgen(n_elements(baryy)*max(channels))/max(channels)

    xcenBegin = xcenBegin*fitxBegin[1]+fitxBegin[0]
    xcenEnd = xcenEnd*fitxEnd[1]+fitxEnd[0]

    ycenBegin = ycenBegin*fityBegin[1]+fityBegin[0] 
    ycenEnd = ycenEnd*fityEnd[1]+fityEnd[0]

    cassure = (idx-2)*max(channels)+idxfine*factor
    print, ada[cassure]

    baryx = convol(baryx, [0.25, 0.5, 0.25], /edge_truncate)
    baryy = convol(baryy, [0.25, 0.5, 0.25], /edge_truncate)

    xcenReal = interpolate(baryx, findgen(n_elements(ada))/(max(channels)))
    ycenReal = interpolate(baryy, findgen(n_elements(ada))/(max(channels)))

    xcen = [xcenReal[0:(idx-2)*max(channels)-1], xcenBegin[(idx-2)*max(channels):cassure-1], xcenEnd[cassure:(idx+1)*max(channels)-1], xcenReal[(idx+1)*max(channels):n_elements(xcenReal)-1]]
    ycen = [ycenReal[0:(idx-2)*max(channels)-1], ycenBegin[(idx-2)*max(channels):cassure-1], ycenEnd[cassure:(idx+1)*max(channels)-1], ycenReal[(idx+1)*max(channels):n_elements(ycenReal)-1]]

  endif else begin

     baryx = convol(baryx, [0.25, 0.5, 0.25], /edge_truncate)
     baryy = convol(baryy, [0.25, 0.5, 0.25], /edge_truncate)

     xcen = interpolate(baryx, findgen(n_elements(ada))/(max(channels)))
     ycen = interpolate(baryy, findgen(n_elements(ada))/(max(channels)))
     
  endelse

  dx = xcen - median(xcen)
  dy = ycen - median(ycen)
  k=0L

  minx = min([min(xcen), min(baryx)])
  miny = min([min(ycen), min(baryy)])

  if p eq 1 then begin
    window, 2
    !P.multi = [0,1,2]
    plot, findgen(n_elements(baryx))+1, baryx-minx, xtickinterval=1, xminor=1, psym=1, title='X error' ;all_xcen-min(all_xcen)
    oplot, findgen(n_elements(baryx)*max(channels))/max(channels)+1, xcen-minx ;new_xcen-min(all_xcen)
    plot, findgen(n_elements(baryy))+1, baryy-miny, xtickinterval=1, xminor=1, psym=1, title='Yerror'  ;all_ycen-min(all_ycen)
    oplot, findgen(n_elements(baryx)*max(channels))/max(channels)+1, ycen-miny ;new_ycen-min(all_ycen)

    !P.multi = 0

    window, 1, xsize=xsize, ysize=ysize
  endif

do_correct_continuum:
;  window, 0, xsize=1400, ysize=800
;  window, 1, xsize=512, ysize=512
;  tvset, 0

  if keyword_set(correct_continuum) then begin

    !P.multi = [0,2,1]
    loadct, 39, /silent

    continuum_binsize=4

    dist_circle, cir, sz[2], cal_cen[0], cal_cen[1]
    hist = histogram(cir, binsize=continuum_binsize, reverse_indices=r)

    max = sqrt(cal_cen[0]^2+cal_cen[1]^2)
    max = max > sqrt((cal_cen[0]-xsize)^2+cal_cen[1]^2)
    max = max > sqrt((cal_cen[0]-xsize)^2+(cal_cen[1]-ysize)^2)
    max = max > sqrt(cal_cen[0]^2+(cal_cen[1]-ysize)^2)
    max -= (max mod continuum_binsize)
    max /= continuum_binsize

    arrayx = (abs(findgen(xsize)-cal_cen[0]) # replicate(1,ysize))/continuum_binsize
    arrayy = (replicate(1,xsize) # abs(findgen(ysize)-cal_cen[1]))/continuum_binsize
    array = sqrt(arrayx^2+arrayy^2)
    array2 = array^2
    array3 = array^3
    array4 = array^4

    all_fits = fltarr(2, zsize, max(cycles))

    for i=0, zsize-1 do begin

      print, 'Processing channel ' + strtrim(i,2)

      test = where(channels eq i+1)
      if test[0] ne -1 then begin
        temp = fltarr(n_elements(test), xsize, ysize)

        for j=0, n_elements(test)-1 do begin
          temp[j,*,*] = readada(ada[test[j]], xsize, ysize, silent=1)
        endfor

        radial_flux = fltarr(n_elements(test), max)

;        filename = current_wd[0] + '/channel-'+strtrim(i,2)+'-radialflux.fits'

;        if file_test(filename, /read) then begin
;          radial_flux = readfits(filename, /silent)
;        endif else begin

        temp=reform(temp, n_elements(test), xsize*ysize)

        for j=0, max-1 do begin
          idx = r[r[j]:r[j+1]-1]
          for k=0, n_elements(test)-1 do begin
            radial_flux[k,j] = mean(temp[k,idx])
;              radial_flux[k,j] = mean((reform(temp[k,*,*]))[idx])
          endfor
        endfor
 
;          writefits, filename, radial_flux

;        endelse

        if keyword_set(p) then begin
           plot, radial_flux[0,*], yrange=[0, max(radial_flux)], title='Before correction', /nodata
           for k=0, n_elements(test)-1 do oplot, radial_flux[k,*], color=40+(250-40)*(float(k)/n_elements(test))

           plot, radial_flux[0,*], /nodata, title='After correction'
        endif

        for k=1, n_elements(test)-1 do begin
          error = reform(radial_flux[k,*] - radial_flux[0,*])
          fit = ladfit(findgen(max*0.8), error[0:max*0.8-1])
          yfit = fit[0]+fit[1]*findgen(max)
;          fit = poly_fit(findgen(max), error, 4, yfit=yfit)

;          tvset, 1
;          if k eq 1 then begin
;            plot, yfit, color=40+(250-40)*(float(k)/n_elements(test))
;          endif else begin
;            oplot, yfit, color=40+(250-40)*(float(k)/n_elements(test))
;          endelse
;          tvset, 0

          if keyword_set(p) then oplot, radial_flux[k,*] - yfit, color=40+(250-40)*(float(k)/n_elements(test))

          all_fits[*,i,k] = fit
;          print, 'Fit for ' + ada[test[k]], fit
        endfor

      endif else continue

;    neb[i, *,*] = temp
    endfor
  wait, 2
  !P.multi=0

  loadct, 0, /silent

  endif
  
  lagx = findgen(xsize)-(xsize/2)
  lagy = findgen(ysize)-(ysize/2)
  lambda = fltarr(zsize, xsize, ysize)
  flux = 0

  for i=0, n_elements(ada)-1 do begin
    eff = where((adtch eq channels[i]) and (adtcy eq cycles[i]))

    if keyword_set(correct_continuum) then begin
      yfit = all_fits[*,channels[i]-1,cycles[i]-1]
      map = fit[0] + fit[1]*array ; + fit[2]*array2 + fit[3]*array3 + fit[4]*array4
    endif

;    mapped_cube = fltarr(zsize, xsize, ysize)

;    print, adteff[eff[0]]

    lambda_temp = readAdaToLambda(ada[i], channels[i]-1, zsize, prb, adteff[eff[0]], np, -dx[i], -dy[i], continuum_map=map)
    flux = [flux, np]
;    lambda_temp = readAdaToLambda(ada[i], channels[i]-1, zsize, prb, np, 0, 0)
    nphotons[i] = np

    print, '  ' + ada[i] + ': ' + strtrim(np,2) + ' photons, correction: ' + strtrim(dx[i],2) + ', ' + strtrim(dy[i],2)
    if keyword_set(correct_continuum) then print, '    yfit: ', yfit

    lambda = lambda + lambda_temp

 ;   for j=0, zsize-1 do begin
 ;     lambda[j,*,*] = lambda[j,*,*] + interpolate(lambda_temp[j,*,*], findgen(xsize)+dx[i], findgen(ysize)+dy[i], /grid, cubic=-0.6)
 ;   endfor
    if ((i mod 10) eq 0) and (p eq 1) then begin
      tvset, 0
      tvscl, rescale(total(lambda_temp,1), 0.0001, 0.9999)
      tvset, 1
      tvscl, rescale(total(lambda,1), 0.001, 0.999)
    endif

;    print, dx[i], dy[i]

  endfor

  flux = flux[1:n_elements(flux)-1]

  popd

  return, lambda

end