%source code by:Bijan Binaee <bij@n@binaee.com>
function solenoidII(a)
    step = 0.6;
    max = 3;
    if nargin < 1
      a = 0.7;
    end
    [x,y,z]=meshgrid(-max:step:max,-max:step:max,-max:step:max);
    %z=x
    U = zeros(size(x));
    V = zeros(size(x));
    W = zeros(size(x));
    s = 0;
    for i = 1 : numel(x)
        [ U(i) , V(i) , W(i) ] = H_feild (x(i) ,y(i) ,z(i),a);
        disp (ceil(i./numel(x) * 100));
        H = sqrt ( U(i)^2 + V(i)^2 + W(i) ^ 2 );
        if ( H > 3 )
            U(i) = U(i) ./ H * 3;
            V(i) = V(i) ./ H * 3;
            W(i) = W(i) ./ H * 3;
        end
        s = s + H;
    end
    disp(s);
    [r,theta] = meshgrid(0:a/2:a,linspace(0,2*pi,20));
    X = r.*cos(theta);
    Y = r.*sin(theta);
    Z = X - X;
    figure;
    surf(X,Y,Z);
    hold on
    quiver3(x,y,z,U,V,W,2);
    hold off
    figure;
    [sx,sy,sz] = meshgrid(-2.0:0.5:2.0,-2.0:0.5:2.0,0);
    streamline(stream3(x,y,z,U,V,W,sx,sy,sz));
    view(3);
    xlabel('X-axis','fontsize',14);
    ylabel('Y-axis','fontsize',14);
    uicontrol('style','slider',...
         'position',[150 0 300 20],'Callback',@(hObject, event) updateplot(hObject, event,x,y,z),...
         'value',a, 'min',0, 'max',4);
end

function [u ,v, w] = H_feild (x,y,z,r0)
    r = sqrt ( x .^ 2 + y .^ 2 + z .^ 2 );
    fun = @(t) z.*r0.*cos(t)./((sqrt(r.^2-2.*r0.*(x.*cos(t)+y.*sin(t))+r0.^2)).^3);
    u = integral( fun , 0 , 2 .* pi );
    fun = @(t) z.*r0.*sin(t)./((sqrt(r.^2-2.*r0.*(x.*cos(t)+y.*sin(t))+r0.^2)).^3);
    v = integral( fun , 0 , 2 .* pi );
    fun = @(t) r0.^2./((sqrt(r.^2-2.*r0.*(x.*cos(t)+y.*sin(t))+r0.^2)).^3);
    w = integral( fun , 0 , 2 .* pi );
    fun = @(t) x.*r0.*cos(t)./((sqrt(r.^2-2.*r0.*(x.*cos(t)+y.*sin(t))+r0.^2)).^3);
    h1 = integral( fun , 0 , 2 .* pi );
    fun = @(t) y.*r0.*sin(t)./((sqrt(r.^2-2.*r0.*(x.*cos(t)+y.*sin(t))+r0.^2)).^3);
    h2 = integral( fun , 0 , 2 .* pi );
    w = w - h1 - h2;
    if ( isnan(w) )
        u = 0;
        v = 0;
        w = 0;
        return;
    end

%{
    mag = sqrt ( u^2 + v^2 + w^2);
    u = u ./ mag;
    v = v ./ mag;
    w = w ./ mag;
%}
end

function updateplot(hObject,~,x,y,z)
    a = round (get(hObject,'Value')) - 0.1;
    U = zeros(size(x));
    V = zeros(size(x));
    W = zeros(size(x));
    for i = 1 : numel(x)
        [ U(i) , V(i) , W(i) ] = H_feild (x(i) ,y(i) ,z(i),a);
        disp (ceil(i./numel(x) * 100));
    end
    figure;
    [r,theta] = meshgrid(0:a,linspace(0,2*pi,20));
    X = r.*cos(theta);
    Y = r.*sin(theta);
    Z = X - X;
    surf(X,Y,Z);
    hold on;
    quiver3(x,y,z,U,V,W,2);
    disp(a);
    drawnow;
end