!
!  this is an example of the *mapread* program (adjust if necessary)
!
!  *mapread* calls 2 routines (load and fields) from *maplib*
!
!  subroutine *load* is called once to read the fields from the map file
!  subroutine *fields* returns fields at a point for each call
!


c=======================================================================
      program mapread
c-----
      implicit none
c
      integer*4 mxs
      parameter (mxs=3000000)
      real*4 field(mxs), x, y, z, compx, compy, compz
      integer*2 ierr
      character*80 dummy, messg
c
      real*4 x0, y0, z0, dx, dy, dz, xmin, ymin, zmin, xmax, ymax,
     +       zmax, angle
      integer*4 nx, ny, nz, ida, idb, idc, ima, imb, imc
      common /size/ x0, y0, z0, dx, dy, dz, xmin, ymin, zmin, xmax,
     +              ymax, zmax, angle, nx, ny, nz, ida, idb, idc, ima,
     +              imb, imc
c-----
c
      write(6, '(x,a,$)') 'MAP filename [<RET>=exit]: '
      read(5, '(a)', end=999) dummy
      if(dummy .eq. ' ') goto 999
      call load(dummy, field, mxs, messg, ierr)  ! open & read MAP file.
      if(ierr .ne. 0) then
        write(6, *) messg
        stop
      endif
c
c display the map header ?
c
      write(6, *)
      if(angle .ne. 0.0) then
        write(6, '(5x,a,f7.1,a,i1,a)')
     +  'Rotation angle      -> ', angle*90.0/atan2(1.0,0.0),
     +  ' (rot.', nint(atan2(1.0,0.0)/abs(angle))-1, 'x)'
      else
        write(6, '(5x,a,f7.1,a)')
     +  'Rotation angle      -> ', 0.0, ' (rot.all)'
      endif
      write(6, '(5X,a,3(i7,3x))')
     +'Reflection indices  -> ', ida, idb, idc 
      write(6, *)
      write(6, '(5x,a,3f10.3)') 'Reference point  -> ',x0, y0, z0
      write(6, *)
      write(6, '(5x,a,3f10.3)') 'Starting point   -> ',xmin, ymin, zmin
      write(6, '(5x,a,3f10.3)') 'End point        -> ',xmax, ymax, zmax
      write(6, '(5x,a,3f10.3)') 'Step size        -> ',dx, dy, dz
c
      write(6, *)
100   write(6, '(x,a,$)') '(x,y,z)_cm point coordinates [<RET>=exit]: '
      read(5, '(a)', end=999) dummy
      if(dummy .eq. ' ') goto 999
      read(dummy, *, err=100) x, y, z

      call fields(field, mxs, x, y, z, compx, compy, compz, messg, ierr)
      if(ierr .ne. 0) then
        write(6, *) messg
      else
        if(ima .eq. 1) then
          write(6, *) 'Bx,By,Bz [Gauss]: ', compx, compy, compz
        elseif(ima.eq.4) then
          write(6, *) 'Ex,Ey,Ez [V/cm]: ', compx, compy, compz
        else
          write(6, *) 'Unknown field components'
        endif
      endif
      goto 100
c
999   call exit
      end
