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

Transforming a matlab code!

  1. Jan 24, 2010 #1
    Hi fellows!!
    I need a bit of help with transforming a matlab code that simulate magnetic field of two concentric solenoid TO a plot that simulate the magnetic fields of high voltage transmission lines (specificly two and four lines), how can it be modified? here is the code:

    % GKLmagnet2.m: magnetic profile for TWO ring magnets.
    % This code simulates the magnetic field on the z-axis of
    % two concentric solenoidal superconducting magnets of specified
    % size for use in the 140GHz Gyroklystron (GKL) experiment, MIT.
    % With two coils, it's possible to have one large region that
    % meets the tolerance spec, otherwise two fields are created which
    % may not meet spec. If the field is not continuous (one field)
    % to within spec, the varialbe 'flag' is set to unity (otherwise 0),
    % and a warning is displayed.
    %
    % run GKLmagnet2
    %
    % Last edit by Colin Joye, 7/10/03
    clear all;
    opengl autoselect;
    % -------------- USER INPUT -------------------------
    Bmax = 5.12; %[Telsa] maximum magnetic field in uniform region
    Tol = 0.3; %[%] percent homogeneity of uniform region

    % Magnet size. The magnet is a finite-length finite-thickness hollow cylinder
    % modeled as a single turn rectangular cross-section current density.
    % The magnet is centered at z=0, extending 'Length'/2 above and below z=0.
    RadIn = 15; %[cm] inner radius of coil
    RadOut = 18; %[cm] outer radius of coil
    Length = 6; %[cm] coil length
    Separ = 18.0; %[cm] center separation distance > Length
    zaxis = 30; %[cm] length above z=0 to view the field profile

    % -------------- End USER INPUT ---------------------

    % extra user constants:
    transparency = 0.25; % controls transparency of slices & patches. '= 1.0' is opaque.
    cyl_frac = 0.74; % azimuthal fraction of cylinder to plot. '= 1.0' for full cyl.

    % check if length is larger than separation.
    if(Separ<Length),
    disp(['Error, Separation must be greater than Length']);
    break;
    end
    flag = 0;

    % Define z-axis vector for field profile.
    z = zaxis*[-1:1/2e4:1];
    %z = [0 1 2 5 10 15 20 25 30 31]; %Kreischer's Table II.

    % Draw cylinders
    sep = Separ/2;
    L = Length;
    [Xi,Yi,Zi] = cylinder(RadIn,100);
    [Xo,Yo,Zo] = cylinder(RadOut,100);
    s = floor( cyl_frac*size(Xi,2) )+1;
    cylXi = (Zi(:,1:s)-0.5)*L;
    cylYi = -Xi(:,1:s);
    cylZi = -Yi(:,1:s);
    cylXo = (Zo(:,1:s)-0.5)*L;
    cylYo = -Xo(:,1:s);
    cylZo = -Yo(:,1:s);

    % Applying Biot-Savart law to a finite-length finite-thickness cylinder
    % with a rectangular cross-section of uniform current density, we get,
    % using an analytic (exact) result from Maple v8.0:
    r1 = RadIn;
    r2 = RadOut;
    Bz = zeros(1,length(z));
    for M=1:2, % loop calculates for both magnets.
    r1p = r1 + sqrt(r1^2 + (L/2 + (z + sep)).^2);
    r1n = r1 + sqrt(r1^2 + (L/2 - (z + sep)).^2);
    r2p = r2 + sqrt(r2^2 + (L/2 + (z + sep)).^2);
    r2n = r2 + sqrt(r2^2 + (L/2 - (z + sep)).^2);
    C0 = (z + sep) .* log( r1n./r1p .* r2p./r2n );
    C1 = L/2 * log( r2n .* r2p ./ r1n ./ r1p );
    A0 = 1 / (r2 - r1) / L; % multiplicative value
    Hz = A0 * (C0 + C1);
    Bz = Bz + Hz * Bmax / max(Hz); % normalizing for desired value of Bmax.
    sep = -sep; % reverse sign for other magnet, reverse again on exit.
    end
    % Renormalize again, since addition of field messes up first normalization.
    Bz = Bz*Bmax/max(Bz);

    % Find uniformity zones.
    zmid = find(z==0);
    u = find( abs(Bz(zmid:end) - Bmax)/Bmax <= (Tol/100) );
    UniformLength = 2*(z(u(end)+zmid)-z(u(1)+zmid));
    if(u(1)~=1), %check if the field is still continuous.
    disp(['Warning, field is broken into two regions']);
    flag = 1;
    end

    % ---------------- Figures -------------------
    figure(1)
    clf(1)
    hold on;
    surf(cylXi-sep,cylYi,cylZi); % inner cylinder
    surf(cylXo-sep,cylYo,cylZo); % outer cylinder
    surf(cylXi+sep,cylYi,cylZi); % inner cylinder
    surf(cylXo+sep,cylYo,cylZo); % outer cylinder
    colormap autumn;
    shading flat;

    % Plot dressups
    a = 1.2*RadOut;
    b = 1.1*zaxis;
    axis equal;
    axis([-b b -a a -a a]);
    view([-15 15]);

    % define scaling variable for making field plot look better.
    scale = floor(a/max(Bz)); % scale line and grid to view field better
    if(scale<=1) % make sure you dont have scale = 0, 3, 6, 7, 9, etc.
    scale=1;
    elseif(mod(scale,2)~=0 & mod(scale,5)~=0)
    scale=scale-1;
    end

    % Add transparent slices. Transparency controlled by "alpha"
    % Horizontal plane slice
    h = patch([-b -b b b],[-a a a -a],[0 0 0 0],[0.2 1 0.5]);
    alpha(h,transparency);
    set(h,'EdgeAlpha',0);
    % Vertical plane slice
    h = patch([-b -b b b],[0 0 0 0],[-a a a -a],[0.2 1 0.5]);
    alpha(h,transparency);
    set(h,'EdgeAlpha',0);
    % Z=0 plane slice
    %h=patch([0 0 0 0],[-a -a a a],[-a a a -a],[0.7 1 0]);
    %alpha(h,.3);
    %set(h,'EdgeAlpha',0);

    % solidify cylinders
    map = colormap;
    px = cylXi(:,end)'*[1 0 0 1;0 1 1 0];
    py = [cylYi(:,end)' cylYo(:,end)'];
    pz = [cylZi(:,end)' cylZo(:,end)'];
    h = patch(px-sep, py, pz,[0.7 0.8 0.9]);
    set(h,'EdgeAlpha',transparency);
    h = patch(px+sep, py, pz,[0.7 0.8 0.9]);
    set(h,'EdgeAlpha',transparency);
    px = cylXi(:,1)'*[1 0 0 1;0 1 1 0];
    py = [cylYi(:,1)' cylYo(:,1)'];
    pz = [cylZi(:,1)' cylZo(:,1)'];
    h = patch(px-sep, py, pz,[0.7 0.8 0.9]);
    set(h,'EdgeAlpha',transparency);
    h = patch(px+sep, py, pz,[0.7 0.8 0.9]);
    set(h,'EdgeAlpha',transparency);

    % Put grid on field line patch
    ytick = get(gca,'Ztick');
    ytick = (ytick + max(ytick))/2; % remove negative ticks and interpolate.
    ztick = get(gca,'Xtick');
    Ly = length(ytick);
    Lz = length(ztick);
    y0 = zeros(1,Ly);
    y1 = b*ones(1,Ly);
    z0 = zeros(1,Lz);
    z1 = a*ones(1,Lz);
    line([1;1]*ztick,[z0;z0],[z0;z1],'Color',0.4*[1 1 1],'LineStyle',:); %vert lines
    line([-y1;y1],[y0;y0],[1;1]*ytick,'Color',0.4*[1 1 1],'LineStyle',:);%horiz lines
    % Add axis line and numbers on new grid.
    % Vertical:
    line(-[b b],[0 0],[0 a],'Color',[0 0.3 1]); % vert axis line for numbers
    line(-[1;1.03]*y1,[y0;y0],[1;1]*ytick,'Color',[0 0.3 1]); % draw tick lines.
    h = text(-1.03*y1',y0',ytick',num2str(ytick'/scale),'HorizontalAlignment','right');
    set(h,'Color',[0 0.3 1]);
    pz = length(get(h,'String'))-3;
    h = text(-1.1*b-pz,0,mean(ytick),'Field [T] (blue)','HorizontalAlignment','center');
    set(h,'Color',[0 0.3 1],'Rotation',90);
    % Horizontal:
    line(-[b zaxis],[0 0],[0 0],'Color',[0 0.3 1]); % horiz axis line for numbers
    line( [zaxis b],[0 0],[0 0],'Color',[0 0.3 1]); % horiz axis line for numbers
    line([1;1]*ztick,[z0;z0],-[0;0.03]*z1,'Color',[0 0.3 1]); % draw tick lines.
    h = text(ztick',z0',-0.055*z1',num2str(ztick'),'HorizontalAlignment','center');
    set(h,'Color',[0 0.3 1]);
    h = line([-zaxis zaxis],[0,0],[0,0]); % z-axis line.
    set(h,'Color','red','LineWidth',2.5);

    % Add plot of field line
    h = plot3(z,zeros(1,length(Bz)),scale*Bz,'b');
    set(h,'Color','blue','Linewidth',2);
    h = text(0,0,0.1*Bmax,'Uniform Region');
    set(h,'Rotation',90,'Color',[0 0.4 0.4])

    % Add patch showing uniform region within 'Tol' percent.
    uz0 = u(1)+zmid;
    uz1 = u(end)+zmid;
    px = [z(uz0) z(uz1) z(uz1) z(uz0)];
    py = [0 0 0 0];
    pz = scale*[0 0 Bz(u(1)+zmid) Bz(u(end)+zmid)];
    h = patch(px, py, pz,[1 0.2 0.5]);
    alpha(h,1.5*transparency);
    set(h,'EdgeAlpha',0);
    px = -[z(uz0) z(uz1) z(uz1) z(uz0)];
    h = patch(px, py, pz,[1 0.2 0.5]);
    alpha(h,1.5*transparency);
    set(h,'EdgeAlpha',0);
    %h = plot3(z(u),zeros(1,length(u)),max(Bz)*ones(1,length(u)),'m');
    %set(h,'LineWidth',2');
    h = text(0,0,0.1*Bmax,'Uniform Region');
    set(h,'Rotation',90,'Color',[0 0.4 0.3])

    % Title and Labels
    title(['Magnet strucure w/ Field: ' num2str(Bmax) ' T max; R1 = ' num2str(r1) ...
    ' cm; R2 = ' num2str(r2) ' cm; L = ' num2str(L) ' cm; S = ' num2str(Separ) ...
    ' cm; ' num2str(Tol) '% Uniform length: ' num2str(UniformLength) ' cm']);
    xlabel('Z-axis [cm]');
    ylabel('X-axis [cm]');
    zlabel('Y-axis [cm]');

    hold off;

    % ------------- END GKLmagnet2.m ---------------------------

    PS: It dooes what i need but i need the magnetic fields of those high voltage transsmision lines, and i dont know to much of programing, sooo i will apreciate a loooot if you can help me out =)
     
  2. jcsd
  3. Jan 30, 2010 #2
    I could use some help you know...
     
  4. Jan 30, 2010 #3

    uart

    User Avatar
    Science Advisor

    Ok then. Hang around the Computer Lab while other students from your class are working on this project. After they leave, look though all the waste paper bins in the hope they've disgarded some hard copies of partially debugged code. Then you'll have some code that's a little bit closer to what you're trying to achieve to try and hack.
     
    Last edited: Jan 31, 2010
  5. Jan 31, 2010 #4

    uart

    User Avatar
    Science Advisor

    Now for a serious reply.

    You haven't got any help because you've posted a silly question. No one here is likely to be writting code by "transforming" code from a completely different problem, it's just not how it's done.

    You should be starting here with the relevent equations for the actual problem at hand. If you don't know what these are then that's what your first question should be in relation to.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Transforming a matlab code!
  1. Help with code (Replies: 1)

  2. Want MATLAB code (Replies: 4)

  3. Inductor color code (Replies: 13)

Loading...