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     Integrate a region of device
c
      subroutine integrate
      implicit none
c
c     Asks user what quantity to look at, determine integration region,
c     then compute 1 or 2D integral.
c
c     Programmer: John Faricelli
c     Date: 30-May-86
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'user.cmn'
      include 'bincodes.cmn'
c
c     locals
c
      logical luser, error
      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
c        Load Grapes data, if needed...
         if( binfiletype .eq. DT_GRAPES ) then
            call pm_select_grapes_data(error)
            if( error ) return
         endif
c        Load SPICE3 data, if needed...
         if( binfiletype .eq. DT_SPICE3 ) then
            call pm_select_spice3_data(error)
            if( error ) return
         endif
         if( data_is_trigrid ) then
            call int_region_tri(arrdat(arrind(menu_to_arr(1))),
     1                   tri_npts,
     2                   does_data_jump_at_oxide(1),
     3                   is_data_in_oxide(1),
     4                   data_units_int_2d(1),
     5                   data_units_int_1d(1))
         else if( is_data_1d(1) ) then
            luser = data_nickname(1) .eq. 'ASCII'
            call integrate_1d(arr2d(1,menu_to_arr(1)),
     1                   COMP_ARRAYSIZE,
     2                   does_data_jump_at_oxide(1),
     3                   is_data_in_oxide(1),
     4                   data_units_int_1d(1),luser)
         else
            call int_region(arrdat(arrind(menu_to_arr(1))),
     1                   ni, nj,
     2                   does_data_jump_at_oxide(1),
     3                   is_data_in_oxide(1),
     3                   data_units_int_2d(1),
     4                   data_units_int_1d(1))
         endif
         return
      endif
c
200   call clear_screen
      write(luttyo,errfmt)
      errfmt = '('' '')'
      call write_data_menu('Integrate region:')
      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 = '('' '',''Item is not present in data file'')'
         call beep_terminal
      else if( data_is_trigrid ) then
            call int_region_tri(arrdat(arrind(menu_to_arr(answer))),
     1                   tri_npts,
     2                   does_data_jump_at_oxide(answer),
     3                   is_data_in_oxide(answer),
     4                   data_units_int_2d(answer),
     5                   data_units_int_1d(answer))
      else if( is_data_1d(answer) ) then
         luser = data_nickname(answer) .eq. 'ASCII'
         call integrate_1d(arr2d(1,menu_to_arr(answer)), COMP_ARRAYSIZE,
     1                   does_data_jump_at_oxide(answer),
     2                   is_data_in_oxide(answer),
     3                   data_units_int_1d(answer),luser)
      else
         call int_region(arrdat(arrind(menu_to_arr(answer))),
     1                   ni, nj,
     2                   does_data_jump_at_oxide(answer),
     3                   is_data_in_oxide(answer),
     3                   data_units_int_2d(answer),
     4                   data_units_int_1d(answer))
      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 int_region(data_array, max_h, max_v,
     1                      discontinuity, value_in_oxide,
     2                      units_2dint, units_1dint)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'units.cmn'
c
      integer max_h, max_v
      logical discontinuity, value_in_oxide
      real data_array(max_h, max_v)
      character*(*) units_1dint, units_2dint
c
c     locals
c
      character*1 answer
      character*6 third_dir
      character*1  type
      integer hl_index, hh_index, vl_index, vh_index, h_index, v_index
      integer i, ulen, string_len, ilen
      real int_h_min, int_h_max, int_v_min, int_v_max, h_coord, v_coord
      real drain_contact_length
      real value
      logical do_1d, do_2d, do_v_int, do_h_int, novalue_in_oxide
c
      ulen = string_len(pm_length_unit)
      novalue_in_oxide = .NOT. value_in_oxide
c
      call clear_screen
c
c     Print out coordinates of landmarks in the device.
c     It would also be nice to have the position of the source and drain
c     junctions, but this is really hard to guess from the doping
c     profile in all cases -- such as a depletion or counterdoped channel,
c     where there strictly is no junction. Any suggestions?
c
c     Find edges of gate. Mesh always has points on exact edges of gate.
      if( .NOT. three_d ) then
         do 1 i=1,nx
            if( xyz(i,1) .eq. 0.0 ) nxg1 = i
            if( xyz(i,1) .eq. xl  ) nxg2 = i
1        continue
      endif
c
      if( binfiletype .eq. DT_MINIMOS .or.
     1    binfiletype .eq. DT_PROMIS ) then
         write(luttyo,5)
5        format(' ','Caution: Integration uses approximate area ',
     1              'weights'/
     1          ' ','for non-planar simulations and may be inaccurate'/)
      endif
c
      if( binfiletype .eq. DT_MINIMOS ) then
      if( plane_type .eq. XY_PLANE .or. plane_type .eq. XZ_PLANE) then
         drain_contact_length  = xyz(nx,1) - xyz(nxsoxm+1,1)
c
         write(luttyo,10) xyz(1,1)*pm_scale_cm2plot,
     1                    xyz(nxsp-1,1)*pm_scale_cm2plot
10       format(' ','Source contact extends from  X = ',1pg12.5,
     1         ' to ',1pg12.5)
         write(luttyo,11) xyz(nxsp,1)*pm_scale_cm2plot,
     1                    xyz(nxg1-1,1)*pm_scale_cm2plot
11       format(' ','Source oxide extends from    X = ',1pg12.5,
     1         ' to ',1pg12.5/)
c
         write(luttyo,12) xyz(nxg1,1)*pm_scale_cm2plot,
     1                    xyz(nxg2,1)*pm_scale_cm2plot
12       format(' ','Gate metal extends from      X = ',1pg12.5,
     1         ' to ',1pg12.5/)
c
         write(luttyo,13) xyz(nxg2+1,1)*pm_scale_cm2plot,
     1                    xyz(nxsoxm,1)*pm_scale_cm2plot
13       format(' ','Drain oxide extends from     X = ',1pg12.5,
     1         ' to ',1pg12.5)
         write(luttyo,14) xyz(nxsoxm+1,1)*pm_scale_cm2plot,
     1                    xyz(nx,1)*pm_scale_cm2plot,
     1                    drain_contact_length*pm_scale_cm2plot,
     1                    pm_length_unit(1:ulen)
14       format(' ','Drain contact extends from   X = ',1pg12.5,
     1         ' to ',1pg12.5,' (length: ',1pg12.5,' ',a,')')
      endif
c
      if( plane_type .eq. YZ_PLANE .or. plane_type .eq. XZ_PLANE) then
         write(luttyo,16) xyz(1,3)*pm_scale_cm2plot,
     1                    xyz(nz,3)*pm_scale_cm2plot
16       format(' ','In the width direction:'/
     1          ' ','Device starts (thin oxide) at Z = ',1pg12.5/
     2          ' ','Edge of gate            is at Z =   0.000'/
     3          ' ','Device ends (field oxide)  at Z = ',1pg12.5)
      endif
      endif
c
c     Ask user for region to integrate
c
      do_1d = .FALSE.
      do_2d = .FALSE.
c
20    call get_char(type,
     1     'Perform 1D (line) or 2D (area) integration (1 or 2)? ',0)
      if( type .eq. '1' ) then
         do_1d = .TRUE.
      else if (type .eq. '2' ) then
         do_2d = .TRUE.
      else
         write(luttyo,25)
25       format(' ','Illegal entry: integration type must be 1 or 2')
         goto 20
      endif
c
      if( do_2d ) then
         call get_int_limits(int_h_min, int_h_max, hl_index, hh_index,
     1        'H', novalue_in_oxide)
         call get_int_limits(int_v_min, int_v_max, vl_index, vh_index,
     1        'V', novalue_in_oxide)
         call do_2d_integrate(data_array, max_h, max_v,
     1        int_h_min, int_h_max, int_v_min, int_v_max,
     2        hl_index, hh_index, vl_index, vh_index,
     3        novalue_in_oxide, value)
      endif
c
      if( do_1d ) then
         do_h_int = .FALSE.
         do_v_int = .FALSE.
30       call get_char(type,
     1      'Integrate along horizontal or vertical line (H or V)? ',0)
         call upcase(type)
         if( type .eq. 'H' ) then
            do_h_int = .TRUE.
         else if (type .eq. 'V' ) then
            do_v_int = .TRUE.
         else
            write(luttyo,40)
40          format(' ','Illegal entry: integration type must be H or V')
            goto 30
         endif
c
         if( do_v_int ) then
            call get_cutline_coord(h_coord, h_index, 'H',
     1           novalue_in_oxide)
            call get_int_limits(int_v_min, int_v_max, vl_index,
     1           vh_index, 'V', novalue_in_oxide)
            write(luttyo,100) charh,
     1                        xyz(h_index,idisti)*pm_scale_cm2plot
100         format(/' ','Integration done along mesh line at ',a1,' = ',
     1           f12.4)
            write(luttyo,101) charv,
     1                        xyz(vl_index,idistj)*pm_scale_cm2plot,
     2                        xyz(vh_index,idistj)*pm_scale_cm2plot
101         format(' ','with limits from mesh lines at ',a1,' = ',
     1           f12.4,' to ',f12.4/)
            call do_1d_integrate(data_array, max_h, max_v,
     1           int_v_min, int_v_max, vl_index, vh_index, h_index,
     2           novalue_in_oxide, 'V', value)
         else
            call get_cutline_coord(v_coord ,v_index, 'V',
     1           novalue_in_oxide)
            call get_int_limits(int_h_min, int_h_max, hl_index,
     1            hh_index, 'H', novalue_in_oxide)
            write(luttyo,200) charv,
     1                        xyz(v_index,idistj)*pm_scale_cm2plot
200         format(/' ','Integration done along mesh line at ',a1,' = ',
     1           f12.4)
            write(luttyo,201) charh,
     1                        xyz(hl_index,idisti)*pm_scale_cm2plot,
     1                        xyz(hh_index,idisti)*pm_scale_cm2plot
201         format(' ','with limits from mesh lines at ',a1,' = ',
     1           f12.4,' to ',f12.4/)
            call do_1d_integrate(data_array, max_h, max_v,
     1           int_h_min, int_h_max, hl_index, hh_index, v_index,
     2           novalue_in_oxide, 'H', value)
         endif
      endif
c
      if( plane_type .eq. XY_PLANE ) then
         third_dir = 'width '
      else if( plane_type .eq. YZ_PLANE ) then
         third_dir = 'length'
      else
         third_dir = 'depth '
      endif
c
      if( do_2d ) then
         ilen = string_len(units_2dint)
         write(luttyo,500) value, units_2dint(1:ilen), third_dir
500      format(' ','Integrated value is: ',1pe15.6,1x,a,
     1          ' (scaled to 1 cm in ',a,')'/)
      else
         ilen = string_len(units_1dint)
         write(luttyo,501) value, units_1dint(1:ilen), third_dir
501      format(' ','Integrated value is: ',1pe15.6,1x,a,
     1          ' (scaled to 1 cm in ',a,')'/)
      endif
c
      call get_char(answer,'Hit return to continue...',0)
c
      return
      end
c
      subroutine do_2d_integrate(data_array, max_h, max_v,
     1     int_h_min, int_h_max, int_v_min, int_v_max,
     2     index_hl, index_hh, index_vl, index_vh,
     3     novalue_in_oxide, value)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'io.cmn'
      include 'units.cmn'
c
c     Finally does the actual integration
c
      logical novalue_in_oxide
      integer index_hl, index_hh, index_vl, index_vh
      integer max_h, max_v
      real data_array(max_h, max_v)
      real int_h_min, int_h_max, int_v_min, int_v_max
      real value
c
c     locals
      integer i,j
      real dy
c
      write(luttyo,100) charh, xyz(index_hl,idisti)*pm_scale_cm2plot,
     1                  xyz(index_hh,idisti)*pm_scale_cm2plot
100   format(/' ','Actual integration from ',a1,' = ',1pg12.5,' to ',
     1             1pg12.5)
      write(luttyo,200) charv, xyz(index_vl,idistj)*pm_scale_cm2plot,
     1                  xyz(index_vh,idistj)*pm_scale_cm2plot
200   format( ' ','               and from ',a1,' = ',1pg12.5,' to ',
     1             1pg12.5/)
c
c     Double sum
c
c     Note well: this code has been generalized for 3D, but
c     does not apply the correct area weights for nonplanar regions.
c
      value = 0.0
      do 1000 j=index_vl, index_vh
         dy = dxyzm(j,idistj)
c        Use correct area weight at oxide/semiconductor
c        interface. This is most critical when integrating
c        currents.
c        We do not take non-planar boundaries into account yet.
         if( plane_type .ne. XZ_PLANE .and.
     1       novalue_in_oxide .and. xyz(j,idistj) .eq. 0.0 ) then
            dy = (xyz(j+1,idistj) - xyz(j,idistj))/2.0
         endif
         do 1000 i=index_hl, index_hh
            value = value + dxyzm(i,idisti)*dy*data_array(i,j)
1000  continue
      return
      end
c
      subroutine do_1d_integrate(data_array, max_h, max_v,
     1     lower, upper, ilower, iupper, iline,
     2     novalue_in_oxide, type, value)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'mesh.cmn'
c
c     Finally does the actual integration
c
      character*1 type
      logical novalue_in_oxide
      integer ilower, iupper, iline
      integer max_h, max_v
      real data_array(max_h, max_v)
      real upper, lower
      real value
c
c     locals
      integer i,j
      real dy
c
      value = 0.0
      if( type .eq. 'V' ) then
         do 1000 j=ilower, iupper
c           Correct integration weight at integration limits
c           Does not take nonplanar boundary into account.
            if( j .eq. ilower ) then
               dy = (xyz(j+1,idistj) - xyz(j,idistj))/2.0
            else if( j .eq. iupper ) then
               dy = (xyz(j,idistj) - xyz(j-1,idistj))/2.0
            else
               dy = dxyzm(j,idistj)
            endif
            value = value + dy*data_array(iline,j)
1000     continue
      else if( type .eq. 'H' ) then
         do 2000 i=ilower, iupper
c           Correct integration weight at integration limits
c           Does not take nonplanar boundary into account.
            if( i .eq. ilower ) then
               dy = (xyz(i+1,idisti) - xyz(i,idisti))/2.0
            else if( i .eq. iupper ) then
               dy = (xyz(i,idisti) - xyz(i-1,idisti))/2.0
            else
               dy = dxyzm(i,idisti)
            endif
            value = value + dy*data_array(i,iline)
2000     continue
      else
         write(luttyo,3000)
3000     format(' ','Internal error: Illegal 1D integration type')
      endif
c
      return
      end

c
      subroutine integrate_1d(data_array, array_size,
     1     discontinuity, value_in_oxide, units_1dint, luser)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'units.cmn'
      include 'mesh.cmn'
c
c     Integration for pure 1D quantities
c
c     For now, discontinuity information is ignored
c
      logical discontinuity, value_in_oxide, luser
      integer array_size
      real data_array(array_size)
      character*(*) units_1dint
c
c     locals
      logical novalue_in_oxide
      integer j, ilen, string_len
      real dy, value, interp, interp2, lower_bound, upper_bound,
     1     lower, upper, real_vec(2)
      character*20 format(2)
      character*1 answer
c
      novalue_in_oxide = .NOT. value_in_oxide
c
      write(luttyo,100) pm_scale_cm2plot*xyz(1,idisti),
     1                  pm_scale_cm2plot*xyz(ni,idisti)
100   format(/' ','Minimum  = ',1pg12.5,'; Maximum  = ',1pg12.5)
c
200   continue
      lower_bound = xyz(1 ,idisti)
      upper_bound = xyz(ni,idisti)
      real_vec(1) = pm_scale_cm2plot*lower_bound
      real_vec(2) = pm_scale_cm2plot*upper_bound
      format(1) = '1pe10.3'
      format(2) = '1pe10.3'
      call get_reals_with_default(real_vec, 2,
     1     'Enter min/max integration bounds: ', format, 0)
      lower = real_vec(1)*pm_scale_plot2cm
      if( lower .lt. lower_bound ) then
         lower = lower_bound
         write(luttyo,203)
203      format(' ',
     1   'Warning: lower bound less than minimum: was set to minimum')
      endif
      upper = real_vec(2)*pm_scale_plot2cm
      if( upper .gt. upper_bound ) then
         upper = upper_bound
         write(luttyo,204)
204      format(' ',
     1  'Warning: upper bound greater than maximum: was set to maximum')
      endif
c
      if( lower .gt. upper ) then
         write(luttyo,201)
201      format(' ','Error: minimum must be less than maximum')
         goto 200
      endif
c
      value = 0.0
      do 1000 j=1, ni-1
c
c        Scan coordinates, performing integration only when we are
c        within the integration limits.
c
c        Check for duplicate points
         if( xyz(j,idisti) .eq. xyz(j+1,idisti) ) goto 1000
c
         if( xyz(j,idisti) .le. lower .and.
     1       upper .le. xyz(j+1,idisti) ) then
c           Integration limits reside within a single interval
            interp = data_array(j) +
     1               (data_array(j+1) - data_array(j))*
     1               (lower - xyz(j,idisti))/
     1               (xyz(j+1,idisti) - xyz(j,idisti))
            interp2 = data_array(j) +
     1                (data_array(j+1) - data_array(j))*
     1                (upper - xyz(j,idisti))/
     1                (xyz(j+1,idisti) - xyz(j,idisti))
            dy = 0.5*(upper - lower)
            value = value + dy*(interp+interp2)
         else if( xyz(j,idisti)   .lt. lower .and.
     1            xyz(j+1,idisti) .ge. lower ) then
c           Integrate partial interval...
            interp = data_array(j) +
     1               (data_array(j+1) - data_array(j))*
     1               (lower - xyz(j,idisti))/
     1               (xyz(j+1,idisti) - xyz(j,idisti))
            dy = 0.5*(xyz(j+1,idisti) - lower)
            value = value + dy*(interp+data_array(j+1))
         else if( xyz(j,idisti)   .lt. upper .and.
     1            xyz(j+1,idisti) .ge. upper ) then
c           Integrate partial interval...
            interp = data_array(j) +
     1               (data_array(j+1) - data_array(j))*
     1               (upper - xyz(j,idisti))/
     1               (xyz(j+1,idisti) - xyz(j,idisti))
            dy = 0.5*(upper - xyz(j,idisti))
            value = value + dy*(interp+data_array(j))
         else if( lower .le. xyz(j,idisti)  .and.
     1            xyz(j+1,idisti) .le. upper ) then
c           Integrate full interval
            dy = 0.5*(xyz(j+1,idisti) - xyz(j,idisti))
            value = value + dy*(data_array(j)+data_array(j+1))
         endif
1000  continue
c
      ilen = string_len(units_1dint)
      if( luser ) then
         write(luttyo,101) value*pm_scale_cm2plot, units_1dint(1:ilen)
      else
         write(luttyo,101) value, units_1dint(1:ilen)
      endif
101   format(' ','Integrated value is: ',1pe15.6,1x,a)
c
      call get_char(answer,'Hit return to continue...',0)
c
      return
      end
c
      subroutine get_int_limits(lower, upper, ilower, iupper, 
     1           type, novalue_in_oxide)
      implicit none
c
c     Get integration limits in horizontal or vertical in device, depending
c     on 'type' ('H' or 'V'). Lower and upper are the lower and upper
c     integration limits (real). Ilower and iupper are the mesh
c     line indicies whose position is closest to lower and upper.
c        Also check to make sure that limits are inside
c     device bounds, upper is greater than lower, and round
c     iupper and ilower to the closest mesh point.
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'io.cmn'
      include 'units.cmn'
c
      real lower, upper
      integer ilower, iupper
      character*1 type
      logical novalue_in_oxide
c
c     locals
      real real_vec(2), lower_bound, upper_bound
      character*20 format(2)
      integer getindex
c
      if( type .eq. 'H' ) then
         write(luttyo,100) charh, pm_scale_cm2plot*xyz(1,idisti),
     1                     charh, pm_scale_cm2plot*xyz(ni,idisti)
100      format(/' ','Minimum ',a1,' = ',1pg12.5,'; Maximum ',a1,' = ',
     1           1pg12.5)
      else if( type .eq. 'V' ) then
         write(luttyo,100) charv, pm_scale_cm2plot*xyz(1,idistj),
     1                     charv, pm_scale_cm2plot*xyz(nj,idistj)
      else
         write(luttyo,900)
900      format(' ','Internal error: Illegal int limits type')
         return
      endif
c
200   continue
      if( type .eq. 'H' ) then
         lower_bound = xyz(1 ,idisti)
         upper_bound = xyz(ni,idisti)
         real_vec(1) = pm_scale_cm2plot*lower_bound
         real_vec(2) = pm_scale_cm2plot*upper_bound
      else
         lower_bound = xyz(1 ,idistj)
         upper_bound = xyz(nj,idistj)
         real_vec(1) = pm_scale_cm2plot*lower_bound
         real_vec(2) = pm_scale_cm2plot*upper_bound
      endif
      format(1) = '1pe10.3'
      format(2) = '1pe10.3'
      call get_reals_with_default(real_vec, 2,
     1     'Enter min/max integration bounds: ', format, 0)
      lower = real_vec(1)*pm_scale_plot2cm
      if( lower .lt. lower_bound ) then
         lower = lower_bound
         write(luttyo,203)
203      format(' ',
     1   'Warning: lower bound less than minimum: was set to minimum')
      endif
      upper = real_vec(2)*pm_scale_plot2cm
      if( upper .gt. upper_bound ) then
         upper = upper_bound
         write(luttyo,204)
204      format(' ',
     1  'Warning: upper bound greater than maximum: was set to maximum')
      endif
c
      if( lower .gt. upper ) then
         write(luttyo,201)
201      format(' ','Error: minimum must be less than maximum')
         goto 200
      endif
c
      if( type .eq. 'V' .and. plane_type .ne. XZ_PLANE .and.
     1    novalue_in_oxide .and. 
     1   (lower .lt. 0.0 .or. upper .lt. 0.0) ) then
         write(luttyo,202)
202      format(' ','Error: Integration in oxide (y < 0) not defined ',
     1              'for this variable')
         goto 200
      endif
c
c     Get equivalent index
      if( type .eq. 'H' ) then
         ilower = getindex(lower, xyz(1,idisti), ni)
      else
         ilower = getindex(lower, xyz(1,idistj), nj)
      endif
      if( ilower .eq. 0 ) then
         write(luttyo,1000) lower*pm_scale_cm2plot
1000     format(' ','Value ',1pe10.3,' not in bounds')
         goto 200
      endif
      if( type .eq. 'H' ) then
         iupper = getindex(upper, xyz(1,idisti), ni)
      else
         iupper = getindex(upper, xyz(1,idistj), nj)
      endif
      if( iupper .eq. 0 ) then
         write(luttyo,1000) upper*pm_scale_cm2plot
         goto 200
      endif
c
c     Adjust index to get closest to specified value.
c     (getindex() function always gets mesh index of left (lower)
c     side of interval.
      if( type .eq. 'H' ) then
c        Round lower limit
         if( (xyz(ilower+1,idisti) - lower) .lt.
     1       (lower - xyz(ilower,idisti)) ) then
            ilower = ilower + 1
         endif
c        Round upper limit
         if( iupper .ne. ni ) then
            if( (xyz(iupper+1,idisti) - upper) .lt.
     1          (upper - xyz(iupper,idisti)) ) then
               iupper = iupper + 1
            endif
         endif
      else
c        Round lower limit
         if( (xyz(ilower+1,idistj) - lower) .lt.
     1       (lower - xyz(ilower,idistj)) ) then
            ilower = ilower + 1
         endif
c
c        Round upper limit
         if( iupper .ne. nj ) then
            if( (xyz(iupper+1,idistj) - upper) .lt.
     1          (upper - xyz(iupper,idistj)) ) then
               iupper = iupper + 1
            endif
         endif
      endif
c
      return
      end
c
      subroutine get_cutline_coord(line, iline, type, novalue_in_oxide)
      implicit none
c
c     Get a valid horizontal or vertical coordinate in device along which a
c     1d integration may be performed. Line is the coordinate
c     value (real). Iline is the mesh line index whose position
c     is closest to line.
c        Also check to make sure that line is inside
c     device bounds, and round line to the closest mesh point.
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'mesh.cmn'
      include 'units.cmn'
c
      real line
      integer iline
      character*1 type
      logical novalue_in_oxide
c
c     locals
      real real_vec(2)
      integer getindex
c
      if( type .eq. 'H' ) then
         write(luttyo,100) charh, pm_scale_cm2plot*xyz(1,idisti),
     1                     charh, pm_scale_cm2plot*xyz(ni,idisti)
100      format(/' ','Minimum ',a1,' = ',1pg12.5,'; Maximum ',a1,' = ',
     1                1pg12.5)
      else if( type .eq. 'V' ) then
         write(luttyo,100) charv, pm_scale_cm2plot*xyz(1,idistj),
     1                     charv, pm_scale_cm2plot*xyz(nj,idistj)
      else
         write(luttyo,900)
900      format(' ','Internal error: Illegal int limits type')
         return
      endif
c
200   continue
      call get_reals(real_vec,1,
     1  'Enter integration line coordinate in microns: ', 0)
      line = real_vec(1)*pm_scale_plot2cm
      if( type .eq. 'V' .and. plane_type .ne. XZ_PLANE .and.
     1    novalue_in_oxide .and. line .lt. 0.0 ) then
         write(luttyo,202)
202      format(' ','Error: Integration in oxide (y < 0) not defined ',
     1              'for this variable')
         goto 200
      endif
c
c     Get equivalent index
      if( type .eq. 'H' ) then
         iline = getindex(line, xyz(1,idisti), ni)
      else
         iline = getindex(line, xyz(1,idistj), nj)
      endif
      if( iline .eq. 0 ) then
         write(luttyo,1000) line
1000     format(' ','Value ',1pe10.3,' not in device bounds')
         goto 200
      endif
c
c     Adjust index to get closest to specified value.
c     (getindex() function always gets mesh index of left (lower)
c     side of interval.
      if( type .eq. 'H' ) then
         if( iline .ne. ni ) then
            if( (xyz(iline+1,idisti) - line) .lt.
     1          (line - xyz(iline,idisti)) ) then
               iline = iline + 1
            endif
         endif
      else
         if( iline .ne. nj ) then
            if( (xyz(iline+1,idistj) - line) .lt.
     1          (line - xyz(iline,idistj)) ) then
               iline = iline + 1
            endif
         endif
      endif
c
      return
      end
c
      subroutine int_region_tri(data_array, array_size,
     1                      discontinuity, value_in_oxide,
     2                      units_2dint, units_1dint)
      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
      logical discontinuity, value_in_oxide
      real data_array(array_size)
      character*(*) units_1dint, units_2dint
c
c     locals
c
      character*1 answer
      character*6 third_dir
      character*1  type
      integer ulen, string_len, ilen
      real h_min, h_max, v_min, v_max, h_coord, v_coord
      real value
      logical do_1d, do_2d, do_v_int, do_h_int, novalue_in_oxide
c
      ulen = string_len(pm_length_unit)
      novalue_in_oxide = .NOT. value_in_oxide
c
      call clear_screen
c
c     Print out coordinates of landmarks in the device.
c     It would also be nice to have the position of the source and drain
c     junctions, but this is really hard to guess from the doping
c     profile in all cases -- such as a depletion or counterdoped channel,
c     where there strictly is no junction. Any suggestions?
c
c
      write(luttyo,5)
5     format(' ','Caution: Integration takes entire area assigned ',
     1           'to each node'/
     1       ' ','inside the integration bounds'/)
c
      write(luttyo,10) xyzlimits(1,1)*pm_scale_cm2plot,
     1                 xyzlimits(2,1)*pm_scale_cm2plot
10       format(' ','Structure extends in X from = ',1pg12.5,
     1         ' to ',1pg12.5)
         write(luttyo,11) xyzlimits(1,2)*pm_scale_cm2plot,
     1                    xyzlimits(2,2)*pm_scale_cm2plot
11       format(' ','Structure extends in Y from = ',1pg12.5,
     1         ' to ',1pg12.5/)
c
c     Ask user for region to integrate
c
      do_1d = .FALSE.
      do_2d = .FALSE.
c
20    call get_char(type,
     1     'Perform 1D (line) or 2D (area) integration (1 or 2)? ',0)
      if( type .eq. '1' ) then
         do_1d = .TRUE.
      else if (type .eq. '2' ) then
         do_2d = .TRUE.
      else
         write(luttyo,25)
25       format(' ','Illegal entry: integration type must be 1 or 2')
         goto 20
      endif
c
      if( do_2d ) then
         call get_int_limits_tri(h_min, h_max, 'X', novalue_in_oxide)
         call get_int_limits_tri(v_min, v_max, 'Y', novalue_in_oxide)
         call do_2d_integrate_tri(data_array, array_size, h_min, h_max,
     1        v_min, v_max,
     3        novalue_in_oxide, value)
      endif
c
      if( do_1d ) then
         do_h_int = .FALSE.
         do_v_int = .FALSE.
30       call get_char(type,
     1      'Integrate along horizontal or vertical line (H or V)? ',0)
         call upcase(type)
         if( type .eq. 'H' ) then
            do_h_int = .TRUE.
         else if (type .eq. 'V' ) then
            do_v_int = .TRUE.
         else
            write(luttyo,40)
40          format(' ','Illegal entry: integration type must be H or V')
            goto 30
         endif
c
         if( do_v_int ) then
            call get_cutline_coord_tri(h_coord, charh, novalue_in_oxide)
            call get_int_limits_tri(v_min, v_max,
     1                              charv, novalue_in_oxide)
            write(luttyo,100) charh,
     1                        h_coord*pm_scale_cm2plot
100         format(/' ','Integration done at ',a1,' = ',
     1           f12.4)
            write(luttyo,101) charv,
     1                        v_min*pm_scale_cm2plot,
     2                        v_max*pm_scale_cm2plot
101         format(' ','with limits from ',a1,' = ',
     1           f12.4,' to ',f12.4/)
            call do_1d_integrate_tri(data_array, array_size,
     1           h_coord, v_min, v_max,
     2           novalue_in_oxide, charh, value)
         else
            call get_cutline_coord_tri(v_coord ,charv, novalue_in_oxide)
            call get_int_limits_tri(h_min, h_max,
     1                              charh, novalue_in_oxide)
            write(luttyo,200) charv,
     1                        v_coord*pm_scale_cm2plot
200         format(/' ','Integration done at ',a1,' = ',
     1           f12.4)
            write(luttyo,201) charh,
     1                        h_min*pm_scale_cm2plot,
     1                        h_max*pm_scale_cm2plot
201         format(' ','with limits from ',a1,' = ',
     1           f12.4,' to ',f12.4/)
            call do_1d_integrate_tri(data_array, array_size,
     1           v_coord, h_min, h_max,
     2           novalue_in_oxide, charv, value)
         endif
      endif
c
      if( plane_type .eq. XY_PLANE ) then
         third_dir = 'width '
      else if( plane_type .eq. YZ_PLANE ) then
         third_dir = 'length'
      else
         third_dir = 'depth '
      endif
c
      if( do_2d ) then
         ilen = string_len(units_2dint)
         write(luttyo,500) value, units_2dint(1:ilen), third_dir
500      format(' ','Integrated value is: ',1pe15.6,1x,a,
     1          ' (scaled to 1 cm in ',a,')'/)
      else
         ilen = string_len(units_1dint)
         write(luttyo,501) value, units_1dint(1:ilen), third_dir
501      format(' ','Integrated value is: ',1pe15.6,1x,a,
     1          ' (scaled to 1 cm in ',a,')'/)
      endif
c
      call get_char(answer,'Hit return to continue...',0)
c
      return
      end
c
      subroutine get_int_limits_tri(dmin, dmax, direction,
     1                              novalue_in_oxide)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'units.cmn'
      include 'mesh.cmn'
c
      real dmin, dmax
      character*1 direction
      logical novalue_in_oxide
c
      real lower, upper, real_vec(2)
c
      if( direction .eq. charh ) then
         lower = xyzlimits(1,1)
         upper = xyzlimits(2,1)
      else
         lower = xyzlimits(1,2)
         upper = xyzlimits(2,2)
      endif
c
      write(luttyo,100) direction,
     1                  lower*pm_scale_cm2plot, upper*pm_scale_cm2plot
100   format(' ','Limits in ',a,' are: ',1pe12.5,' to ',1pe12.5)
101   call get_reals(real_vec,2,'Enter lower and upper limits: ',0)
      if( real_vec(1) .ge. real_vec(2) ) then
         write(luttyo,102)
102      format(' ','Error: invalid limits')
         goto 101
      endif
      dmin = real_vec(1)*pm_scale_plot2cm
      if( dmin .lt. lower ) dmin = lower
      dmax = real_vec(2)*pm_scale_plot2cm
      if( dmax .gt. upper ) dmax = upper
c
      return
      end
c
      subroutine get_cutline_coord_tri(coord, direction,
     1                                 novalue_in_oxide)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'units.cmn'
      include 'io.cmn'
      include 'mesh.cmn'
c
      real coord
      character*1 direction
      logical novalue_in_oxide
c
c
      real lower, upper, real_vec(1)
c
      if( direction .eq. charh ) then
         lower = xyzlimits(1,1)
         upper = xyzlimits(2,1)
      else
         lower = xyzlimits(1,2)
         upper = xyzlimits(2,2)
      endif
c
      write(luttyo,100) lower*pm_scale_cm2plot, upper*pm_scale_cm2plot
100   format(' ','Limits are: ',1pe12.5,' to ',1pe12.5)
101   call get_reals(real_vec,1,'Enter cut line coordinate: ',0)
      real_vec(1) = real_vec(1)*pm_scale_plot2cm
      if( real_vec(1) .lt. lower .or. real_vec(1) .gt. upper ) then
         write(luttyo,102)
102      format(' ','Error: cut line outside structure')
         goto 101
      endif
      coord = real_vec(1)
c
      return
      end
c
      subroutine do_1d_integrate_tri(data_array, array_size,
     1           cut_coord, dmin, dmax,
     2           novalue_in_oxide, direction, value)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'mesh.cmn'
c
      integer array_size
      real data_array(array_size)
      real cut_coord, dmin, dmax, value
      logical novalue_in_oxide
      character*1 direction
c
      integer TMPSIZE
      parameter (TMPSIZE = 4*COMP_ARRAYSIZE)
      real tmp(TMPSIZE,2)
      integer i, np
      logical discontinuity, value_in_oxide, along_oxide, data_is_log,
     1        lerr
c
      discontinuity = .FALSE.
      data_is_log   = .FALSE.
      along_oxide   = .FALSE.
      value_in_oxide = .NOT. novalue_in_oxide
c
      call extract_cut_line_tri(data_array, array_size,
     1                          discontinuity,
     2                          data_is_log, value_in_oxide,
     3                          cut_coord, direction,
     4                          along_oxide,
     5                          tmp, TMPSIZE, 2, np, lerr)
c
c     tmp(*,1) has values
c     tmp(*,2) has coordinates
c
      if( lerr ) then
         write(luttyo,100)
100      format(' ','Error during extraction of cut line')
         value = 0.0
      else
         value = 0.0
         do 200 i=1,np
            if( tmp(i,2) .ge. dmin .and. tmp(i,2) .le. dmax ) then
            if( i .eq. 1 ) then
              value = value +
     1           0.5*tmp(1,1)*(tmp(2,2) - tmp(1,2))
            else if( i .eq. np ) then
              value = value +
     1           0.5*tmp(np,1)*(tmp(np,2) - tmp(np-1,2))
            else
              value = value +
     1           0.5*tmp(i,1)*(tmp(i+1,2) - tmp(i-1,2))
            endif
            endif
200      continue
      endif
c
      return
      end
c
      subroutine do_2d_integrate_tri(data_array, array_size,
     1        h_min, h_max,
     2        v_min, v_max,
     3        novalue_in_oxide, value)
      implicit none
c
      include 'const.h'
      include 'mesh.cmn'
c
      integer array_size
      real data_array(array_size)
      real h_min, h_max, v_min, v_max, value
      logical novalue_in_oxide
c
      integer i, j, ipt
      real x1, y1, x2, y2, x3, y3, xmin, xmax, ymin, ymax, area(3)
c
      value = 0.0
c
      do 1000 i=1,tri_ne
c        Skip if triangle is oxide and data has no value there
         if ( novalue_in_oxide .and. imat(tri2reg(i)) .lt. 0 ) goto 1000
c
c        Skip if triangle outside integration region
         ipt = nop(1,i)
         x1 = cord(1,ipt)
         y1 = cord(2,ipt)
         ipt = nop(2,i)
         x2 = cord(1,ipt)
         y2 = cord(2,ipt)
         ipt = nop(3,i)
         x3 = cord(1,ipt)
         y3 = cord(2,ipt)
         xmin = MIN(x1,x2,x3)
         xmax = MAX(x1,x2,x3)
         ymin = MIN(y1,y2,y3)
         ymax = MAX(y1,y2,y3)
         if( xmax .lt. h_min .or. h_max .lt. xmin .or.
     1       ymax .lt. v_min .or. v_max .lt. ymin ) goto 1000
c
c        Caution: get_tri_area may alter nop array
c                 that's why we re-read coordinates from nop
         call get_tri_area(i,area)
         do 1010 j=1,3
            ipt = nop(j,i)
            x1 = cord(1,ipt)
            y1 = cord(2,ipt)
c           If point is inside bounding box, add its portion of
c           triangle area. THIS IS A VERY CRUDE integration.
c           Note: in regions where data varies very fast,
c           (e.g. inversion layer, or at segregating boundaries) the
c           integral may be overestimated if the grid spacing is large.
            if( h_min .le. x1 .and. x1 .le. h_max .and.
     1          v_min .le. y1 .and. y1 .le. v_max ) then
               value = value + area(j)*data_array(ipt)
            endif
1010     continue
c
1000  continue
c
      return
      end
c
      subroutine get_tri_area(itri,area)
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'mesh.cmn'
      include 'io.cmn'
c
      integer itri
      real    area(3)
c
c     Patterned after PISCES GEOM subroutine
c
      logical swpflg
      integer i,j,k,swpcnt
      double precision xij,yij,xjk,yjk,xki,yki,disq,djsq,dksq,
     +     di,dj,dk,s,s4,r,rsq,den,tani,tanj,tank,hi,hj,hk,ai,aj,ak,
     +     si,sj,sk,xi,yi,xj,yj,xk,yk
c
      area(1) = 0.0
      area(2) = 0.0
      area(3) = 0.0
c
c     Set the swap flag and counter
      swpcnt = 1
      swpflg =.true.
c
10    i=nop(1,itri)
      j=nop(2,itri)
      k=nop(3,itri)
c
      xi=cord(1,i)
      yi=cord(2,i)
      xj=cord(1,j)
      yj=cord(2,j)
      xk=cord(1,k)
      yk=cord(2,k)
c
      xij=xi-xj
      yij=yi-yj
      xjk=xj-xk
      yjk=yj-yk
      xki=xk-xi
      yki=yk-yi
c
c     Compute coupling coefficients (cotangents of opposite angles)
c     and get denominator term
      den=xki*yjk-yki*xjk
c
c     Triangle area is -1/2 of this value
      s=-.5d0*den
      if (s.le.0.d0) then
         call plt$error_message('get_tri_area',
     1                          'swapping triangle nodes')
        if (swpcnt.gt.6) then
           call plt$error_message('get_tri_area',
     1     'Error in obtuse triangle area calculation')
           return
        endif

        swpcnt = swpcnt+1
        if (swpflg) then
c
c      Swap second and third nodes
          nop(2, itri) = k
          nop(3, itri) = j
        else
c
c      Swap first and second nodes
          nop(1, itri) = j
          nop(2, itri) = i
        endif
        swpflg = .not.swpflg
        goto 10
      endif
c
c     We want 1/2 of cotangent so invert and scale
      den=.5d0/den
c
c     compute area weighting
c
c***    2. calculate the distances between the nodes
c
      disq=xjk*xjk+yjk*yjk
      djsq=xki*xki+yki*yki
      dksq=xij*xij+yij*yij
      di=dsqrt(disq)
      dj=dsqrt(djsq)
      dk=dsqrt(dksq)
      if (di.eq.0.or.dj.eq.0.or.dk.eq.0) then
         call plt$error_message('get_tri_area',
     1                          'triangle has zero length side')
         return
      endif
c
c***    3. calculate the area, radius for outer circle and tangents
c
      s4=s*4.d0
      r=di*dj*dk/s4
      rsq=r*r
      den=djsq+dksq-disq
      if(den.eq.0.0d0) goto 400
      tani=s4/den
      den=dksq+disq-djsq
      if(den.eq.0.0d0) goto 400
      tanj=s4/den
      den=disq+djsq-dksq
      if(den.eq.0.0d0) goto 400
      tank=s4/den
c
c***    4. check for obtuse angle
c
      if(tani.lt.0.0d0.or.tanj.lt.0.0d0.or.tank.lt.0.0d0) goto 500
c
c***    5. acute triangle case
c
400   continue
c
      ai=0.25d0*di*dsqrt(dabs(rsq-0.25d0*disq))
      aj=0.25d0*dj*dsqrt(dabs(rsq-0.25d0*djsq))
      ak=0.25d0*dk*dsqrt(dabs(rsq-0.25d0*dksq))
      si=aj+ak
      sj=ak+ai
      sk=ai+aj
      goto 901
c
c***    6. obtuse triangle case
c
500   continue
      if(tani.ge.0.0d0) goto 600
      hi=0.
      hj=0.5d0*dj*tank
      hk=0.5d0*dk*tanj
      sj=0.25d0*dj*hj
      sk=0.25d0*dk*hk
      si=s-(sj+sk)
      goto 900
c
600   if(tanj.ge.0.0d0) goto 700
      hi=0.5d0*di*tank
      hj=0.
      hk=0.5d0*dk*tani
      sk=0.25d0*dk*hk
      si=0.25d0*di*hi
      sj=s-(sk+si)
      goto 900
c
700   if(tank.ge.0.0d0) then
           call plt$error_message('get_tri_area',
     1     'Error in obtuse triangle area calculation')
         return
      endif
      hi=0.5d0*di*tanj
      hj=0.5d0*dj*tani
      hk=0.
      si=0.25d0*di*hi
      sj=0.25d0*dj*hj
      sk=s-(si+sj)
c
900   continue
c
c     Store area contributions
901   area(1)=si
      area(2)=sj
      area(3)=sk
c
      return
      end
