source: trunk/SRC/ToBeReviewed/TRIANGULATION/section.pro @ 495

;+
; @file_comments
;
; @categories
;
; @param FIELD
;
; @param RES
;
; @param GLAMAXE
;
; @param GPHIAXE
;
; @keyword ENDPOINTS
;
; @keyword BOXZOOM
;
; @keyword TYPE
;
; @keyword WDEPTH
;
; @keyword DIREC
;
; @keyword SHOWBUILD
;
; @keyword ONLYBOX
;
; @keyword _EXTRA
; Used to pass keywords
;
; @returns
;
; @uses
; <pro>common</pro>
;
; @restrictions
;
; @examples
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS=endpoints $
             , BOXZOOM=boxzoom, TYPE=type, WDEPTH=wdepth $
             , DIREC=direc, SHOWBUILD=showbuild, ONLYBOX=onlybox $
             , _EXTRA=ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
  ENDIF
;--------------------------------------------------------------
;---------------------
;---------------------
; definition of boxzoom in function of endpoints, then redefinition of the domain
  boxzoom2d = [min([endpoints[0], endpoints[2]], max = ma02), ma02 $
               , min([endpoints[1], endpoints[3]], max = ma13), ma13]
;
  minprof = 0.
  profdefault = 200.
  if n_elements(type) EQ 0 then type = 'nothing'
  Case N_Elements(Boxzoom) OF
    0:localbox = [boxzoom2d, minprof, profdefault]
    1:localbox = [boxzoom2d, minprof, boxzoom[0]]
    2:localbox = [boxzoom2d, boxzoom]
    4:if strpos(type, 'z') NE -1 THEN $
      localbox = [boxzoom2d, minprof, profdefault] ELSE localbox = boxzoom2d
    5:localbox = [boxzoom2d, minprof, boxzoom[4]]
    6:localbox = [boxzoom2d, boxzoom[4:5]]
    Else:BEGIN
      ras = report('Bad definition of the box')
      stop
    END
  ENDCASE
  nelbox = n_elements(localbox)
;
  if keyword_set(wdepth) then grillechoice = [vargrid, 'W'] $
  ELSE grillechoice = vargrid
  domdef, localbox, GRIDTYPE = grillechoice, /findalways, _extra = ex
  grille, -1, -1, -1, -1, nx, ny
; if less than 10 points where found, we apply domdef over the whole domain
; -> problem... why 10 points as a test value???
; how can we find a good test value?
  IF nx * ny LE 10 THEN domdef, [0, jpi-1, 0, jpj-1, localbox[4:5]], /xindex, /yindex, GRIDTYPE = grillechoice, _extra = ex
; We redefine lon1, ... in case findalways has been used in domdef
  lon1 = min([endpoints[0], endpoints[2]], max = lon2)
  lat1 = min([endpoints[1], endpoints[3]], max = lat2)
; we extend the box along the z axis -> i that way the plot will be drawn
; until its bottom part.
  if strpos(type, 'z') NE -1 THEN BEGIN
; We keep yranges (axis z) before changing the boxzoom.
    !y.range = [localbox[nelbox-1], localbox[nelbox-2]]
    if vargrid EQ 'W' OR keyword_set(wdepth) then BEGIN
      firstzw = 0 > (firstzw-1)
      lastzw = (lastzw+1) < (jpk-1)
      nzw = lastzw - firstzw + 1
    ENDIF ELSE BEGIN
      firstzt = 0 > (firstzt-1)
      lastzt = (lastzt+1) < (jpk-1)
      nzt = lastzt - firstzt + 1
    ENDELSE
    IF NOT keyword_set(key_forgetold) THEN BEGIN
     @updateold
   ENDIF
  ENDIF
  !y.range = [localbox[nelbox-1], localbox[nelbox-2]]
; We recuperate the grid on the boxzoom.
;
  IF NOT keyword_set(onlybox) THEN saveboxparam, 'boxparam4section.dat'
  grille, -1, -1, -1, -1, nx, ny, nz $
          , firstx, firsty, firstz, lastx, lasty, lastz
; the extend the box in the x and y direction in order to maximise
; the number of triangles used to build the section
  firstx = 0 > (firstx - 1)
  lastx = (lastx + 1) < (jpi -1)
  firsty = 0 > (firsty - 1)
  lasty = (lasty + 1) < (jpj -1)

  domdef, firstx, lastx, firsty, lasty, firstz, lastz $
          , /index, gridtype = vargrid

  IF keyword_set(onlybox) THEN return
;
  grille, mask, glam, gphi, gdep, nx, ny, nz $
          , firstx, firsty, firstz, lastx, lasty, lastz

;---------------------
; We define the triangulation which will allows us to determinate the section.
; We recalculate it because it must be defined on the Earth and on oceans.
; Following the direction of the section -rather longitude or rather latitude-
; we define the way to triangulate.
  if strpos(type, 'x') NE -1 then BEGIN
    downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2]
    tri = definetri(nx, ny, (downward)[*])
  ENDIF ELSE tri = definetri(nx, ny)
; If we have an irregular grid that is periodic, then it is possible that
; some of the triangle have a very large size (neighborg points on the
; sphere but far away when doing the projection) and should not be
; taken into account..
  IF keyword_set(key_irregular) AND keyword_set(key_periodic) THEN BEGIN
    glamtri = glam[tri]
    glamtri = abs(glamtri - shift(glamtri, 1, 0))
    good = temporary(glamtri) LT (10.*max(glam)/nx)
    good = where(total(temporary(good), 1) EQ 3)
    tri = (temporary(tri))[*, temporary(good)]
  ENDIF
;---------------------
; Equation of the line on which we do the section.
;---------------------
  abc = linearequation(endpoints[0:1], endpoints[2:3])
  glamtri = glam[tri]
  gphitri = gphi[tri]
; Which points of the triangulation are above and below the line?
  if abc[1] NE 0 THEN $
    test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $
  ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0])

  zero123 = total(test, 1)
; to keep: triangles of the triangulation which are over the line.
  tokeep1 = where(zero123 EQ 1)
  tokeep2 = where(temporary(zero123) EQ 2)
  tokeep = [tokeep1, tokeep2]

  test = test[*, tokeep]
  tri = tri[*, tokeep]
; Which summit of the triangle is alone in a side of the line?
  single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1)
  single1 = single1-(single1/3)*3
  single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0)
  single2 = single2-(single2/3)*3

  undefine, tokeep
  undefine, tokeep1
  undefine, tokeep2
  undefine, test

  single = [temporary(single1), temporary(single2)]
; points1 the point (of the triangle) alone in a side of the line.
; point2 the other point of the triangle in the other side of the line.
  point1 = [single, single]
  point2 = [single EQ 0, 1 + (single LE 1)]

  undefine,  single

  ntri = (size(tri))[2]
  index = [lindgen(ntri), lindgen(ntri)]

  points1 = tri[point1, index]
  points2 = tri[point2, temporary(index)]
; points : complex containing couples of points in a side and the other
; side of the line. We have to delete duplicates.
  points = dcomplex(points1, points2)
  points = points[uniq(points, sort(points))]
  symetrique = dcomplex(imaginary(points), double(points))
  points = points[where(points-shift(temporary(symetrique), 1) NE 0)]
; points1 coordinates of the point of the triangle which is alone in a side of the line.
; point2 coordinates of the other point of the triangle in the other side of the line.
  points1 = complex(glam[   double(points)], gphi[   double(points)])
  points2 = complex(glam[imaginary(points)], gphi[imaginary(points)])
; droites equations of line whose we look for the intersection wit the section.
  droites = linearequation(points1, points2)
  inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) )

; Geographic coordinates of points we look for on the section.
  glamaxe = float(inter)
  gphiaxe = imaginary(inter)
; We arrange them in the growing order between boundaries of the section.
  if strpos(type, 'x') NE -1 then BEGIN
    sort = sort(glamaxe)
    glamaxe = glamaxe[sort]
    inbox = where(glamaxe GE lon1 AND glamaxe LE lon2)
    glamaxe = glamaxe[inbox]
    sort = sort[inbox]
    gphiaxe = gphiaxe[sort]
  ENDIF ELSE BEGIN
    sort = sort(gphiaxe)
    gphiaxe = gphiaxe[sort]
    inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2)
    gphiaxe = gphiaxe[inbox]
    sort = sort[inbox]
    glamaxe = glamaxe[sort]
  ENDELSE
  points = points[sort]
  points1 = points1[sort]
  points2 = points2[sort]
  inter = inter[sort]
  poids = abs(points2-inter)/abs(points2-points1)

;---------------------
  array = litchamp(field)
  array = fitintobox(array)
  if array[0] EQ -1 THEN BEGIN
    res = -1
    return
  ENDIF
  if n_elements(valmask) EQ 0 THEN valmask = 1e20
  taille = size(array)
  if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN
    direc = 't'
    array = grossemoyenne(array, 't')
    taille = size(array)
    jpt = 1
  ENDIF
  case 1 of
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
    taille[0] EQ 2:BEGIN
      value1 = array[double(points)]
      terre = where(value1 GT valmask/10)
      if terre[0] NE -1 then value1[terre] = !values.f_nan
      value2 = array[imaginary(points)]
      terre = where(value2 GT valmask/10)
      if terre[0] NE -1 then value2[terre] = !values.f_nan
      res = poids*value1+(1-poids)*value2
    END
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
    taille[0] EQ 3 AND jpt EQ 1:BEGIN
      npoints = n_elements(points)
      index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz)
      value1 = array[index]
      terre = where(value1 GT valmask/10)
      if terre[0] NE -1 then value1[terre] = !values.f_nan
      index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz)
      value2 = array[index]
      terre = where(value2 GT valmask/10)
      if terre[0] NE -1 then value2[terre] = !values.f_nan
      poids = poids#replicate(1, nz)
      res = poids*value1+(1-poids)*value2
; average following z ?
      if strpos(type, 'z') EQ -1 then begin
        nan = where(finite(res) EQ 0)
        if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt]
        weight = replicate(1, npoints)#e3
        if nan[0] NE -1 then weight[nan] = !values.f_nan
        totalweight = total(weight, 2, /nan)
        zero = where(totalweight EQ 0)
        if zero[0] NE -1 then totalweight[zero] = !values.f_nan
        res = total(res*weight, 2, /nan)/totalweight
        direc = 'z'+string(byte(testvar(var = toto)))
      endif
    END
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
    taille[0] EQ 3 AND jpt NE 1:BEGIN
      npoints = n_elements(points)
      index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt)
      value1 = array[index]
      terre = where(value1 GT valmask/10)
      if terre[0] NE -1 then value1[terre] = !values.f_nan
      index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt)
      value2 = array[index]
      terre = where(value2 GT valmask/10)
      if terre[0] NE -1 then value2[terre] = !values.f_nan
      poids = poids#replicate(1, jpt)
      res = poids*value1+(1-poids)*value2
    END
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
    taille[0] EQ 4:BEGIN
      npoints = n_elements(points)
      index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt)
      index = reform(index, npoints, nz, jpt, /over)
      value1 = array[index]
      terre = where(value1 GT valmask/10)
      if terre[0] NE -1 then value1[terre] = !values.f_nan
      index = imaginary(points)#replicate(1, nz*jpt)+replicate(nx*ny, npoints)#lindgen(nz*jpt)
      index = reform(index, npoints, nz, jpt, /over)
      value2 = array[index]
      terre = where(value2 GT valmask/10)
      if terre[0] NE -1 then value2[terre] = !values.f_nan
      poids = poids#replicate(1, nz*jpt)
      poids = reform(poids, npoints, nz, jpt, /over)
      res = poids*value1+(1-poids)*value2
; average following z ?
      if strpos(type, 'z') EQ -1 then begin
        nan = where(finite(res) EQ 0)
        if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt]
        weight = replicate(1, npoints)#e3
        weight = weight[*]#replicate(1, jpt)
        weight = reform(weight, npoints, nz, jpt, /over)
        if nan[0] NE -1 then weight[nan] = !values.f_nan
        totalweight = total(weight, 2, /nan)
        zero = where(totalweight EQ 0)
        if zero[0] NE -1 then totalweight[zero] = !values.f_nan
        res = total(res*weight, 2, /nan)/totalweight
        direc = 'z'+string(byte(testvar(var = toto)))
      endif
    END
  endcase
;---------------------

  terre = where(finite(res) EQ 0)
  if terre[0] NE -1 then res[terre] = valmask

  if n_elements(showbuild) then BEGIN
    winsave = !window
    psave = !p
    xsave = !x
    ysave = !y
    plt, findgen(nx, ny), /nodata, /nofill, /rempli, title = '', subtitle = '' $
         , coast_thick = 2, window = showbuild
    !p.title = ''
    !p.subtitle = ''

    plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50
    plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50, psym = 2, thick = 2

    FOR i = 0, n_elements(points1)-1 DO $
      plots, [float(points1[i]), float(points2[i])] $
             , [imaginary(points1[i]), imaginary(points2[i])], color = 150

    plots, float(points1), imaginary(points1), color = 150, psym = 1
    plots, float(points2), imaginary(points2), color = 150, psym = 1
    plots, float(inter), imaginary(inter), color = 250, psym = 1

;  ?? bug ??    IF terre[0] NE -1 THEN plots, float(terre[inter]), imaginary(terre[inter]), color = 0, psym = 1

;      dummy = ''
;      read, dummy,  prompt = 'press return to continue'
    IF !d.name EQ 'PS' THEN erase ELSE wset, winsave
    !p = psave
    !x = xsave
    !y = ysave
  ENDIF

  restoreboxparam, 'boxparam4section.dat'

;---------------------
  return
end