Code covered by the BSD License  

Highlights from
Bloch Equation Vector Diagram Simulation Toolbox

image thumbnail

Bloch Equation Vector Diagram Simulation Toolbox

by

 

01 Feb 2012 (Updated )

Two functions that carry out pulses and evolution time periods for an array of spins.

[Spinsout Sim]=pulse(Spins,Sim)
%% Pulse Function
% Ohio Advanced EPR Laboratory
% Robert McCarrick
% Vector Simulation Package, Version 2.0
% Last Modified 1/23/2012

function [Spinsout Sim]=pulse(Spins,Sim)

%% Checks for variables and falls back on a default value if missing

if (isfield(Spins,'degrees') == 0);
    error('No spin array (.degrees) specified in structure.');
end

if (isfield(Sim,'angle') == 0);
    error('No angle (.angle) specified in structure.');
end

if (isfield(Sim,'direction') == 0);
    Sim.direction='-y';
end

if (isfield(Sim,'framerate') == 0);
    Sim.framerate=50;
end

if (isfield(Sim,'nframes') == 0);
    Sim.nframes=round(Sim.angle/2);
end

if (isfield(Sim,'video') == 0);
    Sim.video='n';
end

if strcmp(Sim.video,'gif') | strcmp(Sim.video,'avi');
    if (isfield(Sim,'filename') == 0);
        error('Please specify filename for animation.');
    end
end

if (isfield(Sim,'numloops') == 0);
    Sim.numloops=inf;
end

if (isfield(Sim,'figuresize') == 0);
    Sim.figuresize= [900 500];
end

%%  Sets the figure size and position

figurehandle = get(0,'CurrentFigure');
screensize = get(0,'ScreenSize');
currentposition = get(figurehandle,'Position');
desiredposition = [50, screensize(4)-Sim.figuresize(2)-150, Sim.figuresize(1), Sim.figuresize(2)];

    if isempty(figurehandle) == 1
        figure('Position',[50, screensize(4)-Sim.figuresize(2)-150, Sim.figuresize(1), Sim.figuresize(2)]);
    elseif isequal(desiredposition,currentposition) == 0
        set(gcf,'Position',[50, screensize(4)-Sim.figuresize(2)-150, Sim.figuresize(1), Sim.figuresize(2)]);
    else
    end

%% Sets the rotation axis based on the input

angle = Sim.angle/360*2*pi;
theta=angle/(Sim.nframes);

if (Sim.direction == 'x');
    rotmat = [1 0 0;0 cos(-theta) -sin(-theta);0 sin(-theta) cos(-theta)];
elseif (Sim.direction == 'y');
    rotmat = [cos(-theta) 0 sin(-theta);0 1 0;-sin(-theta) 0 cos(-theta)];
elseif (Sim.direction == '-x');
    rotmat = [1 0 0;0 cos(theta) -sin(theta);0 sin(theta) cos(theta)];
elseif (Sim.direction == '-y');
    rotmat = [cos(theta) 0 sin(theta);0 1 0;-sin(theta) 0 cos(theta)];
end

%% Gets the number of spins in the array

mdegrees = size(Spins.degrees);
starts = zeros(mdegrees(1),3);

%% Converts spins to rectangular coordinates

coordinates = zeros(mdegrees(1),3);

for i=1:mdegrees(1)
    x = sin(Spins.degrees(i,2)/360*2*pi)*cos(Spins.degrees(i,1)/360*2*pi)*Spins.degrees(i,3);
    y = sin(Spins.degrees(i,2)/360*2*pi)*sin(Spins.degrees(i,1)/360*2*pi)*Spins.degrees(i,3);
    z = cos(Spins.degrees(i,2)/360*2*pi)*Spins.degrees(i,3);
    coordinates(i,:)=[x,y,z];
end

%% Creates time axes for the plots

time = 0:Sim.angle/Sim.nframes:Sim.angle;
xsum = zeros(1,numel(time));
ysum = zeros(1,numel(time));
zsum = zeros(1,numel(time));

if strcmp(Sim.video,'n');
%% Performs the incremental rotation and plots each frame with no animation output

for j=1:Sim.nframes+1
    if j ~= 1
        coordinates = coordinates*rotmat;
    end
    ends=coordinates;
	xsum(1,j)=round(sum(ends(:,1))/mdegrees(1)*100)/100;
	ysum(1,j)=round(sum(ends(:,2))/mdegrees(1)*100)/100;
    zsum(1,j)=round(sum(ends(:,3))/mdegrees(1)*100)/100;
    subplot('Position',[0.03 0.03 0.4 0.9]);
    quiver3(starts(:,1),starts(:,2),starts(:,3),coordinates(:,1),coordinates(:,2),coordinates(:,3));
    hold on
    plot3([1 0],[0 0],[0 0],'k')
    plot3([0 0],[1 0],[0 0],'k')
    plot3([0 0],[0 0],[1 0],'k')
    text(1.1,0,0,'X_R')
    text(0,1.1,0,'Y_R')
    text(0,0,1.1,'Z_R')
    hold off
    set(gca,'XDir','rev','YDir','rev','GridLineStyle','none','XTick',[],'YTick',[],'ZTick',[]);
    title('Vector Diagram Simulation')
    axis ([-1 1 -1 1 -1 1]);
    subplot('Position',[0.5 0.76 0.47 0.2]);
    plot(time,xsum)
    axis([min(time) max(time) -1 1])
    ylabel('M_X');
    xlabel('degree rotation');
    subplot('Position',[0.5 0.43 0.47 0.2]);
    plot(time,ysum)
    axis([min(time) max(time) -1 1])
    ylabel('M_Y');
    xlabel('degree rotation');
    subplot('Position',[0.5 0.1 0.47 0.2]);
    plot(time,zsum)
    axis([min(time) max(time) -1 1])
    ylabel('M_Z');
    xlabel('degree rotation');
    pause(1/Sim.framerate);
end

% Converts the rectangular coordinates back to spherical coordinates

r = zeros(1,mdegrees(1));

for k=1:mdegrees(1)
    r(k)=((ends(k,1)-0)^2+(ends(k,2)-0)^2+(ends(k,3)-0)^2)^0.5;
end

angles = zeros(mdegrees(1),2);

for l=1:mdegrees(1)
    angles(l,:)= [atan2(ends(l,2),ends(l,1))*180/pi acos(ends(l,3)/r(l))*180/pi];
end

Spinsout = [angles r'];

elseif strcmp(Sim.video,'gif');
%% Performs the incremental rotation and plots each frame with animated GIF output

for j=1:Sim.nframes+1
    if j ~= 1
        coordinates = coordinates*rotmat;
    end
    ends=coordinates;
	xsum(1,j)=round(sum(ends(:,1))/mdegrees(1)*100)/100;
	ysum(1,j)=round(sum(ends(:,2))/mdegrees(1)*100)/100;
    zsum(1,j)=round(sum(ends(:,3))/mdegrees(1)*100)/100;
    subplot('Position',[0.03 0.03 0.4 0.9]);
    quiver3(starts(:,1),starts(:,2),starts(:,3),coordinates(:,1),coordinates(:,2),coordinates(:,3));
    hold on
    plot3([1 0],[0 0],[0 0],'k')
    plot3([0 0],[1 0],[0 0],'k')
    plot3([0 0],[0 0],[1 0],'k')
    text(1.1,0,0,'X_R')
    text(0,1.1,0,'Y_R')
    text(0,0,1.1,'Z_R')
    hold off
    set(gca,'XDir','rev','YDir','rev','GridLineStyle','none','XTick',[],'YTick',[],'ZTick',[]);
    title('Vector Diagram Simulation')
    axis ([-1 1 -1 1 -1 1]);
    subplot('Position',[0.5 0.76 0.47 0.2]);
    plot(time,xsum)
    axis([min(time) max(time) -1 1])
    ylabel('M_X');
    xlabel('degree rotation');
    subplot('Position',[0.5 0.43 0.47 0.2]);
    plot(time,ysum)
    axis([min(time) max(time) -1 1])
    ylabel('M_Y');
    xlabel('degree rotation');
    subplot('Position',[0.5 0.1 0.47 0.2]);
    plot(time,zsum)
    axis([min(time) max(time) -1 1])
    ylabel('M_Z');
    xlabel('degree rotation');
    frame = getframe(gcf);
    image = frame2im(frame);
    [imind,cm] = rgb2ind(image,256);
    if exist(Sim.filename,'file') == 0;
        imwrite(imind,cm,Sim.filename,'gif','Loopcount',Sim.numloops,'DelayTime',1/Sim.framerate);
        Sim.gifout = 1;
    else
        if (isfield(Sim,'gifout') == 0);
            error('File already exists!');
        else
            imwrite(imind,cm,Sim.filename,'gif','WriteMode','append','DelayTime',1/Sim.framerate);
        end
    end
end

% Converts the rectangular coordinates back to spherical coordinates

r = zeros(1,mdegrees(1));

for k=1:mdegrees(1)
    r(k)=((ends(k,1)-0)^2+(ends(k,2)-0)^2+(ends(k,3)-0)^2)^0.5;
end

angles = zeros(mdegrees(1),2);

for l=1:mdegrees(1)
    angles(l,:)= [atan2(ends(l,2),ends(l,1))*180/pi acos(ends(l,3)/r(l))*180/pi];
end

Spinsout = [angles r'];

elseif strcmp(Sim.video,'avi');
%% Performs the incremental rotation and plots each frame with AVI output

    if exist(Sim.filename,'file') == 0;

        aviobj = VideoWriter(Sim.filename);
        open(aviobj)

        for j=1:Sim.nframes+1
            if j ~= 1
                coordinates = coordinates*rotmat;
            end
            ends=coordinates;
        	xsum(1,j)=round(sum(ends(:,1))/mdegrees(1)*100)/100;
        	ysum(1,j)=round(sum(ends(:,2))/mdegrees(1)*100)/100;
            zsum(1,j)=round(sum(ends(:,3))/mdegrees(1)*100)/100;
            subplot('Position',[0.03 0.03 0.4 0.9]);
            quiver3(starts(:,1),starts(:,2),starts(:,3),coordinates(:,1),coordinates(:,2),coordinates(:,3));
            hold on
            plot3([1 0],[0 0],[0 0],'k')
            plot3([0 0],[1 0],[0 0],'k')
            plot3([0 0],[0 0],[1 0],'k')
            text(1.1,0,0,'X_R')
            text(0,1.1,0,'Y_R')
            text(0,0,1.1,'Z_R')
            hold off
            set(gca,'XDir','rev','YDir','rev','GridLineStyle','none','XTick',[],'YTick',[],'ZTick',[]);
            title('Vector Diagram Simulation')
            axis ([-1 1 -1 1 -1 1]);
            subplot('Position',[0.5 0.76 0.47 0.2]);
            plot(time,xsum)
            axis([min(time) max(time) -1 1])
            ylabel('M_X');
            xlabel('degree rotation');
            subplot('Position',[0.5 0.43 0.47 0.2]);
            plot(time,ysum)
            axis([min(time) max(time) -1 1])
            ylabel('M_Y');
            xlabel('degree rotation');
            subplot('Position',[0.5 0.1 0.47 0.2]);
            plot(time,zsum)
            axis([min(time) max(time) -1 1])
            ylabel('M_Z');
            xlabel('degree rotation');
            aviframes(j) = getframe(gcf);
        end

        writeVideo(aviobj,aviframes)
        close(aviobj)
 
        % Converts the rectangular coordinates back to spherical coordinates
        
        r = zeros(1,mdegrees(1));
        
        for k=1:mdegrees(1)
            r(k)=((ends(k,1)-0)^2+(ends(k,2)-0)^2+(ends(k,3)-0)^2)^0.5;
        end

        angles = zeros(mdegrees(1),2);
        
        for l=1:mdegrees(1)
            angles(l,:)= [atan2(ends(l,2),ends(l,1))*180/pi acos(ends(l,3)/r(l))*180/pi];
        end

        Sim.aviout = aviframes;
        Spinsout = [angles r'];

    else

        
        if (isfield(Sim,'aviout') == 0);
            error('File already exists!');
        end
        
        aviobj = VideoWriter(Sim.filename);
        open(aviobj)
        aviframes = Sim.aviout;
        framenumber = size(Sim.aviout);

        for j=1:Sim.nframes+1
            if j ~= 1
                coordinates = coordinates*rotmat;
            end
            ends=coordinates;
            xsum(1,j)=round(sum(ends(:,1))/mdegrees(1)*100)/100;
            ysum(1,j)=round(sum(ends(:,2))/mdegrees(1)*100)/100;
            zsum(1,j)=round(sum(ends(:,3))/mdegrees(1)*100)/100;
            subplot('Position',[0.03 0.03 0.4 0.9]);
            quiver3(starts(:,1),starts(:,2),starts(:,3),coordinates(:,1),coordinates(:,2),coordinates(:,3));
            hold on
            plot3([1 0],[0 0],[0 0],'k')
            plot3([0 0],[1 0],[0 0],'k')
            plot3([0 0],[0 0],[1 0],'k')
            text(1.1,0,0,'X_R')
            text(0,1.1,0,'Y_R')
            text(0,0,1.1,'Z_R')
            hold off
            set(gca,'XDir','rev','YDir','rev','GridLineStyle','none','XTick',[],'YTick',[],'ZTick',[]);
            title('Vector Diagram Simulation')
            axis ([-1 1 -1 1 -1 1]);
            subplot('Position',[0.5 0.76 0.47 0.2]);
            plot(time,xsum)
            axis([min(time) max(time) -1 1])
            ylabel('M_X');
            xlabel('degree rotation');
            subplot('Position',[0.5 0.43 0.47 0.2]);
            plot(time,ysum)
            axis([min(time) max(time) -1 1])
            ylabel('M_Y');
            xlabel('degree rotation');
            subplot('Position',[0.5 0.1 0.47 0.2]);
            plot(time,zsum)
            axis([min(time) max(time) -1 1])
            ylabel('M_Z');
            xlabel('degree rotation');
            aviframes(j+framenumber(2)) = getframe(gcf);
        end
    
        writeVideo(aviobj,aviframes)
        close(aviobj)

        % Converts the rectangular coordinates back to spherical coordinates
        
        r = zeros(1,mdegrees(1));
        
        for k=1:mdegrees(1)
            r(k)=((ends(k,1)-0)^2+(ends(k,2)-0)^2+(ends(k,3)-0)^2)^0.5;
        end

        angles = zeros(mdegrees(1),2);
        
        for l=1:mdegrees(1)
            angles(l,:)= [atan2(ends(l,2),ends(l,1))*180/pi acos(ends(l,3)/r(l))*180/pi];
        end

        Sim.aviout = aviframes;
        Spinsout = [angles r'];

    end
end

Contact us