Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab: help with biot savart law

  1. Apr 22, 2012 #1
    Hi, if anyone good with matlab and knows biot savart law then i hope you can help. I have the following program:

    Code (Text):

     clear all
     Img = imread('littlecircle.png');
     Img = Img(:,:,1);
     Img = double(Img);
     w = size(Img,1);               % width size
     h = size(Img,2);               % height size
     [Ix,Iy] = gradient(Img);       %gradient of image
     i=1;     %iteration for magnetic field loop
     b=0;     %initialize b to zero
     % Magnetic Field
     for pxRow = 1:h % fixed pixel row
     for pxCol = 1:w % fixed pixel column
     for r = 1:h % row of distant pixel
     for c = 1:w % column of distant pixel
     R(c,r) = sqrt((r-pxRow)^2 + (c-pxCol)^2);                               % pythagoras theorem to get distance to each pixel
    [COLOR="Red"] O(c,r) = atan(Iy(c,r)./Ix(c,r));                         % direction of current
     If(c,r) = sqrt((Ix(c,r)).^2 + (Iy(c,r)).^2);                            % magnitude of current
     Rxs(c,r) = R(c,r)./norm(R(c,r));                                        % unit vector from x to s                    
     b(c,r) = If(c,r).*O(c,r).*(Rxs(c,r)./(R(c,r)).^2);                      % b field = If(s)O(s) x Rxs/R.^2  BIOT SAVART LAW[/COLOR]Rxs/R.^2  BIOT SAVART LAW
     B(i) = {b}; % filling a cell array with results. read below
     i = i+1;
    all lines in red are currently not calculating what I want them to. Can anyone help? Thanks!
  2. jcsd
  3. Apr 22, 2012 #2
  4. Apr 24, 2012 #3
    Debugging hints/help:

    They are not calculating what you want them to?
    What ARE they calculating? and
    what do YOU expect them to be?

    Are you sure your starting numbers are correct (Ix,Iy) ?

    What does the function 'norm' do for you?
    You are taking the norm of scalar R(c,r)...what is that? Its length? Its magnitude? Is it always a positive number?
    By the way, you just calculated R(c,r) above as the sqrt() of something...so, R(c,r) IS a scalar and IS positive; also, it is just a scalar...not a matrix.
    So....what do you expect Rxs(c,r) = R(c,r)./norm(R(c,r)) to be? You call it a unit vector...but I don't think it is a vector at all...it is just a number.

    It seems to me that Rxs(c,r) will always be 1.0 ?
  5. Apr 24, 2012 #4
    I am struggling with getting to grips with vectors in matlab i'll admit. but here's what i need to do:

    1. read in an image
    2. compute image gradient vectors [Ix,Iy], and then magnitude sqrt(Ix.^2+Iy.^2) = If(c,r)
    3. compute the edge orientation of the image by rotating the gradient vectors 90 degrees clockwise or anti clockwise (O(c,r))
    4. take each pixel in the image, and find its distance to all other pixels (R(c,r))
    5. compute the unit vectors of the pixel distances Rxs(c,r)
    6. Compute b field with Biot-Savart law B = If.*O cross (Rxs/R.^2)

    and yes my understanding of how matlab handles vectors isn't helping the issue.
  6. Apr 24, 2012 #5
    I hear you...but you did not answer ANY of my questions. :-(
  7. Apr 25, 2012 #6
    sorry i'll have a go -

    only R(c,r) calculates what I want. it calculates the distance one pixel to all others and stores into an array.

    O(c,r) is supposed to be the current direction, ie taking the gradient vectors Ix(c,r) and Iy(c,r) and rotating them by 90 degrees. When I quiver plot these rotated vectors it looks as expected, however I am told that the actual result of these are just scalars, ie numbers - but they need to be vectors.

    Rxs(c,r) is supposed to be the unit vector of R(c,r). Matlab computes unit vectors with the norm function:

    Rxs(c,r) = R(c,r)/norm(R(c,r)) give me the unit vectors, and they are indeed all supposed = 1

    Where I could be going wrong, as you say is that none of the above are vectors in which I have a problem. I am new to matlab and when you do your sums on the blackboard you take the nature of vectors for granted and just punch them in the calculator. but i'm having difficulty representing them here in matlab.

    The clue "you are taking the norm of a scalar R(c,r)" as you say.....well how do I turn R(c,r) into a vector? If i knew this im sure the rest of the program would fall into place. I have been relying on matlabs quiver function (vector plot) to show me the direction of these numbers, and they do indeed have direction, hence why i'm assuming I'm returning vectors.....argh!
  8. Apr 25, 2012 #7
    Yeap...you are struggling.

    When on the blackboard you take the nature of vectors granted?
    What do vectors look like?

    Recall that a vector has a magnitude AND a direction (angle).

    So, step back a bit and take a good look at your own code as if you had not written it...just because you had something in mind and had certain intent, that does not mean that it went through...take a look at your calculations and you will see that you are just doing a bunch of scalar operations...nothing there is a vector...

    That is not necessarily bad, i.e., keeping track of your own magnitudes and your own angles separately is o.k....if you know that you are doing it!.

    Other than that, you say that O(c,r) is a current direction (angle)? I guess it is, after all, it is the arctan() of something; so it is an angle, but just an angle, not a vector...are these the angles of your vectors? Where are the corresponding magnitudes? If()?

    And Rxs() is not giving you unit vectors, it is just giving you units!

    Alternatively...if you want to take vectors for granted in matlab, than you first need to find out how to handle vectors in matlab...really, vectors, you know, those things that do have magnitude and angle...can matlab handle this quantity in a single entity, variable? I think it might but I have not used matlab in about 20 years...so, take a few well invested minutes and do a google search or read the user's manual...I guess I could do it, too; but it is best if you learn this by yourself.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook