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     Subroutines to calculate an analytical function for input into Postmini
c
      subroutine create_analytical
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'io.cmn'
      include 'analytical.cmn'
c
c     Calculate an data set (1D/2D)
c
c     locals
c
c     Menu code items
      character*80 string
      character*40 str
      real         real_vec(3)
      integer      i, item_number, int_vec(1),
     1             imax, imin, string_len, ilen
      logical      eof, do_eval, delayed_message, read_error, lerr,
     1             use_last_input
c
c     Intialize
      if( analytical_nx .eq. 0 ) then
         analytical_dim    = 1
         analytical_type   = ANALYTICAL_TYPE_X
         analytical_xstart = 0.0
         analytical_xend   = 1.0
         analytical_xstep  = 0.1
         analytical_pstart = 0.0
         analytical_pend   = 1.0
         analytical_pstep  = 0.1
         analytical_ystart = 0.0
         analytical_yend   = 1.0
         analytical_ystep  = 0.1
         analytical_nx     =
     1      INT((analytical_xend - analytical_xstart)/analytical_xstep
     1           + 0.5) + 1
         analytical_np     =
     1      INT((analytical_pend - analytical_pstart)/analytical_pstep
     1           + 0.5) + 1
         analytical_ny     =
     1      INT((analytical_yend - analytical_ystart)/analytical_ystep
     1           + 0.5) + 1
         analytical_expression = 'x'
         analytical_sweep = .FALSE.
      endif
c
      do_eval = .FALSE.
      use_last_input = .FALSE.
c
100   call clear_screen
      call display_string(1, 1,
     1     '  Analytical expression',41,0)
      call display_string(3 ,1 ,' 1) y=f(x), y=f(x,p), z=f(x,y) :',32,0)
      call display_string(4 ,1 ,' 2) X start, end, step values  :',32,0)
      call display_string(5 ,1 ,' 3) P start, end, step values  :',32,0)
      call display_string(6 ,1 ,' 4) Y start, end, step values  :',32,0)
      call display_string(7 ,1 ,' 5) Function :',14,0)
c
      if( analytical_dim .eq. 1 ) then
         if( analytical_type .eq. ANALYTICAL_TYPE_X ) then
            str = 'y=f(x)'
         else
            str = 'y=f(x,p)'
         endif
      else if( analytical_dim .eq. 2 ) then
         str = 'z=f(x,y)'
      endif
      call display_string (3,34,str,40,1)
c
      call display_real (4,34,analytical_xstart,1)
      call display_real (4,45,analytical_xend  ,1)
      call display_real (4,56,analytical_xstep ,1)
c
      if( analytical_dim .eq. 2 .or. (analytical_dim .eq. 1 .and.
     1    analytical_type .eq. ANALYTICAL_TYPE_X) ) then
         str = 'NA'
         call display_string ( 5,34,str,40,1)
      else
         call display_real ( 5,34,analytical_pstart,1)
         call display_real ( 5,45,analytical_pend  ,1)
         call display_real ( 5,56,analytical_pstep ,1)
      endif
c
      if( analytical_dim .eq. 1 ) then
         str = 'NA'
         call display_string (6,34,str,40,1)
      else
         call display_real ( 6,34,analytical_ystart,1)
         call display_real ( 6,45,analytical_yend  ,1)
         call display_real ( 6,56,analytical_ystep ,1)
      endif
c
      call display_string( 7,27,analytical_expression,52,1)
c
c     Display any messaged from before window refresh
      if( delayed_message ) then
         call display_message(string)
         delayed_message = .FALSE.
      endif
c
200   call display_prompt('Enter item number, ' //
     1   'or Repaint, Accept (evaluate function), Quit ')
      if( .NOT. use_last_input ) then
         call read_from_input_area(string, eof)
      endif
      if( eof ) goto 200
      call upcase(string)
      if( string(1:1) .eq. 'A') then
         do_eval = .TRUE.
         goto 300
      else if( string(1:1) .eq. 'Q') then
         do_eval = .FALSE.
         goto 300
      elseif (string(1:1) .eq. 'R' ) then
         goto 100
      else
         read(string,*,err=1000) item_number
         if ( item_number .eq. 1 ) then
c           Function type
            call get_ints_with_prompt(int_vec, 1,
     1         'Enter 1 for y=f(x), 2 for y=f(x,p), 3 for z=f(x,y)')
            if( int_vec(1) .lt. 1 .or. int_vec(1) .gt. 3 ) then
               call display_message('Invalid type, must be 1, 2 or 3')
            else
               analytical_type = int_vec(1)
               if( analytical_type .eq. ANALYTICAL_TYPE_X ) then
                  analytical_dim = 1
                  str = 'y=f(x)'
               else if( analytical_type .eq. ANALYTICAL_TYPE_XP ) then
                  analytical_dim = 1
                  str = 'y=f(x,p)'
               else
                  analytical_dim = 2
                  str = 'z=f(x,y)'
               endif
               call display_string (3,34,str,40,1)
c
               if( analytical_type .eq. ANALYTICAL_TYPE_X ) then
                  str = 'NA'
                  call display_string ( 5,34,str,40,1)
                  call display_string ( 6,34,str,40,1)
                  analytical_sweep = .FALSE.
               else if( analytical_type .eq. ANALYTICAL_TYPE_XP ) then
                  call display_real ( 5,34,analytical_pstart,1)
                  call display_real ( 5,45,analytical_pend  ,1)
                  call display_real ( 5,56,analytical_pstep ,1)
                  str = 'NA'
                  call display_string ( 6,34,str,40,1)
                  analytical_sweep = .TRUE.
               else
                  str = 'NA'
                  call display_string ( 5,34,str,40,1)
                  call display_real ( 6,34,analytical_ystart,1)
                  call display_real ( 6,45,analytical_yend  ,1)
                  call display_real ( 6,56,analytical_ystep ,1)
                  analytical_sweep = .FALSE.
               endif
            endif
         else if ( item_number .eq. 2 ) then
c           X start,end,step
            call get_reals_with_prompt(real_vec, 3,
     1         'Enter start, end, step values for x (e.g. 0 1 .1)')
            if( real_vec(2) .le. real_vec(1) ) then
               call display_message(
     1            'End value must be larger than start value')
            else if( real_vec(3) .gt. (real_vec(2) - real_vec(1))) then
               call display_message(
     1            'Step value must be <= (end-start)')
            else
               analytical_xstart = real_vec(1)
               analytical_xend   = real_vec(2)
               analytical_xstep  = real_vec(3)
               call display_real (4,34,analytical_xstart,1)
               call display_real (4,45,analytical_xend  ,1)
               call display_real (4,56,analytical_xstep ,1)
            endif
         else if ( item_number .eq. 3 ) then
c           P start, end, step
            if( analytical_type .ne. ANALYTICAL_TYPE_XP ) then
               call display_message(
     1              'Not applicable for this type function')
            else
               call get_reals_with_prompt(real_vec, 3,
     1            'Enter start, end, step values for p (e.g. 0 1 .1)')
               if( real_vec(2) .le. real_vec(1) ) then
                  call display_message(
     1               'End value must be larger than start value')
               else if(real_vec(3) .gt. (real_vec(2)-real_vec(1))) then
                  call display_message(
     1               'Step value must be <= (end-start)')
               else
                  analytical_pstart = real_vec(1)
                  analytical_pend   = real_vec(2)
                  analytical_pstep  = real_vec(3)
                  call display_real (5,34,analytical_pstart,1)
                  call display_real (5,45,analytical_pend  ,1)
                  call display_real (5,56,analytical_pstep ,1)
               endif
            endif
         else if ( item_number .eq. 4 ) then
c           Y start, end, step
            if( analytical_dim .eq. 1 ) then
               call display_message('Invalid for f(x) or f(x,p)')
            else
               call get_reals_with_prompt(real_vec, 3,
     1            'Enter start, end, step values for y (e.g. 0 1 .1)')
               if( real_vec(2) .le. real_vec(1) ) then
                  call display_message(
     1               'End value must be larger than start value')
               else if(real_vec(3) .gt. (real_vec(2)-real_vec(1))) then
                  call display_message(
     1               'Step value must be <= (end-start)')
               else
                  analytical_ystart = real_vec(1)
                  analytical_yend   = real_vec(2)
                  analytical_ystep  = real_vec(3)
                  call display_real (6,34,analytical_ystart,1)
                  call display_real (6,45,analytical_yend  ,1)
                  call display_real (6,56,analytical_ystep ,1)
               endif
            endif
         else if ( item_number .eq. 5 ) then
c           Expression
            if( analytical_type .eq. ANALYTICAL_TYPE_X ) then
               call get_string_with_prompt(analytical_expression,
     1         'Enter expression to generate f(x)'//
     2         '(e.g. 0.1 + x**2 + sin(x) )')
            else if( analytical_type .eq. ANALYTICAL_TYPE_XP ) then
               call get_string_with_prompt(analytical_expression,
     1         'Enter expression to generate f(x,p)'//
     2         '(e.g. 0.1 + p**2 + sin(x) )')
            else
               call get_string_with_prompt(analytical_expression,
     1         'Enter expression to generate f(x,y)'//
     2         '(e.g. 0.1 + y**2 + sin(x) )')
            endif
            ilen = string_len(analytical_expression)
            if( ilen .eq. 0 ) then
               call display_message('Must have an expression')
            else
               call display_string( 7,27,analytical_expression,52,1)
            endif
         else
            call display_message('Item number out of range')
         endif
      endif
      goto 200
c
1000  call display_message('Illegal input')
      goto 200
c
300   if( do_eval ) then
         if( analytical_dim .eq. 1 ) then
            call do_eval_analytical
         else
            call do_eval_analytical_2d
         endif
      else
         have_analytical_data = .FALSE.
      endif
c
      return
      end
c
      subroutine do_eval_analytical
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'mesh.cmn'
      include 'analytical.cmn'
      include 'gui.cmn'
      include 'evaluate.cmn'
c
c     evaluate a 1D dataset, using parameters in analytical.cmn
c
c     locals
c
      integer i, j, k, ioerr, ilen, itype,
     1        ierr, ind, nvars, npts, ioff
      integer string_len
c
      integer istat, pm_put_error_msg
      character*8 var_names(2)
c
c     GUI error level: 0=success, 1=warning, 2=error
      gui_error_level = 0
      call pm_clear_error_msg
c
      prev_binfiletype = binfiletype
      binfiletype = DT_ANALYTICAL
      binfilename = 'analytical'
      datafile_record = 1
c     get stem of file name (for hardcopy)
      call get_filestem(binfilename, binfilestem)
c     Initialize x/y name and units
      analytical_x_name  = ' '
      analytical_x_units = ' '
      analytical_y_name  = 'User defined'
      analytical_y_units = ' '
c
      if( .NOT. lgui ) then
         write(luttyo,10)
10       format(' ','Evaluating function...',$)
      endif
c
      if( analytical_array_size .eq. 0 ) then
         analytical_array_size = COMP_ARRAYSIZE
         analytical_xval_ptr = pm_malloc(4*analytical_array_size)
         analytical_yval_ptr = pm_malloc(4*analytical_array_size)
         analytical_zval_ptr = pm_malloc(4*analytical_array_size)
      endif
c
      call zero_real_array(analytical_xval,analytical_array_size,1)
      call zero_real_array(analytical_yval,analytical_array_size,1)
      call zero_real_array(analytical_zval,analytical_array_size,1)
c
      have_analytical_data = .FALSE.
c
c     Evaluate function
c
      if( analytical_type .eq. ANALYTICAL_TYPE_X ) then
c        Generate y=f(x) or y=f(x,p) data
         analytical_nx = INT((analytical_xend - analytical_xstart)/
     1                        analytical_xstep + 0.5) + 1
         npts = analytical_nx
         if( npts .gt. COMP_ARRAYSIZE ) goto 10040
c
         do i=1,analytical_nx
            analytical_xval(i) = analytical_xstart +
     1                           analytical_xstep*float(i-1)
            xyz(i,1) = analytical_xval(i)
         end do
         var_names(1) = 'x'
         call ev_expression_f(analytical_expression,var_names,1,
     1                        analytical_nx,
     2                        analytical_xval,COMP_ARRAYSIZE,
     3                        xyz(1,2),ierr,ind)
         if( ierr .ne. EVAL_SUCCESS ) goto 10050
      else
c
c        Generate y=f(x,p) data
         analytical_nx = INT((analytical_xend - analytical_xstart)/
     1                        analytical_xstep + 0.5) + 1
         analytical_np = INT((analytical_pend - analytical_pstart)/
     1                        analytical_pstep + 0.5) + 1
         npts = analytical_nx*analytical_np
         if( npts .gt. COMP_ARRAYSIZE ) goto 10040
c
         do j=1,analytical_np
            ioff = (j-1)*analytical_nx
            do i=1,analytical_nx
               analytical_xval(i+ioff) =
     1             analytical_xstart + analytical_xstep*float(i-1)
               xyz(i+ioff,1)      = analytical_xval(i+ioff)
               work_vec(i+ioff,1) = analytical_xval(i+ioff)
               work_vec(i+ioff,2) = analytical_pstart +
     1                              analytical_pstep*float(j-1)
            end do
         end do
         var_names(1) = 'x'
         var_names(2) = 'p'
         call ev_expression_f(analytical_expression,var_names,2,
     1                        npts,
     2                        work_vec,COMP_ARRAYSIZE,
     3                        xyz(1,2),ierr,ind)
         if( ierr .ne. EVAL_SUCCESS ) goto 10050
      endif
      do i=1,npts
         analytical_yval(i) = xyz(i,2)
      end do
c
      have_analytical_data = .TRUE.
c
c     Now, copy data to Postmini main storage area
c
c     Clear out arrays
      do 200 j=1,3
         do 200 i=1,COMP_ARRAYSIZE
            xyz(i,j) = 0.0
            edxyz(i,j) = 0.0
            dxyzm(i,j) = 0.0
200   continue
      do 201 j=1,MENUSIZE
         in_binary_file(j) = .FALSE.
201   continue
      call zero_real_array(arr2d,COMP_ARRAYSIZE,MENUSIZE)
c
      do 203 i=1,npts
c        Scaled data goes in slot 1
         xyz(i,1)   = analytical_xval(i)*MU2CM
         arr2d(i,1) = analytical_yval(i)
c        Unscaled data goes in slots 2 and 3 (so we change scaling
c        later without re-reading file)
         arr2d(i,2) = analytical_xval(i)*MU2CM
         arr2d(i,3) = analytical_yval(i)
203   continue
c
      nx = npts
c     Fake second and third dimension
      ny = 1
      xyz(1,2) = 0.0
c
      nz = 1
      nt = ny
      nyoxm = 0
      nytot= ny
      nxsp = 1
      nxsoxm = 1
      plane_index = 1
      prev_plane_type = plane_type
      plane_type = XY_PLANE
      xyz(1,3) = 0.0
      nxy=nx*ny
      nxytot=nx*nytot
c
c     This variable must be set by all read routines...
      data_is_trigrid = .FALSE.
c
      in_binary_file(1) = .TRUE.
c
c     Load generic mesh pointers, depending on plane type
      call set_3d_index_limits
c
c     Set up boundary limits information
      call coord_limits
c
      nxm = nx - 1
      nym = ny - 1
      nytotm = nytot - 1
c
c     inverse mesh spacing...
c
      do 204 i=1,nxm
         if( (xyz(i+1,1) - xyz(i,1)) .eq. 0.0 ) then
            edxyz(i,1) = PM_MAX_FLOAT
         else
            edxyz(i,1) = 1.0/(xyz(i+1,1) - xyz(i,1))
         endif
204   continue
      do 205 i=1,nytotm
205      edxyz(i,2) = 1.0/(xyz(i+1,2) - xyz(i,2))
      do 206 i=1,nz-1
206      edxyz(i,3) = 1.0/(xyz(i+1,3) - xyz(i,3))
c
c        average mesh spacing...
c
      dxyzm(1,1)  = (xyz(2,1)  - xyz(1,1))/2.0
      dxyzm(nx,1) = (xyz(nx,1) - xyz(nx-1,1))/2.0
      do 207 i=2,nxm
207     dxyzm(i,1) = (xyz(i+1,1) - xyz(i-1,1))/2.0
c
      dxyzm(1,2)     = (xyz(2,2)  - xyz(1,2))/2.0
      dxyzm(nytot,2) = (xyz(nytot,2) - xyz(nytot-1,2))/2.0
      do 208 i=2,nytotm
208      dxyzm(i,2) = (xyz(i+1,2) - xyz(i-1,2))/2.0
c
      dxyzm(1,3)     = (xyz(2,3)  - xyz(1,3))/2.0
      dxyzm(nz,3) = (xyz(nz,3) - xyz(nz-1,3))/2.0
      do 209 i=2,nz-1
209      dxyzm(i,3) = (xyz(i+1,3) - xyz(i-1,3))/2.0
c
      itype = 1
      call load_analytical_specific(itype)
c
      if( lfirst_read ) then
         lfirst_read = .FALSE.
         prev_plane_type = plane_type
         prev_binfilename = binfilename
         prev_binfiletype = binfiletype
      endif
c     Force plot to rescale each time we evaluate function
      oned_lfirst = .TRUE.
      twod_lfirst = .TRUE.
      threed_lfirst = .TRUE.
c
      datafile_loaded = .TRUE.
c
      goto 20000
c
c     errors of commission/omission
c
10040 write(msgbuf,10041) COMP_ARRAYSIZE
10041 format('Error: too many points. Maximum = ',i10)
      istat = pm_put_error_msg(msgbuf)
      have_analytical_data = .FALSE.
      gui_error_level = 2
      goto 20000
c
10050 write(msgbuf,10051)
10051 format('Error: while evaluating expression')
      istat = pm_put_error_msg(msgbuf)
      call ev_lookup_error_message(ierr,msgbuf)
      istat = pm_put_error_msg(msgbuf)
      if( ind .gt. 0 ) then
         istat = pm_put_error_msg(analytical_expression)
         msgbuf = ' '
         msgbuf(ind:ind) = '^'
         istat = pm_put_error_msg(msgbuf)
      endif
      have_analytical_data = .FALSE.
      gui_error_level = 2
      goto 20000
c
20000 continue
      return
      end
c
      subroutine load_analytical_specific(itype)
      implicit none
c
c     Load descriptive arrays with ASCII specific data
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'data.cmn'
      include 'analytical.cmn'
c
      integer itype   ! 1== 1-D, 2==2-D data
c
c     Data varies linear or logarithmic
      logical LOG, LIN
      parameter (LOG = .TRUE., LIN = .FALSE.)
c     Data has value in oxide, or not
      logical OX, NO_OX
      parameter (OX = .TRUE., NO_OX = .FALSE.)
c     Data jumps at oxide interfaces
      logical JUMP, NO_JUMP
      parameter (JUMP = .TRUE., NO_JUMP = .FALSE.)
c     Data is 1D or 2D
      logical ONED, TWOD
      parameter (ONED = .TRUE., TWOD = .FALSE.)
c
c     Number of menu items
      max_menu_num = 1
c     Array index of net doping
      net_dop_array_index = 0
c
      if( itype .eq. 1 ) then
         call load_desc_arr(1, analytical_y_name, 'ANALYTICAL',
     1                      analytical_y_units,' ',' ','.dat',
     2                      LIN,NO_OX,NO_JUMP,ONED,1)
      else
         call load_desc_arr(1, 'User defined', 'ANALYTICAL',
     1                      ' ',' ',' ','.dat',
     2                      LIN,NO_OX,NO_JUMP,TWOD,1)
      endif
c
      return
      end
c
      subroutine do_eval_analytical_2d
      implicit none
c
      include 'const.h'
      include 'sysdep.cmn'
      include 'plot.cmn'
      include 'data.cmn'
      include 'io.cmn'
      include 'mesh.cmn'
      include 'analytical.cmn'
      include 'evaluate.cmn'
      include 'gui.cmn'
c
c     evaluate an analytical 2D expression, using parameters in analytical.cmn
c
c     locals
c
      real, allocatable :: tmp(:,:)
c
      integer i, j, ncols, ioerr, ilen, itype, ierr, istat, ind
      integer string_len
      character*8 var_names(2)
      integer pm_put_error_msg
c
      have_analytical_data = .FALSE.
      prev_binfiletype = binfiletype
      binfiletype = DT_ANALYTICAL
      binfilename = 'analytical'
      datafile_record = 1
c     Initialize x/y name and units
      analytical_x_name  = ' '
      analytical_x_units = ' '
      analytical_y_name  = 'User defined'
      analytical_y_units = ' '
c
      write(luttyo,10)
10    format(' ','Evaluating function ...',$)
c
c     clear out arrays
      do 45 j=1,3
         do 45 i=1,COMP_ARRAYSIZE
            xyz(i,j) = 0.0
            edxyz(i,j) = 0.0
            dxyzm(i,j) = 0.0
45    continue
      do 51 j=1,MENUSIZE
         in_binary_file(j) = .FALSE.
51       continue
c
      if( analytical_array_size .eq. 0 ) then
         analytical_array_size = COMP_ARRAYSIZE
         analytical_xval_ptr = pm_malloc(4*analytical_array_size)
         analytical_yval_ptr = pm_malloc(4*analytical_array_size)
         analytical_zval_ptr = pm_malloc(4*analytical_array_size)
      endif
c
      call zero_real_array(analytical_xval,analytical_array_size,1)
      call zero_real_array(analytical_yval,analytical_array_size,1)
      call zero_real_array(analytical_zval,analytical_array_size,1)
c
      have_analytical_data = .FALSE.
c
c     Generate x data
      analytical_nx = INT((analytical_xend - analytical_xstart)/
     1                     analytical_xstep + 0.5) + 1
      do i=1,analytical_nx
         analytical_xval(i) = analytical_xstart +
     1                        analytical_xstep*float(i-1)
         xyz(i,1) = MU2CM*analytical_xval(i)
      end do
c
c     Generate y data
      analytical_ny = INT((analytical_yend - analytical_ystart)/
     1                     analytical_ystep + 0.5) + 1
      do i=1,analytical_ny
         analytical_yval(i) = analytical_ystart +
     1                        analytical_ystep*float(i-1)
         xyz(i,2) = MU2CM*analytical_yval(i)
      end do
c
      if( analytical_nx*analytical_ny .gt. ARR2DSIZE ) goto 10030
c
c     Generate z=f(x,y) data
      var_names(1) = 'x'
      var_names(2) = 'y'
c
c     We allocate a (large) temporary to hold the x,y pairs
c     We could avoid this at the cost of calling the expression evaluator
c     in a doubly nested loop...
      allocate(tmp(analytical_nx*analytical_ny,2),stat=istat)
      if( istat .ne. 0 ) goto 10040
      do j=1,analytical_ny
         do i=1,analytical_nx
            tmp(i+(j-1)*analytical_nx,1) = analytical_xval(i)
            tmp(i+(j-1)*analytical_nx,2) = analytical_yval(j)
         end do
      end do
c
      call ev_expression_f(analytical_expression,var_names,2,
     1                     analytical_nx*analytical_ny,
     2                     tmp,analytical_nx*analytical_ny,
     3                     work,ierr,ind)
      if( ierr .ne. EVAL_SUCCESS ) goto 10050
c
      deallocate(tmp,stat=istat)
c
1500  continue
      write(luttyo,150)
150   format(' ','done')
c
      have_analytical_data = .TRUE.
c
c     Now, copy data to Postmini main storage area
c
      nx = analytical_nx
      ny = analytical_ny
      call allocate_rect_arrays(nx,ny,1,ierr)
         if( ierr .ne. 0 ) then
c           Failure during array allocation
            write(luttyo,10080)
10080       format(' ','Internal error: could not allocate data array')
            have_analytical_data = .FALSE.
            return
         endif
c        Zero array (not really necessary)
         call zero_real_array(arrdat,nx*ny,1)
c
c        Copy data to the arrdat array
         do i=1,nx*ny
            arrdat(i) = work(i)
         end do
c
         nz = 1
         nt = ny
         nyoxm = 0
         nytot= ny
         nxsp = 1
         nxsoxm = 1
         plane_index = 1
         prev_plane_type = plane_type
         plane_type = XY_PLANE
         xyz(1,3) = 0.0
         nxy=nx*ny
         nxytot=nx*nytot
c
c        This variable must be set by all read routines...
         data_is_trigrid = .FALSE.
c
         in_binary_file(1) = .TRUE.
c
c        Load generic mesh pointers, depending on plane type
         call set_3d_index_limits
c
c        Set up boundary limits information
         call coord_limits
c
         nxm = nx - 1
         nym = ny - 1
         nytotm = nytot - 1
c
c        inverse mesh spacing...
c
         do 204 i=1,nxm
            if( (xyz(i+1,1) - xyz(i,1)) .eq. 0.0 ) then
               edxyz(i,1) = PM_MAX_FLOAT
            else
               edxyz(i,1) = 1.0/(xyz(i+1,1) - xyz(i,1))
            endif
204      continue
         do 205 i=1,nytotm
205         edxyz(i,2) = 1.0/(xyz(i+1,2) - xyz(i,2))
         do 206 i=1,nz-1
206         edxyz(i,3) = 1.0/(xyz(i+1,3) - xyz(i,3))
c
c        average mesh spacing...
c
         dxyzm(1,1)  = (xyz(2,1)  - xyz(1,1))/2.0
         dxyzm(nx,1) = (xyz(nx,1) - xyz(nx-1,1))/2.0
         do 207 i=2,nxm
207        dxyzm(i,1) = (xyz(i+1,1) - xyz(i-1,1))/2.0
c
         dxyzm(1,2)     = (xyz(2,2)  - xyz(1,2))/2.0
         dxyzm(nytot,2) = (xyz(nytot,2) - xyz(nytot-1,2))/2.0
         do 208 i=2,nytotm
208         dxyzm(i,2) = (xyz(i+1,2) - xyz(i-1,2))/2.0
c
         dxyzm(1,3)     = (xyz(2,3)  - xyz(1,3))/2.0
         dxyzm(nz,3) = (xyz(nz,3) - xyz(nz-1,3))/2.0
         do 209 i=2,nz-1
209         dxyzm(i,3) = (xyz(i+1,3) - xyz(i-1,3))/2.0
c
         itype = 2
         call load_analytical_specific(itype)
c
         if( lfirst_read ) then
            lfirst_read = .FALSE.
            prev_plane_type = plane_type
            prev_binfilename = binfilename
            prev_binfiletype = binfiletype
         endif
c        Force plot to rescale each time we evaluate the function
         oned_lfirst = .TRUE.
         twod_lfirst = .TRUE.
         threed_lfirst = .TRUE.
c
         datafile_loaded = .TRUE.
c
      goto 20000
c
c     errors of commission/omission
c
10030 write(luttyo,10031) ARR2DSIZE
10031 format(' Error: analytical function has too many points, max = ',
     1       i8)
      have_analytical_data = .FALSE.
      gui_error_level = 2
      goto 20000
c
10040 write(luttyo,10041) 
10041 format(' Error: could not allocate temporary storage ',
     1       'during function evaluation')
      have_analytical_data = .FALSE.
      gui_error_level = 2
      goto 20000
c
10050 write(luttyo,10051) 
10051 format('Error: while evaluating expression')
      istat = pm_put_error_msg(msgbuf)
      call ev_lookup_error_message(ierr,msgbuf)
      istat = pm_put_error_msg(msgbuf)
      if( ind .gt. 0 ) then
         istat = pm_put_error_msg(analytical_expression)
         msgbuf = ' '
         msgbuf(ind:ind) = '^'
         istat = pm_put_error_msg(msgbuf)
      endif
      have_analytical_data = .FALSE.
      gui_error_level = 2
      goto 20000
c
20000 continue
      return
      end

