;
; Modifications:
; 13_Dec-2007, Kim.  Added misc_struct to return values function return
; 1-jun-2017, ras,  fixes half pixel offset for even number of pixels,
; 14-Mar-2018, Kim. Only set plot params if show is set
;
;---


pro mem_njit_prep, vis, u, v, flux = flux, $
  ferr = ferr_, imsiz = imsiz_, misc = misc, $
  svis = svis_, tol = tol_, $
  xyint = xyint, show=show, caller=caller, $
  nf = nf, ptrs=ptrs, otol=tol, oferr=ferr, osvis=svis, $
  chi2_str=chi2_str, $ ;mmap=chi2mmap, chi2=chi2, cherr=cherr,$
  info=info, visobsk=visobsk, $
  entparm=entparm, m_=m, tflux=tflux, tferr=tferr, tb2flux=tb2flux,$
  dispset=dispset

  ;set parameters
  if keyword_set(tol_) then tol = tol_ > 1.e-6 < 1. else tol = 3.e-2
  if keyword_set(svis_) then svis = svis_ else svis = tol * abs(vis)
  if keyword_set(ferr_) then ferr = ferr_ else ferr = tol * flux
  if keyword_set(imsiz_) then m = long(imsiz_) else m = 128l
  if keyword_set(show) then show = 1 else show = 0

  ;backup current display setting
  if show then begin
    device, get_decomposed = dc0
    device, decomposed = 0
    tvlct, r0, g0, b0, /get
    multi0 = !p.multi
    ch0 = !p.charsize
    f0 = !p.font
    xm0 = !x.margin
    ym0 = !y.margin
    win0 = !d.window
    dispset = {dc0: dc0, r0: r0, g0:g0, b0:b0, multi:multi0, charsize:ch0, $
      font:f0, xmargin:xm0, ymargin:ym0, window:win0}
  endif else dispset={dummy:0}

  good = where(finite(vis) and svis gt 0.)
  vis = vis[good]
  svis = svis[good]
  u = u[good]
  v = v[good]

  if keyword_set(xyint) then begin ; *090105*
    tot_uv_siz = 206265. / xyint ; [radian]     ; *090105*
  endif else begin            ; *090105*
    tot_uv_siz = 4.0 * max([abs(u), abs(v)]) ; *090105*
    xyint = 206265. / tot_uv_siz ; [arcsec]     ; *090105*
  endelse                     ; *090105*
  uvint = tot_uv_siz / float(m) ; *090105*

  uv = complex(u, v)
  wgt = 1. / svis^2

  m2 = m / 2l
  f = 1.0   ;frequency array
  nf = n_elements(f)
  nf1 = nf - 1l
  tb2flux = 1.
  xyint = xyint * f[nf1] / f
  uvint = uvint * f / f[nf1]   ; *090105*
  spwgt = 0.
  spset = 0b

  pvisobs = ptrarr(nf)
  ptmap = ptrarr(nf)
  pmmap = ptrarr(nf)
  nul = fltarr(nf)
  chi2 = fltarr(nf)
  chi2mmap = fltarr(nf)
  info = {iter : 0, rchi2 : nul, rflux : nul, a : nul, b : nul, $
    perror : ptr_new([0l]), control : ptr_new(0), tol : tol, init : 1b}

  mmap = fltarr(m, m, nf)
  for k = 0l, nf1 do mmap[*, *, k] = flux[k] / tb2flux[k] / float(m)^2
  for k = 0l, nf1 do pmmap[k] = ptr_new(mmap[*, *, k] )
  tflux = flux / tb2flux
  tferr = tflux * tol
  tferr = ferr / tb2flux > tferr
  for k = 0l, nf1 do ptmap[k] = ptr_new(*pmmap[k])

  for k = 0l, nf1 do begin
    visobsk = ssmem_uvgrid(uv[*, *, k], vis[*, *, k] / tb2flux[k], $
      wgt[*, *, k] * tb2flux[k]^2, m, uvint[k])  ;*090105*
    chi2[k] = n_elements(visobsk.iu)
    pvisobs[k] = ptr_new(visobsk)
    dummy = ssmem_grd_chi2(*pmmap[k], visobsk, chi2mem = chi2memk)
    chi2mmap[k] = chi2memk
  endfor
  cherr = chi2 * tol > sqrt(2. * chi2)
  chi2_str = {chi2: chi2, chi2mmap:chi2mmap, cherr:cherr}

  ;check validity
  entparm = {pmmap : pmmap, spset : spset, spwgt : spwgt, $
    f : f, nf : nf, fqwgt : 0.}

  ; why block the caller = 'objref'?, this should now work with the
  ; change in the window call in ssmem_show_uv
  ;	if (show eq 1 and caller ne 'OBJREF') then begin
  if (show eq 1) then begin
    chk = ssmem_show_uv(m, xyint[nf1], uv[*, *, nf1], *pvisobs[nf1], caller = caller)
    if not chk then begin
      ptr_free, pmmap, ptmap, pvisobs

    endif
  endif
  ptrs = [pmmap, ptmap, pvisobs]
end



function mem_njit, vis, u, v, flux = flux, $
  ferr = ferr_, imsiz = imsiz_, misc = misc, $
  svis = svis_, tol = tol_, $
  xyint = xyint, show=show, caller=caller  ; *090105*



  mem_njit_prep, vis, u, v, flux = flux, $
    ferr = ferr_, imsiz = imsiz_, misc = misc, $
    svis = svis_, tol = tol_, $
    xyint = xyint, show=show, caller=caller,$
    nf = nf, ptrs=ptrs, otol=tol,oferr=ferr,osvis=svis,$
    chi2_str=chi2_str, $;chi2mmap=chi2mmap, chi2=chi2, cherr=cherr, $
    info=info, visobsk=visobsk,$
    entparm=entparm, m_=m, tflux=tflux, tferr=tferr, tb2flux=tb2flux,$
    dispset=dispset

  ;ptrs = [pmmap, ptmap, pvisobs]
  pmmap = ptrs[0]
  ptmap = ptrs[1]
  pvisobs = ptrs[2]

  ; chi2_str = {chi2: chi2, chi2mmap:chi2mmap, cherr=cherr}
  chi2 = chi2_str.chi2
  chi2mmap = chi2_str.chi2mmap
  cherr= chi2_str.cherr

  ;entparm = {pmmap : pmmap, spset : spset, spwgt : spwgt, $
  ;     f : f, nf : nf, fqwgt : 0.}
  f = entparm.f
  chi2chk = bytarr(nf)
  nf1 = nf - 1l
  for k = 0l, nf1 do chi2chk[k] = ssmem_chi2chk(*pvisobs[k], *ptmap[k], tol)

  bad = where(chi2mmap le chi2, nbad)
  if nbad gt 0l then begin
    error = [3l, nbad, bad, 0l] ;error code 3, invalid input
    *info.perror = error
    print, '%MEM_NJIT: Error. Too large data error.'
    goto, skip
  endif

  bad = where(chi2chk, nbad)
  if nbad gt 0l then begin
    error = [3l, nbad, bad, 0l] ;error code 3, invalid input
    *info.perror = error
    print, '%MEM_NJIT: Error. Too small map size.'
    goto, skip
  endif


  if show then begin
    control = ssmem_control(init = {f : f, m : m}, /single, $
      caller=caller )
    *info.control = control.tlb
  endif



  ;main algorithm
  ssmem_optimiz, ptmap, chi2, cherr, tflux, tferr, pvisobs, entparm, info, $
    show, caller = caller

  ;destroy control panel
  if (show eq 1) and not obj_valid( caller)  then status = ssmem_control(control.tlb, /destroy)

  skip:

  nmap = fltarr(m, m, nf)
  for k = 0l, nf1 do begin
    nmap[*, *, k] = *ptmap[k]
  endfor

  ptr_free, pmmap, ptmap, pvisobs

  ;organize outputs:

  error = *info.perror
  ;error information. [code, nbad, [bad....]].
  case error[0] of
    0l: print, '%MEM_NJIT: Successful end.'
    1l: print, '%MEM_NJIT: Error. Maximization failure.'
    2l: print, '%MEM_NJIT: Error. Aborted by user.'
    3l: print, '%MEM_NJIT: Error. Inconsistent constraint.'
    4l: print, '%MEM_NJIT: Error. Optimization failure.'
  endcase

  i0 = 1l
  nsgl = error[i0]
  if nsgl gt 0l then begin
    print, '%MEM_NJIT: Inconsistency detected.'
    i0 = i0 + 2l
  endif else i0 = i0 + 1l
  nbad = error[i0]
  if nbad gt 0l then begin
    print, '%MEM_NJIT: Constraint not satisfied.'
  endif

  misc = $
    ['%MEM_NJIT: iterat   = ' + string(info.iter), $
    '%MEM_NJIT: chi2_nrm = ' + string(info.rchi2), $
    '%MEM_NJIT: flux_nrm = ' + string(info.rflux), $
    '%MEM_NJIT: alpha    = ' + string(info.a), $
    '%MEM_NJIT: beta     = ' + string(info.b)]

  ptr_free, info.perror

  ;restore previous display setting
  ;    dispset = {dc0: dc0, r0: r0, g0:g0, b0:b0, multi:multi0, charsize:ch0, $
  ; 	font:f0, xmargin:xm0, ymargin:ym0, window:win0}

  if show then begin
    wset, dispset.window
    !y.margin = dispset.ymargin
    !x.margin = dispset.xmargin
    !p.font = dispset.font
    !p.charsize = dispset.charsize
    !p.multi = dispset.multi
    tvlct, dispset.r0, dispset.g0, dispset.b0
    device, decomposed = dispset.dc0
  endif

  ;x=(findgen(m)-m/2)*xyint[0] & y=x   ; *090105*
  ; 1-jun-2017, ras,  fixes half pixel offset for even number of pixels,
  x = reform(  ( pixel_coord( [ m, 1] ) )[ 0, *] ) * xyint[0]
  y = x
  xtit='[arcsec]' &ytit=xtit          ; *090105*

  entropy = total(nmap * alog(nmap/(flux[0]/tb2flux[0]/float(m)^2 * exp(1))))
  hsi_chi2 = info.rchi2[0] * chi2[0]
  lambda = double(entropy/hsi_chi2)

  misc_struct = str_subset(info, ['iter', 'rchi2', 'rflux', 'a', 'b'])

  return, {map:nmap, x:x, y:y, xtit:xtit, ytit:ytit, $
    flux:flux, ferr:ferr, imsiz:m, misc:misc, misc_struct:misc_struct, $
    svis:svis, tol:tol, xyint:xyint, $
    entropy:entropy, chi2:hsi_chi2, lambda:lambda,visobsk:visobsk} ; *090105*

  ;   return, nmap ; *090105*

end