% calculate the gradients waveforms for the k-space trajectory given by
% k_x = alpha*t*cos(omega*t);
% k_y = alpha*t*sin(omega*t);

clear all
close all

%% set the fixed parameters
dt = 0.001; % time sampling interval
L = 1; % field of view (valid for both spatial directions)
N = 64; % number of voxels for each direction
gamma = .01; % the physical constant
%% define a good value for omega:
% according to Nyquist criterion
omega_Nyquist = 2/(dt*N);
% according to constraint on gradient amplitude
Gmax = 400000;
omega_Gmax = Gmax*L*gamma/(pi*N);
% choose value which fullfils both constraints:
omega = min([omega_Nyquist,omega_Gmax]) ;

%% now set alpha, the other parameter for the trajectory
alpha = omega/(2*pi*L);
%% calculate and plot gradients
Nth = 2*pi/(omega*dt); % number of sample points in one "circle"
% total number of samples along trajectory is thus Nth*N/2 
t = 0:dt:(Nth*N/2)*dt; % the time samples 
G_x = alpha*2*pi/gamma*(cos(omega*t)-omega*t.*sin(omega*t));
G_y = alpha*2*pi/gamma*(sin(omega*t)+omega*t.*cos(omega*t));
figure;
plot(t,G_x);hold on;
plot(t,G_y);hold on;
plot(t,alpha*2*pi/gamma*omega*t);hold off;xlabel('t')
legend('Gx','Gy','2\alpha\omega\pi t/\gamma'); title('The gradient functions')
figure;plot((alpha*t.*exp(1i*omega*t)),'-o');axis image;xlabel('k_x');ylabel('k_y');title('k-space trajectory')