Asked by fang-chu
on 6 Jul 2013

clc %parameters B1=10; gamma = 267.513*10^6; %rad/(s*T) Bzi = -5; %in HZ Bzj = 5; %in HZ a = 1000000; G = 42.576; %Gyromagnetic ratio (MHz/T) MHZ=1*10^6 B_zi = Bzi/(a*G); %in tesla B_zj = Bzj/(a*G); %in tesla mag_initial = [0 0 1; 0 0 1]; dt = 10;

%functions % equation magnetizaion = [m_x1 m_y1 m_z1; m_x2 m_y2 m_z2]; M = mean(mag); %average of mag m_x = M(1); m_y = M(2); euler = complex(cos(theta),sin(theta)); m = complex(m_x,m_y); B0 = euler*m; B_x = real(B0); B_y = imag(B0); B= [B_x+B1 B_y B_zi; B_x+B1 B_y B_zj];

%time evolution for n = 1:10 mag = 0; mag = mag + n; [t mag] = ode45('bloch', [0,dt/2,dt], [mag_initial, B]); end

I try to get a series of matrix- mag by having the output to be the input of the matrix at the next time step.

My odefun is function dmdt = bloch(gamma,mag,B) %outputs and inputs mag = reshape(mag,length(mag)/3,3); %length(array) finds the number of elements along the largest dimension of an array dmdt = gamma*cross(mag,B); return

*No products are associated with this question.*

Answer by Jan Simon
on 7 Jul 2013

Accepted answer

Your code is very strange. I cannot understand its intention, but it does not seem to be a task for an ODE integrator. It is not clear, how the function `bloch` could obtain the 3rd input "B" or "gamma".

`euler = complex(cos(theta),sin(theta));` fails due to a missing definition of "theta".

This looks confused:

mag = 0; mag = mag + n;

What about this instead:

mag = n

The initial values of an integration have to be a vector, not a matrix.

Summary: I cannot understand the intention of the code. It should throw a lot of error, when you try to run it. Without knowing what you want to achieve, suggesting an improvement is not possible.

Opportunities for recent engineering grads.

## 1 Comment

## Jan Simon (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/81300#comment_158659

Do you see, that omitting the code formatting makes your code unreadable? Do you think of improving this as explained in the "? Help" link?