c-----------------------------------------------------------------------
!!!!!!!!!!!!!!!!!! PLEASE DO NOT CHANGE THESE ROUTINES !!!!!!!!!!!!!!!!!
c-----------------------------------------------------------------------
Contact: Vjeran Vrankovic
       ! PSI Villigen, Switzerland
       ! WMHA-C30, tel. +41 (0)56 310-3591
       ! vjeran.vrankovic@psi.ch



c=======================================================================
      subroutine load(filename, field, mxs, messg, ierr)
c-----
      implicit none

      integer*4 mxs
      real*4 field(mxs), pi
      integer*4 isa, isb, isc, np, ix, iy, iz, ib
      integer*2 ierr
      character*80 filename, title, messg
      character*18 dates
      character*8 axa, axb, axc
      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-----
      open(21, file=filename, status='old', form='unformatted', err=901)
c
      read(21, err=902) title                                 ! comment.
      read(21, err=903) dates                           ! creation date.
      read(21, err=904) axa, axb, axc, angle            ! axes sequence,
                             ! rotation angle in the 1st (X-Y) quadrant,
                         ! (90DEG-dipol 45-quad 30-sext ... 0-solenoid).
      pi = atan2(0.0, -1.0)
      angle = angle / 180.0 * pi
      read(21, err=905) x0, y0, z0, isa, isb, isc, ima, imb, imc,
     +                  ida, idb, idc
                             ! starting grid point, coordinate sequence,
                   ! component sequence (1=Bx 2=By 3=Bz 4=Ex 5=Ey 6=Ez),
                       ! plane reflections (0=NONE 1=MIRROR -1=INVERSE).
      read(21, err=906) dx, dy, dz, nx, ny, nz             ! increments,
                                                     ! number of points.
      np = nx * ny * nz
      if(3*np .gt. mxs) then
        goto 907
      else
        do iz = 1, nz
          read(21, err=908) (((field(3*(nx*(ny*(iz-1)+iy-1)+ix-1)+ib),
     +    ib=1,3), ix=1,nx), iy=1,ny)                          ! fields.
        enddo
      endif
c
      xmax = x0 + (nx - 1) * dx
      if(abs(angle-0.5*pi) .lt. 1.0e-5) then
        ymax = y0 + (ny - 1) * dy
      else
        ymax = y0 + (nx - 1) * dx
      endif
      zmax = z0 + (nz - 1) * dz
c
      if(ida .eq. 0) then
        xmin = x0
      else
        xmin = x0 - (nx - 1) * dx
      endif
      if(idb .eq. 0) then
        ymin = y0
      else
        if(abs(angle-0.5*pi) .lt. 1.0e-5) then
          ymin = y0 - (ny - 1) * dy
        else
          ymin = y0 - (nx - 1) * dx
        endif
      endif
      if(idc .eq. 0) then
        zmin = z0
      else
        zmin = z0 - (nz - 1) * dz
      endif
c
      ierr = 0
999   close(21)
      return
c
901   ierr = 1
      messg = 'Cannot open the MAP file'
      goto 999
902   ierr = 2
      messg = 'Reading error: TITLE'
      goto 999
903   ierr = 3
      messg = 'Reading error: DATE'
      goto 999
904   ierr = 4
      messg = 'Reading error: AX~A/B/C, ANGLE'
      goto 999
905   ierr = 5
      messg = 'Reading error: X/Y/Z~MIN, IS~A/B/C, IM~A/B/C, ID~A/B/C'
      goto 999
906   ierr = 6
      messg = 'Reading error: D~X/Y/Z, N~X/Y/Z'
      goto 999
907   ierr = 7
      write(messg,*) 'Change the program parameter MXS to', 3*np
      goto 999
908   ierr = 8
      messg = 'Reading error: FIELDs'
      goto 999
c
      end



c=======================================================================
      subroutine fields(field,mxs,xp,yp,zp,cx,cy,cz,messg,ierr)
c-----
      implicit none

      integer*4 mxs
      real*4 field(mxs), xp, yp, zp, cx, cy, cz, xa, ya, atem, atemang,
     +       phi, cosphi, sinphi, rad, x, y, z, hx0, hx1, hy0, hy1, hz0,
     +       hz1, br, bp, cosatem, sinatem,
     +       h000, h001, h010, h011, h100, h101, h110, h111,
     +       cx000, cx001, cx010, cx011, cx100, cx101, cx110, cx111,
     +       cy000, cy001, cy010, cy011, cy100, cy101, cy110, cy111,
     +       cz000, cz001, cz010, cz011, cz100, cz101, cz110, cz111
      integer*4 i, j, k, ijk000, ijk001, ijk010, ijk011, ijk100, ijk101,
     +          ijk110, ijk111
      integer*2 ierr, nsek, nnsek
      character*80 messg
      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-----
      if(zp .lt. zmin .or. zp .gt. zmax) then        ! outside Z-limits.
        goto 903
      elseif(xp .lt. xmin) then                    ! outside Xmin-limit.
        goto 901
      elseif(yp .lt. ymin) then                    ! outside Ymin-limit.
        goto 902
      else
c
        xa = abs(xp - x0)
        ya = abs(yp - y0)
        if(xa .ne. 0.0 .or. ya .ne. 0.0) then
          atem = atan2(ya, xa)
        else
          atem = 0.0
        endif
        if(atem .gt. abs(angle)) then                        ! rotation.
          if(angle .ne. 0.0) then
            atemang = 0.5 * atem / abs(angle)
            nsek = int(atemang)
            nnsek = nint(atemang)
            phi = abs(atem - 2.0 * nnsek * abs(angle))
            if(phi .lt. 1.0e-06) phi = 0.0
          else
            nsek = 0
            nnsek = 0
            phi = 0.0
          endif
          cosphi = cos(phi)
          sinphi = sin(phi)
          rad = sqrt(xa * xa + ya * ya)
          xa = rad * cosphi
          ya = rad * sinphi
        endif
        x = x0 + xa
        y = y0 + ya
        z = z0 + abs(zp - z0)
        if(x .gt. xmax) then                         ! outside X-limits.
          goto 901
        elseif(y .gt. ymax) then                     ! outside Y-limits.
          goto 902
        endif
c
        i = int((x - x0) / dx) + 1
        j = int((y - y0) / dy) + 1
        k = int((z - z0) / dz) + 1
        ijk000 = 3 * (nx * (ny * (k - 1) + j - 1) + i - 1)
        ijk001 = 3 * (nx * (ny * (k - 1 + 1) + j - 1) + i - 1)
        ijk010 = 3 * (nx * (ny * (k - 1) + j - 1 + 1) + i - 1)
        ijk011 = 3 * (nx * (ny * (k - 1 + 1) + j - 1 + 1) + i - 1)
        ijk100 = 3 * (nx * (ny * (k - 1) + j - 1) + i - 1 + 1)
        ijk101 = 3 * (nx * (ny * (k - 1 + 1) + j - 1) + i - 1 + 1)
        ijk110 = 3 * (nx * (ny * (k - 1) + j - 1 + 1) + i - 1 + 1)
        ijk111 = 3 * (nx * (ny * (k - 1 + 1) + j - 1 + 1) + i - 1 + 1)
        cx000 = field(ijk000 + 1)
        cx001 = field(ijk001 + 1)
        cx010 = field(ijk010 + 1)
        cx011 = field(ijk011 + 1)
        cx100 = field(ijk100 + 1)
        cx101 = field(ijk101 + 1)
        cx110 = field(ijk110 + 1)
        cx111 = field(ijk111 + 1)
        cy000 = field(ijk000 + 2)
        cy001 = field(ijk001 + 2)
        cy010 = field(ijk010 + 2)
        cy011 = field(ijk011 + 2)
        cy100 = field(ijk100 + 2)
        cy101 = field(ijk101 + 2)
        cy110 = field(ijk110 + 2)
        cy111 = field(ijk111 + 2)
        cz000 = field(ijk000 + 3)
        cz001 = field(ijk001 + 3)
        cz010 = field(ijk010 + 3)
        cz011 = field(ijk011 + 3)
        cz100 = field(ijk100 + 3)
        cz101 = field(ijk101 + 3)
        cz110 = field(ijk110 + 3)
        cz111 = field(ijk111 + 3)
c
        hx1 = (x - x0 - (i - 1) * dx) / dx       ! linear interpolation.
        hy1 = (y - y0 - (j - 1) * dy) / dy
        hz1 = (z - z0 - (k - 1) * dz) / dz
        hx0 = 1.0 - hx1
        hy0 = 1.0 - hy1
        hz0 = 1.0 - hz1
        h000 = hx0 * hy0 * hz0
        if(abs(h000) .gt. 0.0 .and. cx000 .gt. 9.0e05 .and.
     +  cy000 .gt. 9.0e05 .and. cz000 .gt. 9.0e05) goto 904
        h001 = hx0 * hy0 * hz1
        if(abs(h001) .gt. 0.0 .and. cx001 .gt. 9.0e05 .and.
     +  cy001 .gt. 9.0e05 .and. cz001 .gt. 9.0e05) goto 904
        h010 = hx0 * hy1 * hz0
        if(abs(h010) .gt. 0.0 .and. cx010 .gt. 9.0e05 .and.
     +  cy010 .gt. 9.0e05 .and. cz010 .gt. 9.0e05) goto 904
        h011 = hx0 * hy1 * hz1
        if(abs(h011) .gt. 0.0 .and. cx011 .gt. 9.0e05 .and.
     +  cy011 .gt. 9.0e05 .and. cz011 .gt. 9.0e05) goto 904
        h100 = hx1 * hy0 * hz0
        if(abs(h100) .gt. 0.0 .and. cx100 .gt. 9.0e05 .and.
     +  cy100 .gt. 9.0e05 .and. cz100 .gt. 9.0e05) goto 904
        h101 = hx1 * hy0 * hz1
        if(abs(h101) .gt. 0.0 .and. cx101 .gt. 9.0e05 .and.
     +  cy101 .gt. 9.0e05 .and. cz101 .gt. 9.0e05) goto 904
        h110 = hx1 * hy1 * hz0
        if(abs(h110) .gt. 0.0 .and. cx110 .gt. 9.0e05 .and.
     +  cy110 .gt. 9.0e05 .and. cz110 .gt. 9.0e05) goto 904
        h111 = hx1 * hy1 * hz1
        if(abs(h111) .gt. 0.0 .and. cx111 .gt. 9.0e05 .and.
     +  cy111 .gt. 9.0e05 .and. cz111 .gt. 9.0e05) goto 904
        cx = cx000 * h000 + cx001 * h001 + cx010 * h010 + cx011 * h011 +
     +       cx100 * h100 + cx101 * h101 + cx110 * h110 + cx111 * h111
        cy = cy000 * h000 + cy001 * h001 + cy010 * h010 + cy011 * h011 +
     +       cy100 * h100 + cy101 * h101 + cy110 * h110 + cy111 * h111
        cz = cz000 * h000 + cz001 * h001 + cz010 * h010 + cz011 * h011 +
     +       cz100 * h100 + cz101 * h101 + cz110 * h110 + cz111 * h111
c
        if(atem.gt.abs(angle)) then                      ! backrotation.
          br=(cx*cosphi+cy*sinphi)*(-1)**nsek
          bp=(cy*cosphi-cx*sinphi)*(-1)**nnsek
          cosatem=cos(atem)
          sinatem=sin(atem)
          cx=br*cosatem-bp*sinatem
          cy=bp*cosatem+br*sinatem
        endif
c
        if(xp .lt. x0) then                             ! YZ-reflection.
          if(ida .eq. 1) then
            cx = -cx
          else
            cy = -cy
            cz = -cz
          endif
        endif
        if(yp .lt. y0) then                             ! ZX-reflection.
          if(idb .eq. 1) then
            cy = -cy
          else
            cx = -cx
            cz = -cz
          endif
        endif
        if(zp .lt. z0) then                             ! XY-reflection.
          if(idc .eq. 1) then
            cz = -cz
          else
            cx = -cx
            cy = -cy
          endif
        endif
      endif
c
      ierr = 0
999   return
c
901   ierr = 1
      messg = 'Point is outside X-limits'
      goto 999
902   ierr = 2
      messg = 'Point is outside Y-limits'
      goto 999
903   ierr = 3
      messg = 'Point is outside Z-limits'
      goto 999
904   ierr = 4
      messg = 'Non measured point - probably iron'
      goto 999
c
      end
