% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creates MagVenture_MC_B35.ccd
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


close all
%clear all
% %%%%%%%%%%%%
% Coil Constants -- coil defining parameters
% %%%%%%%%%%%%
% all intial dimensions in m or degree


ccd_filename = 'MagVenture_MC_B35.ccd';
VersionStr = '1.0';

Mirrored = 1;      % Two identical sides of the coil
nl = 3;            % number different coils/layers
a=[0 0 0];         % angle between coil and plane
Roff = [0.026 0.026 0.026];% offset of coil centre to zero-point in coil plane in x-direction
Zoff = 0.00710;      % Offset of internal coil-centre to isocentre on coil surface (in negative z-direction)
ri = [0.0135 0.0135 0.0135]; % inner radii of coils
ro = [0.0255 0.0255 0.0255]; % outer radii of coils
zoff = [0.003 0 -0.003];     % z-offset of coils from coil centre (in positive z-direction)
t = [0.003 0.003 0.003];         % thickness of each layer together
n = [4 4 4];         % Number of Windings for each layer
StepRadius = 0.002;     % steps for dipole radii in mm, lower values --> more dipoles
StepAngle = 0.3;     % steps for dipole angle in rad, lower values --> more dipoles


xCoG =[]; yCoG =[]; zCoG =[];
xDir =[]; yDir =[]; zDir =[];


for C = 1:nl
    % Calculate positions and weightings for each dipole moment

    % define radii inside the innermost ring in 3mm steps
    r = linspace(0,ri(C),int16(ri(C)/StepRadius));
    % and in between the Rings
    r(end:end+n(C)) = linspace(ri(C),ro(C),n(C)+1);

    % Number of Magnetic Dipoles for each C
    nPhi = double(int16(1000*r(1:end-1)/StepAngle));
    nPhi(1) = int16(1000*mean(r(1:2))/StepAngle); 

    % Number of Turns for each radius
    nTurns =  ones(1,length(r)-1)*n(C);
    nTurns(end:-1:end-n(C)+1) = (1:n(C)) - 0.5;
    %nTurns =  ones(1,length(r))*n(C);
    %nTurns(end:-1:end-n(C)+1) = (1:n(C)) - 0.5;


    % Calculate number of elements

    nElements = sum(nPhi);  
    offsetPhi = 0.001;      % Rotation of the elements

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CoG, Positions and area of elements
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    counter = 0;

    for i = 1:length(r)-1
      deltaPhi = 2*pi/nPhi(i);

      rLower = r(i);
      rUpper = r(i + 1);

      rCoG = 2*sqrt( 2*(1-cos(deltaPhi)))*(rLower^2 + rUpper^2 + rLower*rUpper) ...
        /(3*deltaPhi*(rLower + rUpper));  % rCoG: radius of the center of gravity

      areaHelp = deltaPhi*(rUpper^2 - rLower^2)/2;

      for j = 1:nPhi(i)
        phiCoG = (2*j-1)*deltaPhi/2;
        counter = counter +1;
        xCoGL(counter) = rCoG*cos(phiCoG + offsetPhi) + Roff(C);
        yCoGL(counter) = rCoG*sin(phiCoG + offsetPhi);
        areaOfElement(counter) = areaHelp;
        nTurnsElement(counter) = nTurns(i);
        zCoGL(counter) = zoff(C);
      end;

    end

    xCoGLtemp = xCoGL;

    xCoGL =  - zCoGL * sind(a(C)) + xCoGLtemp*cosd(a(C));
    zCoGL = zCoGL * cosd(a(C)) + xCoGLtemp*sind(a(C));

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %       Direction and weighting of Elements
    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    zDirL = areaOfElement.*nTurnsElement;

    % no angulation --> carry out now

    xDirL =  -zDirL*sind(a(C));
    zDirL =  zDirL*cosd(a(C));
    yDirL = zeros(1,length(areaOfElement));




    %       Split winding in layers of approx 3mm thickness
    nt = max(1,round(t(C)/3));
    xCoGL = repmat(xCoGL,[1 nt]);
    yCoGL = repmat(yCoGL,[1 nt]);
    zCoGL = repmat(zCoGL,[1 nt]);
    xDirL = repmat(xDirL,[1 nt])/nt;
    yDirL = repmat(yDirL,[1 nt])/nt;
    zDirL = repmat(zDirL,[1 nt])/nt;

    offnt = linspace(-t(C)/2,t(C)/2,2*nt+1);
    xoffnt = -sind(a(C))*offnt(2:2:end-1);
    zoffnt = cosd(a(C))*offnt(2:2:end-1);

    for j = 1:nt
        xCoGL((j-1)*nElements+1:(j)*nElements) = xCoGL((j-1)*nElements+1:(j)*nElements) + xoffnt(j);
        zCoGL((j-1)*nElements+1:(j)*nElements) = zCoGL((j-1)*nElements+1:(j)*nElements) + zoffnt(j);
    end

    xCoG = [xCoG xCoGL];  
    yCoG = [yCoG yCoGL];  
    zCoG = [zCoG zCoGL];  
    xDir = [xDir xDirL];
    yDir = [yDir yDirL];
    zDir = [zDir zDirL];

    clear xCoGL yCoGL zCoGL xDirL yDirL zDirL nTurnsElement areaOfElement;
end

zCoG = zCoG - Zoff;

% Create a Mirrored part of the coil if exists
if  Mirrored
    xCoG = [xCoG (-xCoG)]; 
    yCoG = [yCoG (-yCoG)]; 
    zCoG = [zCoG (zCoG)]; 
    xDir = [xDir xDir];
    yDir = [yDir yDir];
    zDir = [zDir -zDir];
end

% xCoG = xCoG/1000; yCoG = yCoG/1000; zCoG = zCoG/1000; % Convert from mm to meter


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Grafik 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;
plot3(xCoG,yCoG,zCoG,'Marker','.','LineStyle','none');
axis equal;

figure;
quiver3(xCoG,yCoG,zCoG,xDir,yDir,zDir,2,'.');
axis equal;
axis([-0.01+min(xCoG) 0.01+max(xCoG) -0.01+min(yCoG) 0.01+max(yCoG) -0.01+min(zCoG) 0.01+max(zCoG)]);



% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Datei schreiben
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



fid = fopen(ccd_filename,'w');

strhlp = sprintf(['# ',ccd_filename,' version ',VersionStr,' number of elements\n']);
fwrite(fid, strhlp, 'char');

strhlp = sprintf('%.0f\n',length(zCoG));
fwrite(fid, strhlp, 'char');

strhlp = sprintf('# centers and weighted directions of the elements (magnetic dipoles)\n');
fwrite(fid, strhlp, 'char');


for i = 1:length(zCoG)
  strhlp = sprintf('%.15e %.15e %.15e %.15e %.15e %.15e\n', xCoG(i), yCoG(i), ...
		   zCoG(i), xDir(i), yDir(i), zDir(i));
  fwrite(fid, strhlp, 'char');
end;


fclose(fid);