c
c     Copyright (c) 1999,2000, Compaq Computer Corporation
c
c     See the file "license.terms" for information on usage and redistribution
c     of this file, and for a DISCLAIMER OF ALL WARRANTIES.
c
c *
c * PLOT3D.C
c * Graphics Utility layer -- 3D 
c *
c * Routines:
c *
c *     surf_set_view()
c *	surf_set_window_3d (umin, umax, vmin, vmax, wmin, wmax)
c *     surf_init_xform_3d
c *     surf_calculate_ndc_transform_3d
c *     surf_set_parallel_transform_3d
c *     surf_tr_3d_2d(start)
c *	surf_xform_line_3d
c *     surf_distance_3d
c *     surf_get_min
c *     surf_get_max
c *     surf_set_limits
c *     surf_scale_array
c *
      subroutine surf_set_view
      implicit none
c
c     surf_set_view() calculates the viewing transformation matrix
c
      include 'const.h'
      include 'surf.cmn'
      include 'plot3d.cmn'
c
      call surf_set_window_3d (surf_view%umin,surf_view%umax,
     +                         surf_view%vmin,surf_view%vmax,
     +                         surf_view%wmin,surf_view%wmax)
      call surf_init_xform_3d
      call surf_calculate_ndc_transform_3d
c     Port only supports parallel projection
      call surf_set_parallel_transform_3d
      return
      end
c
      subroutine surf_set_window_3d (umin, umax, vmin, vmax, wmin, wmax)
      implicit none
c
      include 'const.h'
      include 'surf.cmn'
      include 'plot3d.cmn'
c
      real	umin, umax
      real	vmin, vmax
      real	wmin, wmax
c      
      u%delta = umax - umin
c      
      if (u%delta .LT. 0) then
          u%delta = (- u%delta)
          u%min = umax
          u%max = umin
      else
          u%min = umin
          u%max = umax
      endif
c      
      if (u%delta .LE. TOLERANCE) then
          call surf_error('Delta u is too small, quitting',0,1)
      endif
      
      v%delta = vmax - vmin
      
      if (v%delta .LT. 0) then
          v%delta = (- v%delta)
          v%min = vmax
          v%max = vmin
      else
          v%min = vmin
          v%max = vmax
      endif
      
      if (v%delta .LE. TOLERANCE) then
          call surf_error('Delta v is too small, quitting.',0,1)
      endif
      
      w%delta = wmax - wmin
      
      if (w%delta .LT. 0) then
          w%delta = (- w%delta)
          w%min = wmax
          w%max = wmin
      else
          w%min = wmin
          w%max = wmax
      endif
      
      if (w%delta .LE. TOLERANCE) then
          call surf_error('Delta w is too small, quitting.',0,1)
      endif
c
      return
      end
c
      subroutine surf_init_xform_3d
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'surf.cmn'
      include 'plot3d.cmn'
      include 'io.cmn'
c      
      eye_point(1) = surf_eye%xposition
      eye_point(2) = surf_eye%yposition
      eye_point(3) = surf_eye%zposition
c      
      vp_normal(1) = surf_view%xnormal
      vp_normal(2) = surf_view%ynormal
      vp_normal(3) = surf_view%znormal
c      
      vp_up(1) = surf_view%xupvec
      vp_up(2) = surf_view%yupvec
      vp_up(3) = surf_view%zupvec
c      
      vp_distance = surf_view%distance
      
      if (vp_distance .LT. TOLERANCE) then
          write(luttyo,100) vp_distance
100       format(' ','Viewplane distance cannot be negative or zero.'/
     +           ' ','You specified ',1pe12.3)
          vp_distance = 2*TOLERANCE
      endif
c
      return
      end
c      
      subroutine surf_calculate_ndc_transform_3d
      implicit none
c      
c     calculate_ndc_transform() checks for the case that the viewplane 
c     normal vector is parallel to the viewplane up vector
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'surf.cmn'
      include 'plot3d.cmn'
      include 'io.cmn'
c
      real a
      real scale(4)
      real vpn(4), vup(4)
      real offset(4)
      real surf_distance_3d
c      
c     calculate ndc_xform_3d from the window and perspective data
c
c     Only support parallel projection
c     start with a unity xform
      call surf_copy_transform_3d (surf_unity_xform_3d,
     1                             surf_ndc_xform_3d)
c      
c     make the origin the eyepoint
      offset(1) = (- eye_point(1))
      offset(2) = (- eye_point(2))
      offset(3) = (- eye_point(3))
c      
      call surf_translate_transform_3d (offset, surf_ndc_xform_3d)
c       
c     rotate about X until VPN and W are in the same plane
      a = surf_distance_3d (vp_normal(2), vp_normal(3), 0.0)
c      
      if (ABS(a) .LT. TOLERANCE) then
c        don't do anything, already in the plane
         continue
      else 
         call surf_rotate_x_transform_3d(vp_normal(3)/a,vp_normal(2)/a,
     +                                   surf_ndc_xform_3d)
      endif
c      
c     rotate about Y until VPN and W are collinear
      call surf_transform_ray_3d (vp_normal, surf_ndc_xform_3d, vpn)
c      
      a = surf_distance_3d (vpn(1), vpn(3), 0.0)
c      
      if (ABS(a) .LT. TOLERANCE) then 
         continue  ! don't do anything, already collinear
      else 
          call surf_rotate_y_transform_3d(vpn(3)/a, vpn(1)/a,
     +         surf_ndc_xform_3d)
      endif
c      
c     rotate about z until VUP and V are coplaner
      call surf_transform_ray_3d (vp_up, surf_ndc_xform_3d, vup)
c      
      a = surf_distance_3d (vup(1), vup(2), 0.0)
c      
      if (ABS(a) .LT. TOLERANCE) then
         write(luttyo,100) vp_up(1), vp_up(2), vp_up(3),
     +                     vup(1), vup(2), vup(3)
100      format(' ','Viewplane UP vector too close to the viewplane ',
     +              'NORMAL vector.'/
     +          ' ','Original UP vector [',1p3e12.3,']'/
     +          ' ','Transformed UP vector [',1p3e12.3,']')
          call surf_error('SURF terminating',0,1)
      endif
c      
      call surf_rotate_z_transform_3d(vup(2)/a, vup(1)/a,
     1                                surf_ndc_xform_3d)
c      
c     Going from a right to a left coordinate system.
      call surf_right_left(surf_ndc_xform_3d)
c
c     Only support parallel projection:
c     No shear needed: the direction of sight is the same as the w direction
c     Might need to modify it if we have oblique projection ?
c     Next step is transforming view volume into the unit cube
      offset(1) = - u%min
      offset(2) = - v%min
      offset(3) = - w%min
      call surf_translate_transform_3d (offset, surf_ndc_xform_3d)
      scale(1) = 1.0 / (u%max - u%min)
      scale(2) = 1.0 / (v%max - v%min)
      scale(3) = 1.0 / (w%max - w%min)  
      call surf_scale_transform_3d (scale, surf_ndc_xform_3d)  
      frontplane_ndc = w%min
      backplane_ndc  = w%max
c
c     everything applied but perspective.
c     we will do depth clipping before applying perspective.
c
      return
      end
c
      subroutine surf_set_parallel_transform_3d
      implicit none
c
c     parallel perspective transformation
c
      include 'const.h'
      include 'surf.cmn'
      include 'plot3d.cmn'
c      
      call surf_copy_transform_3d (surf_unity_xform_3d,
     +                             surf_perspective_xform_3d)
c      
      surf_perspective_xform_3d(3,3) = 0.0  ! no more z coordinate
c
      return
      end
c      
      subroutine surf_tr_3d_2d(x3d, y3d, z3d, x2d, y2d, z2d, n)
      implicit none
c      
c     A routine to transform  3D points into NDC points doing the ndc 
c     transformation and the perspective transformation. The input to this 
c     function have been changed so that it can take an array of 3D points
c     coordinates and change it to an array of 2D points by the active 
c     transformations thus eliminating the repetitive calling of this function.
c      
      integer n
      real    x3d(n),y3d(n),z3d(n)
      real    x2d(n),y2d(n),z2d(n)
c
      include 'const.h'
      include 'surf.cmn'
      include 'plot3d.cmn'
c      
      real pt_3d(4), pt_2d(4)
      integer i
      
      do i=1,n
         pt_3d(1)  = x3d(i)    
         pt_3d(2)  = y3d(i)    
         pt_3d(3)  = z3d(i)    
c
c        We're only supporting parallel projection. That projection
c        only zero's the z coordinate. So omit call. jvf
         call surf_transform_point_3d(pt_3d, surf_ndc_xform_3d, pt_2d)
CC
CC         call transform_point_3d(pt_temp, surf_perspective_xform_3d,
CCc     +                           pt_2d)
c      
          x2d(i) = pt_2d(1)
          y2d(i) = pt_2d(2)
          z2d(i) = 0.0
      end do
c
      return
      end
c
      subroutine surf_xform_line_3d(xstart, ystart, zstart,
     +                              xend, yend, zend,
     +                              xstart_2d, ystart_2d,
     +                              xend_2d, yend_2d)
      implicit none
c
c     Transforms a point in 3-D space to one in the view plane.
c
c     (1) change to NDC
c     (2) perform z (w -- depth) clipping - leave the other
c         two axes to the 2D package, after perspective xform.
c     (3) apply perspective transformation
c
      include 'const.h'
      include 'surf.cmn'
      include 'plot3d.cmn'
c
      real     xstart,ystart,zstart,xend,yend,zend
      real     xstart_2d,ystart_2d,xend_2d,yend_2d
c      
      real     start(4),end(4)
      real s(4), e(4)    ! start and end points, clipped
      real ss(4), ee(4)  ! start and end points, after perspective
c      
      start(1)  = xstart
      start(2)  = ystart
      start(3)  = zstart
      end(1)    = xend
      end(2)    = yend
      end(3)    = zend
c      
c     Change to NDC
c     Starting point
      call surf_transform_point_3d (start, surf_ndc_xform_3d, s)
c     End point
      call surf_transform_point_3d (end, surf_ndc_xform_3d, e)
c      
CC  Depth clipping omitted, since we don't use it anyway (front/backplane_ndc
CC  are set very large). jvf
CCc     Depth clipping
CCc     First check for equal-Z lines ( s(3) is z coordinate)
CC      if (s(3) .EQ. e(3)) then
CC         if (s(3) .LT. frontplane_ndc .OR. s(3) .GT. backplane_ndc) then
CCc           outside range, do nothing
CC            return
CC         endif      
CC      else
CC         if (s(3) .LE. e(3)) then
CCc           s is in front of e
CC            flip = .FALSE.
CC            do i=1,4
CC               f(i) = s(i)
CC               b(i) = e(i)
CC            end do
CC         else
CCc           s is behind e
CC            flip = .TRUE.
CC            do i=1,4
CC               f(i) = e(i)
CC               b(i) = s(i)
CC            end do
CC         endif
CCc      
CC         if ( (f(3) .GT. backplane_ndc) .OR.
CC     +        (b(3) .LT. frontplane_ndc) ) then
CC            return
CC         endif
CCc      
CC         if (f(3) .LT. frontplane_ndc) then
CCc           Must compute clipping intersection
CC            temp = (frontplane_ndc - f(3)) / (b(3) - f(3))
CC            f(1) = f(1) + (b(1) - f(1)) * temp
CC            f(2) = f(2) + (b(2) - f(2)) * temp
CC            f(3) = frontplane_ndc
CC         endif
CC         if (b(3) .GT. backplane_ndc) then
CCc           Must compute clipping
CC            temp = (b(3) - backplane_ndc) / (b(3) - f(3))
CC            b(1) = b(1) - (b(1) - f(1)) * temp 
CC            b(2) = b(2) - (b(2) - f(2)) * temp
CC            b(3) = backplane_ndc
CC         endif
CC         if( flip ) then
CC            do i=1,4
CC               s(i) = b(i)
CC               e(i) = f(i)
CC            end do
CC         else
CC            do i=1,4
CC               s(i) = f(i)
CC               e(i) = b(i)
CC            end do
CC         endif
CC      endif
CCc     End of ELSE part of IF s(3) == e(3)
c      
c     we now have the two points, clipped to the frontplane and backplane
c     now transform into 2D space
      call surf_transform_point_3d (s, surf_perspective_xform_3d, ss)
      call surf_transform_point_3d (e, surf_perspective_xform_3d, ee)
c      
      xstart_2d  = ss(1)
      ystart_2d  = ss(2)
      xend_2d    = ee(1)
      yend_2d    = ee(2)
c
      return
      end
c
      real function surf_distance_3d (dx, dy, dz)
      implicit none
c
c     surf_distance_3d returns the euclidean distance for the
c     three offsets in x, y and z
c
      real   dx, dy, dz
c
      surf_distance_3d = SQRT(dx*dx + dy*dy + dz*dz)
      return
      end
c      
      real function surf_get_min(x,n)
      implicit none
c
c     surf_get_min() gets the minimum value of an array
c      
      integer  n
      real x(n)
c
      real min
      integer i
c         
      min  = x(1)
      do i=2,n
          if(x(i) .LT. min)  min = x(i)
      end do
c
      surf_get_min = min
      return
      end
c      
      real function surf_get_max(x,n)
      implicit none
c
c     surf_get_max() gets the maximum value of an array
c      
      integer  n
      real x(n)
c      
      real max
      integer i
c         
      max  = x(1)
      do i=2,n
          if (x(i) .GT. max)  max = x(i)
      end do
      surf_get_max = max
      return
      end
c      
      subroutine surf_set_limits
      implicit none
c
c     surf_set_limits() set the minimum and maximum of the axis.
c
      include 'const.h'
      include 'surf.cmn'
c
      real round_up_mant, round_down_mant, surf_get_min, surf_get_max
      
      if (.NOT. surf_xaxis%min_flag) then
          surf_xaxis%min = surf_get_min(surf_curve%xposition,
     +                                  surf_curve%npoints)
          surf_xaxis%min = round_down_mant(surf_xaxis%min)
      endif
      if (.NOT. surf_xaxis%max_flag) then
          surf_xaxis%max = surf_get_max(surf_curve%xposition,
     +                                  surf_curve%npoints)
          surf_xaxis%max = round_up_mant(surf_xaxis%max)
      endif
      if (.NOT. surf_yaxis%min_flag) then
          surf_yaxis%min = surf_get_min(surf_curve%yposition,
     +                                  surf_curve%npoints)
          surf_yaxis%min = round_down_mant(surf_yaxis%min)
      endif
      if (.NOT. surf_yaxis%max_flag) then
          surf_yaxis%max = surf_get_max(surf_curve%yposition,
     +                                  surf_curve%npoints)
          surf_yaxis%max = round_up_mant(surf_yaxis%max)
      endif
      if (.NOT. surf_zaxis%min_flag) then
          surf_zaxis%min = surf_get_min(surf_curve%zposition,
     +                                  surf_curve%npoints)
          surf_zaxis%min = round_down_mant(surf_zaxis%min)
      endif
      if (.NOT. surf_zaxis%max_flag) then
          surf_zaxis%max = surf_get_max(surf_curve%zposition,
     +                                  surf_curve%npoints)
          surf_zaxis%max = round_up_mant(surf_zaxis%max)
      endif
c
c     Now scale the data between 0.0 and 1 (min=0.0,max=1.0)
c
      call surf_scale_array(0.0, 1.0, surf_xaxis%min, surf_xaxis%max, 
     +                 surf_curve%xposition,surf_curve%scaled_xposition,
     +                 surf_curve%npoints,surf_xaxis%log)
      call surf_scale_array(0.0, 1.0, surf_yaxis%min, surf_yaxis%max,
     +                 surf_curve%yposition,surf_curve%scaled_yposition,
     +                 surf_curve%npoints,surf_yaxis%log)
      call surf_scale_array(0.0, 1.0, surf_zaxis%min, surf_zaxis%max,
     +                 surf_curve%zposition,surf_curve%scaled_zposition,
     +                 surf_curve%npoints,surf_zaxis%log)
c
      return
      end
c
      subroutine surf_scale_array(wmin,wmax,min,max,array,scaled_array,
     +                            n,log)
      implicit none
c
c     function that scaled the values of a certain array to a range specified
c     by the wmin,wmax arguments and stores the results in another array.
c
      real wmin,wmax   ! window to scale to
      real min,max     ! array values to map to window
      integer  n       ! number of data points
      real array(n),scaled_array(n) ! array and scaled array
      logical log
c
      real ratio
      integer  i
      
      if (max .eq. min) then
         ratio = wmax - wmin
      else
         if (log) then
            ratio = (wmax - wmin)/(ALOG10(max) - ALOG10(min))
         else 
            ratio = (wmax - wmin)/(max - min)
         endif
      endif
c
      do i=1,n
         if (log) then
            scaled_array(i) = (ALOG10(array(i)) - ALOG10(min))*ratio
     +                        + wmin
         else
            scaled_array(i) = (array(i) - min)* ratio + wmin
         endif
      end do
c
      return
      end
