240 lines
8.1 KiB
Matlab
240 lines
8.1 KiB
Matlab
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);
|
|
|
|
|
|
% 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);
|
|
|
|
%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
|
|
|
|
|
|
% 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);
|
|
|
|
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
|
|
|
|
|
|
% 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
|