source: trunk/TRIANGULATION/section.pro @ 2

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:
;
; PURPOSE:
;
; CATEGORY:
;
; CALLING SEQUENCE:
; 
; INPUTS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------

PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints, BOITE = boite, TYPE = type, DIREC = direc
@common
; 
;---------------------
   array = litchamp(field)
;---------------------
; definition de boite en fonction de endpoints
; puis redefinition du domaine
   boite2d = [min([endpoints[0], endpoints[2]]), max([endpoints[0], endpoints[2]]), min([endpoints[1], endpoints[3]]), max([endpoints[1], endpoints[3]])]
;
   minprof = 0
   profdefault = 200
   if n_elements(type) EQ 0 then type = 'nothing'
   Case N_Elements(Boite) OF
      1:boite=[boite2d, minprof, boite[0]]
      2:boite=[boite2d, boite[0],boite[1]]
      Else:$ 
       if strpos(type, 'z') NE -1 THEN boite=[boite2d, minprof, profdefault] $
      ELSE boite=[boite2d, prof1, prof2]
   ENDCASE
   if strpos(type, 'z') NE -1 THEN BEGIN
      !y.range = [boite[5], boite[4]] ;on garde les yranges (axe z) avant de changer la boite.
      profmax = boite[5]
      vraiboite = boite
      if vargrid EQ 'W' then gdep = gdepw ELSE gdep = gdept
; check with vertical grid limits (nearest level)
      gwork = gdep
; check the increse or decrese of the z axis
      IF gwork[1] LE gwork[0] THEN gwork = reverse(gdep, 1)
      niveauprof = where(gwork ge boite[5]) & niveauprof = niveauprof[0]
      if niveauprof NE -1 then boite[5] = gwork[niveauprof]+1
   ENDIF
   domdef, boite, GRILLE=[vargrid], /findalways, _extra = ex
; on redefinit lon1, ... au cas ou findalways ait ete utilise ds domdef
   lon1 = min([endpoints[0], endpoints[2]])
   lon2 = max([endpoints[0], endpoints[2]])
   lat1 = min([endpoints[1], endpoints[3]])
   lat2 = max([endpoints[1], endpoints[3]])
; on recupere la grille sur la boite
   grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz
;---------------------
; on definit la triangulation qui va nous permetre de determiner la
; section. on la recalcule car elle doit etre definie sur la terre
; aussi bien que sur la mer. suivant le sens de la section -plutot
; longitude ou plutot latitude- on definit la facon de trianguler
   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)
;---------------------
; equation de la droite suivant laquelle on fait la section
;---------------------
   abc = linearequation(endpoints[0:1], endpoints[2:3])
   glamtri = glam[tri]
   gphitri = gphi[tri]
; quels sont les points de la triangulation qui sont au dessus et au
; dessous de la droite ?
   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 de la triangulation qui sont a cheval sur la droite.
   tokeep1 = where(zero123 EQ 1)
   tokeep2 = where(temporary(zero123) EQ 2)
   tokeep = [tokeep1, tokeep2]

   test = test[*, tokeep]
   tri = tri[*, tokeep]
; quel est le sommet du triangle qui est seul d''un cote de la droite?
   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 le point du triangles qui est seul d''un cote de la droite.
; point2 l''autre point du triangle de l''autre cote de la droite
   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 : complexe contenant les couples de points de part et
; d''autre de la droite. Ils faut supprimer les doublons.
   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 les coordnnees du point du triangles qui est seul d''un cote de la droite.
; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite
   points1 = complex(glam[   double(points)], gphi[   double(points)])
   points2 = complex(glam[imaginary(points)], gphi[imaginary(points)])
; droites les equations des droites dont on cherche l''intersection
; avec la section.
   droites = linearequation(points1, points2)
   inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) )

; les ccordonnes geographiques des points que l''on cherche sur la section.
   glamaxe = float(inter)
   gphiaxe = imaginary(inter)
; on les range ds l''ordre croissant entre les bornes de la 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 = 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
; moyenne suivant z ?
         if strpos(type, 'z') EQ -1 then begin
            nan = where(finite(res) EQ 0)
            if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt]
            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
; moyenne suivant z ?
         if strpos(type, 'z') EQ -1 then begin
            nan = where(finite(res) EQ 0)
            if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt]
            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


;    plt, tmask[*,*,0],/nocouleur, /rempli, title = '', subtitle = ''
;    !p.title=''
;    !p.subtitle=''
;    dessinetri
;    plot,[endpoints[0],endpoints[2]],[endpoints[1],endpoints[3]],/noerase,color=50

;    plot,float(points1),imaginary(points1),/noerase,color=150,psym=1
;    plot,float(points2),imaginary(points2),/noerase,color=150,psym=1
;    plot,float(inter),imaginary(inter),/noerase,color=250,psym=1

;    plot,[float(points1[0])],[imaginary(points1[0])],/noerase,color=150,psym=1
;    plot,[float(points2[0])],[imaginary(points2[0])],/noerase,color=150,psym=1
;    plot,[float(inter[0])],[imaginary(inter[0])],/noerase,color=250,psym=1


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