Final version of field_import.m
This commit is contained in:
shiroka
parent
76c7f4062e
commit
38e0151aef
@ -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);
|
||||
|
||||
% 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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user