;*************************************************
; tdpack_3.ncl
;************************************************
;
; This file is loaded by default in NCL V6.2.0 and newer
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"

; 
; Function for generating random data.
;
function dsrnd1(ifrst,nextn)
begin
  MPLIER = 16807
  MODLUS = 2147483647
  MOBYMP = 127773
  MOMDMP = 2836
  JSEED  = 123456789

  if (ifrst .eq. 0) then
    nextn = JSEED
    ifrst = 1
  end if

  hvlue = nextn / MOBYMP
  lvlue = nextn % MOBYMP
  testv = MPLIER*lvlue - MOMDMP*hvlue

  if (testv .gt. 0) then
    nextn = testv
  else
    nextn = testv + MODLUS
 end if

  return((1.*nextn)/(1.*MODLUS))
end


begin
;N       = 1331
;N       = 133
 N       = 12
;NEAREST =  500
;NEAREST =  50
 NEAREST =  N
 MTRI    = 150000
 NFARTHER = N - NEAREST

;
; Create our input and work arrays.
;
  x = new(N,float)
  y = new(N,float)
  z = new(N,float)
  rtri = new((/MTRI,10/),float)
  rtwk = new((/2,MTRI/),float)

;
; Fill up the dummy input arrays.
;
  ifrst = 0
  nextn = 0
  do i = 0,N-1
    x(i) = dsrnd1(ifrst,nextn)    
    y(i) = dsrnd1(ifrst,nextn)
    z(i) = dsrnd1(ifrst,nextn)    
  end do
  print(x)
  print(y)
  print(z)

;
;  Specify the reference point from which we want to find the NEAREST
;  nearest points.
;
  px = 0.5
  py = 0.5
  pz = 0.5

  wks = gsn_open_wks("pdf","tdpack")     ; send graphics to PNG file

;
; Set some TDPACK parameters and initialize. These four are viewport
; specifiers.
;
  tdsetp("VPB", 0.09)
  tdsetp("VPT", 0.99)
  tdsetp("VPL", 0.11)
  tdsetp("VPR", 1.00)

  tdinit((/4.6, 3.0, 3.3/), (/0.5, 0.5, 0.5/), (/0.5, 0.5, 2.7/), 0.)

;
;  Set up some colors using the standard Tdpack entry for that.
;
  tdclrs(wks, 1, 0., 0.8, 8, 37, 8)

;
;  Define style indices for shades of gray, green, and red.
;
  tdstrs(1,  8, 37,   8,  37, 1, 1, 0, 0.05, 0.05, 0.)
  tdstrs(3,  8, 37,  68,  97, 1, 1, 0, 0.05, 0.05, 0.)
  tdstrs(4,  8, 37,  98, 127, 1, 1, 0, 0.05, 0.05, 0.)

;
;  Store the indices of the nearest points in npts and the complement
;  of that set (with respect to the entire input dataset) in mpts.
;
  npts = new(N,integer)
  mpts = new(N,integer)

  npts(0) = shgetnp(px,py,pz,x,y,z,0)
  do i=0,N-1
      npts(i) = shgetnp(px,py,pz,x,y,z,1)
  end do
  print(npts)
  npts(N-1) = 4

;
;  Plot the near points in green.
;
  ntri = 0
  dotsize = 0.035

  ; Assuming N is the number of points
  letters = (/"B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L"/)  ; Adjust length as needed

  ; Set the font height
  font_height = 0.02

  ; Loop through the points
  ntri = 0
  dotsize = 0.03
  do i = 0, N-1
    print("i = " + i)
    print("npts = " + npts(i))
    tdmtri(-5, (/x(npts(i)), y(npts(i)), z(npts(i))/), dotsize, \
           rtri, ntri, 4, (/0.,0.,0./),(/1.,1.,1./))
 ; Combine the coordinates into a single array for tdprpt
    uvwi3d = (/x(i), y(i), z(i)/)

    ; Project the 3D coordinates to 2D
    xy_proj = tdprpt(uvwi3d)

;   ; Overlay the corresponding letter at the projected coordinates
;   gsn_text(wks, "NCL", tostring(letters(i)), xy_proj(0), xy_proj(1), font_height)
   ; Draw the corresponding letter at the projected coordinates using gtext
    gsn_text_ndc(wks, tostring(letters(i)), xy_proj(0), xy_proj(1), font_height)
  end do

;
;  Mark the reference point in red.
;
 ;tdmtri(-5,(/px,py,pz/),1.2*dotsize,rtri,ntri,3,(/0.,0.,0./),(/1.,1.,1./))

;
;  Draw
;
  itwk = tdotri(rtri, ntri, rtwk, 0)
  tddtri(wks,rtri, ntri, itwk)

;
;  Draw a box around the perimeter.
; removed top line
; tdgrds(wks,(/0., 1., 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)

; back wall vertical right and bottom
; tdgrds(wks,(/0., 1., 0./), (/1., 0.01, 1.0/), (/-1., -1., -1./),12,1)
; tdgrds(wks,(/0., 1., 0./), (/1., 0.20, 1.0/), (/-1., -1., -1./),12,1)
; tdgrds(wks,(/0., 1., 0./), (/1., 0.40, 1.0/), (/-1., -1., -1./),12,1)
; tdgrds(wks,(/0., 1., 0./), (/1., 0.60, 1.0/), (/-1., -1., -1./),12,1)
; tdgrds(wks,(/0., 1., 0./), (/1., 0.80, 1.0/), (/-1., -1., -1./),12,1)
; by changing 1 to 0, removed near vertical and horizontal lines
  tdgrds(wks,(/0., 0., 0./), (/1., 0.001, 1.0/), (/-1., -1., -1./),12,1)
  tdgrds(wks,(/0., 0., 0./), (/1., 0.20, 1.0/), (/-1., -1., -1./),12,1)
  tdgrds(wks,(/0., 0., 0./), (/1., 0.40, 1.0/), (/-1., -1., -1./),12,1)
  tdgrds(wks,(/0., 0., 0./), (/1., 0.60, 1.0/), (/-1., -1., -1./),12,1)
  tdgrds(wks,(/0., 0., 0./), (/1., 0.80, 1.0/), (/-1., -1., -1./),12,1)
;
  tdgrds(wks,(/0., 0., 0./), (/1., 1.00, 1.0/), (/-1., -1., -1./),12,1)
;
; tdgrds(wks,(/0.8, 1., 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)
; back wall left vertical
  tdgrds(wks,(/0.8, 0., 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0.6, 0., 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0.4, 0., 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0.2, 0., 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0.001, 0., 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)
; front and top vertical
; tdgrds(wks,(/0., 0.5, 0./), (/1., 0., 1./), (/-1., -1., -1./),12,0)
; tdgrds(wks,(/0., 1., 0./), (/1., 0., 0.0/), (/-1., -1., -1./),12,1)
; this is box within 
; tdgrds(wks,(/0., 1., 0./), (/1., 0., 1./), (/0., 0., 0./),12,0)
; back wall left horizontal
  tdgrds(wks,(/0.0, 0., 0.8/), (/1., 0., 1./), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0.0, 0., 0.6/), (/1., 0., 1./), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0.0, 0., 0.4/), (/1., 0., 1./), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0.0, 0., 0.2/), (/1., 0., 1./), (/-1., -1., -1./),12,0)
; back wall right horizontal

  tdgrds(wks,(/0., 1., 0./), (/0.0, 0., 0.2/), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0., 1., 0./), (/0.0, 0., 0.4/), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0., 1., 0./), (/0.0, 0., 0.6/), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0., 1., 0./), (/0.0, 0., 0.8/), (/-1., -1., -1./),12,0)
; bottom grid x direction
  tdgrds(wks,(/0., 1., 0./), (/0.8, 0., 0.0/), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0., 1., 0./), (/0.6, 0., 0.0/), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0., 1., 0./), (/0.4, 0., 0.0/), (/-1., -1., -1./),12,0)
  tdgrds(wks,(/0., 1., 0./), (/0.2, 0., 0.0/), (/-1., -1., -1./),12,0)


  
  ; Add vertical lines from points to x-y plane
; Define resources for line appearance
  line_res = True
  line_res@gsLineThicknessF = 2.0  ; Increase the line thickness (default is 1.0)
  do i = 0, N-1
    
  uvw1 = (/x(i), y(i), z(i)/)  ; Starting point in 3D
  uvw2 = (/x(i), y(i), 0.0/)   ; End point at the x-y plane
  tdline(wks, uvw1, uvw2)      ; Draw the vertical line

  end do
  frame(wks) 
end