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     Minimum/maximum derivatives utilities
c
      subroutine pm_minmax
      implicit none
c
c     Asks user what array to process
c
c     Programmer: John Faricelli
c     Date: 13-Jun-85
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'bincodes.cmn'
c
c     locals
c
      integer answer, string_len
      character*10 buffer
      character*80 errfmt
c
      errfmt = '('' '')'
c
c     Short circuit this menu if there is only one data item to select from
c
      if( max_menu_num .eq. 1 ) then
         if( is_data_1d(1) ) then
            call clear_screen
            call beep_terminal
            write(luttyo,100)
100         format(' ','Sorry, 1D data not supported yet'/
     1             ' ','Hit return to continue...')
            read(luttyi,106,end=300,err=300) buffer
         else if( data_is_trigrid ) then
            call min_max_array_tri(arrdat(arrind(menu_to_arr(1))),
     1                             tri_npts)
         else
            call min_max_array(arrdat(arrind(menu_to_arr(1))), ni, nj)
         endif
         return
      endif
c
200   call clear_screen
      write(luttyo,errfmt)
      errfmt = '('' '')'
      call write_data_menu('Min/max computation')
      write(luttyo,105)
105   format(/' ','Request? ',$)
      read(luttyi,106,err=300,end=300) buffer
106   format(a)
      call upcase(buffer)
      if( string_len(buffer) .eq. 0 ) then
         answer = EXIT_MENU
      else if( buffer(1:1) .eq. 'Q' .or.  buffer(1:1) .eq. 'E' ) then
         answer = EXIT_MENU
      else
         read(buffer,*,err=300) answer
      endif
c
c    Error check:
c       1) Exit menu
c       2) Answer in range
c       3) Is requested data in binary file?
c
      if( answer .eq. EXIT_MENU ) then
         goto 1000
      else if ( answer .lt. 1 .or. answer .gt. max_menu_num ) then
         errfmt = '('' '',''Invalid input: Please try again...'')'
         call beep_terminal
      else if( .NOT. in_binary_file(answer) ) then
         errfmt = '('' '',''Data not available from binary file'')'
         call beep_terminal
      else if( data_is_trigrid ) then
         call min_max_array_tri(arrdat(arrind(menu_to_arr(answer))),
     1                          tri_npts)
      else if( is_data_1d(answer) ) then
         errfmt = '('' '',''Min/max not implemented for 1D data yet'')'
         call beep_terminal
      else
         call min_max_array(arrdat(arrind(menu_to_arr(answer))), ni, nj)
      endif
      goto 200
c
c     illegal input
300   errfmt = '('' '',''Invalid input: Please try again...'')'
      call beep_terminal
      goto 200
c
1000  continue
      return
      end
c
      subroutine min_max_array(data_array, isize, jsize)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'units.cmn'
c
      integer isize, jsize
      real data_array(isize, jsize)
c
c     locals
c
      character*1 answer
      character*27 prompt
      character*10 format_array(2)
      logical do_zero
      integer i,j, min_i, min_j, max_i, max_j, v_start, getindex
      integer ulen, string_len
      integer min_h, max_h, min_v, max_v
      real min_value, max_value
      real low_h, hi_h, low_v, hi_v, low_temp, hi_temp, real_vec(2),
     1     d1, d2
c
      ulen = string_len(pm_length_unit)
      min_value = PM_MAX_FLOAT
      max_value = -PM_MAX_FLOAT
c
c     Get search bounds
      write(luttyo,50)
50    format(/' ','Note: Some field quantities can be discontinuous'/
     1        ' ','      Field quantities at Y=0.0 (oxide/semiconductor'
     2            ,' boundary)'/
     3        ' ','      are taken in the semiconductor')
c
c     Horizontal limits
c
      prompt = 'Enter ? limits in ' // pm_length_unit(1:ulen) // ': '
      prompt(7:7) = charh
110   real_vec(1) = pm_scale_cm2plot*xyz(1,idisti)
      real_vec(2) = pm_scale_cm2plot*xyz(ni,idisti)
      format_array(1) = '1pg12.5'
      format_array(2) = '1pg12.5'
      call get_reals_with_default(real_vec, 2, prompt, format_array, 0)
      low_temp = real_vec(1)*pm_scale_plot2cm
      hi_temp  = real_vec(2)*pm_scale_plot2cm
      if( hi_temp .le. low_temp ) then
         write(luttyo,111)
111      format(' ','Error: max must be larger than min')
         goto 110
      else if( low_temp .lt. xyz(1, idisti) .or.
     1         hi_temp  .gt. xyz(ni,idisti) ) then
         write(luttyo,112)
112      format(' ','Error: Bounds are outside of the device')
         goto 110
      endif
      low_h = low_temp
      hi_h  = hi_temp
      if( low_h .le. xyz(1,idisti) ) then
         min_h = 1
      else if( low_h .ge. xyz(ni,idisti) ) then
         min_h = ni
      else
         min_h = getindex(low_h, xyz(1,idisti), ni)
c        Find mesh index which is closest to min_h
         if( min_h .ne. ni ) then
            d1 = low_h - xyz(min_h,idisti)
            d2 = xyz(min_h+1,idisti) - low_h
            if( d1 .gt. d2 ) min_h = min_h + 1
         endif
      endif
      if( hi_h .le. xyz(1,idisti) ) then
         max_h = 1
      else if( hi_h .ge. xyz(ni,idisti) ) then
         max_h = ni
      else
         max_h = getindex(hi_h, xyz(1,idisti), ni)
c        Find mesh index which is closest to max_h
         if( max_h .ne. ni ) then
            d1 = hi_h - xyz(max_h,idisti)
            d2 = xyz(max_h+1,idisti) - hi_h
            if( d1 .gt. d2 ) max_h = max_h + 1
         endif
      endif
c
c     Vertical limits
c
      prompt(7:7) = charv
120   real_vec(1) = pm_scale_cm2plot*xyz(1,idistj)
      real_vec(2) = pm_scale_cm2plot*xyz(nj,idistj)
      format_array(1) = '1pg12.5'
      format_array(2) = '1pg12.5'
      call get_reals_with_default(real_vec, 2, prompt, format_array, 0)
      low_temp = real_vec(1)*pm_scale_plot2cm
      hi_temp  = real_vec(2)*pm_scale_plot2cm
      if( hi_temp .le. low_temp ) then
         write(luttyo,121)
121      format(' ','error: max must be larger than min')
         goto 120
      else if( low_temp .lt. xyz(1,idistj) .or.
     1         hi_temp  .gt. xyz(nj,idistj) ) then
         write(luttyo,122)
122      format(' ','Error: Bounds are outside of the device')
         goto 120
      endif
      hi_v = hi_temp
      low_v = low_temp
c
c     Find closest mesh lines
      if( low_v .le. xyz(1,idistj) ) then
         min_v = 1
      else if( low_v .ge. xyz(nj,idistj) ) then
         min_v = nj
      else
         min_v = getindex(low_v, xyz(1,idistj), nj)
c        Find mesh index which is closest to min_v
         if( min_v .ne. nj ) then
            d1 = low_v - xyz(min_v,idistj)
            d2 = xyz(min_v+1,idistj) - low_v
            if( d1 .gt. d2 ) min_v = min_v + 1
         endif
      endif
      if( hi_v .le. xyz(1,idistj) ) then
         max_v = 1
      else if( hi_v .ge. xyz(nj,idistj) ) then
         max_v = nj
      else
         max_v = getindex(hi_v, xyz(1,idistj), nj)
c        Find mesh index which is closest to max_v
         if( max_v .ne. nj ) then
            d1 = hi_v - xyz(max_v,idistj)
            d2 = xyz(max_v+1,idistj) - hi_v
            if( d1 .gt. d2 ) max_v = max_v + 1
         endif
      endif
c
      write(luttyo,130) charh, pm_scale_cm2plot*xyz(min_h,idisti),
     1                         pm_scale_cm2plot*xyz(max_h,idisti),
     2                  charv, pm_scale_cm2plot*xyz(min_v,idistj),
     3                         pm_scale_cm2plot*xyz(max_v,idistj)
130   format(' ','Min/Max taken from mesh lines at:'/
     1       ' ','   ',a1,'= ',1pg12.5,' to ',1pg12.5/
     2       ' ','   ',a1,'= ',1pg12.5,' to ',1pg12.5)
c
      call get_char(answer,
     1     'Exclude zero values in search for minimum? (Y/N) [N] ',0)
      call upcase(answer)
      if( answer .eq. 'Y' ) then
         do_zero = .FALSE.
      else
         do_zero = .TRUE.
      endif
c
c     find minimum/maximum value of array
c
      do 10 j=min_v,max_v
         do 10 i=min_h,max_h
            if( min_value .gt. data_array(i,j) ) then
               if( data_array(i,j) .ne. 0.0 .or. do_zero ) then
                  min_value = data_array(i,j)
                  min_i = i
                  min_j = j
               endif
            endif
            if( max_value .lt. data_array(i,j) ) then
               max_value = data_array(i,j)
               max_i = i
               max_j = j
            endif
10          continue
c
      write(luttyo,200) min_value, pm_scale_cm2plot*xyz(min_i,idisti),
     1                             pm_scale_cm2plot*xyz(min_j,idistj)
200   format(' ','Minimum value ',1pe12.3,' found at (',1pg12.5,',',
     1       1pg12.5,')')
      write(luttyo,210) max_value, pm_scale_cm2plot*xyz(max_i,idisti),
     1                             pm_scale_cm2plot*xyz(max_j,idistj)
210   format(' ','Maximum value ',1pe12.3,' found at (',1pg12.5,',',
     1       1pg12.5,')')
      call get_char(answer,
     1     'Do you want derivatives at min/max? (Y/N) ',0)
      call upcase(answer)
      if( answer .eq. 'Y' ) then
         write(luttyo,225)
225      format(/' ','WARNING: Taking derivatives of simulation ',
     1           'quantities is prone to error.'/
     2           ' ','         Quantitative use of derivative values ',
     3           'should be done with caution.'/
     4           ' ','         Second derivative values, '
     5           'especially of field quantities,'/
     6           ' ','         should only be used qualitatively.')
         write(luttyo,230)
230      format(/' ','At minimum point:')
         call derivatives(data_array,isize,jsize,v_start,min_i,min_j,0)
c        pause
         call get_char(answer,'Hit return for more... ',0)
         write(luttyo,240)
240      format(/' ','At maximum point:')
         call derivatives(data_array,isize,jsize,v_start,max_i,max_j,0)
c        pause
         call get_char(answer,'Hit return for more... ',0)
         write(luttyo,250)
250      format(/' ','At maximum point (with array values inverted):')
         call derivatives(data_array,isize,jsize,v_start,max_i,max_j,1)
      endif
c
c     pause
      call get_char(answer,'Hit return for more... ',0)
      return
      end
c
      subroutine derivatives(array,isize,jsize,v_start,pt_i,pt_j,iflgs)
      implicit none
c
c     computes 1st and 2nd derivatives at point array(pt_i,pt_j)
c
c     Programmer: John Faricelli
c     Date: 13-Jun-85
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
c
      integer isize, jsize, v_start, pt_i, pt_j, iflgs
      real array(isize,jsize)
c
c     locals
      logical can_do_h, can_do_v
      real denomh, denomv
      real dfdh, dfdv, df2dh2, df2dv2, df2dhdv, dfdhjp, dfdhjm
      real term1, term2, term3
      real array_im_jm, array_im_j, array_im_jp
      real array_i_jm, array_i_j, array_i_jp
      real array_ip_jm, array_ip_j, array_ip_jp
      real deltah_i, deltah_im, deltah_sum
      real deltav_j, deltav_jm, deltav_sum
      real invert
c
c     check: do we have two neighboring points?
      if( pt_i .gt. 1 .and. pt_i .lt. ni ) then
         can_do_h = .TRUE.
      else
         can_do_h = .FALSE.
         write(luttyo,100) charh
100      format(' ','Point is on left or right boundary;'
     1              ' can''t compute ',a1,' derivatives')
      endif
      if( pt_j .gt. v_start .and. pt_j .lt. nj ) then
         can_do_v = .TRUE.
      else
         can_do_v = .FALSE.
         write(luttyo,110) charv
110   format(' ','Point is on top or bottom boundary;',
     1      ' can''t compute ',a1,' derivatives')
      endif
c
      if( can_do_h ) then
c        some constants...
         if( iflgs .eq. 1 ) then
            array_im_j  = invert(array(pt_i-1,pt_j))
            array_i_j   = invert(array(pt_i,pt_j))
            array_ip_j  = invert(array(pt_i+1,pt_j))
         else
            array_im_j  = array(pt_i-1,pt_j)
            array_i_j   = array(pt_i,pt_j)
            array_ip_j  = array(pt_i+1,pt_j)
         endif
         deltah_im = xyz(pt_i,idisti) - xyz(pt_i - 1,idisti)
         deltah_i = xyz(pt_i+1,idisti) - xyz(pt_i,idisti)
         deltah_sum = deltah_im + deltah_i
         denomh = 1.0/(deltah_i*deltah_im*(deltah_sum))
c
c        do second derivative
         df2dh2 = (2.0*( array_im_j *deltah_i -
     1                   array_i_j  *deltah_sum +
     1                   array_ip_j *deltah_im ) )*denomh
c        do first derivative
         term1 = ( array_ip_j*deltah_im ) /
     1           ( deltah_i*deltah_sum )
         term2 = ( array_im_j*deltah_i ) /
     1           ( deltah_im*deltah_sum )
         term3 = ( array_i_j*( deltah_i - deltah_im ) )/
     1           ( deltah_i*deltah_im )
         dfdh = (term1 - term2) + term3
         write(luttyo,200) charh, dfdh, charh, df2dh2
200      format(' ','First partial  wrt ',a1,' = ',1pe12.3/
     1          ' ','Second partial wrt ',a1,' = ',1pe12.3)
      endif
c
c     do y derivatives
c
      if( can_do_v ) then
         if( iflgs .eq. 1 ) then
            array_i_jm  = invert(array(pt_i,pt_j-1))
            array_i_j   = invert(array(pt_i,pt_j))
            array_i_jp  = invert(array(pt_i,pt_j+1))
         else
            array_i_jm  = array(pt_i,pt_j-1)
            array_i_j   = array(pt_i,pt_j)
            array_i_jp  = array(pt_i,pt_j+1)
         endif
c        some constants...
         deltav_jm = xyz(pt_j,idistj) - xyz(pt_j - 1,idistj)
         deltav_j = xyz(pt_j+1,idistj) - xyz(pt_j,idistj)
         deltav_sum = deltav_jm + deltav_j
         denomv = 1.0/(deltav_j*deltav_jm*(deltav_sum))
c
c        do second derivative
         df2dv2 = (2.0*( array_i_jm*deltav_j -
     1                   array_i_j*(deltav_sum) +
     1                   array_i_jp*deltav_jm ) ) * denomv
c        do first derivative
         term1 = ( array_i_jp*deltav_jm ) /
     1           ( deltav_j*deltav_sum )
         term2 = ( array_i_jm*deltav_j ) /
     1           ( deltav_jm*deltav_sum )
         term3 = ( array_i_j*( deltav_j - deltav_jm ) )/
     1           ( deltav_j*deltav_jm )
         dfdv = (term1 - term2) + term3
         write(luttyo,210) charv, dfdv, charv, df2dv2
210      format(' ','First partial  wrt ',a1,' = ',1pe12.3/
     1          ' ','Second partial wrt ',a1,' = ',1pe12.3)
      endif
c
c     do mixed partial
      if( can_do_h .and. can_do_v ) then
         if( iflgs .eq. 1 ) then
            array_im_jm  = invert(array(pt_i-1,pt_j-1))
            array_im_jp  = invert(array(pt_i-1,pt_j+1))
            array_ip_jm  = invert(array(pt_i+1,pt_j-1))
            array_ip_jp  = invert(array(pt_i+1,pt_j+1))
         else
            array_im_jm  = array(pt_i-1,pt_j-1)
            array_im_jp  = array(pt_i-1,pt_j+1)
            array_ip_jm  = array(pt_i+1,pt_j-1)
            array_ip_jp  = array(pt_i+1,pt_j+1)
         endif
c
c     compute first x derivative at j-1 and j+1;
c     then use these derivative values to compute first derivative
c     in y (partial x, partial y). In theory, the order of
c     differentiation does not matter...
c
c        do first derivative at j+1
         term1 = ( array_ip_jp*deltah_im ) /
     1           ( deltah_i*deltah_sum )
         term2 = ( array_im_jp*deltah_i ) /
     1           ( deltah_im*deltah_sum )
         term3 = ( array_i_jp*( deltah_i - deltah_im ) )/
     1           ( deltah_i*deltah_im )
         dfdhjp = (term1 - term2) + term3
c
c        do first derivative at j-1
         term1 = ( array_ip_jm*deltah_im ) /
     1           ( deltah_i*deltah_sum )
         term2 = ( array_im_jm*deltah_i ) /
     1           ( deltah_im*deltah_sum )
         term3 = ( array_i_jm*( deltah_i - deltah_im ) )/
     1           ( deltah_i*deltah_im )
         dfdhjm = (term1 - term2) + term3
c
c        do first derivative in y direction
c
         term1 = ( dfdhjp*deltav_jm ) /
     1           ( deltav_j*deltav_sum )
         term2 = ( dfdhjm*deltav_j ) /
     1           ( deltav_jm*deltav_sum )
         term3 = ( dfdh*( deltav_j - deltav_jm ) )/
     1           ( deltav_j*deltav_jm )
         df2dhdv = (term1 - term2) + term3
         write(luttyo,220) df2dhdv
220      format(' ','Mixed second partial = ',1pe12.3)
      else
         write(luttyo,300)
300      format(' ','Point is on boundary;'
     1              ' can''t compute mixed partial')
      endif
c
      return      
      end
c
      real function invert(value)
      implicit none
c
c     invert value, with test for zero divide
c
      real value
c
      include 'sysdep.cmn'
      include 'io.cmn'
c
      if( value .eq. 0.0 ) then
         write(luttyo,100)
100      format(' ','(invert) warning: tried to divide by zero ',
     1          '...result set to PM_MAX_FLOAT')
         invert = PM_MAX_FLOAT  ! large single float number
      else
         invert = 1.0/value
      endif
c
      return
      end
c
      subroutine min_max_array_tri(data_array, array_size)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'units.cmn'
c
      integer array_size
      real data_array(array_size)
c
c     locals
c
      character*1 answer
      character*27 prompt
      character*10 format_array(2)
      logical do_zero
      integer i
      real min_h, max_h, min_v, max_v
      real min_value, max_value
      real low_h, hi_h, low_v, hi_v, low_temp, hi_temp, real_vec(2)
c
      min_value = PM_MAX_FLOAT
      max_value = -PM_MAX_FLOAT
c
c     Get search bounds
c
c     Horizontal limits
c
      prompt = 'Enter ? limits: '
      prompt(7:7) = charh
110   real_vec(1) = pm_scale_cm2plot*xyzlimits(1,idisti)
      real_vec(2) = pm_scale_cm2plot*xyzlimits(2,idisti)
      format_array(1) = '1pg12.5'
      format_array(2) = '1pg12.5'
      call get_reals_with_default(real_vec, 2, prompt, format_array, 0)
      low_temp = real_vec(1)*pm_scale_plot2cm
      hi_temp  = real_vec(2)*pm_scale_plot2cm
      if( hi_temp .le. low_temp ) then
         write(luttyo,111)
111      format(' ','Error: max must be larger than min')
         goto 110
      else if( low_temp .lt. xyzlimits(1, idisti) .or.
     1         hi_temp  .gt. xyzlimits(2, idisti) ) then
         write(luttyo,112)
112      format(' ','Error: Bounds are outside of the device')
         goto 110
      endif
      low_h = low_temp
      hi_h  = hi_temp
c
c     Vertical limits
c
      prompt(7:7) = charv
120   real_vec(1) = pm_scale_cm2plot*xyzlimits(1,idistj)
      real_vec(2) = pm_scale_cm2plot*xyzlimits(2,idistj)
      format_array(1) = '1pg12.5'
      format_array(2) = '1pg12.5'
      call get_reals_with_default(real_vec, 2, prompt, format_array, 0)
      low_temp = real_vec(1)*pm_scale_plot2cm
      hi_temp  = real_vec(2)*pm_scale_plot2cm
      if( hi_temp .le. low_temp ) then
         write(luttyo,121)
121      format(' ','error: max must be larger than min')
         goto 120
      else if( low_temp .lt. xyzlimits(1,idistj) .or.
     1         hi_temp  .gt. xyzlimits(2,idistj) ) then
         write(luttyo,122)
122      format(' ','Error: Bounds are outside of the device')
         goto 120
      endif
      hi_v = hi_temp
      low_v = low_temp
c
      call get_char(answer,
     1     'Exclude zero values in search for minimum? (Y/N) [N] ',0)
      call upcase(answer)
      if( answer .eq. 'Y' ) then
         do_zero = .FALSE.
      else
         do_zero = .TRUE.
      endif
c
c     find minimum/maximum value of array
c
      do 10 i=1,tri_npts
         if( min_value .gt. data_array(i) ) then
            if( data_array(i) .ne. 0.0 .or. do_zero ) then
               min_value = data_array(i)
               min_h = cord(1,i)
               min_v = cord(2,i)
               endif
            endif
         if( max_value .lt. data_array(i) ) then
            max_value = data_array(i)
            max_h = cord(1,i)
            max_v = cord(2,i)
         endif
10    continue
c
      write(luttyo,200) min_value, pm_scale_cm2plot*min_h,
     1                             pm_scale_cm2plot*min_v
200   format(' ','Minimum value ',1pe12.3,' found at (',1pg12.5,',',
     1       1pg12.5,')')
      write(luttyo,210) max_value, pm_scale_cm2plot*max_h,
     1                             pm_scale_cm2plot*max_v
210   format(' ','Maximum value ',1pe12.3,' found at (',1pg12.5,',',
     1       1pg12.5,')')
c
      call get_char(answer,'Hit return for more... ',0)
c
      return
      end
