From 38e0151aeffdcb88b055e9a64a9a757bfcffcd98 Mon Sep 17 00:00:00 2001 From: shiroka Date: Thu, 8 Jan 2009 14:15:16 +0000 Subject: [PATCH] Final version of field_import.m --- geant4/spin_rot/run/field_import.m | 249 ++++++++++++++++++++++++++--- 1 file changed, 228 insertions(+), 21 deletions(-) diff --git a/geant4/spin_rot/run/field_import.m b/geant4/spin_rot/run/field_import.m index 8d4ecaf..38bf732 100644 --- a/geant4/spin_rot/run/field_import.m +++ b/geant4/spin_rot/run/field_import.m @@ -1,32 +1,239 @@ -[x y z bx by bz] = textread('sep_0nl.lis','%f %f %f %f %f %f', 'headerlines', 13); -[nx ny nz] = textread('sep_0nl.lis','%*s %*s %*s %*s %u %u %u', 1, 'headerlines', 10); +function field_import(varargin) +%FIELD_IMPORT Imports B-field maps calculated by Opera in G4. +% B-field maps calculated by Opera in the first (positive) octant +% are extended to all octants and modified for importing in Geant4. +% Notice that the variable export order in Opera is x->y->z, while +% Geant4 expects variables to vary in the opposite order z->y->x. +% Notice that the original By values have been multiplied by -1. +% USAGE: field_import or field_import(B_map). +% OUTPUT: Extended B-map with name B_map_G4.dat +% +% See also ELMAT, FLIPDIM, CAT, SUB2IND, PERMUTE. - +% Copyright (c) 2008 by Toni Shiroka +% $Revision: 1.0 $ $Date: 2008/12/31 13:45:13 $ + +error(nargchk(0,1,nargin)); +wd=cd; +if nargin == 0, + %[filename, pathname] = uigetfile('*.lis', 'Choose a B-field map to read');%, 10, 24); + [filename, pathname] = uigetfile({'*.lis';'*.dat';'*.txt';'*.*'}, 'Choose a B-field map to read'); + fnam=filename; direc=pathname; + if filename==0, + disp(' '); disp('No file loaded ...'); + disp(' '); cd(wd); return + end +end + +if nargin == 1, + fnam = char(varargin); + if ~exist(fnam,'file'), + disp(' '); disp(['File ',fnam,' does not exist. No file loaded ...']); + disp(' '); cd(wd); return + end +end + +%Read variables, step sizes and number of points +[x y z bx by bz] = textread(fnam,'%f %f %f %f %f %f', 'headerlines', 13); +[units sx sy sz] = textread(fnam,'%*s %s %*s %f %f %f', 1, 'headerlines', 9); +[nx ny nz] = textread(fnam,'%*s %*s %*s %*s %u %u %u', 1, 'headerlines', 10); + +disp(' '); disp(['File ',fnam,' read successfully. Processing data ...']); + +% Prepare length units for later use +% Use length units as mm, cm or m, but NOT dm! +units = char(units); +str_s = findstr(units,'['); str_e = findstr(units,']'); +l_unit = units(str_s+1:str_e-1); + + +% Recreate B field in the first octant Bx = reshape(bx,nx,ny,nz); By = reshape(by,nx,ny,nz); Bz = reshape(bz,nx,ny,nz); -x1 = reshape(x,nx,ny,nz); -y1 = reshape(y,nx,ny,nz); -z1 = reshape(z,nx,ny,nz); -i=1:32; -test = [x1(i);y1(i);z1(i)]' +% Extend B-field map in negative x direction: +Bx_extx = -flipdim(Bx,1); % Notice the NEGATIVE sign! +By_extx = flipdim(By,1); +Bz_extx = flipdim(Bz,1); -%min(x) max(x) +%The (1:end-1) argument is to avoid repeating the origin Bx_extx(end) twice! +Bx_X = cat(1, Bx_extx(1:end-1,:,:), Bx); % Extend Bx in X direction +By_X = cat(1, By_extx(1:end-1,:,:), By); % Extend By in X direction +Bz_X = cat(1, Bz_extx(1:end-1,:,:), Bz); % Extend Bz in X direction -sx = max(diff(x)); % read step size along x -sy = max(diff(y)); % read step size along y -sz = max(diff(z)); % read step size along z -x2 = -x(end:-1:1); -x3 = -x(end:-1:1); -x4 = -x(end:-1:1); +% Extend B-field map in negative y direction: +Bx_exty = flipdim(Bx_X,2); +By_exty = flipdim(By_X,2); % Notice the ABSENCE of negative sign due to -1 in orig. file! +Bz_exty = flipdim(Bz_X,2); -x2 = reshape(x2,nx,ny,nz); -x3 = reshape(x3,nx,ny,nz); -x4 = reshape(x4,nx,ny,nz); +Bx_XY = cat(2, Bx_exty(:,1:end-1,:), Bx_X); % Extend Bx_X in Y direction +By_XY = cat(2, By_exty(:,1:end-1,:), By_X); % Extend By_X in Y direction +Bz_XY = cat(2, Bz_exty(:,1:end-1,:), Bz_X); % Extend Bz_X in Y direction -%p = [x(:)'; y(:)'; z(:)']; -%ind2sub -%[i j] = ind2sub(3,x1); \ No newline at end of file + +% Extend B-field map in negative z direction: +Bx_extz = -flipdim(Bx_XY,3); % Notice the NEGATIVE sign! +By_extz = flipdim(By_XY,3); % Notice the ABSENCE of negative sign due to -1 in orig. file! +Bz_extz = flipdim(Bz_XY,3); + +Bx_XYZ = cat(3, Bx_extz(:,:,1:end-1), Bx_XY); % Extend Bx_XY in Z direction +By_XYZ = cat(3, By_extz(:,:,1:end-1), By_XY); % Extend By_XY in Z direction +Bz_XYZ = cat(3, Bz_extz(:,:,1:end-1), Bz_XY); % Extend Bz_XY in Z direction + + +% The final B-field map is stored in Bx_XYZ, By_XYZ, Bz_XYZ. +% The coordinates are reconstructed using the number of points and the step size. + + + + +% Test plots: By is the main magnetic field component + +k = 45; +figure(1) % plot original By field map (xy projection) +Byo_plot = squeeze(By(:,:,k)); +[c,h] = contourf(Byo_plot); clabel(c,h); +axis('square') + + +figure(2) % plot extended By field map (xy projection) +Bye_plot = squeeze(By_XYZ(:,:,nz + k)); % nz +/- k +[x_plot, y_plot] = meshgrid(-x(end):sx:x(end), -y(end):sy:y(end)); +[c,h] = contourf(x_plot, y_plot, Bye_plot); clabel(c,h); +axis('square') +xlabel('x-axis (cm)','FontSize',14) +ylabel('y-axis (cm)','FontSize',14) +title(['xy section of the extended B_y-field map at z = ',num2str(sz*k),' cm'],'FontSize',16) + +%for kk = 1:4:2*nz; +% Bye_plot = squeeze(By_XYZ(:,:,kk)); +% [c,h] = contourf(Bye_plot); clabel(c,h); +% axis('square') +% hold on +% quiver(squeeze(Bx_XYZ(:,:,nz+k)),squeeze(By_XYZ(:,:,nz+k))); +% drawnow +%end + +p = 3; +figure(3) % plot Bz original field map (yz projection) +Bzo_plot = squeeze(Bz(p,:,:)); +[c,h] = contourf(Bzo_plot); clabel(c,h); +axis('equal','tight') + +figure(4) % plot Bz extended field map (yz projection) +Bze_plot = squeeze(Bz_XYZ(nx+p,:,:)); +[y_plot, z_plot] = meshgrid(-y(end):sy:y(end), -z(end):sz:z(end)); +[c,h] = contourf(y_plot, z_plot, Bze_plot'); clabel(c,h); +view(90,90) +axis('equal','tight') +xlabel('y-axis (cm)','FontSize',14) +ylabel('z-axis (cm)','FontSize',14) +title(['yz section of the extended B_z-field map at x = ',num2str(sx*p),' cm'],'FontSize',16) + + +%for pp = 1:2*nx, +% Bze_plot = squeeze(Bz_XYZ(pp,:,:)); +% [c,h] = contourf(Bze_plot); clabel(c,h); pause(1) +% axis('equal','tight') +% hold on +% quiver(squeeze(By_XYZ(nx+p,:,:)),squeeze(Bz_XYZ(nx+p,:,:))); +% drawnow +%end + + +figure(5) % plot By extended field map (yz projection) +Bye_plot = squeeze(By_XYZ(nx+p,:,:)); +%[c,h] = contourf(Bye_plot); clabel(c,h); +[c,h] = contourf(y_plot, z_plot, Bye_plot'); clabel(c,h); +view(90,90) +axis('equal','tight') +xlabel('y-axis (cm)','FontSize',14) +ylabel('z-axis (cm)','FontSize',14) +title(['yz section of the extended B_y-field map at x = ',num2str(sx*p),' cm'],'FontSize',16) + + +figure(6) % plot By extended field map along the z-axis +Bye_zplot = squeeze(By_XYZ(nx,ny,:)); +plot(-z(end):sz:z(end),Bye_zplot); +xlabel('z-axis (cm)','FontSize',14) +ylabel('Field value (arb. units)','FontSize',14) +title('Extended B_y-field map along {\itz} axis','FontSize',16) + + + + + + +%Prepare extended coord. and field maps for writing into the final file. + +% Extend x, y, z coordinates to negative values. +% Use ndgrid for reverse order, or meshgrid for mixed order! +% Variation order: ndgrd: z -> y -> x, ndgrid: x -> y -> z +[x_ext, y_ext, z_ext] = ndgrd(-x(end):sx:x(end), -y(end):sy:y(end), -z(end):sz:z(end)); + +%Change also the B order according to the same coord. order, i.e.: z -> y -> x +Bx_ext = permute(Bx_XYZ, [3 2 1]); +By_ext = permute(By_XYZ, [3 2 1]); +Bz_ext = permute(Bz_XYZ, [3 2 1]); + +%Export coordinates and field values +data = [x_ext(:) y_ext(:) z_ext(:) Bx_ext(:) By_ext(:) Bz_ext(:)]; +%Export field values only +%data = [Bx_ext(:) By_ext(:) Bz_ext(:)]; + +imax = find(data(:,1)==0 & data(:,2)==0 & data(:,3)==0); % Center of magnet +M = data(imax,5); % Consider as max. the By field value at magnet center! +%M = max(max(data(:,4:6))); % consider only the last three columns (Bx, By and Bz) +norm_fact = 1/M; % normalization factor + +% The filename of the extended field map +fname = [strtok(fnam,'.'),'_ext.map']; + +% Prepare data for writing +header1 = sprintf('%d\t%d\t%d\t%s\t%.8f', 2*nx-1, 2*ny-1, 2*nz-1, l_unit, norm_fact); +header2 = sprintf('%s\t%s\t%s\t%s\t%s', '%nx', 'ny', 'nz', 'length_unit', 'norm_fact'); + +disp(' ') +disp('Data being written to the extended field map. Please wait ...') +%Write data to the extended field map file +fid = fopen(fname,'w'); +fprintf(fid,'%s\n',header1); +fprintf(fid,'%s\n',header2); +fprintf(fid,'%s\n','%'); +fprintf(fid,'%2.5g\t%2.5g\t%2.5g\t%2.5g\t%2.5g\t%2.5g\n',data'); +fclose(fid); + +% This works only in more advanced matlab versions! +% dlmwrite(filename,header,'delimiter','') +% dlmwrite(filename,data,'delimiter','\t','precision','%-2.5g','-append'); +% try also f, or e instead of g, or fprintf + +disp(['Data successfully written to ',fname,'.']) +disp(' ') + + +return + + +% Initial check of correctness +%[Bx_extx(:,2,13),Bx(end:-1:1,2,13)] +%[By_extx(:,2,13),By(end:-1:1,2,13)] +%[Bz_extx(:,2,13),Bz(end:-1:1,2,13)] + + +% Final Check of correctness +% Extended data +%data(237000,:) +%4.0 2.0 1.5 -2.46 128.60 -3.06 +%data(237002,:) +%4.0 2.0 2.5 -2.42 129.63 -5.13 + +% Original data! +%datao = [x(:) y(:) z(:) bx(:) by(:) bz(:)]; +%datao(744,:) +%4.0 2.0 1.5 -2.46 128.60 -3.06 +%datao(1194,:) +%4.0 2.0 2.5 -2.42 129.63 -5.13