C                  INTERNATIONAL AVS CENTER
C        (This disclaimer must remain at the top of all files)
C
C WARRANTY DISCLAIMER
C
C This module and the files associated with it are distributed free of charge.
C It is placed in the public domain and permission is granted for anyone to use,
C duplicate, modify, and redistribute it unless otherwise noted.  Some modules
C may be copyrighted.  You agree to abide by the conditions also included in
C the AVS Licensing Agreement, version 1.0, located in the main module
C directory located at the International AVS Center ftp site and to include
C the AVS Licensing Agreement when you distribute any files downloaded from
C that site.
C
C The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module provide absolutely
C NO WARRANTY OF ANY KIND with respect to this software.  The entire risk as to
C the quality and performance of this software is with the user.  IN NO EVENT
C WILL The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module BE LIABLE TO
C ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING,
C WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY
C SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES.
C
C This AVS module and associated files are public domain software unless
C otherwise noted.  Permission is hereby granted to do whatever you like with
C it, subject to the conditions that may exist in copyrighted materials. Should
C you wish to make a contribution toward the improvement, modification, or
C general performance of this module, please send us your comments:  why you
C liked or disliked it, how you use it, and most important, how it helps your
C work. We will receive your comments at avs@ncsc.org.
C
C Please send AVS module bug reports to avs@ncsc.org.
C
C**********************************************************************
C*
C*      file name:      linecont.f
C*
C*      description:    generates iso-line contour geometory.
C*
C*                      Copyright (c) 1991 by
C*                      KUBOTA COMPUTER INC.
C*                      All Rights Reserved
C*
C**********************************************************************/


C/**********************************************************************
C*
C*      function name:  AVSinit_modules
C*
C*      description:    description function
C*
C**********************************************************************/

      	subroutine AVSinit_modules
	include '/usr/avs/include/avs.inc'
	integer iport, op1,op2,op3, iparm

	external contour

        call AVSset_module_name('line_contour', 'mapper')

C Create an input port (name, type) 
	op1 = AVScreate_input_port('line_cont_input',
     &     'field 2D 1-vector real ',REQUIRED)

	op2 = AVScreate_input_port('input field',
     &     'colormap',REQUIRED)

C Create an output port for the result (name, type) 
	op3 = AVScreate_output_port('Geometry', 'geom')

C Add one paramter: the filename of the polygon object 
        ip1 = AVSadd_parameter('On/Off','boolean',1,0,1)
	ip2 = AVSadd_parameter
     &  ('min ', 'float',-0.5,float_unbound,float_unbound) 
	ip3 = AVSadd_parameter
     &   ('max ', 'float',0.5,float_unbound,float_unbound) 
	ip4 = AVSadd_parameter('line_num ', 'integer',2,5,55)
	ip10 = AVSadd_parameter('just', 'oneshot',' ',' ',' ') 

C Tell avs what subroutine to call to do the compute 
        call AVSautofree_output(op3)
	call AVSset_compute_proc(contour)

	return
	end

C/**********************************************************************
C*
C*      function name:  contour
C*
C*      description:    computation function
C*
C*      input:          field 2D scalar float
C*                      colormap
C*
C*      output:         geom
C*
C*      return:         0/1
C*
C**********************************************************************/


C
 	integer function contour(f,nx,ny,coordflag,ncoord,xyz
     &     ,size,lower,upper,hue,sat,val,alpha
     &     ,output
     &     ,ionoff,pmin,pmax,num_line
     &     ,ijust)

	include '/usr/avs/include/avs.inc'
	include '/usr/avs/include/geom.inc'

C input port data for " field 2D 1-vector real irregular"
        dimension f(nx,ny),xyz(*)
	integer*4 nx,ny,coordflag,ncoord

C input port data for " colormap " (optional)
        integer size
	real lower, upper
	real hue(256), sat(256), val(256), alpha(256)

C output port data for " geometry"
	integer output

C option data for min_value,max_value, number of contour line.
	real*4 pmin,pmax
	integer num_line
C In C this is a ptr to a GEOMobj struct.
	integer		obj
C dimension for verts & color 
        dimension verts(6),colors(3,2)
        dimension color_tbl(3,256)
C memory alloc
        real*4 c(1)
        integer pc,oc,dim(2),dim2(3),ialcf
        real*4 w(1)
        integer pw,ow
C
C The first time through, the filname will be NULL.  Return a silent error
C indicating that downstream modules should not be executed.
C
      if( coodflag .eq. irregular ) then
          call AVSwarning(' Only irregular data')
          goto 9999
      endif
      if( AVSparameter_changed('just') ) then
          call minmax(f,nx*ny,ppmin,ppmax)
          halfval = 0.4 * ( ppmin + ppmax)
          call AVSmodify_parameter('max ',
     &       IOR(AVS_VALUE,IOR(AVS_MAXVAL,AVS_MINVAL)),
     &      ppmax,ppmax-abs(halfval),ppmax+abs(halfval))
          call AVSmodify_parameter('min ',
     &       IOR(AVS_VALUE,IOR(AVS_MAXVAL,AVS_MINVAL)),
     &      ppmin,ppmin-abs(halfval),ppmin+abs(halfval))
      endif
c     write(6,*) ' end ione param '
      if( ionoff .eq. 0) then
          write(6,*) ' NULL obj '
          obj =  GEOM_create_obj(GEOM_POLYTRI,GEOM_NULL)
          output = geom_init_edit_list(output)
          call geom_edit_geometry(output,'line_cont',obj)
          call geom_destroy_obj(obj)
          contour = 1
          return
      endif
C
	obj = geom_create_obj(GEOM_POLYTRI,GEOM_NULL)

      do L = 1,256
          call hsv_to_rgb(color_tbl(1,L),color_tbl(2,L),
     &                    color_tbl(3,L),hue(L),sat(L),val(L))
      enddo
      call falloc(nx*ny*6,4,0,w,pw,ow)
      do L = 1,num_line
        target = pmin + float(L-1)*(pmax-pmin)/float(num_line-1)
        indexc = 255*max(0.,min(1.0,(target - pmin)/(pmax-pmin))) + 1
        colors(1,1) = color_tbl(1,indexc)
        colors(2,1) = color_tbl(2,indexc)
        colors(3,1) = color_tbl(3,indexc)
        colors(1,2) = colors(1,1)
        colors(2,2) = colors(2,1)
        colors(3,2) = colors(3,1)
            call con(target,xyz,f,nx,ny,w(ow+1),numl)
c       write(6,*) ' numl = ',numl
        do LL = 1,numl
	  do ii = 1,6
	    verts(ii) = w(ow+ii+6*(LL-1))
	  enddo
          call geom_add_disjoint_line(obj,verts,colors,2,GEOM_COPY_DATA)
        enddo
      enddo

   	output = geom_init_edit_list(output)
	call geom_edit_geometry(output,'line_cont',obj)
   	call geom_destroy_obj(obj)
C
        call free(pw)
	contour = 1
	return

 9999   contour = 0
	return
	end

	subroutine hsv_to_rgb(r,g,b,ht,s,v)
	real f, p, q, t
	real r, g, b
	real ht, s, v
	real h
C Make sure not to trash the input colormap
	h = ht
	if(v.eq.0) then
	 r=0
	 g=0
	 b=0
	 goto 100
	elseif(s.eq.0) then
	 r = v
	 g = v
	 b = v
	 goto 100
	else
	 h = h * 6.0
         if(h.ge.6.0) then
            h = 0.0
         endif

         i = h
         f = h - i
         p = v*(1.0-s)
         q = v*(1.0-s*f)
         t = v*(1.0-s*(1.0-f))
      endif

      if(i.eq.0) then
          r = v
          g = t
          b = p
      elseif(i.eq.1) then
         r = q
         g = v
         b = p
      elseif(i.eq.2) then
         r = p
         g = v
         b = t
      elseif(i.eq.3) then
         r = p
         g = q
         b = v
      elseif(i.eq.4) then
         r = t
         g = p
         b = v
      elseif(i.eq.5) then
         r = v
         g = p
         b = q
      endif
  100 continue
      return
      end

      subroutine minmax(f,imax,ppmin,ppmax)
      real*4 f(*)
      ppmin = f(1)
      ppmax = f(1)
      do i = 2,imax
          ppmin = min(f(i),ppmin)
          ppmax = max(f(i),ppmax)
      enddo
      return
      end
          
  
C
C   Contur search sample program    1989.9.12   
C
C *********************************************************************
C
      subroutine con(al,xyz,t,nx,ny,line,numl)
C
C     In case of triangle element
C
C  al  value of line
C  xp  x cordinate value
C  yp  y cordinate value
C   t  scalar value
C  line /* return value */   line data
C  numl /* return value */   line number
C
c     include 'sample3d.h'
	  real*4 al,xyz(nx,ny,3),t(nx,ny),line(6,*)
	  real*4 w1,w2,w3
	  real*4 xyz12(6)
	  real*4 t1,t2,t3
	  integer*4 numl

	  numl = 0
      do i = 1, nx - 1
        do j = 1, ny - 1
	  w1 = t(i,j) - al
	  w2 = t(i+1,j) - al
	  w3 = t(i,j+1) - al
	  if(w1*w2 .lt. 0.)then
		w4 = t(i+1,j) - t(i,j)
		do L = 1,3
		  xyz12(L) = xyz(i+1,j,L) - w2/w4*(xyz(i+1,j,L)-xyz(i,j,L))
		enddo
		if(w1*w3 .lt. 0.)then
		  w4 = t(i,j+1) - t(i,j)
		  do L = 1,3
		    xyz12(L+3) = xyz(i,j,L) - w1/w4*(xyz(i,j+1,L)-xyz(i,j,L))
		  enddo
		else
		  w4 = t(i,j+1) - t(i+1,j)
		  do L = 1,3
		    xyz12(L+3) = xyz(i+1,j,L) - w2/w4*(xyz(i,j+1,L)-xyz(i+1,j,L))
		  enddo
		endif
		numl = numl + 1
		do L = 1,6
		  line(L,numl) = xyz12(L)
		enddo
	  else
	    if(w1*w3 .lt. 0.)then
              w4 = t(i,j+1) - t(i,j)
	      do L = 1,3
	        xyz12(L) = xyz(i,j,L) - w1/w4*(xyz(i,j+1,L)-xyz(i,j,L))
	      enddo
	      w4 = t(i,j+1) - t(i+1,j)
	      do L = 1,3
	        xyz12(L+3) = xyz(i+1,j,L) - w2/w4*(xyz(i,j+1,L)-xyz(i+1,j,L))
	      enddo
              numl = numl + 1
	      do L = 1,6
		    line(L,numl) = xyz12(L)
	      enddo
            endif 
	  endif
	  w1 = t(i+1,j+1) - al
	  w2 = t(i+1,j) - al
	  w3 = t(i,j+1) - al
	  if(w1*w2 .lt. 0.)then
            w4 = t(i+1,j) - t(i+1,j+1)
	    do L = 1,3
	      xyz12(L) = xyz(i+1,j,L) - w2/w4*(xyz(i+1,j,L)-xyz(i+1,j+1,L))
	    enddo
	    if(w1*w3 .lt. 0.)then
		w4 = t(i,j+1) - t(i+1,j+1)
		do L = 1,3
		  xyz12(L+3) = xyz(i+1,j+1,L) - w1/w4*(xyz(i,j+1,L)-xyz(i+1,j+1,L))
		enddo
		else
		w4 = t(i,j+1) - t(i+1,j)
		do L = 1,3
		  xyz12(L+3) = xyz(i+1,j,L) - w2/w4*(xyz(i,j+1,L)-xyz(i+1,j,L))
		enddo
		endif
		numl = numl + 1
		do L = 1,6
		  line(L,numl) = xyz12(L)
		enddo
	  else
	    if(w1*w3 .lt. 0.)then
	      w4 = t(i,j+1) - t(i+1,j+1)
	      do L = 1,3
	        xyz12(L) = xyz(i+1,j+1,L) - w1/w4*(xyz(i,j+1,L)-xyz(i+1,j+1,L))
	      enddo
	      w4 = t(i,j+1) - t(i+1,j)
	      do L = 1,3
	        xyz12(L+3) = xyz(i+1,j,L) - w2/w4*(xyz(i,j+1,L)-xyz(i+1,j,L))
	      enddo
	      numl = numl + 1
	      do L = 1,6
	        line(L,numl) = xyz12(L)
	      enddo
	    endif 
	  endif
        enddo
      enddo
	  return
	  end