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